From 0f9eb6550b5ee8af3dc0d435f6d692025b985728 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 2 Oct 2026 19:41:17 +0200 Subject: [PATCH 1/3] Preserve Cartesian vacuum query coordinates and expose unavailable stencils --- gpec/gpout.f | 19 ++- tests/.gitignore | 1 + tests/CMakeLists.txt | 79 +++++++++++ tests/README.md | 43 ++++++ tests/expect_rejection.py | 13 ++ tests/fixtures/vac.in | 82 ++++++++++++ tests/generate_query_fixture.py | 36 +++++ tests/test_cartesian_postprocessors.f90 | 50 +++++++ tests/test_vacuum_query.f90 | 138 +++++++++++++++++++ vacuum/vacuum_ma.f | 169 ++++++++++++++++++++++-- 10 files changed, 617 insertions(+), 13 deletions(-) create mode 100644 tests/.gitignore create mode 100644 tests/CMakeLists.txt create mode 100644 tests/README.md create mode 100644 tests/expect_rejection.py create mode 100644 tests/fixtures/vac.in create mode 100644 tests/generate_query_fixture.py create mode 100644 tests/test_cartesian_postprocessors.f90 create mode 100644 tests/test_vacuum_query.f90 diff --git a/gpec/gpout.f b/gpec/gpout.f index fa5ccf325..19a36e5d2 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -4320,6 +4320,7 @@ END SUBROUTINE gpout_xbtangent c write perturbed rzphi components on rzphi grid. c----------------------------------------------------------------------- SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) + USE vacuum_mod, ONLY: check_cartesian_postprocessors c----------------------------------------------------------------------- c declaration. c----------------------------------------------------------------------- @@ -4508,7 +4509,9 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) $ "Computing external vacuum fields by surface currents" CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vbr,vbz,vbp) + $ nr,nz,vgdl,vgdr,vgdz,vbr,vbz,vbp, + $ preserve_query_coordinates=.TRUE., + $ query_r=gdr,query_z=gdz) IF (helicity<0) THEN vbr=CONJG(vbr) vbz=CONJG(vbz) @@ -4534,7 +4537,9 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) bnomn=bnomn-bnimn CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp) + $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp, + $ preserve_query_coordinates=.TRUE., + $ query_r=gdr,query_z=gdz) IF (helicity<0) THEN vpbr=CONJG(vpbr) @@ -4578,6 +4583,8 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) ENDDO ENDDO + CALL check_cartesian_postprocessors(btr,btz,btp, + $ divzero_flag,chebyshev_flag) IF (divzero_flag) THEN CALL gpeq_rzpdiv(nr,nz,gdr,gdz,btr,btz,btp) ENDIF @@ -4742,7 +4749,9 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) bnomn=bnomn-bnimn CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp) + $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp, + $ preserve_query_coordinates=.TRUE., + $ query_r=gdr,query_z=gdz) IF (helicity<0) THEN vpbr=CONJG(vpbr) vpbz=CONJG(vpbz) @@ -4755,7 +4764,9 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) $ "Computing vacuum fields by external surface currents" CALL gpvacuum_bnormal(psilim,bnimn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vvbr,vvbz,vvbp) + $ nr,nz,vgdl,vgdr,vgdz,vvbr,vvbz,vvbp, + $ preserve_query_coordinates=.TRUE., + $ query_r=gdr,query_z=gdz) IF (helicity<0) THEN vvbr=CONJG(vvbr) vvbz=CONJG(vvbz) diff --git a/tests/.gitignore b/tests/.gitignore new file mode 100644 index 000000000..84c048a73 --- /dev/null +++ b/tests/.gitignore @@ -0,0 +1 @@ +/build/ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 000000000..422b3f138 --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,79 @@ +cmake_minimum_required(VERSION 3.20) +project(gpec_cartesian_query_tests LANGUAGES C Fortran) +set(CMAKE_BUILD_TYPE Release CACHE STRING "Build configuration") +get_filename_component(GPEC_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/.." ABSOLUTE) +if(NOT CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") + message(FATAL_ERROR "This standalone native test adapter supports GNU Fortran") +endif() +find_package(BLAS REQUIRED) +find_package(LAPACK REQUIRED) +find_package(Python3 REQUIRED COMPONENTS Interpreter) +find_path(NETCDF_INCLUDE netcdf.mod REQUIRED) +find_library(NETCDFF_LIBRARY netcdff REQUIRED) +find_library(NETCDF_LIBRARY netcdf REQUIRED) +file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/version.inc" + " CHARACTER(len=*), PARAMETER :: version='cartesian-query-tests'\n") + +function(gpec_library name) + file(GLOB src "${GPEC_SOURCE}/${name}/*.f" "${GPEC_SOURCE}/${name}/*.F" + "${GPEC_SOURCE}/${name}/*.f90" "${GPEC_SOURCE}/${name}/*.F90" + "${GPEC_SOURCE}/${name}/*.c") + list(FILTER src EXCLUDE REGEX "/(pentrc\\.F90|slayer\\.f)$") + add_library(g_${name} STATIC ${src}) + set_target_properties(g_${name} PROPERTIES + Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/modules/${name}") + target_include_directories(g_${name} PUBLIC + "${CMAKE_CURRENT_BINARY_DIR}/modules/${name}") + target_include_directories(g_${name} PRIVATE "${NETCDF_INCLUDE}" + "${CMAKE_CURRENT_BINARY_DIR}") + target_compile_options(g_${name} PRIVATE + $<$:-ffixed-line-length-none;-fallow-argument-mismatch;-frecursive;-fopenmp>) + target_link_libraries(g_${name} PUBLIC ${ARGN} "${NETCDFF_LIBRARY}" + "${NETCDF_LIBRARY}" ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + target_link_options(g_${name} PUBLIC -fopenmp) +endfunction() +gpec_library(equil) +gpec_library(lsode g_equil) +gpec_library(vacuum g_equil g_lsode) +gpec_library(pentrc g_equil g_lsode) +gpec_library(coil g_equil) +gpec_library(slayer g_equil g_lsode) +foreach(code dcon gpec) + file(GLOB src "${GPEC_SOURCE}/${code}/*.f" "${GPEC_SOURCE}/${code}/*.F") + add_executable(${code} ${src}) + set_target_properties(${code} PROPERTIES + Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/modules/${code}") + target_include_directories(${code} PRIVATE "${NETCDF_INCLUDE}" + "${CMAKE_CURRENT_BINARY_DIR}" "${CMAKE_CURRENT_BINARY_DIR}/modules/${code}") + target_compile_options(${code} PRIVATE -ffixed-line-length-none + -fallow-argument-mismatch -frecursive -fopenmp) + target_link_libraries(${code} PRIVATE g_pentrc g_vacuum g_slayer g_coil) +endforeach() + +set(QUERY_CASE "${CMAKE_CURRENT_BINARY_DIR}/query_case") +file(MAKE_DIRECTORY "${QUERY_CASE}") +configure_file(fixtures/vac.in "${QUERY_CASE}/vac.in" COPYONLY) +execute_process(COMMAND "${Python3_EXECUTABLE}" + "${CMAKE_CURRENT_SOURCE_DIR}/generate_query_fixture.py" "${QUERY_CASE}" + RESULT_VARIABLE FIXTURE_STATUS) +if(NOT FIXTURE_STATUS EQUAL 0) + message(FATAL_ERROR "Synthetic native query fixture generation failed") +endif() +enable_testing() +add_executable(test_vacuum_query test_vacuum_query.f90) +target_link_libraries(test_vacuum_query PRIVATE g_vacuum) +target_compile_options(test_vacuum_query PRIVATE -ffree-line-length-none -fcheck=all) +add_test(NAME vacuum_query_exact COMMAND test_vacuum_query) +set_tests_properties(vacuum_query_exact PROPERTIES + WORKING_DIRECTORY "${QUERY_CASE}" TIMEOUT 90) +add_executable(test_cartesian_postprocessors test_cartesian_postprocessors.f90) +target_link_libraries(test_cartesian_postprocessors PRIVATE g_vacuum) +add_test(NAME cartesian_postprocessors_finite_and_raw COMMAND test_cartesian_postprocessors) +foreach(mode chebyshev divzero) + add_test(NAME cartesian_postprocessors_reject_${mode} + COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/expect_rejection.py" + "$" "${mode}") +endforeach() +set_tests_properties(cartesian_postprocessors_finite_and_raw + cartesian_postprocessors_reject_chebyshev cartesian_postprocessors_reject_divzero + PROPERTIES TIMEOUT 15) diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 000000000..67097c782 --- /dev/null +++ b/tests/README.md @@ -0,0 +1,43 @@ +# Cartesian vacuum query tests + +This standalone adapter builds the checkout's native libraries and the GPEC and +DCON executables, then runs the Cartesian export regression tests. It uses no +external case or private path. The upstream installation makefiles remain +unchanged. + +Requirements are GNU C/Fortran, CMake 3.20 or newer, BLAS/LAPACK, NetCDF C and +Fortran development libraries, and Python 3 (standard library only). CMake's +normal prefix/library/include search paths select those dependencies. Run from +this directory with the local Fortran workflow: + +```sh +FO_JOBS=2 OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 fo +``` + +For installations without `fo`, the standalone adapter can also be run using +the standard CMake/CTest commands: + +```sh +cmake -S tests -B tests/build -DCMAKE_BUILD_TYPE=Release +cmake --build tests/build --parallel 2 +OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 ctest --test-dir tests/build --output-on-failure +``` + +Run the CMake commands from the checkout root. All generated inputs and native +diagnostics stay in the chosen build directory. The synthetic fixture is a +circular clockwise surface with 128 unique nodes, five poloidal modes, signed complex +normal-field coefficients and n=3. It is a bounded numerical source, not a +device equilibrium or a physical forcing convergence study. + +`vacuum_query_exact` independently evaluates native CHI at double query +coordinates, then checks pickup and MSCFLD forwarding, source/stencil masks, +and identical absent/false legacy optional behavior. It includes a resolved +query inside the old snapping band and unresolved source/crossing queries. +The independent reference bypasses both coordinate construction and pickup. +The remaining tests permit finite fields and raw masked exports, and require +the specific native rejection diagnostics for each optional postprocessor. +They do not treat an arbitrary crash as a successful rejection. + +These tests establish the query/output contract. They do not establish source +quadrature convergence, accuracy near a singular sheet, arbitrary malformed +derivative settings, global plasma response or physical forcing validity. diff --git a/tests/expect_rejection.py b/tests/expect_rejection.py new file mode 100644 index 000000000..bfd712f8f --- /dev/null +++ b/tests/expect_rejection.py @@ -0,0 +1,13 @@ +"""Require the native guard's specific diagnostic, rather than any failure.""" +import subprocess +import sys + +executable, mode = sys.argv[1:] +assert mode in ("chebyshev", "divzero") +result = subprocess.run([executable, mode], capture_output=True, text=True, timeout=10) +output = result.stdout + result.stderr +expected = f"disable {mode}_flag for raw masked exports" +if result.returncode == 0 or expected not in output: + print(output) + raise SystemExit(f"Expected native {mode} rejection was not observed") +print(f"PASS native {mode} rejection with specific diagnostic") diff --git a/tests/fixtures/vac.in b/tests/fixtures/vac.in new file mode 100644 index 000000000..7cc8612d0 --- /dev/null +++ b/tests/fixtures/vac.in @@ -0,0 +1,82 @@ +&MODES + mth=128 + xiin(1:9)=0 0 0 0 0 0 0 1 0 + lsymz=.true. + leqarcw=1 + lzio=0 + lgato=0 + lrgato=0 +/ +&DEBUGS + checkd=.false. + check1=.false. + check2=.false. + checke=.false. + checks=.false. + wall=.false. + lkplt=0 +/ +&VACDAT + ishape=6 + nsing=128 + epsq=1e-5 + noutv=37 + idgt=6 + idot=0 + idsk=0 + delg=15.01 + delfac=0.001 + cn0=1 +/ +&SHAPE + ipshp=0 + xpl=100 + apl=1 + a=0.2413 + b=170 + bpl=1 + dpl=0 + r=1 + abulg=0.932 + bbulg=17 + tbulg=0.02 + qain=2.5 +/ +&DIAGNS + lkdis=.false. + ieig=0 + iloop=0 + lpsub=1 + nloop=128 + nloopr=0 + nphil=3 + nphse=1 + xofsl=0 + ntloop=32 + aloop=0.01 + bloop=1.6 + dloop=0.5 + rloop=1 + deloop=0.001 + mx=21 + mz=21 + nph=0 + nxlpin=2 + nzlpin=2 + epslp=0.02 + xlpmin=0.7 + xlpmax=2.7 + zlpmin=-1.5 + zlpmax=1.5 + linterior=2 +/ +&SPRK + nminus=0 + nplus=0 + mphi=16 + lwrt11=0 + civ=0 + lff=0 + ff=1.6 + fv=1.6 +/ diff --git a/tests/generate_query_fixture.py b/tests/generate_query_fixture.py new file mode 100644 index 000000000..7d96a7d86 --- /dev/null +++ b/tests/generate_query_fixture.py @@ -0,0 +1,36 @@ +"""Write a bounded synthetic circular source in native VACIN5 format.""" +import math +from pathlib import Path +import sys + +out = Path(sys.argv[1]) +out.mkdir(parents=True, exist_ok=True) +nodes, low, high, n = 128, -2, 2, 3 +center, radius, zcenter = 6.2, 0.62, 0.0 +theta = [i / nodes for i in range(nodes + 1)] +R = [center + radius * math.cos(2 * math.pi * t) for t in theta[:-1]] +# The native vacuum contour convention is clockwise in the (R,Z) plane. +Z = [zcenter - radius * math.sin(2 * math.pi * t) for t in theta[:-1]] +R.append(R[0]) +Z.append(Z[0]) +coefficients = {-1: 0.04 - 0.03j, 0: 0.02 + 0.07j, 2: -0.015 + 0.01j} +lines = ["scalars", "", "2", "2", f"{center-radius:.17e}", + f"{center+radius:.17e}", f"{zcenter-radius:.17e}", + f"{zcenter+radius:.17e}", str(nodes), str(low), str(high), str(n), "2.5"] + + +def block(label, values): + lines.extend([label, ""]) + for i in range(0, len(values), 4): + lines.append(" ".join(f"{v:.17e}" for v in values[i:i+4])) + + +block("Poloidal Coordinate Theta", theta) +block("Normalized Arclength", theta) +block("Major Radius", R) +block("Vertical Coordinate", Z) +block("Toroidal Angle Shift", [0.0] * (nodes + 1)) +block("Normal Field Real", [coefficients.get(m, 0j).real for m in range(low, high + 1)]) +block("Normal Field Imaginary", [coefficients.get(m, 0j).imag for m in range(low, high + 1)]) +(out / "vacin5").write_text("\n".join(lines) + "\n") +(out / "query_contract.in").write_text(f"{nodes} {high-low+1}\n") diff --git a/tests/test_cartesian_postprocessors.f90 b/tests/test_cartesian_postprocessors.f90 new file mode 100644 index 000000000..c17152ca7 --- /dev/null +++ b/tests/test_cartesian_postprocessors.f90 @@ -0,0 +1,50 @@ +! Exercise the production guard with manufactured finite and unavailable fields. +program test_cartesian_postprocessors + use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use vacuum_mod, only: check_cartesian_postprocessors + implicit none + complex(dp) :: br(3,4), bz(3,4), bp(3,4) + real(dp) :: unavailable + character(24) :: mode + integer :: component + + call get_command_argument(1, mode) + unavailable = ieee_value(0.0_dp, ieee_quiet_nan) + br = cmplx(1.0_dp, 2.0_dp, dp) + bz = cmplx(-3.0_dp, 4.0_dp, dp) + bp = cmplx(5.0_dp, -6.0_dp, dp) + call check_cartesian_postprocessors(br, bz, bp, .true., .true.) + select case (trim(mode)) + case ('chebyshev') + bz(2,3) = cmplx(-3.0_dp, unavailable, dp) + call check_cartesian_postprocessors(br, bz, bp, .false., .true.) + print *, 'UNEXPECTED: Chebyshev accepted an unavailable total field' + case ('divzero') + bp(2,3) = cmplx(unavailable, -6.0_dp, dp) + call check_cartesian_postprocessors(br, bz, bp, .true., .false.) + print *, 'UNEXPECTED: divzero accepted an unavailable total field' + case default + do component = 1, 6 + br = cmplx(1.0_dp, 2.0_dp, dp) + bz = cmplx(-3.0_dp, 4.0_dp, dp) + bp = cmplx(5.0_dp, -6.0_dp, dp) + select case (component) + case (1) + br(2,3) = cmplx(unavailable, 2.0_dp, dp) + case (2) + br(2,3) = cmplx(1.0_dp, unavailable, dp) + case (3) + bz(2,3) = cmplx(unavailable, 4.0_dp, dp) + case (4) + bz(2,3) = cmplx(-3.0_dp, unavailable, dp) + case (5) + bp(2,3) = cmplx(unavailable, -6.0_dp, dp) + case (6) + bp(2,3) = cmplx(5.0_dp, unavailable, dp) + end select + call check_cartesian_postprocessors(br, bz, bp, .false., .false.) + end do + print *, 'PASS finite postprocessor inputs and raw masked exports' + end select +end program test_cartesian_postprocessors diff --git a/tests/test_vacuum_query.f90 b/tests/test_vacuum_query.f90 new file mode 100644 index 000000000..e38e81646 --- /dev/null +++ b/tests/test_vacuum_query.f90 @@ -0,0 +1,138 @@ +! Native query oracle; supply query_contract.in, vac.in and vacin5 externally. +! The reference calls CHI directly at doubles, without pickup or loops. +program test_vacuum_query + use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: ieee_arithmetic, only: ieee_is_nan, ieee_is_finite + use vacuum_mod, only: defglo, ent33, make_bltobp, diaplt, chi, pickup, mscfld + use vglobal_mod, only: global_alloc, global_dealloc, kernelsign, ntsin0, & + nths0, nfm, mtot, ndimlp, xobp, zobp, xloop, zloop, ldcon, lgpec, ieig, & + iotty, outpest, inmode, outmod, farwal, mth, lmax, lmin, vacmat, vacmtiu, & + xpla, zpla, plrad, deloop, epslp, grri, nths, nloop, nloopr, xplap, & + zplap, bnlr, bnli, n, nxlpin, nzlpin + implicit none + integer :: surface_nodes, modes, iu, peak, side, l, k + real(dp) :: qr(0:1,0:1), qz(0:1,0:1), rr(0:1,0:1), zz(0:1,0:1) + real(dp) :: lr(0:1,0:1), lz(0:1,0:1), delta, band, center_r, center_z + integer :: flags(0:1,0:1), legacy_flags(0:1,0:1) + complex(dp) :: br(0:1,0:1), bz(0:1,0:1), bp(0:1,0:1) + complex(dp) :: legacy_br(0:1,0:1), legacy_bz(0:1,0:1), legacy_bp(0:1,0:1) + complex(dp) :: reference(3,2) + complex(dp), allocatable :: matrix_ref(:,:), matrix_out(:,:) + real(dp), allocatable :: re(:,:), im(:,:), region(:), cc(:,:), ss(:,:) + real(dp) :: field_error, field_scale + + open(newunit=iu,file='query_contract.in',status='old') + read(iu,*) surface_nodes,modes + close(iu) + if(surface_nodes<16 .or. surface_nodes>512) error stop 'Bounded surface size required' + if(modes<1 .or. modes>40) error stop 'Bounded mode count required' + kernelsign=1;ntsin0=surface_nodes+1;nths0=surface_nodes;nfm=modes;mtot=modes + call global_alloc(nths0,nfm,mtot,ntsin0) + ndimlp=4 + allocate(xobp(4),zobp(4),xloop(4),zloop(4)) + call defglo(surface_nodes) + farwal=.true.;ldcon=0;lgpec=1;ieig=0 + open(iotty,file='mscvac.out',status='replace') + open(outpest,file='pestotv',status='replace',form='formatted') + open(inmode,file='vac.in',status='old',form='formatted') + open(outmod,file='modovmc',status='replace',form='formatted') + call ent33 + if(mth/=surface_nodes .or. lmax(1)-lmin(1)+1/=modes) error stop 'Source basis mismatch' + call make_bltobp + call diaplt + allocate(matrix_ref(modes,modes),matrix_out(modes,modes)) + matrix_ref=cmplx(vacmat,vacmtiu,dp) + open(newunit=iu,file='query_native_matrix.dat',status='replace') + do k=1,modes + do l=1,modes + write(iu,'(2i6,2es26.17)') l,k,real(matrix_ref(l,k)),aimag(matrix_ref(l,k)) + end do + end do + close(iu) + plrad=.5_dp*(maxval(xpla(:mth))-minval(xpla(:mth))) + delta=plrad*deloop;band=plrad*epslp + if(delta<=0 .or. band<=8*delta) error stop 'Oracle requires a resolvable legacy snap band' + peak=maxloc(xpla(:mth),dim=1) + center_r=.5_dp*(maxval(xpla(:mth))+minval(xpla(:mth))) + center_z=.5_dp*(maxval(zpla(:mth))+minval(zpla(:mth))) + qr(0,0)=xpla(peak)-.5_dp*band;qz(0,0)=zpla(peak) + qr(0,1)=center_r+.13_dp*plrad;qz(0,1)=center_z+.07_dp*plrad + qr(1,0)=xpla(peak);qz(1,0)=zpla(peak) + qr(1,1)=xpla(peak)-.5_dp*delta;qz(1,1)=zpla(peak) + + ! The reference independently evaluates only the two resolved interior queries. + allocate(re(5,4),im(5,4),region(4),cc(nths,modes),ss(nths,modes)) + cc=0;ss=0 + do l=1,modes + cc(:mth,l)=grri(:mth,l) + ss(:mth,l)=grri(:mth,modes+l) + end do + nloop=2;nloopr=0 + xloop(1:2)=qr(0,:);zloop(1:2)=qz(0,:) + region=-1;re=0;im=0 + do side=1,5 + xobp(1:2)=xloop(1:2);zobp(1:2)=zloop(1:2) + select case(side) + case(1) + zobp(1:2)=zobp(1:2)+delta + case(2) + zobp(1:2)=zobp(1:2)-delta + case(3) + xobp(1:2)=xobp(1:2)+delta + case(4) + xobp(1:2)=xobp(1:2)-delta + end select + call chi(xpla,zpla,xplap,zplap,-1,cc,ss,mth,1,re,im,side,bnlr,bnli,region) + end do + do k=1,2 + reference(1,k)=cmplx(re(3,k)-re(4,k),im(3,k)-im(4,k),dp)/(2*delta) + reference(2,k)=cmplx(re(1,k)-re(2,k),im(1,k)-im(2,k),dp)/(2*delta) + reference(3,k)=cmplx(n*im(5,k),-n*re(5,k),dp)/xloop(k) + end do + nxlpin=2;nzlpin=2 + call pickup(bnlr,bnli,1,1,flags,rr,zz,br,bz,bp, & + preserve_query_coordinates=.true.,query_r=qr,query_z=qz) + call check_exact_query() + close(iotty);close(outpest);close(inmode);close(outmod) + call global_dealloc + call cleanup + + ! Exercise actual MSCFLD forwarding; its legacy GPEC path does not write WV. + farwal=.true. + call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & + 1,1,flags,rr,zz,br,bz,bp,preserve_query_coordinates=.true.,query_r=qr,query_z=qz) + call check_exact_query() + + farwal=.true. + call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & + 1,1,legacy_flags,lr,lz,legacy_br,legacy_bz,legacy_bp) + farwal=.true. + call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & + 1,1,flags,rr,zz,br,bz,bp,preserve_query_coordinates=.false.) + if(any(flags/=legacy_flags) .or. any(rr/=lr) .or. any(zz/=lz)) & + error stop 'Legacy optional contract changed coordinates or flags' + if(any(br/=legacy_br) .or. any(bz/=legacy_bz) .or. any(bp/=legacy_bp)) & + error stop 'Legacy optional contract changed field values' + print *, 'PASS exact native query coordinates, direct CHI fields, masks and legacy contract' + print *, 'DIRECT_CHI_RELATIVE_FIELD_ERROR',field_error +contains + subroutine check_exact_query() + integer :: j + if(any(rr/=qr) .or. any(zz/=qz)) error stop 'Query coordinates relocated or downcast' + if(any(flags(0,:)/=1)) error stop 'Resolved interior query misclassified' + if(any(flags(1,:)/=-2)) error stop 'Source or crossing-stencil query not masked' + if(.not.all(ieee_is_nan(real(br(1,:))))) error stop 'Invalid BR not unavailable' + if(.not.all(ieee_is_nan(aimag(br(1,:))))) error stop 'Invalid imaginary BR not unavailable' + if(.not.all(ieee_is_nan(real(bz(1,:))))) error stop 'Invalid BZ not unavailable' + if(.not.all(ieee_is_nan(aimag(bz(1,:))))) error stop 'Invalid imaginary BZ not unavailable' + if(.not.all(ieee_is_nan(real(bp(1,:))))) error stop 'Invalid Bphi not unavailable' + if(.not.all(ieee_is_nan(aimag(bp(1,:))))) error stop 'Invalid imaginary Bphi not unavailable' + field_error=0;field_scale=max(1e-20_dp,maxval(abs(reference))) + do j=0,1 + field_error=max(field_error, & + maxval(abs([br(0,j),bz(0,j),bp(0,j)]-reference(:,j+1)))/field_scale) + end do + if(.not.ieee_is_finite(field_error)) error stop 'Nonfinite resolved field' + if(field_error>2e-11_dp) error stop 'Pickup field differs from independent direct CHI query' + end subroutine check_exact_query +end program test_vacuum_query diff --git a/vacuum/vacuum_ma.f b/vacuum/vacuum_ma.f index b4b3da596..3484b8cc7 100644 --- a/vacuum/vacuum_ma.f +++ b/vacuum/vacuum_ma.f @@ -170,7 +170,8 @@ subroutine mscvac(wv,mpert,mtheta,mthvac,complex_flag, c declarations. c----------------------------------------------------------------------- subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, - $ lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp,op_ahgfile) + $ lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp,op_ahgfile, + $ preserve_query_coordinates,query_r,query_z) USE vglobal_mod implicit real(r8) (a-h,o-z) implicit integer (i-n) @@ -186,6 +187,8 @@ subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, complex(r8), parameter :: ifac=(0,1) dimension xi(nfm), xii(nfm), xilnq(nfm), xiilnq(nfm) character(128), intent(in), optional :: op_ahgfile + logical, intent(in), optional :: preserve_query_coordinates + real(r8), intent(in), optional :: query_r(0:,0:),query_z(0:,0:) if (present(op_ahgfile)) then ahgfile = trim(op_ahgfile) @@ -287,7 +290,8 @@ subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, 99 continue call make_bltobp call diaplt - call pickup(bnlr,bnli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) + call pickup(bnlr,bnli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, + $ preserve_query_coordinates,query_r,query_z) c----------------------------------------------------------------------- c termination. c----------------------------------------------------------------------- @@ -621,14 +625,21 @@ subroutine diaplt c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- - subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) + subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, + $ preserve_query_coordinates,query_r,query_z) USE vglobal_mod + USE, INTRINSIC :: ieee_arithmetic, ONLY: + $ ieee_value,ieee_quiet_nan implicit real(r8) (a-h,o-z) implicit integer (i-n) INTEGER, DIMENSION(0:lx,0:lz) :: vgdl REAL(r8), DIMENSION(0:lx,0:lz) :: vgdx,vgdz COMPLEX(r8), DIMENSION(0:lx,0:lz) :: vbx,vbz,vbp + logical, intent(in), optional :: preserve_query_coordinates + real(r8), intent(in), optional :: query_r(0:,0:),query_z(0:,0:) + logical :: keep_query,query_valid(ndimlp) + real(r8) :: clearance,unavailable DIMENSION blr(*), bli(*) ! CHARACTER(130), DIMENSION(10) :: string @@ -643,6 +654,28 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) $ btr(ndimlp),bti(ndimlp), bpr(ndimlp),bpi(ndimlp), $ bphir(ndimlp),bphii(ndimlp) + keep_query = .false. + IF (PRESENT(preserve_query_coordinates)) + $ keep_query = preserve_query_coordinates + IF (keep_query) THEN + IF (.NOT. PRESENT(query_r)) + $ ERROR STOP 'Exact pickup requires query_r' + IF (.NOT. PRESENT(query_z)) + $ ERROR STOP 'Exact pickup requires query_z' + IF (SIZE(query_r,1) /= lx+1 .OR. + $ SIZE(query_r,2) /= lz+1) + $ ERROR STOP 'Exact pickup query_r shape mismatch' + IF (SIZE(query_z,1) /= lx+1 .OR. + $ SIZE(query_z,2) /= lz+1) + $ ERROR STOP 'Exact pickup query_z shape mismatch' + IF (nxlpin /= lx+1 .OR. nzlpin /= lz+1) + $ ERROR STOP 'Exact pickup grid dimensions mismatch' + IF (nloopr /= 0) + $ ERROR STOP 'Exact pickup supports the Cartesian grid only' + ENDIF + query_valid = .true. + unavailable = ieee_value(0.0_r8,ieee_quiet_nan) + ndlp = mth / ntloop jmax1 = lmax(1) - lmin(1) + 1 @@ -680,7 +713,19 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) zmax = max(zmxp,zmxw) endif - CALL loops + IF (keep_query) THEN +c Bypass both legacy surface snapping and REAL(4) grid generation. + nloop = (lx+1)*(lz+1) + DO i = 0, lx + DO j = 0, lz + indxl = i*(lz+1)+j+1 + xloop(indxl) = query_r(i,j) + zloop(indxl) = query_z(i,j) + ENDDO + ENDDO + ELSE + CALL loops + ENDIF nobs = nloop + 3*nloopr @@ -694,6 +739,11 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) ns = mth delx = plrad * deloop delz = plrad * deloop + IF (keep_query) THEN + IF (delx == 0.0_r8 .OR. delz == 0.0_r8) + $ ERROR STOP 'Exact pickup requires a nonzero derivative step' + ENDIF + clearance = 2.0_r8*MAX(ABS(delx),ABS(delz)) igdl = 8 isgchi = -1 @@ -703,6 +753,23 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) xloops(i) = xloop(i) zloops(i) = zloop(i) + IF (keep_query) THEN +c -2: no resolved same-side finite-difference field here. + query_valid(i) = pickup_query_clear(xloop(i),zloop(i), + $ mth,xpla,zpla,clearance) + IF (xloop(i) <= 0.0_r8) query_valid(i) = .false. + IF (.NOT. farwal) THEN + IF (query_valid(i)) query_valid(i) = + $ pickup_query_clear(xloop(i),zloop(i), + $ mw,xwal,zwal,clearance) + ENDIF + IF (.NOT. query_valid(i)) THEN + igdl(i) = -2 + rgdl(i) = 1.0_r8 + GO TO 239 + ENDIF + ENDIF + fintjj = 0.0 DO jj = 1, mth dxjj = xloop(i) - xinf(jj) @@ -731,6 +798,16 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) 238 CONTINUE + IF (keep_query) THEN + IF (igdl(i) == 8) THEN + query_valid(i) = .false. + igdl(i) = -2 + rgdl(i) = 1.0_r8 + GO TO 239 + ENDIF + ENDIF + + IF (.NOT. keep_query) THEN DO j = 1, mth IF ( (xinf(j)-xloop(i))**2 + (zinf(j)-zloop(i))**2 $ < (epslp*plrad)**2 ) THEN ! points on surface for these. @@ -758,6 +835,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) bphii(i) = - n * zchipr(i) / xloops(i) END IF END DO + ENDIF 239 CONTINUE @@ -820,7 +898,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) isg = -1 call chi ( xpla,zpla,xplap,zplap,isg, chiwc,chiws, ns,1, - $ cwrkr,cwrki,nsew, blr,bli,rgdl ) + $ cwrkr,cwrki,nsew, blr,bli,rgdl,query_valid ) do i = 1, nobs chir(nsew,i) = cwrkr(nsew,i) @@ -843,7 +921,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) isg = 1 call chi ( xwal,zwal,xwalp,zwalp,isg,chiwc,chiws, ns,0, - $ cwrkr,cwrki,nsew, blr,bli,rwall ) + $ cwrkr,cwrki,nsew, blr,bli,rwall,query_valid ) do i = 1, nobs chir(nsew,i) = chir(nsew,i) + cwrkr(nsew,i) @@ -854,7 +932,14 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) END DO DO i = 1, nobs - IF (igdl(i) .NE. -1) THEN + IF (.NOT. query_valid(i)) THEN + bxr(i) = unavailable + bxi(i) = unavailable + bzr(i) = unavailable + bzi(i) = unavailable + bphir(i) = unavailable + bphii(i) = unavailable + ELSE IF (igdl(i) .NE. -1) THEN bxr(i) = ( chir(3,i) - chir(4,i) ) / (2.0*delx) bxi(i) = ( chii(3,i) - chii(4,i) ) / (2.0*delx) bzr(i) = ( chir(1,i) - chir(2,i) ) / (2.0*delz) @@ -883,6 +968,67 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) return end c----------------------------------------------------------------------- +c Nonlocal postprocessors cannot consume unavailable query fields. +c Raw exports retain the mask when both operations are disabled. +c----------------------------------------------------------------------- + subroutine check_cartesian_postprocessors(br,bz,bp,divzero, + $ chebyshev) + USE local_mod, ONLY: r8 + USE, INTRINSIC :: ieee_arithmetic, ONLY: ieee_is_finite + IMPLICIT NONE + complex(r8), intent(in) :: br(:,:),bz(:,:),bp(:,:) + logical, intent(in) :: divzero,chebyshev + logical :: finite + + finite = ALL(ieee_is_finite(REAL(br))) .AND. + $ ALL(ieee_is_finite(AIMAG(br))) .AND. + $ ALL(ieee_is_finite(REAL(bz))) .AND. + $ ALL(ieee_is_finite(AIMAG(bz))) .AND. + $ ALL(ieee_is_finite(REAL(bp))) .AND. + $ ALL(ieee_is_finite(AIMAG(bp))) + IF (finite) RETURN + IF (divzero) ERROR STOP 'Unresolved Cartesian total field: ' + $ //'disable divzero_flag for raw masked exports' + IF (chebyshev) ERROR STOP 'Unresolved Cartesian total field: ' + $ //'disable chebyshev_flag for raw masked exports' + end subroutine check_cartesian_postprocessors +c----------------------------------------------------------------------- +c Clearance from a closed discretized source contour. The margin +c excludes source crossings by the central finite-difference stencil. +c----------------------------------------------------------------------- + pure logical function pickup_query_clear(x,z,nsrc,xsrc,zsrc, + $ margin) result(clear) + USE local_mod, ONLY: r8 + USE, INTRINSIC :: ieee_arithmetic, ONLY: ieee_is_finite + IMPLICIT NONE + integer, intent(in) :: nsrc + real(r8), intent(in) :: x,z,xsrc(:),zsrc(:),margin + real(r8) :: dx,dz,len2,fraction,distance2,roundoff,limit + integer :: j,k + + clear = .false. + IF (.NOT. ieee_is_finite(x)) RETURN + IF (.NOT. ieee_is_finite(z)) RETURN + IF (nsrc < 3) RETURN + IF (SIZE(xsrc) < nsrc .OR. SIZE(zsrc) < nsrc) RETURN + roundoff = 64.0_r8*EPSILON(1.0_r8)*MAX(1.0_r8,ABS(x), + $ ABS(z),MAXVAL(ABS(xsrc(:nsrc))),MAXVAL(ABS(zsrc(:nsrc)))) + limit = (margin+roundoff)**2 + DO j = 1, nsrc + k = MOD(j,nsrc)+1 + dx = xsrc(k)-xsrc(j) + dz = zsrc(k)-zsrc(j) + len2 = dx**2+dz**2 + fraction = 0.0_r8 + IF (len2 > 0.0_r8) fraction = MAX(0.0_r8,MIN(1.0_r8, + $ ((x-xsrc(j))*dx+(z-zsrc(j))*dz)/len2)) + distance2 = (x-xsrc(j)-fraction*dx)**2+ + $ (z-zsrc(j)-fraction*dz)**2 + IF (distance2 <= limit) RETURN + ENDDO + clear = .true. + end function pickup_query_clear +c----------------------------------------------------------------------- c subprogram 9. loops. c grid loops. c----------------------------------------------------------------------- @@ -967,7 +1113,7 @@ subroutine loops c declarations. c----------------------------------------------------------------------- subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, - $ chir,chii,nsew,blr,bli,rgdl) + $ chir,chii,nsew,blr,bli,rgdl,query_valid) USE vglobal_mod implicit real(r8) (a-h,o-z) implicit integer (i-n) @@ -975,6 +1121,7 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, DIMENSION blr(*),bli(*),xsce(*),zsce(*),xscp(*),zscp(*) DIMENSION creal(nths,nfm), cimag(nths,nfm) DIMENSION chir(5,ndimlp), chii(5,ndimlp), rgdl(ndimlp) + logical, intent(in), optional :: query_valid(:) REAL nq factpi = twopi @@ -989,6 +1136,10 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, do io = 1, nobs + IF (PRESENT(query_valid)) THEN + IF (.NOT. query_valid(io)) CYCLE + ENDIF + xs = xobp(io) zs = zobp(io) @@ -1031,4 +1182,4 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, return end - end module vacuum_mod \ No newline at end of file + end module vacuum_mod From c394973462b16166ad0496c0f0209c08f617d580 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 8 Oct 2026 07:06:41 +0200 Subject: [PATCH 2/3] VACUUM: keep new pickup diagnostics within fixed-form columns --- vacuum/vacuum_ma.f | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/vacuum/vacuum_ma.f b/vacuum/vacuum_ma.f index 3484b8cc7..aa9884f54 100644 --- a/vacuum/vacuum_ma.f +++ b/vacuum/vacuum_ma.f @@ -714,7 +714,8 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, endif IF (keep_query) THEN -c Bypass both legacy surface snapping and REAL(4) grid generation. +c Bypass both legacy surface snapping and REAL(4) grid +c generation. nloop = (lx+1)*(lz+1) DO i = 0, lx DO j = 0, lz @@ -741,7 +742,8 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, delz = plrad * deloop IF (keep_query) THEN IF (delx == 0.0_r8 .OR. delz == 0.0_r8) - $ ERROR STOP 'Exact pickup requires a nonzero derivative step' + $ ERROR STOP 'Exact pickup requires a nonzero ' + $ //'derivative step' ENDIF clearance = 2.0_r8*MAX(ABS(delx),ABS(delz)) @@ -994,7 +996,8 @@ subroutine check_cartesian_postprocessors(br,bz,bp,divzero, end subroutine check_cartesian_postprocessors c----------------------------------------------------------------------- c Clearance from a closed discretized source contour. The margin -c excludes source crossings by the central finite-difference stencil. +c excludes source crossings by the central finite-difference +c stencil. c----------------------------------------------------------------------- pure logical function pickup_query_clear(x,z,nsrc,xsrc,zsrc, $ margin) result(clear) From e338152093fe0a7bc1cccf0cec3598b5c5b4927d Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 8 Oct 2026 08:43:58 +0200 Subject: [PATCH 3/3] VACUUM: retain double precision in observer-grid construction --- gpec/gpout.f | 19 +-- tests/.gitignore | 1 - tests/CMakeLists.txt | 79 ----------- tests/README.md | 43 ------ tests/expect_rejection.py | 13 -- tests/fixtures/vac.in | 82 ----------- tests/generate_query_fixture.py | 36 ----- tests/test_cartesian_postprocessors.f90 | 50 ------- tests/test_vacuum_query.f90 | 138 ------------------ vacuum/vacuum_ma.f | 180 ++---------------------- 10 files changed, 17 insertions(+), 624 deletions(-) delete mode 100644 tests/.gitignore delete mode 100644 tests/CMakeLists.txt delete mode 100644 tests/README.md delete mode 100644 tests/expect_rejection.py delete mode 100644 tests/fixtures/vac.in delete mode 100644 tests/generate_query_fixture.py delete mode 100644 tests/test_cartesian_postprocessors.f90 delete mode 100644 tests/test_vacuum_query.f90 diff --git a/gpec/gpout.f b/gpec/gpout.f index 19a36e5d2..fa5ccf325 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -4320,7 +4320,6 @@ END SUBROUTINE gpout_xbtangent c write perturbed rzphi components on rzphi grid. c----------------------------------------------------------------------- SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) - USE vacuum_mod, ONLY: check_cartesian_postprocessors c----------------------------------------------------------------------- c declaration. c----------------------------------------------------------------------- @@ -4509,9 +4508,7 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) $ "Computing external vacuum fields by surface currents" CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vbr,vbz,vbp, - $ preserve_query_coordinates=.TRUE., - $ query_r=gdr,query_z=gdz) + $ nr,nz,vgdl,vgdr,vgdz,vbr,vbz,vbp) IF (helicity<0) THEN vbr=CONJG(vbr) vbz=CONJG(vbz) @@ -4537,9 +4534,7 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) bnomn=bnomn-bnimn CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp, - $ preserve_query_coordinates=.TRUE., - $ query_r=gdr,query_z=gdz) + $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp) IF (helicity<0) THEN vpbr=CONJG(vpbr) @@ -4583,8 +4578,6 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) ENDDO ENDDO - CALL check_cartesian_postprocessors(btr,btz,btp, - $ divzero_flag,chebyshev_flag) IF (divzero_flag) THEN CALL gpeq_rzpdiv(nr,nz,gdr,gdz,btr,btz,btp) ENDIF @@ -4749,9 +4742,7 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) bnomn=bnomn-bnimn CALL gpvacuum_bnormal(psilim,bnomn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp, - $ preserve_query_coordinates=.TRUE., - $ query_r=gdr,query_z=gdz) + $ nr,nz,vgdl,vgdr,vgdz,vpbr,vpbz,vpbp) IF (helicity<0) THEN vpbr=CONJG(vpbr) vpbz=CONJG(vpbz) @@ -4764,9 +4755,7 @@ SUBROUTINE gpout_xbrzphi(egnum,xspmn,nr,nz,bnimn,bnomn) $ "Computing vacuum fields by external surface currents" CALL gpvacuum_bnormal(psilim,bnimn,nr,nz) CALL mscfld(wv,mpert,mthsurf,mthsurf,complex_flag, - $ nr,nz,vgdl,vgdr,vgdz,vvbr,vvbz,vvbp, - $ preserve_query_coordinates=.TRUE., - $ query_r=gdr,query_z=gdz) + $ nr,nz,vgdl,vgdr,vgdz,vvbr,vvbz,vvbp) IF (helicity<0) THEN vvbr=CONJG(vvbr) vvbz=CONJG(vvbz) diff --git a/tests/.gitignore b/tests/.gitignore deleted file mode 100644 index 84c048a73..000000000 --- a/tests/.gitignore +++ /dev/null @@ -1 +0,0 @@ -/build/ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt deleted file mode 100644 index 422b3f138..000000000 --- a/tests/CMakeLists.txt +++ /dev/null @@ -1,79 +0,0 @@ -cmake_minimum_required(VERSION 3.20) -project(gpec_cartesian_query_tests LANGUAGES C Fortran) -set(CMAKE_BUILD_TYPE Release CACHE STRING "Build configuration") -get_filename_component(GPEC_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/.." ABSOLUTE) -if(NOT CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") - message(FATAL_ERROR "This standalone native test adapter supports GNU Fortran") -endif() -find_package(BLAS REQUIRED) -find_package(LAPACK REQUIRED) -find_package(Python3 REQUIRED COMPONENTS Interpreter) -find_path(NETCDF_INCLUDE netcdf.mod REQUIRED) -find_library(NETCDFF_LIBRARY netcdff REQUIRED) -find_library(NETCDF_LIBRARY netcdf REQUIRED) -file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/version.inc" - " CHARACTER(len=*), PARAMETER :: version='cartesian-query-tests'\n") - -function(gpec_library name) - file(GLOB src "${GPEC_SOURCE}/${name}/*.f" "${GPEC_SOURCE}/${name}/*.F" - "${GPEC_SOURCE}/${name}/*.f90" "${GPEC_SOURCE}/${name}/*.F90" - "${GPEC_SOURCE}/${name}/*.c") - list(FILTER src EXCLUDE REGEX "/(pentrc\\.F90|slayer\\.f)$") - add_library(g_${name} STATIC ${src}) - set_target_properties(g_${name} PROPERTIES - Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/modules/${name}") - target_include_directories(g_${name} PUBLIC - "${CMAKE_CURRENT_BINARY_DIR}/modules/${name}") - target_include_directories(g_${name} PRIVATE "${NETCDF_INCLUDE}" - "${CMAKE_CURRENT_BINARY_DIR}") - target_compile_options(g_${name} PRIVATE - $<$:-ffixed-line-length-none;-fallow-argument-mismatch;-frecursive;-fopenmp>) - target_link_libraries(g_${name} PUBLIC ${ARGN} "${NETCDFF_LIBRARY}" - "${NETCDF_LIBRARY}" ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) - target_link_options(g_${name} PUBLIC -fopenmp) -endfunction() -gpec_library(equil) -gpec_library(lsode g_equil) -gpec_library(vacuum g_equil g_lsode) -gpec_library(pentrc g_equil g_lsode) -gpec_library(coil g_equil) -gpec_library(slayer g_equil g_lsode) -foreach(code dcon gpec) - file(GLOB src "${GPEC_SOURCE}/${code}/*.f" "${GPEC_SOURCE}/${code}/*.F") - add_executable(${code} ${src}) - set_target_properties(${code} PROPERTIES - Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/modules/${code}") - target_include_directories(${code} PRIVATE "${NETCDF_INCLUDE}" - "${CMAKE_CURRENT_BINARY_DIR}" "${CMAKE_CURRENT_BINARY_DIR}/modules/${code}") - target_compile_options(${code} PRIVATE -ffixed-line-length-none - -fallow-argument-mismatch -frecursive -fopenmp) - target_link_libraries(${code} PRIVATE g_pentrc g_vacuum g_slayer g_coil) -endforeach() - -set(QUERY_CASE "${CMAKE_CURRENT_BINARY_DIR}/query_case") -file(MAKE_DIRECTORY "${QUERY_CASE}") -configure_file(fixtures/vac.in "${QUERY_CASE}/vac.in" COPYONLY) -execute_process(COMMAND "${Python3_EXECUTABLE}" - "${CMAKE_CURRENT_SOURCE_DIR}/generate_query_fixture.py" "${QUERY_CASE}" - RESULT_VARIABLE FIXTURE_STATUS) -if(NOT FIXTURE_STATUS EQUAL 0) - message(FATAL_ERROR "Synthetic native query fixture generation failed") -endif() -enable_testing() -add_executable(test_vacuum_query test_vacuum_query.f90) -target_link_libraries(test_vacuum_query PRIVATE g_vacuum) -target_compile_options(test_vacuum_query PRIVATE -ffree-line-length-none -fcheck=all) -add_test(NAME vacuum_query_exact COMMAND test_vacuum_query) -set_tests_properties(vacuum_query_exact PROPERTIES - WORKING_DIRECTORY "${QUERY_CASE}" TIMEOUT 90) -add_executable(test_cartesian_postprocessors test_cartesian_postprocessors.f90) -target_link_libraries(test_cartesian_postprocessors PRIVATE g_vacuum) -add_test(NAME cartesian_postprocessors_finite_and_raw COMMAND test_cartesian_postprocessors) -foreach(mode chebyshev divzero) - add_test(NAME cartesian_postprocessors_reject_${mode} - COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/expect_rejection.py" - "$" "${mode}") -endforeach() -set_tests_properties(cartesian_postprocessors_finite_and_raw - cartesian_postprocessors_reject_chebyshev cartesian_postprocessors_reject_divzero - PROPERTIES TIMEOUT 15) diff --git a/tests/README.md b/tests/README.md deleted file mode 100644 index 67097c782..000000000 --- a/tests/README.md +++ /dev/null @@ -1,43 +0,0 @@ -# Cartesian vacuum query tests - -This standalone adapter builds the checkout's native libraries and the GPEC and -DCON executables, then runs the Cartesian export regression tests. It uses no -external case or private path. The upstream installation makefiles remain -unchanged. - -Requirements are GNU C/Fortran, CMake 3.20 or newer, BLAS/LAPACK, NetCDF C and -Fortran development libraries, and Python 3 (standard library only). CMake's -normal prefix/library/include search paths select those dependencies. Run from -this directory with the local Fortran workflow: - -```sh -FO_JOBS=2 OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 fo -``` - -For installations without `fo`, the standalone adapter can also be run using -the standard CMake/CTest commands: - -```sh -cmake -S tests -B tests/build -DCMAKE_BUILD_TYPE=Release -cmake --build tests/build --parallel 2 -OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 ctest --test-dir tests/build --output-on-failure -``` - -Run the CMake commands from the checkout root. All generated inputs and native -diagnostics stay in the chosen build directory. The synthetic fixture is a -circular clockwise surface with 128 unique nodes, five poloidal modes, signed complex -normal-field coefficients and n=3. It is a bounded numerical source, not a -device equilibrium or a physical forcing convergence study. - -`vacuum_query_exact` independently evaluates native CHI at double query -coordinates, then checks pickup and MSCFLD forwarding, source/stencil masks, -and identical absent/false legacy optional behavior. It includes a resolved -query inside the old snapping band and unresolved source/crossing queries. -The independent reference bypasses both coordinate construction and pickup. -The remaining tests permit finite fields and raw masked exports, and require -the specific native rejection diagnostics for each optional postprocessor. -They do not treat an arbitrary crash as a successful rejection. - -These tests establish the query/output contract. They do not establish source -quadrature convergence, accuracy near a singular sheet, arbitrary malformed -derivative settings, global plasma response or physical forcing validity. diff --git a/tests/expect_rejection.py b/tests/expect_rejection.py deleted file mode 100644 index bfd712f8f..000000000 --- a/tests/expect_rejection.py +++ /dev/null @@ -1,13 +0,0 @@ -"""Require the native guard's specific diagnostic, rather than any failure.""" -import subprocess -import sys - -executable, mode = sys.argv[1:] -assert mode in ("chebyshev", "divzero") -result = subprocess.run([executable, mode], capture_output=True, text=True, timeout=10) -output = result.stdout + result.stderr -expected = f"disable {mode}_flag for raw masked exports" -if result.returncode == 0 or expected not in output: - print(output) - raise SystemExit(f"Expected native {mode} rejection was not observed") -print(f"PASS native {mode} rejection with specific diagnostic") diff --git a/tests/fixtures/vac.in b/tests/fixtures/vac.in deleted file mode 100644 index 7cc8612d0..000000000 --- a/tests/fixtures/vac.in +++ /dev/null @@ -1,82 +0,0 @@ -&MODES - mth=128 - xiin(1:9)=0 0 0 0 0 0 0 1 0 - lsymz=.true. - leqarcw=1 - lzio=0 - lgato=0 - lrgato=0 -/ -&DEBUGS - checkd=.false. - check1=.false. - check2=.false. - checke=.false. - checks=.false. - wall=.false. - lkplt=0 -/ -&VACDAT - ishape=6 - nsing=128 - epsq=1e-5 - noutv=37 - idgt=6 - idot=0 - idsk=0 - delg=15.01 - delfac=0.001 - cn0=1 -/ -&SHAPE - ipshp=0 - xpl=100 - apl=1 - a=0.2413 - b=170 - bpl=1 - dpl=0 - r=1 - abulg=0.932 - bbulg=17 - tbulg=0.02 - qain=2.5 -/ -&DIAGNS - lkdis=.false. - ieig=0 - iloop=0 - lpsub=1 - nloop=128 - nloopr=0 - nphil=3 - nphse=1 - xofsl=0 - ntloop=32 - aloop=0.01 - bloop=1.6 - dloop=0.5 - rloop=1 - deloop=0.001 - mx=21 - mz=21 - nph=0 - nxlpin=2 - nzlpin=2 - epslp=0.02 - xlpmin=0.7 - xlpmax=2.7 - zlpmin=-1.5 - zlpmax=1.5 - linterior=2 -/ -&SPRK - nminus=0 - nplus=0 - mphi=16 - lwrt11=0 - civ=0 - lff=0 - ff=1.6 - fv=1.6 -/ diff --git a/tests/generate_query_fixture.py b/tests/generate_query_fixture.py deleted file mode 100644 index 7d96a7d86..000000000 --- a/tests/generate_query_fixture.py +++ /dev/null @@ -1,36 +0,0 @@ -"""Write a bounded synthetic circular source in native VACIN5 format.""" -import math -from pathlib import Path -import sys - -out = Path(sys.argv[1]) -out.mkdir(parents=True, exist_ok=True) -nodes, low, high, n = 128, -2, 2, 3 -center, radius, zcenter = 6.2, 0.62, 0.0 -theta = [i / nodes for i in range(nodes + 1)] -R = [center + radius * math.cos(2 * math.pi * t) for t in theta[:-1]] -# The native vacuum contour convention is clockwise in the (R,Z) plane. -Z = [zcenter - radius * math.sin(2 * math.pi * t) for t in theta[:-1]] -R.append(R[0]) -Z.append(Z[0]) -coefficients = {-1: 0.04 - 0.03j, 0: 0.02 + 0.07j, 2: -0.015 + 0.01j} -lines = ["scalars", "", "2", "2", f"{center-radius:.17e}", - f"{center+radius:.17e}", f"{zcenter-radius:.17e}", - f"{zcenter+radius:.17e}", str(nodes), str(low), str(high), str(n), "2.5"] - - -def block(label, values): - lines.extend([label, ""]) - for i in range(0, len(values), 4): - lines.append(" ".join(f"{v:.17e}" for v in values[i:i+4])) - - -block("Poloidal Coordinate Theta", theta) -block("Normalized Arclength", theta) -block("Major Radius", R) -block("Vertical Coordinate", Z) -block("Toroidal Angle Shift", [0.0] * (nodes + 1)) -block("Normal Field Real", [coefficients.get(m, 0j).real for m in range(low, high + 1)]) -block("Normal Field Imaginary", [coefficients.get(m, 0j).imag for m in range(low, high + 1)]) -(out / "vacin5").write_text("\n".join(lines) + "\n") -(out / "query_contract.in").write_text(f"{nodes} {high-low+1}\n") diff --git a/tests/test_cartesian_postprocessors.f90 b/tests/test_cartesian_postprocessors.f90 deleted file mode 100644 index c17152ca7..000000000 --- a/tests/test_cartesian_postprocessors.f90 +++ /dev/null @@ -1,50 +0,0 @@ -! Exercise the production guard with manufactured finite and unavailable fields. -program test_cartesian_postprocessors - use, intrinsic :: iso_fortran_env, only: dp => real64 - use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan - use vacuum_mod, only: check_cartesian_postprocessors - implicit none - complex(dp) :: br(3,4), bz(3,4), bp(3,4) - real(dp) :: unavailable - character(24) :: mode - integer :: component - - call get_command_argument(1, mode) - unavailable = ieee_value(0.0_dp, ieee_quiet_nan) - br = cmplx(1.0_dp, 2.0_dp, dp) - bz = cmplx(-3.0_dp, 4.0_dp, dp) - bp = cmplx(5.0_dp, -6.0_dp, dp) - call check_cartesian_postprocessors(br, bz, bp, .true., .true.) - select case (trim(mode)) - case ('chebyshev') - bz(2,3) = cmplx(-3.0_dp, unavailable, dp) - call check_cartesian_postprocessors(br, bz, bp, .false., .true.) - print *, 'UNEXPECTED: Chebyshev accepted an unavailable total field' - case ('divzero') - bp(2,3) = cmplx(unavailable, -6.0_dp, dp) - call check_cartesian_postprocessors(br, bz, bp, .true., .false.) - print *, 'UNEXPECTED: divzero accepted an unavailable total field' - case default - do component = 1, 6 - br = cmplx(1.0_dp, 2.0_dp, dp) - bz = cmplx(-3.0_dp, 4.0_dp, dp) - bp = cmplx(5.0_dp, -6.0_dp, dp) - select case (component) - case (1) - br(2,3) = cmplx(unavailable, 2.0_dp, dp) - case (2) - br(2,3) = cmplx(1.0_dp, unavailable, dp) - case (3) - bz(2,3) = cmplx(unavailable, 4.0_dp, dp) - case (4) - bz(2,3) = cmplx(-3.0_dp, unavailable, dp) - case (5) - bp(2,3) = cmplx(unavailable, -6.0_dp, dp) - case (6) - bp(2,3) = cmplx(5.0_dp, unavailable, dp) - end select - call check_cartesian_postprocessors(br, bz, bp, .false., .false.) - end do - print *, 'PASS finite postprocessor inputs and raw masked exports' - end select -end program test_cartesian_postprocessors diff --git a/tests/test_vacuum_query.f90 b/tests/test_vacuum_query.f90 deleted file mode 100644 index e38e81646..000000000 --- a/tests/test_vacuum_query.f90 +++ /dev/null @@ -1,138 +0,0 @@ -! Native query oracle; supply query_contract.in, vac.in and vacin5 externally. -! The reference calls CHI directly at doubles, without pickup or loops. -program test_vacuum_query - use, intrinsic :: iso_fortran_env, only: dp => real64 - use, intrinsic :: ieee_arithmetic, only: ieee_is_nan, ieee_is_finite - use vacuum_mod, only: defglo, ent33, make_bltobp, diaplt, chi, pickup, mscfld - use vglobal_mod, only: global_alloc, global_dealloc, kernelsign, ntsin0, & - nths0, nfm, mtot, ndimlp, xobp, zobp, xloop, zloop, ldcon, lgpec, ieig, & - iotty, outpest, inmode, outmod, farwal, mth, lmax, lmin, vacmat, vacmtiu, & - xpla, zpla, plrad, deloop, epslp, grri, nths, nloop, nloopr, xplap, & - zplap, bnlr, bnli, n, nxlpin, nzlpin - implicit none - integer :: surface_nodes, modes, iu, peak, side, l, k - real(dp) :: qr(0:1,0:1), qz(0:1,0:1), rr(0:1,0:1), zz(0:1,0:1) - real(dp) :: lr(0:1,0:1), lz(0:1,0:1), delta, band, center_r, center_z - integer :: flags(0:1,0:1), legacy_flags(0:1,0:1) - complex(dp) :: br(0:1,0:1), bz(0:1,0:1), bp(0:1,0:1) - complex(dp) :: legacy_br(0:1,0:1), legacy_bz(0:1,0:1), legacy_bp(0:1,0:1) - complex(dp) :: reference(3,2) - complex(dp), allocatable :: matrix_ref(:,:), matrix_out(:,:) - real(dp), allocatable :: re(:,:), im(:,:), region(:), cc(:,:), ss(:,:) - real(dp) :: field_error, field_scale - - open(newunit=iu,file='query_contract.in',status='old') - read(iu,*) surface_nodes,modes - close(iu) - if(surface_nodes<16 .or. surface_nodes>512) error stop 'Bounded surface size required' - if(modes<1 .or. modes>40) error stop 'Bounded mode count required' - kernelsign=1;ntsin0=surface_nodes+1;nths0=surface_nodes;nfm=modes;mtot=modes - call global_alloc(nths0,nfm,mtot,ntsin0) - ndimlp=4 - allocate(xobp(4),zobp(4),xloop(4),zloop(4)) - call defglo(surface_nodes) - farwal=.true.;ldcon=0;lgpec=1;ieig=0 - open(iotty,file='mscvac.out',status='replace') - open(outpest,file='pestotv',status='replace',form='formatted') - open(inmode,file='vac.in',status='old',form='formatted') - open(outmod,file='modovmc',status='replace',form='formatted') - call ent33 - if(mth/=surface_nodes .or. lmax(1)-lmin(1)+1/=modes) error stop 'Source basis mismatch' - call make_bltobp - call diaplt - allocate(matrix_ref(modes,modes),matrix_out(modes,modes)) - matrix_ref=cmplx(vacmat,vacmtiu,dp) - open(newunit=iu,file='query_native_matrix.dat',status='replace') - do k=1,modes - do l=1,modes - write(iu,'(2i6,2es26.17)') l,k,real(matrix_ref(l,k)),aimag(matrix_ref(l,k)) - end do - end do - close(iu) - plrad=.5_dp*(maxval(xpla(:mth))-minval(xpla(:mth))) - delta=plrad*deloop;band=plrad*epslp - if(delta<=0 .or. band<=8*delta) error stop 'Oracle requires a resolvable legacy snap band' - peak=maxloc(xpla(:mth),dim=1) - center_r=.5_dp*(maxval(xpla(:mth))+minval(xpla(:mth))) - center_z=.5_dp*(maxval(zpla(:mth))+minval(zpla(:mth))) - qr(0,0)=xpla(peak)-.5_dp*band;qz(0,0)=zpla(peak) - qr(0,1)=center_r+.13_dp*plrad;qz(0,1)=center_z+.07_dp*plrad - qr(1,0)=xpla(peak);qz(1,0)=zpla(peak) - qr(1,1)=xpla(peak)-.5_dp*delta;qz(1,1)=zpla(peak) - - ! The reference independently evaluates only the two resolved interior queries. - allocate(re(5,4),im(5,4),region(4),cc(nths,modes),ss(nths,modes)) - cc=0;ss=0 - do l=1,modes - cc(:mth,l)=grri(:mth,l) - ss(:mth,l)=grri(:mth,modes+l) - end do - nloop=2;nloopr=0 - xloop(1:2)=qr(0,:);zloop(1:2)=qz(0,:) - region=-1;re=0;im=0 - do side=1,5 - xobp(1:2)=xloop(1:2);zobp(1:2)=zloop(1:2) - select case(side) - case(1) - zobp(1:2)=zobp(1:2)+delta - case(2) - zobp(1:2)=zobp(1:2)-delta - case(3) - xobp(1:2)=xobp(1:2)+delta - case(4) - xobp(1:2)=xobp(1:2)-delta - end select - call chi(xpla,zpla,xplap,zplap,-1,cc,ss,mth,1,re,im,side,bnlr,bnli,region) - end do - do k=1,2 - reference(1,k)=cmplx(re(3,k)-re(4,k),im(3,k)-im(4,k),dp)/(2*delta) - reference(2,k)=cmplx(re(1,k)-re(2,k),im(1,k)-im(2,k),dp)/(2*delta) - reference(3,k)=cmplx(n*im(5,k),-n*re(5,k),dp)/xloop(k) - end do - nxlpin=2;nzlpin=2 - call pickup(bnlr,bnli,1,1,flags,rr,zz,br,bz,bp, & - preserve_query_coordinates=.true.,query_r=qr,query_z=qz) - call check_exact_query() - close(iotty);close(outpest);close(inmode);close(outmod) - call global_dealloc - call cleanup - - ! Exercise actual MSCFLD forwarding; its legacy GPEC path does not write WV. - farwal=.true. - call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & - 1,1,flags,rr,zz,br,bz,bp,preserve_query_coordinates=.true.,query_r=qr,query_z=qz) - call check_exact_query() - - farwal=.true. - call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & - 1,1,legacy_flags,lr,lz,legacy_br,legacy_bz,legacy_bp) - farwal=.true. - call mscfld(matrix_out,modes,surface_nodes,surface_nodes,.true., & - 1,1,flags,rr,zz,br,bz,bp,preserve_query_coordinates=.false.) - if(any(flags/=legacy_flags) .or. any(rr/=lr) .or. any(zz/=lz)) & - error stop 'Legacy optional contract changed coordinates or flags' - if(any(br/=legacy_br) .or. any(bz/=legacy_bz) .or. any(bp/=legacy_bp)) & - error stop 'Legacy optional contract changed field values' - print *, 'PASS exact native query coordinates, direct CHI fields, masks and legacy contract' - print *, 'DIRECT_CHI_RELATIVE_FIELD_ERROR',field_error -contains - subroutine check_exact_query() - integer :: j - if(any(rr/=qr) .or. any(zz/=qz)) error stop 'Query coordinates relocated or downcast' - if(any(flags(0,:)/=1)) error stop 'Resolved interior query misclassified' - if(any(flags(1,:)/=-2)) error stop 'Source or crossing-stencil query not masked' - if(.not.all(ieee_is_nan(real(br(1,:))))) error stop 'Invalid BR not unavailable' - if(.not.all(ieee_is_nan(aimag(br(1,:))))) error stop 'Invalid imaginary BR not unavailable' - if(.not.all(ieee_is_nan(real(bz(1,:))))) error stop 'Invalid BZ not unavailable' - if(.not.all(ieee_is_nan(aimag(bz(1,:))))) error stop 'Invalid imaginary BZ not unavailable' - if(.not.all(ieee_is_nan(real(bp(1,:))))) error stop 'Invalid Bphi not unavailable' - if(.not.all(ieee_is_nan(aimag(bp(1,:))))) error stop 'Invalid imaginary Bphi not unavailable' - field_error=0;field_scale=max(1e-20_dp,maxval(abs(reference))) - do j=0,1 - field_error=max(field_error, & - maxval(abs([br(0,j),bz(0,j),bp(0,j)]-reference(:,j+1)))/field_scale) - end do - if(.not.ieee_is_finite(field_error)) error stop 'Nonfinite resolved field' - if(field_error>2e-11_dp) error stop 'Pickup field differs from independent direct CHI query' - end subroutine check_exact_query -end program test_vacuum_query diff --git a/vacuum/vacuum_ma.f b/vacuum/vacuum_ma.f index aa9884f54..97014f3c9 100644 --- a/vacuum/vacuum_ma.f +++ b/vacuum/vacuum_ma.f @@ -170,8 +170,7 @@ subroutine mscvac(wv,mpert,mtheta,mthvac,complex_flag, c declarations. c----------------------------------------------------------------------- subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, - $ lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp,op_ahgfile, - $ preserve_query_coordinates,query_r,query_z) + $ lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp,op_ahgfile) USE vglobal_mod implicit real(r8) (a-h,o-z) implicit integer (i-n) @@ -187,8 +186,6 @@ subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, complex(r8), parameter :: ifac=(0,1) dimension xi(nfm), xii(nfm), xilnq(nfm), xiilnq(nfm) character(128), intent(in), optional :: op_ahgfile - logical, intent(in), optional :: preserve_query_coordinates - real(r8), intent(in), optional :: query_r(0:,0:),query_z(0:,0:) if (present(op_ahgfile)) then ahgfile = trim(op_ahgfile) @@ -290,8 +287,7 @@ subroutine mscfld(wv,mpert,mtheta,mthvac,complex_flag, 99 continue call make_bltobp call diaplt - call pickup(bnlr,bnli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, - $ preserve_query_coordinates,query_r,query_z) + call pickup(bnlr,bnli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) c----------------------------------------------------------------------- c termination. c----------------------------------------------------------------------- @@ -625,21 +621,14 @@ subroutine diaplt c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- - subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, - $ preserve_query_coordinates,query_r,query_z) + subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp) USE vglobal_mod - USE, INTRINSIC :: ieee_arithmetic, ONLY: - $ ieee_value,ieee_quiet_nan implicit real(r8) (a-h,o-z) implicit integer (i-n) INTEGER, DIMENSION(0:lx,0:lz) :: vgdl REAL(r8), DIMENSION(0:lx,0:lz) :: vgdx,vgdz COMPLEX(r8), DIMENSION(0:lx,0:lz) :: vbx,vbz,vbp - logical, intent(in), optional :: preserve_query_coordinates - real(r8), intent(in), optional :: query_r(0:,0:),query_z(0:,0:) - logical :: keep_query,query_valid(ndimlp) - real(r8) :: clearance,unavailable DIMENSION blr(*), bli(*) ! CHARACTER(130), DIMENSION(10) :: string @@ -654,28 +643,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, $ btr(ndimlp),bti(ndimlp), bpr(ndimlp),bpi(ndimlp), $ bphir(ndimlp),bphii(ndimlp) - keep_query = .false. - IF (PRESENT(preserve_query_coordinates)) - $ keep_query = preserve_query_coordinates - IF (keep_query) THEN - IF (.NOT. PRESENT(query_r)) - $ ERROR STOP 'Exact pickup requires query_r' - IF (.NOT. PRESENT(query_z)) - $ ERROR STOP 'Exact pickup requires query_z' - IF (SIZE(query_r,1) /= lx+1 .OR. - $ SIZE(query_r,2) /= lz+1) - $ ERROR STOP 'Exact pickup query_r shape mismatch' - IF (SIZE(query_z,1) /= lx+1 .OR. - $ SIZE(query_z,2) /= lz+1) - $ ERROR STOP 'Exact pickup query_z shape mismatch' - IF (nxlpin /= lx+1 .OR. nzlpin /= lz+1) - $ ERROR STOP 'Exact pickup grid dimensions mismatch' - IF (nloopr /= 0) - $ ERROR STOP 'Exact pickup supports the Cartesian grid only' - ENDIF - query_valid = .true. - unavailable = ieee_value(0.0_r8,ieee_quiet_nan) - ndlp = mth / ntloop jmax1 = lmax(1) - lmin(1) + 1 @@ -713,20 +680,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, zmax = max(zmxp,zmxw) endif - IF (keep_query) THEN -c Bypass both legacy surface snapping and REAL(4) grid -c generation. - nloop = (lx+1)*(lz+1) - DO i = 0, lx - DO j = 0, lz - indxl = i*(lz+1)+j+1 - xloop(indxl) = query_r(i,j) - zloop(indxl) = query_z(i,j) - ENDDO - ENDDO - ELSE - CALL loops - ENDIF + CALL loops nobs = nloop + 3*nloopr @@ -740,12 +694,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, ns = mth delx = plrad * deloop delz = plrad * deloop - IF (keep_query) THEN - IF (delx == 0.0_r8 .OR. delz == 0.0_r8) - $ ERROR STOP 'Exact pickup requires a nonzero ' - $ //'derivative step' - ENDIF - clearance = 2.0_r8*MAX(ABS(delx),ABS(delz)) igdl = 8 isgchi = -1 @@ -755,23 +703,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, xloops(i) = xloop(i) zloops(i) = zloop(i) - IF (keep_query) THEN -c -2: no resolved same-side finite-difference field here. - query_valid(i) = pickup_query_clear(xloop(i),zloop(i), - $ mth,xpla,zpla,clearance) - IF (xloop(i) <= 0.0_r8) query_valid(i) = .false. - IF (.NOT. farwal) THEN - IF (query_valid(i)) query_valid(i) = - $ pickup_query_clear(xloop(i),zloop(i), - $ mw,xwal,zwal,clearance) - ENDIF - IF (.NOT. query_valid(i)) THEN - igdl(i) = -2 - rgdl(i) = 1.0_r8 - GO TO 239 - ENDIF - ENDIF - fintjj = 0.0 DO jj = 1, mth dxjj = xloop(i) - xinf(jj) @@ -800,16 +731,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, 238 CONTINUE - IF (keep_query) THEN - IF (igdl(i) == 8) THEN - query_valid(i) = .false. - igdl(i) = -2 - rgdl(i) = 1.0_r8 - GO TO 239 - ENDIF - ENDIF - - IF (.NOT. keep_query) THEN DO j = 1, mth IF ( (xinf(j)-xloop(i))**2 + (zinf(j)-zloop(i))**2 $ < (epslp*plrad)**2 ) THEN ! points on surface for these. @@ -837,7 +758,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, bphii(i) = - n * zchipr(i) / xloops(i) END IF END DO - ENDIF 239 CONTINUE @@ -900,7 +820,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, isg = -1 call chi ( xpla,zpla,xplap,zplap,isg, chiwc,chiws, ns,1, - $ cwrkr,cwrki,nsew, blr,bli,rgdl,query_valid ) + $ cwrkr,cwrki,nsew, blr,bli,rgdl ) do i = 1, nobs chir(nsew,i) = cwrkr(nsew,i) @@ -923,7 +843,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, isg = 1 call chi ( xwal,zwal,xwalp,zwalp,isg,chiwc,chiws, ns,0, - $ cwrkr,cwrki,nsew, blr,bli,rwall,query_valid ) + $ cwrkr,cwrki,nsew, blr,bli,rwall ) do i = 1, nobs chir(nsew,i) = chir(nsew,i) + cwrkr(nsew,i) @@ -934,14 +854,7 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, END DO DO i = 1, nobs - IF (.NOT. query_valid(i)) THEN - bxr(i) = unavailable - bxi(i) = unavailable - bzr(i) = unavailable - bzi(i) = unavailable - bphir(i) = unavailable - bphii(i) = unavailable - ELSE IF (igdl(i) .NE. -1) THEN + IF (igdl(i) .NE. -1) THEN bxr(i) = ( chir(3,i) - chir(4,i) ) / (2.0*delx) bxi(i) = ( chii(3,i) - chii(4,i) ) / (2.0*delx) bzr(i) = ( chir(1,i) - chir(2,i) ) / (2.0*delz) @@ -970,68 +883,6 @@ subroutine pickup(blr,bli,lx,lz,vgdl,vgdx,vgdz,vbx,vbz,vbp, return end c----------------------------------------------------------------------- -c Nonlocal postprocessors cannot consume unavailable query fields. -c Raw exports retain the mask when both operations are disabled. -c----------------------------------------------------------------------- - subroutine check_cartesian_postprocessors(br,bz,bp,divzero, - $ chebyshev) - USE local_mod, ONLY: r8 - USE, INTRINSIC :: ieee_arithmetic, ONLY: ieee_is_finite - IMPLICIT NONE - complex(r8), intent(in) :: br(:,:),bz(:,:),bp(:,:) - logical, intent(in) :: divzero,chebyshev - logical :: finite - - finite = ALL(ieee_is_finite(REAL(br))) .AND. - $ ALL(ieee_is_finite(AIMAG(br))) .AND. - $ ALL(ieee_is_finite(REAL(bz))) .AND. - $ ALL(ieee_is_finite(AIMAG(bz))) .AND. - $ ALL(ieee_is_finite(REAL(bp))) .AND. - $ ALL(ieee_is_finite(AIMAG(bp))) - IF (finite) RETURN - IF (divzero) ERROR STOP 'Unresolved Cartesian total field: ' - $ //'disable divzero_flag for raw masked exports' - IF (chebyshev) ERROR STOP 'Unresolved Cartesian total field: ' - $ //'disable chebyshev_flag for raw masked exports' - end subroutine check_cartesian_postprocessors -c----------------------------------------------------------------------- -c Clearance from a closed discretized source contour. The margin -c excludes source crossings by the central finite-difference -c stencil. -c----------------------------------------------------------------------- - pure logical function pickup_query_clear(x,z,nsrc,xsrc,zsrc, - $ margin) result(clear) - USE local_mod, ONLY: r8 - USE, INTRINSIC :: ieee_arithmetic, ONLY: ieee_is_finite - IMPLICIT NONE - integer, intent(in) :: nsrc - real(r8), intent(in) :: x,z,xsrc(:),zsrc(:),margin - real(r8) :: dx,dz,len2,fraction,distance2,roundoff,limit - integer :: j,k - - clear = .false. - IF (.NOT. ieee_is_finite(x)) RETURN - IF (.NOT. ieee_is_finite(z)) RETURN - IF (nsrc < 3) RETURN - IF (SIZE(xsrc) < nsrc .OR. SIZE(zsrc) < nsrc) RETURN - roundoff = 64.0_r8*EPSILON(1.0_r8)*MAX(1.0_r8,ABS(x), - $ ABS(z),MAXVAL(ABS(xsrc(:nsrc))),MAXVAL(ABS(zsrc(:nsrc)))) - limit = (margin+roundoff)**2 - DO j = 1, nsrc - k = MOD(j,nsrc)+1 - dx = xsrc(k)-xsrc(j) - dz = zsrc(k)-zsrc(j) - len2 = dx**2+dz**2 - fraction = 0.0_r8 - IF (len2 > 0.0_r8) fraction = MAX(0.0_r8,MIN(1.0_r8, - $ ((x-xsrc(j))*dx+(z-zsrc(j))*dz)/len2)) - distance2 = (x-xsrc(j)-fraction*dx)**2+ - $ (z-zsrc(j)-fraction*dz)**2 - IF (distance2 <= limit) RETURN - ENDDO - clear = .true. - end function pickup_query_clear -c----------------------------------------------------------------------- c subprogram 9. loops. c grid loops. c----------------------------------------------------------------------- @@ -1043,12 +894,12 @@ subroutine loops implicit real(r8) (a-h,o-z) implicit integer (i-n) - REAL, DIMENSION(:,:), ALLOCATABLE :: xloopin, zloopin - REAL, DIMENSION(:), ALLOCATABLE :: sourcemat + REAL(r8), DIMENSION(:,:), ALLOCATABLE :: xloopin, zloopin + REAL(r8), DIMENSION(:), ALLOCATABLE :: sourcemat - dxlin = 1.0 / (nxlpin-1) - dzlin = 1.0 / (nzlpin-1) + dxlin = 1.0_r8 / (nxlpin-1) + dzlin = 1.0_r8 / (nzlpin-1) nxzlin = nxlpin * nzlpin nloop = nxzlin @@ -1116,7 +967,7 @@ subroutine loops c declarations. c----------------------------------------------------------------------- subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, - $ chir,chii,nsew,blr,bli,rgdl,query_valid) + $ chir,chii,nsew,blr,bli,rgdl) USE vglobal_mod implicit real(r8) (a-h,o-z) implicit integer (i-n) @@ -1124,7 +975,6 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, DIMENSION blr(*),bli(*),xsce(*),zsce(*),xscp(*),zscp(*) DIMENSION creal(nths,nfm), cimag(nths,nfm) DIMENSION chir(5,ndimlp), chii(5,ndimlp), rgdl(ndimlp) - logical, intent(in), optional :: query_valid(:) REAL nq factpi = twopi @@ -1139,10 +989,6 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, do io = 1, nobs - IF (PRESENT(query_valid)) THEN - IF (.NOT. query_valid(io)) CYCLE - ENDIF - xs = xobp(io) zs = zobp(io) @@ -1185,4 +1031,4 @@ subroutine chi(xsce,zsce,xscp,zscp,isg,creal,cimag,ns,ip, return end - end module vacuum_mod + end module vacuum_mod \ No newline at end of file