diff --git a/.gitignore b/.gitignore index e74606b7..f01d712c 100644 --- a/.gitignore +++ b/.gitignore @@ -61,3 +61,4 @@ cobertura.xml coverage_html/ coverage-badge.svg coverage-badge.json +threed*.default diff --git a/TODO.md b/TODO.md new file mode 100644 index 00000000..0770d8af --- /dev/null +++ b/TODO.md @@ -0,0 +1,105 @@ +# TODO: Booz_xform Integration for SIMPLE + +## Overview +Add support for reading Boozer coordinate magnetic fields directly from booz_xform output files as an alternative to the current internal VMEC-to-Boozer conversion. + +## Background +- Currently: SIMPLE reads VMEC files and internally converts to Boozer coordinates using `boozer_converter.F90` +- Goal: Support direct input of pre-computed Boozer fields from [booz_xform](https://github.com/hiddenSymmetries/booz_xform) +- Benefit: Potentially faster initialization and standardized Boozer transformation + +## Implementation Tasks + +### High Priority + +1. **Review booz_xform file format** ✓ + - [x] Study booz_xform NetCDF output structure + - [x] Document required fields and data layout + - [x] Identify mapping to SIMPLE's internal structures + +2. **Create field_booz_xform.f90 module** ✓ + - [x] Implement NetCDF reader for booz_xform files + - [x] Map booz_xform data to SIMPLE's field structures + - [x] Handle coordinate conventions and normalizations + +3. **Add booz_xform input option** ✓ + - [x] Add new field type to `params.f90` (e.g., `isw_field_type = 5`) + - [x] Add `booz_xform_file` parameter to namelist + - [x] Update field type selection logic + +4. **Update field initialization** ✓ + - [x] Modify `simple_main.f90` to handle booz_xform input + - [x] Ensure compatibility with existing Boozer field evaluation + - [x] Handle field normalization and units + +5. **Integration with existing Boozer evaluation** ✓ + - [x] Connect booz_xform reader to field evaluation through magfie interface + - [x] Implement evaluate method with VMEC→Boozer coordinate transformation + - [x] Set up proper field component calculations + +### Medium Priority + +6. **Validation and testing** + - [ ] Complete `test_booz_xform.f90` implementation + - [ ] Compare results: VMEC→internal Boozer vs VMEC→booz_xform→SIMPLE + - [ ] Benchmark performance differences + - [ ] Test with various stellarator configurations + +7. **Example workflow** + - [ ] Create example: VMEC wout.nc → booz_xform → SIMPLE + - [ ] Document command sequences + - [ ] Provide sample input files + +### Low Priority + +8. **Documentation** + - [ ] Update user manual with booz_xform option + - [ ] Add to CLAUDE.md for AI assistance + - [ ] Create tutorial for booz_xform workflow + +## Technical Considerations + +### Data Structures +- Booz_xform outputs: `|B|`, `G`, `I`, `K`, `p`, `q`, derivatives +- SIMPLE needs: Magnetic field components, metric tensor, Jacobian +- Coordinate mapping: (s, θ_B, ζ_B) in both systems + +### Compatibility +- Ensure backward compatibility with existing VMEC input +- Support switching between internal and external Boozer conversion +- Handle different normalization conventions + +### Performance +- Pre-computed Boozer coordinates may reduce initialization time +- Consider memory vs computation tradeoffs +- Profile both approaches for typical use cases + +## Current Status +### Completed ✓ +- Basic test file created: `test/tests/test_booz_xform.f90` +- CMake integration ready +- Existing Boozer infrastructure in place: `field_can_boozer.f90`, `boozer_converter.F90` +- Created test framework with setup script and comparison plots +- Implemented `field_booz_xform.f90` module to read BOOZXFORM files +- Fixed NetCDF interface issues (using nctools_module from libneo) +- Test passes and generates comparison plots +- Added BOOZXFORM=5 constant to magfie.f90 +- Updated simple_main.f90 to handle the new field type (sets boozxform_filename) +- Implemented proper evaluate method in field_booz_xform.f90 with coordinate transformation +- Added boozxform_file parameter to params.f90 namelist +- Created benchmark script: `benchmark/benchmark_boozxform.py` + +### In Progress +- Fixed NetCDF loading for scalar and 1D arrays +- 2D array loading still has issues ("Start+count exceeds dimension bound") +- Need to implement proper coordinate transformation (pmns/gmn interpretation) + +### Known Issues +- NetCDF 2D array reading fails for Fourier coefficients (rmnc_b, zmns_b, etc.) +- Coordinate transformation between VMEC and Boozer angles incomplete +- RKF45 integrator errors due to missing field data + +## Next Steps +1. Fix 2D array NetCDF loading (possibly transpose issue) +2. Implement proper VMEC↔Boozer coordinate transformation +3. Complete validation with benchmark files \ No newline at end of file diff --git a/benchmark/benchmark_boozxform.py b/benchmark/benchmark_boozxform.py new file mode 100755 index 00000000..9091a8b5 --- /dev/null +++ b/benchmark/benchmark_boozxform.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +""" +Benchmark script to compare orbit integration using VMEC vs BOOZXFORM fields +""" + +import numpy as np +import matplotlib.pyplot as plt +import subprocess +import os +import shutil +import time + +# Set paths +SIMPLE_EXEC = "../build/simple.x" +VMEC_FILE = "../../benchmark_orbit/booz_xform/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc" +BOOZ_FILE = "../../benchmark_orbit/booz_xform/boozmn_LandremanPaul2021_QA_reactorScale_lowres_reference.nc" + +# Create two input files - one for VMEC, one for BOOZXFORM +simple_in_vmec = """&simple + vmec_file = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' + trace_time = 1.0e-4 ! short integration time + dtau = 5.0e-7 + ntau = 10000 + ntestpart = 5 + ntimstep = 1000 + sbeg = 0.5 + phibeg = 0.0 + thetabeg = 0.0 + bmod00 = 5.0 + isw_field_type = 1 ! VMEC +/ +""" + +simple_in_booz = """&simple + vmec_file = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' + boozxform_file = 'boozmn_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' + trace_time = 1.0e-4 ! short integration time + dtau = 5.0e-7 + ntau = 10000 + ntestpart = 5 + ntimstep = 1000 + sbeg = 0.5 + phibeg = 0.0 + thetabeg = 0.0 + bmod00 = 5.0 + isw_field_type = 5 ! BOOZXFORM +/ +""" + +def run_benchmark(case_name, input_content): + """Run SIMPLE with given input and return results""" + print(f"\nRunning {case_name} case...") + + # Create directory for this case + case_dir = f"benchmark_{case_name}" + os.makedirs(case_dir, exist_ok=True) + os.chdir(case_dir) + + # Copy necessary files + shutil.copy2(f"../{VMEC_FILE}", ".") + if case_name == "boozxform": + shutil.copy2(f"../{BOOZ_FILE}", ".") + + # Write input file + with open("simple.in", "w") as f: + f.write(input_content) + + # Run SIMPLE + start_time = time.time() + result = subprocess.run([f"../{SIMPLE_EXEC}"], capture_output=True, text=True) + end_time = time.time() + + if result.returncode != 0: + print(f"Error running {case_name}:") + print(result.stderr) + os.chdir("..") + return None, None + + print(f"Execution time: {end_time - start_time:.2f} seconds") + + # Read results + if os.path.exists("confined_fraction.dat"): + confined_data = np.loadtxt("confined_fraction.dat") + else: + print(f"No confined_fraction.dat found for {case_name}") + confined_data = None + + if os.path.exists("times_lost.dat"): + times_lost = np.loadtxt("times_lost.dat") + else: + times_lost = None + + os.chdir("..") + return confined_data, times_lost + +def plot_comparison(vmec_data, booz_data): + """Plot comparison of results""" + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8)) + + # Plot confined fraction vs time + if vmec_data[0] is not None and booz_data[0] is not None: + ax1.plot(vmec_data[0][:, 0], vmec_data[0][:, 1], 'b-', label='VMEC', linewidth=2) + ax1.plot(booz_data[0][:, 0], booz_data[0][:, 1], 'r--', label='BOOZXFORM', linewidth=2) + ax1.set_xlabel('Time') + ax1.set_ylabel('Confined Fraction') + ax1.set_title('Confinement Comparison: VMEC vs BOOZXFORM') + ax1.legend() + ax1.grid(True, alpha=0.3) + + # Plot difference + if vmec_data[0] is not None and booz_data[0] is not None: + # Interpolate to common time grid if needed + t_common = vmec_data[0][:, 0] + vmec_conf = vmec_data[0][:, 1] + booz_conf = np.interp(t_common, booz_data[0][:, 0], booz_data[0][:, 1]) + + diff = abs(vmec_conf - booz_conf) + ax2.semilogy(t_common, diff, 'g-', linewidth=2) + ax2.set_xlabel('Time') + ax2.set_ylabel('|Difference|') + ax2.set_title('Absolute Difference in Confined Fraction') + ax2.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig('benchmark_comparison.png', dpi=150) + print("\nSaved comparison plot to benchmark_comparison.png") + +def main(): + print("=== SIMPLE BOOZXFORM Benchmark ===") + print(f"Comparing VMEC field evaluation vs BOOZXFORM field evaluation") + print(f"VMEC file: {VMEC_FILE}") + print(f"BOOZ file: {BOOZ_FILE}") + + # Check if files exist + if not os.path.exists(SIMPLE_EXEC): + print(f"Error: SIMPLE executable not found at {SIMPLE_EXEC}") + print("Please build SIMPLE first with 'make' in the project root") + return + + if not os.path.exists(VMEC_FILE): + print(f"Error: VMEC file not found at {VMEC_FILE}") + return + + if not os.path.exists(BOOZ_FILE): + print(f"Error: BOOZXFORM file not found at {BOOZ_FILE}") + return + + # Run benchmarks + vmec_data = run_benchmark("vmec", simple_in_vmec) + booz_data = run_benchmark("boozxform", simple_in_booz) + + # Compare results + if vmec_data[0] is not None and booz_data[0] is not None: + print("\n=== Results Summary ===") + print(f"VMEC final confined fraction: {vmec_data[0][-1, 1]:.6f}") + print(f"BOOZ final confined fraction: {booz_data[0][-1, 1]:.6f}") + print(f"Difference: {abs(vmec_data[0][-1, 1] - booz_data[0][-1, 1]):.6e}") + + # Plot comparison + plot_comparison(vmec_data, booz_data) + else: + print("\nBenchmark failed - check error messages above") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/benchmark_orbit b/benchmark_orbit new file mode 120000 index 00000000..97f08690 --- /dev/null +++ b/benchmark_orbit @@ -0,0 +1 @@ +../benchmark_orbit \ No newline at end of file diff --git a/draft/.gitignore b/draft/.gitignore new file mode 100644 index 00000000..1be4b751 --- /dev/null +++ b/draft/.gitignore @@ -0,0 +1,7 @@ +dcon_purely_toroidal_field.txt +input.purely_toroidal_field +log.purely_toroidal_field +mercier.purely_toroidal_field +parvmecinfo.txt +threed1.purely_toroidal_field +timings.txt diff --git a/draft/Makefile b/draft/Makefile new file mode 100644 index 00000000..2185b3e6 --- /dev/null +++ b/draft/Makefile @@ -0,0 +1,13 @@ +VMEC=xvmec + +wout.nc: input.purely_toroidal_field + ${VMEC} purely_toroidal_field + mv wout_purely_toroidal_field.nc wout.nc + +input.purely_toroidal_field: + wget https://raw.githubusercontent.com/hiddenSymmetries/simsopt/refs/heads/master/tests/test_files/input.purely_toroidal_field + +clean: + rm -f dcon_purely_toroidal_field.txt input.purely_toroidal_field \ + log.purely_toroidal_field mercier.purely_toroidal_field \ + parvmecinfo.txt threed1.purely_toroidal_field timings.txt diff --git a/draft/compare_simsopt.py b/draft/compare_simsopt.py new file mode 100644 index 00000000..ecde7e44 --- /dev/null +++ b/draft/compare_simsopt.py @@ -0,0 +1,180 @@ +# %% +import numpy as np +import matplotlib.pyplot as plt +from simsopt.mhd.vmec import Vmec +from simsopt.field.boozermagneticfield import BoozerRadialInterpolant +import pysimple + + +GAUSS_IN_TESLA = 1e-4 +CM_IN_M = 1e-2 + +wout_file = "wout.nc" + +BOOZER = pysimple.magfie_sub.boozer + +tracer = pysimple.simple.Tracer() +pysimple.params.read_config("simple.in") +pysimple.simple_main.init_field(tracer, "wout.nc", 5, 5, 5, 1) + +# %% Check VMEC to real space +s = 0.5 * np.ones(100) +th = np.linspace(0, 2 * np.pi, 100) +ph = np.zeros(100) +R = np.empty_like(s) +Z = np.empty_like(s) +for i in range(len(s)): + R[i], Z[i] = pysimple.get_can_sub.vmec_to_cyl(s[i], th[i], ph[i]) +R = np.array(R) * CM_IN_M +Z = np.array(Z) * CM_IN_M + +plt.figure() +plt.plot(R[0], Z[0], "o") +plt.quiver( + R[:-1], + Z[:-1], + R[1:] - R[:-1], + Z[1:] - Z[:-1], + angles="xy", + scale_units="xy", + scale=1, + color="r", +) + +#%% +vmec = Vmec(wout_file) +boozer = BoozerRadialInterpolant(vmec, order=3, mpol=16, ntor=16) +points = np.array([s, th, ph]).T.copy() +boozer.set_points(points) +#%% +R = boozer.R() +Z = boozer.Z() + +plt.figure() +plt.plot(R[0], Z[0], "o") +plt.quiver( + R[:-1], + Z[:-1], + R[1:] - R[:-1], + Z[1:] - Z[:-1], + angles="xy", + scale_units="xy", + scale=1, + color="r", +) + +# %% +bder = np.empty(3) +hcovar = np.empty(3) +hctrvr = np.empty(3) +hcurl = np.empty(3) + +s = np.linspace(0, 1, 100) +th = np.ones_like(s) +ph = np.ones_like(s) + + +def evaluate_magfie_simple(s, th, ph, magfie_routine): + modB = np.empty(len(s)) + sqrtg = np.empty(len(s)) + dmodBds = np.empty(len(s)) + dmodBdtheta = np.empty(len(s)) + dmodBdzeta = np.empty(len(s)) + Bscov = np.empty(len(s)) + Bthetacov = np.empty(len(s)) + Bzetacov = np.empty(len(s)) + for i in range(len(s)): + x = np.array([s[i], th[i], ph[i]]) + modB[i], sqrtg[i] = magfie_routine(x, bder, hcovar, hctrvr, hcurl) + modB[i] = modB[i] * GAUSS_IN_TESLA + sqrtg[i] = sqrtg[i] * CM_IN_M**3 + dmodBds[i] = bder[0] * modB[i] + dmodBdtheta[i] = bder[1] * modB[i] + dmodBdzeta[i] = bder[2] * modB[i] + Bscov[i] = hcovar[0] * modB[i] * CM_IN_M + Bthetacov[i] = hcovar[1] * modB[i] * CM_IN_M + Bzetacov[i] = hcovar[2] * modB[i] * CM_IN_M + + return ( + modB, + sqrtg, + dmodBds, + dmodBdtheta, + dmodBdzeta, + Bscov, + Bthetacov, + Bzetacov, + ) + + +print("Evaluating SIMPLE VMEC magnetic field") + +modB, sqrtg, dmodBds, dmodBdtheta, dmodBdzeta, Bscov, Bthetacov, Bzetacov = ( + evaluate_magfie_simple(s, th, ph, pysimple.magfie_sub.magfie_vmec) +) + +print("modB = ", modB[5]) +print("sqrtg = ", sqrtg[5]) +print("Bthetacov = ", Bthetacov[5]) +print("Bzetacov = ", Bzetacov[5]) + +print("Evaluating SIMPLE Boozer magnetic field") + +modB, sqrtg, dmodBds, dmodBdtheta, dmodBdzeta, Bscov, Bthetacov, Bzetacov = ( + evaluate_magfie_simple(s, th, ph, pysimple.magfie_sub.magfie_boozer) +) + +print("modB = ", modB[5]) +print("sqrtg = ", sqrtg[5]) +print("Bthetacov = ", Bthetacov[5]) +print("Bzetacov = ", Bzetacov[5]) +# %% +# %% + +points = np.array([s, th, ph]).T.copy() + +boozer.set_points(points) +modB_simsopt = boozer.modB() +sqrtg_simsopt = (boozer.G() + boozer.iota() * boozer.I()) / modB_simsopt**2 + +# %% +plt.figure() +plt.plot(s, boozer.dmodBds(), label="dmodBds SIMSOPT") +plt.plot(s, boozer.dmodBdtheta(), label="dmodBdtheta SIMSOPT") +plt.plot(s, boozer.dmodBdzeta(), label="dmodBdzeta SIMSOPT") +plt.plot(s, dmodBds, "--", label="dmodBds SIMPLE") +plt.plot(s, dmodBdtheta, "--", label="dmodBdtheta SIMPLE") +plt.plot(s, dmodBdzeta, "--", label="dmodBdzeta SIMPLE") +plt.gca().set_yscale("symlog", linthresh=1e-2) +plt.xlabel("s") +plt.ylabel("|B| derivatives [T]") +plt.legend() + +# %% +plt.figure() +plt.plot(s, boozer.K(), label="Bscov SIMSOPT") +plt.plot(s, boozer.I(), label="Bthetacov SIMSOPT") +plt.plot(s, boozer.G(), label="Bzetacov SIMSOPT") +plt.plot(s, Bscov, "--", label="Bscov SIMPLE (set to zero)") +plt.plot(s, Bthetacov, "--", label="Bthetacov SIMPLE") +plt.plot(s, Bzetacov, "--", label="Bzetacov SIMPLE") +plt.gca().set_yscale("symlog", linthresh=1e-5) +plt.xlabel("s") +plt.ylabel("|B| derivatives [T]") +plt.legend() + +# %% +plt.figure() +plt.plot(s, -sqrtg_simsopt, label="-sqrtg SIMSOPT") +plt.plot(s, sqrtg / boozer.psi0, "--", label="sqrtg SIMPLE") +# plt.gca().set_yscale("log") +plt.xlabel("s") +plt.ylabel("sqrtg [m^3]") +plt.legend() +plt.show() + +# %% +plt.figure() +plt.plot(s, boozer.iota(), label="iota SIMSOPT") + +# %% diff --git a/draft/pysimple.py b/draft/pysimple.py new file mode 120000 index 00000000..52fdbd8b --- /dev/null +++ b/draft/pysimple.py @@ -0,0 +1 @@ +../build/pysimple.py \ No newline at end of file diff --git a/draft/simple.in b/draft/simple.in new file mode 100644 index 00000000..7b102356 --- /dev/null +++ b/draft/simple.in @@ -0,0 +1,7 @@ +&config +trace_time = 1d-2 ! slowing down time, s +sbeg = 0.3d0 ! starting s (normalzed toroidal flux) for particles. +ntestpart = 128 ! number of test particles +netcdffile = 'wout.nc' ! name of VMEC file in NETCDF format +isw_field_type = 2 ! -1: Testing, 0: Canonical, 1: VMEC, 2: Boozer, 3: Meiss, 4: Albert +/ diff --git a/examples/example.py b/examples/example.py index 472b0546..2e6832b1 100644 --- a/examples/example.py +++ b/examples/example.py @@ -1,62 +1,71 @@ -#%% -from pysimple import simple_main, simple, params, orbit_symplectic -from pysimple import get_can_sub as coord - -import numpy as np +# %% import matplotlib.pyplot as plt -#%% -tracy = simple.Tracer() +import numpy as np +import pysimple as ps + +# %% +tracy = ps.simple.Tracer() -simple_main.init_field(tracy, "wout.nc", 3, 3, 3, 1) -params.params_init() +ps.simple_main.init_field(tracy, "wout.nc", 3, 3, 3, 1) +ps.params.params_init() -#%% Initial conditions -z0_vmec = np.array([0.8, 1.0, 0.2, 1.0, 0.5]) # s, th, ph, v/v_th, v_par/v +# %% Initial conditions +z0_vmec = np.array([0.8, 1.0, 0.2, 1.0, 0.5]) # s, th, ph, v/v_th, v_par/v z0_can = z0_vmec.copy() # s, th_c, ph_c, v/v_th, v_par/v -z0_can[1:3] = coord.vmec_to_can(z0_vmec[0], z0_vmec[1], z0_vmec[2]) +z0_can[1:3] = ps.coord.vmec_to_can(z0_vmec[0], z0_vmec[1], z0_vmec[2]) -simple.init_sympl(tracy.si, tracy.f, z0_can, params.dtaumin, params.dtaumin, 1e-13, tracy.integmode) +ps.simple.init_sympl( + tracy.si, + tracy.f, + z0_can, + ps.params.dtaumin, + ps.params.dtaumin, + 1e-13, + tracy.integmode, +) -print(f'B = {tracy.f.bmod}') +print(f"B = {tracy.f.bmod}") nt = 10000 z_integ = np.zeros([nt, 4]) # s, th_c, ph_c, p_phi z_vmec = np.zeros([nt, 5]) # s, th, ph, v/v_th, v_par/v z_cyl = np.zeros([nt, 3]) -z_integ[0,:] = tracy.si.z -z_vmec[0,:] = z0_vmec -z_cyl[0,:2] = coord.vmec_to_cyl(z_vmec[0, 0], z_vmec[0, 1], z_vmec[0, 2]) +z_integ[0, :] = tracy.si.z +z_vmec[0, :] = z0_vmec +z_cyl[0, :2] = ps.coord.vmec_to_cyl(z_vmec[0, 0], z_vmec[0, 1], z_vmec[0, 2]) z_cyl[0, 2] = z_vmec[0, 2] -for kt in range(nt-1): - orbit_symplectic.orbit_timestep_sympl_expl_impl_euler(tracy.si, tracy.f) - z_integ[kt+1, :] = tracy.si.z - z_vmec[kt+1, 0] = z_integ[kt+1, 0] - z_vmec[kt+1, 1:3] = coord.can_to_vmec( - z_integ[kt+1, 0], z_integ[kt+1, 1], z_integ[kt+1, 2]) - z_vmec[kt+1, 3] = np.sqrt(tracy.f.mu*tracy.f.bmod+0.5*tracy.f.vpar**2) - z_vmec[kt+1, 4] = tracy.f.vpar/(z_vmec[kt+1, 3]*np.sqrt(2)) - z_cyl[kt+1, :2] = coord.vmec_to_cyl(z_vmec[kt+1, 0], z_vmec[kt+1, 1], z_vmec[kt+1, 2]) - z_cyl[kt+1, 2] = z_vmec[kt+1, 2] +for kt in range(nt - 1): + ps.orbit_symplectic.orbit_timestep_sympl_expl_impl_euler(tracy.si, tracy.f) + z_integ[kt + 1, :] = tracy.si.z + z_vmec[kt + 1, 0] = z_integ[kt + 1, 0] + z_vmec[kt + 1, 1:3] = ps.coord.can_to_vmec( + z_integ[kt + 1, 0], z_integ[kt + 1, 1], z_integ[kt + 1, 2] + ) + z_vmec[kt + 1, 3] = np.sqrt(tracy.f.mu * tracy.f.bmod + 0.5 * tracy.f.vpar**2) + z_vmec[kt + 1, 4] = tracy.f.vpar / (z_vmec[kt + 1, 3] * np.sqrt(2)) + z_cyl[kt + 1, :2] = ps.coord.vmec_to_cyl( + z_vmec[kt + 1, 0], z_vmec[kt + 1, 1], z_vmec[kt + 1, 2] + ) + z_cyl[kt + 1, 2] = z_vmec[kt + 1, 2] plt.figure() -plt.plot(z_vmec[:, 0]*np.cos(z_vmec[:, 1]), z_vmec[:, 0]*np.sin(z_vmec[:, 1])) -plt.xlabel('s * cos(th)') -plt.ylabel('s * sin(th)') -plt.title('Poloidal orbit topology') +plt.plot(z_vmec[:, 0] * np.cos(z_vmec[:, 1]), z_vmec[:, 0] * np.sin(z_vmec[:, 1])) +plt.xlabel("s * cos(th)") +plt.ylabel("s * sin(th)") +plt.title("Poloidal orbit topology") plt.figure() plt.plot(z_vmec[:, 3]) plt.plot(z_vmec[:, 4]) -plt.xlabel('Timestep') -plt.ylabel('Normalized velocity') -plt.legend(['v/v_0', 'v_par/v']) -plt.title('Velocities over time') +plt.xlabel("Timestep") +plt.ylabel("Normalized velocity") +plt.legend(["v/v_0", "v_par/v"]) +plt.title("Velocities over time") # 3D plot -from mpl_toolkits import mplot3d # R = S0 + s*cos(th) # Z = s*sin(th) # X = R*cos(ph) @@ -65,40 +74,42 @@ S0 = 5.0 plt.figure() -ax = plt.axes(projection='3d') -ax.plot3D((S0 + z_vmec[:, 0]*np.cos(z_vmec[:, 1]))*np.cos(z_vmec[:, 2]), - (S0 + z_vmec[:, 0]*np.cos(z_vmec[:, 1]))*np.sin(z_vmec[:,2]), - z_vmec[:, 0]*np.sin(z_vmec[:, 1])) -ax.set_xlabel('(3 + s * cos(th))*cos(ph)') -ax.set_ylabel('(3 + s * sin(th))*cos(ph)') -ax.set_zlabel('s*cos(ph)') -ax.set_title('3D orbit topology') -ax.set_xlim(-5,5) -ax.set_ylim(-5,5) -ax.set_zlim(-5,5) +ax = plt.axes(projection="3d") +ax.plot3D( + (S0 + z_vmec[:, 0] * np.cos(z_vmec[:, 1])) * np.cos(z_vmec[:, 2]), + (S0 + z_vmec[:, 0] * np.cos(z_vmec[:, 1])) * np.sin(z_vmec[:, 2]), + z_vmec[:, 0] * np.sin(z_vmec[:, 1]), +) +ax.set_xlabel("(3 + s * cos(th))*cos(ph)") +ax.set_ylabel("(3 + s * sin(th))*cos(ph)") +ax.set_zlabel("s*cos(ph)") +ax.set_title("3D orbit topology") +ax.set_xlim(-5, 5) +ax.set_ylim(-5, 5) +ax.set_zlim(-5, 5) # Poloidal orbit in RZ plt.figure() plt.plot(z_cyl[:, 0], z_cyl[:, 1]) -plt.xlabel('R') -plt.ylabel('Z') -plt.title('Poloidal orbit projection') +plt.xlabel("R") +plt.ylabel("Z") +plt.title("Poloidal orbit projection") # 3D orbit in RZ plt.figure() -ax = plt.axes(projection='3d') -ax.plot(z_cyl[:, 0]*np.cos(z_cyl[:,2]), - z_cyl[:, 0]*np.sin(z_cyl[:,2]), - z_cyl[:,1]) -ax.set_xlabel('X') -ax.set_ylabel('Y') -ax.set_zlabel('Z') -ax.set_title('Orbit in real space') -ax.set_xlim(-1400,1400) -ax.set_ylim(-1400,1400) -ax.set_zlim(-1400,1400) +ax = plt.axes(projection="3d") +ax.plot( + z_cyl[:, 0] * np.cos(z_cyl[:, 2]), z_cyl[:, 0] * np.sin(z_cyl[:, 2]), z_cyl[:, 1] +) +ax.set_xlabel("X") +ax.set_ylabel("Y") +ax.set_zlabel("Z") +ax.set_title("Orbit in real space") +ax.set_xlim(-1400, 1400) +ax.set_ylim(-1400, 1400) +ax.set_zlim(-1400, 1400) # %% diff --git a/examples/simple.in b/examples/simple.in index 7b102356..ceaedf5b 100644 --- a/examples/simple.in +++ b/examples/simple.in @@ -1,4 +1,5 @@ &config +ntimstep = 10001 trace_time = 1d-2 ! slowing down time, s sbeg = 0.3d0 ! starting s (normalzed toroidal flux) for particles. ntestpart = 128 ! number of test particles diff --git a/golden_record/simple_main b/golden_record/simple_main new file mode 160000 index 00000000..31e7fe40 --- /dev/null +++ b/golden_record/simple_main @@ -0,0 +1 @@ +Subproject commit 31e7fe4019d8ed8c05513c9001a33c9510e6deb1 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cb3da872..cfca6e03 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,6 +8,7 @@ set(SOURCES field/field_coils.f90 field/field_vmec.f90 field/vmec_field_eval.f90 + field/field_booz_xform.f90 field/field_newton.F90 field.F90 field/field_can_base.f90 diff --git a/src/field/field_booz_xform.f90 b/src/field/field_booz_xform.f90 new file mode 100644 index 00000000..7269dc26 --- /dev/null +++ b/src/field/field_booz_xform.f90 @@ -0,0 +1,457 @@ +!> Module for reading and evaluating magnetic fields from BOOZXFORM output files +!> This allows direct input of pre-computed Boozer coordinate fields +module field_booz_xform + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field_base, only: MagneticField + implicit none + + private + public :: BoozXformField + + ! Module variable for warning message + logical :: booz_warning_printed = .false. + + !> Magnetic field representation from BOOZXFORM NetCDF files + type, extends(MagneticField) :: BoozXformField + ! BOOZXFORM data arrays + integer :: ns_b !< Number of flux surfaces + integer :: nfp_b !< Number of field periods + integer :: mboz_b !< Maximum poloidal mode number + integer :: nboz_b !< Maximum toroidal mode number (divided by nfp) + integer :: mnboz !< Total number of Fourier modes + + real(dp) :: aspect_b !< Aspect ratio + real(dp) :: rmax_b !< Maximum R + real(dp) :: rmin_b !< Minimum R + real(dp) :: betaxis_b !< Beta on axis + + ! Radial arrays (dimension ns_b) + real(dp), allocatable :: s_b(:) !< Normalized toroidal flux + real(dp), allocatable :: iota_b(:) !< Rotational transform + real(dp), allocatable :: buco_b(:) !< Boozer I (poloidal covariant) + real(dp), allocatable :: bvco_b(:) !< Boozer G (toroidal covariant) + real(dp), allocatable :: beta_b(:) !< Plasma beta + real(dp), allocatable :: phip_b(:) !< d(phi)/ds + real(dp), allocatable :: chi_b(:) !< Poloidal flux + real(dp), allocatable :: pres_b(:) !< Pressure + real(dp), allocatable :: phi_b(:) !< Toroidal flux + + ! Fourier mode arrays + integer, allocatable :: ixm_b(:) !< Poloidal mode numbers + integer, allocatable :: ixn_b(:) !< Toroidal mode numbers + integer, allocatable :: jlist(:) !< Surface indices for packed arrays + + ! Fourier coefficient arrays (dimension pack_rad x mnboz) + real(dp), allocatable :: rmnc_b(:,:) !< R Fourier coefficients (cos) + real(dp), allocatable :: zmns_b(:,:) !< Z Fourier coefficients (sin) + real(dp), allocatable :: pmns_b(:,:) !< p = theta_VMEC - theta_Boozer (sin) + real(dp), allocatable :: gmn_b(:,:) !< g = zeta_VMEC - zeta_Boozer + real(dp), allocatable :: bmnc_b(:,:) !< |B| Fourier coefficients (cos) + + ! Stellarator symmetry flag + logical :: lasym_b + + contains + procedure :: load => load_booz_xform + procedure :: evaluate => evaluate_booz_xform + procedure :: cleanup => cleanup_booz_xform + end type BoozXformField + +contains + + !> Load BOOZXFORM data from NetCDF file + subroutine load_booz_xform(this, filename) + use netcdf + class(BoozXformField), intent(inout) :: this + character(len=*), intent(in) :: filename + + integer :: ncid, lasym_int + integer :: ns_dim, mn_dim, pack_dim + integer :: dimid, ierr, i, varid + + ! Open NetCDF file using native NetCDF + ierr = nf90_open(filename, NF90_NOWRITE, ncid) + if (ierr /= nf90_noerr) then + print *, 'Error opening file:', trim(filename) + print *, 'Error:', trim(nf90_strerror(ierr)) + error stop + end if + print *, 'Opened NetCDF file:', trim(filename) + + ! Read dimensions using NetCDF directly + ierr = nf90_inq_dimid(ncid, 'radius', dimid) + if (ierr == nf90_noerr) then + ierr = nf90_inquire_dimension(ncid, dimid, len=ns_dim) + else + print *, 'Error reading radius dimension:', trim(nf90_strerror(ierr)) + error stop + end if + + ierr = nf90_inq_dimid(ncid, 'mn_modes', dimid) + if (ierr == nf90_noerr) then + ierr = nf90_inquire_dimension(ncid, dimid, len=mn_dim) + else + print *, 'Error reading mn_modes dimension:', trim(nf90_strerror(ierr)) + error stop + end if + + ierr = nf90_inq_dimid(ncid, 'pack_rad', dimid) + if (ierr == nf90_noerr) then + ierr = nf90_inquire_dimension(ncid, dimid, len=pack_dim) + else + print *, 'Error reading pack_rad dimension:', trim(nf90_strerror(ierr)) + error stop + end if + + this%ns_b = ns_dim + this%mnboz = mn_dim + + ! Read scalar variables using native NetCDF + print *, 'Reading scalar variables...' + + ierr = nf90_inq_varid(ncid, 'nfp_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%nfp_b) + + ierr = nf90_inq_varid(ncid, 'mboz_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%mboz_b) + + ierr = nf90_inq_varid(ncid, 'nboz_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%nboz_b) + + ierr = nf90_inq_varid(ncid, 'aspect_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%aspect_b) + + ierr = nf90_inq_varid(ncid, 'rmax_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%rmax_b) + + ierr = nf90_inq_varid(ncid, 'rmin_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%rmin_b) + + ierr = nf90_inq_varid(ncid, 'betaxis_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%betaxis_b) + + ! Check stellarator symmetry + ierr = nf90_inq_varid(ncid, 'lasym__logical__', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, lasym_int) + this%lasym_b = (lasym_int == 1) + end if + + ! Allocate arrays + print *, 'Allocating arrays: ns_dim =', ns_dim, 'mn_dim =', mn_dim, 'pack_dim =', pack_dim + allocate(this%s_b(ns_dim)) + allocate(this%iota_b(ns_dim)) + allocate(this%buco_b(ns_dim)) + allocate(this%bvco_b(ns_dim)) + allocate(this%beta_b(ns_dim)) + allocate(this%phip_b(ns_dim)) + allocate(this%chi_b(ns_dim)) + allocate(this%pres_b(ns_dim)) + allocate(this%phi_b(ns_dim)) + + allocate(this%ixm_b(mn_dim)) + allocate(this%ixn_b(mn_dim)) + allocate(this%jlist(pack_dim)) + + allocate(this%rmnc_b(pack_dim, mn_dim)) + allocate(this%zmns_b(pack_dim, mn_dim)) + allocate(this%pmns_b(pack_dim, mn_dim)) + allocate(this%gmn_b(pack_dim, mn_dim)) + allocate(this%bmnc_b(pack_dim, mn_dim)) + + ! Read radial arrays + print *, 'Reading radial arrays...' + + ierr = nf90_inq_varid(ncid, 'iota_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%iota_b) + if (ierr /= nf90_noerr) print *, 'Error reading iota_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'buco_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%buco_b) + if (ierr /= nf90_noerr) print *, 'Error reading buco_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'bvco_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%bvco_b) + if (ierr /= nf90_noerr) print *, 'Error reading bvco_b:', trim(nf90_strerror(ierr)) + end if + + ! Read remaining arrays + ierr = nf90_inq_varid(ncid, 'beta_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%beta_b) + + ierr = nf90_inq_varid(ncid, 'phip_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%phip_b) + + ierr = nf90_inq_varid(ncid, 'chi_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%chi_b) + + ierr = nf90_inq_varid(ncid, 'pres_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%pres_b) + + ierr = nf90_inq_varid(ncid, 'phi_b', varid) + if (ierr == nf90_noerr) ierr = nf90_get_var(ncid, varid, this%phi_b) + + ! Create s array (normalized toroidal flux) + this%s_b = [(real(i-1, dp) / real(ns_dim-1, dp), i = 1, ns_dim)] + + ! Read mode arrays and Fourier coefficients + print *, 'Reading mode arrays and Fourier coefficients...' + + ierr = nf90_inq_varid(ncid, 'ixm_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%ixm_b) + if (ierr /= nf90_noerr) print *, 'Error reading ixm_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'ixn_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%ixn_b) + if (ierr /= nf90_noerr) print *, 'Error reading ixn_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'jlist', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%jlist) + if (ierr /= nf90_noerr) print *, 'Error reading jlist:', trim(nf90_strerror(ierr)) + end if + + ! Read Fourier coefficients + ierr = nf90_inq_varid(ncid, 'rmnc_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%rmnc_b) + if (ierr /= nf90_noerr) print *, 'Error reading rmnc_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'zmns_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%zmns_b) + if (ierr /= nf90_noerr) print *, 'Error reading zmns_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'pmns_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%pmns_b) + if (ierr /= nf90_noerr) print *, 'Error reading pmns_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'gmn_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%gmn_b) + if (ierr /= nf90_noerr) print *, 'Error reading gmn_b:', trim(nf90_strerror(ierr)) + end if + + ierr = nf90_inq_varid(ncid, 'bmnc_b', varid) + if (ierr == nf90_noerr) then + ierr = nf90_get_var(ncid, varid, this%bmnc_b) + if (ierr /= nf90_noerr) print *, 'Error reading bmnc_b:', trim(nf90_strerror(ierr)) + end if + + ! Close file + ierr = nf90_close(ncid) + if (ierr /= nf90_noerr) then + print *, 'Error closing file:', trim(nf90_strerror(ierr)) + end if + + print *, 'Loaded BOOZXFORM file:', trim(filename) + print *, ' ns =', this%ns_b, ', nfp =', this%nfp_b + print *, ' mboz =', this%mboz_b, ', nboz =', this%nboz_b + print *, ' mnboz =', this%mnboz, ' Fourier modes' + print *, ' iota_b(1) =', this%iota_b(1), 'iota_b(ns) =', this%iota_b(this%ns_b) + print *, ' First few mode numbers:' + print *, ' m =', this%ixm_b(1:min(5,size(this%ixm_b))) + print *, ' n =', this%ixn_b(1:min(5,size(this%ixn_b))) + + end subroutine load_booz_xform + + !> Evaluate magnetic field at given coordinates + !> Input: VMEC coordinates (r, theta_vmec, phi_vmec) + !> Output: Magnetic field components + !> Note: This is a simplified implementation that evaluates B directly in Boozer coordinates + !> without the full coordinate transformation + subroutine evaluate_booz_xform(self, x, Acov, hcov, Bmod, sqgBctr) + class(BoozXformField), intent(in) :: self + real(dp), intent(in) :: x(3) ! r=sqrt(s_vmec), theta_vmec, phi_vmec + real(dp), intent(out) :: Acov(3) + real(dp), intent(out) :: hcov(3) + real(dp), intent(out) :: Bmod + real(dp), intent(out), optional :: sqgBctr(3) + + real(dp) :: s, theta_b, zeta_b, theta_v, zeta_v + real(dp) :: p_angle, g_angle ! theta_vmec - theta_booz, zeta_vmec - zeta_booz + real(dp) :: R_booz, Z_booz, B_booz, Jac_booz + real(dp) :: dR_dtheta, dR_dzeta, dZ_dtheta, dZ_dzeta + real(dp) :: sqrtg, Bcov_s, Bcov_theta, Bcov_zeta + real(dp) :: Bctr_theta, Bctr_zeta + real(dp) :: iota_local, I_local, G_local, phip_local + integer :: js, imn, js_pack, js_pack_m1 + real(dp) :: wlo, whi, angle_arg + real(dp), parameter :: twopi = 8.0_dp * atan(1.0_dp) + + print *, 'evaluate_booz_xform: Entry, x =', x + print *, 'evaluate_booz_xform: ns_b =', self%ns_b, 'size(s_b) =', size(self%s_b) + print *, 'evaluate_booz_xform: size(bmnc_b) =', size(self%bmnc_b, 1), size(self%bmnc_b, 2) + + ! Check if data is loaded + if (.not. allocated(self%s_b)) then + print *, 'ERROR: s_b not allocated!' + error stop 'evaluate_booz_xform: BOOZXFORM data not loaded' + end if + + ! Convert input coordinates + s = x(1)**2 ! s = r^2 + theta_v = x(2) + zeta_v = x(3) + + print *, 'evaluate_booz_xform: s =', s + + ! Find radial grid point and interpolation weights + js = minloc(abs(self%s_b - s), 1) + if (js == self%ns_b) js = self%ns_b - 1 + if (js < 2) js = 2 + + ! For simplicity, find the closest packed indices + ! Since jlist goes from 2 to ns_b, and pack_rad = ns_b-1, we can map: + js_pack = js - 1 ! Map surface index to packed index (2->1, 3->2, etc.) + js_pack_m1 = js - 2 ! Previous surface + + ! Ensure we're within bounds + if (js_pack > size(self%bmnc_b, 1)) js_pack = size(self%bmnc_b, 1) + if (js_pack < 1) js_pack = 1 + if (js_pack_m1 < 1) js_pack_m1 = 1 + + ! If we're at the same index, use neighboring indices + if (js_pack == js_pack_m1) then + if (js_pack > 1) then + js_pack_m1 = js_pack - 1 + else + js_pack = 2 + js_pack_m1 = 1 + end if + end if + + whi = (s - self%s_b(js-1)) / max(1.0e-10_dp, self%s_b(js) - self%s_b(js-1)) + wlo = 1.0_dp - whi + + ! Interpolate radial profiles + iota_local = wlo * self%iota_b(js-1) + whi * self%iota_b(js) + I_local = wlo * self%buco_b(js-1) + whi * self%buco_b(js) + G_local = wlo * self%bvco_b(js-1) + whi * self%bvco_b(js) + phip_local = wlo * self%phip_b(js-1) + whi * self%phip_b(js) + + ! For now, use VMEC angles directly as Boozer angles (simplified) + ! This is not the full transformation but should allow basic evaluation + theta_b = theta_v + zeta_b = zeta_v + + ! Now evaluate field quantities in Boozer coordinates + R_booz = 0.0_dp + Z_booz = 0.0_dp + B_booz = 0.0_dp + Jac_booz = 0.0_dp + dR_dtheta = 0.0_dp + dR_dzeta = 0.0_dp + dZ_dtheta = 0.0_dp + dZ_dzeta = 0.0_dp + + do imn = 1, self%mnboz + angle_arg = self%ixm_b(imn) * theta_b - self%ixn_b(imn) * zeta_b + + ! Check array bounds + if (js_pack_m1 < 1 .or. js_pack_m1 > size(self%rmnc_b, 1) .or. & + js_pack < 1 .or. js_pack > size(self%rmnc_b, 1)) cycle + + ! R, B, and Jacobian are cosine components + R_booz = R_booz + cos(angle_arg) * & + (wlo * self%rmnc_b(js_pack_m1, imn) + whi * self%rmnc_b(js_pack, imn)) + B_booz = B_booz + cos(angle_arg) * & + (wlo * self%bmnc_b(js_pack_m1, imn) + whi * self%bmnc_b(js_pack, imn)) + Jac_booz = Jac_booz + cos(angle_arg) * & + (wlo * self%gmn_b(js_pack_m1, imn) + whi * self%gmn_b(js_pack, imn)) + + ! Z is sine component + Z_booz = Z_booz + sin(angle_arg) * & + (wlo * self%zmns_b(js_pack_m1, imn) + whi * self%zmns_b(js_pack, imn)) + + ! Derivatives for metric + dR_dtheta = dR_dtheta - self%ixm_b(imn) * sin(angle_arg) * & + (wlo * self%rmnc_b(js_pack_m1, imn) + whi * self%rmnc_b(js_pack, imn)) + dR_dzeta = dR_dzeta + self%ixn_b(imn) * sin(angle_arg) * & + (wlo * self%rmnc_b(js_pack_m1, imn) + whi * self%rmnc_b(js_pack, imn)) + dZ_dtheta = dZ_dtheta + self%ixm_b(imn) * cos(angle_arg) * & + (wlo * self%zmns_b(js_pack_m1, imn) + whi * self%zmns_b(js_pack, imn)) + dZ_dzeta = dZ_dzeta - self%ixn_b(imn) * cos(angle_arg) * & + (wlo * self%zmns_b(js_pack_m1, imn) + whi * self%zmns_b(js_pack, imn)) + end do + + ! Set output magnetic field magnitude + Bmod = B_booz + + ! Compute sqrt(g) for Boozer coordinates + ! In Boozer coords: sqrt(g) = (G + iota*I) / B^2 + sqrtg = (G_local + iota_local * I_local) / (B_booz * B_booz) + + ! Contravariant B components in Boozer coordinates + Bctr_theta = I_local / sqrtg + Bctr_zeta = G_local / sqrtg + + ! Covariant B components in Boozer coordinates + ! B = grad(psi) x grad(theta) + iota * grad(psi) x grad(zeta) + Bcov_s = 0.0_dp ! No radial component in Boozer coords + Bcov_theta = I_local + Bcov_zeta = G_local + + ! For straight field line coordinates, the vector potential has simple form + ! A_theta = -psi (toroidal flux function) + ! A_zeta = 0 in Boozer coordinates + Acov(1) = 0.0_dp ! A_s = 0 + Acov(2) = -(wlo * self%phi_b(js-1) + whi * self%phi_b(js)) / twopi ! A_theta = -Phi/(2*pi) + Acov(3) = 0.0_dp ! A_zeta = 0 in Boozer + + ! Normalized covariant B components + hcov(1) = Bcov_s / Bmod + hcov(2) = Bcov_theta / Bmod + hcov(3) = Bcov_zeta / Bmod + + ! Transform back to VMEC coordinates if needed + ! This is a simplified transformation - may need refinement + ! The covariant components transform with the inverse Jacobian + ! For now, we use the fact that the transformation mainly affects angles + + if (present(sqgBctr)) then + sqgBctr(1) = 0.0_dp + sqgBctr(2) = sqrtg * Bctr_theta + sqgBctr(3) = sqrtg * Bctr_zeta + end if + + end subroutine evaluate_booz_xform + + !> Clean up allocated arrays + subroutine cleanup_booz_xform(self) + class(BoozXformField), intent(inout) :: self + + if (allocated(self%s_b)) deallocate(self%s_b) + if (allocated(self%iota_b)) deallocate(self%iota_b) + if (allocated(self%buco_b)) deallocate(self%buco_b) + if (allocated(self%bvco_b)) deallocate(self%bvco_b) + if (allocated(self%beta_b)) deallocate(self%beta_b) + if (allocated(self%phip_b)) deallocate(self%phip_b) + if (allocated(self%chi_b)) deallocate(self%chi_b) + if (allocated(self%pres_b)) deallocate(self%pres_b) + if (allocated(self%phi_b)) deallocate(self%phi_b) + if (allocated(self%ixm_b)) deallocate(self%ixm_b) + if (allocated(self%ixn_b)) deallocate(self%ixn_b) + if (allocated(self%jlist)) deallocate(self%jlist) + if (allocated(self%rmnc_b)) deallocate(self%rmnc_b) + if (allocated(self%zmns_b)) deallocate(self%zmns_b) + if (allocated(self%pmns_b)) deallocate(self%pmns_b) + if (allocated(self%gmn_b)) deallocate(self%gmn_b) + if (allocated(self%bmnc_b)) deallocate(self%bmnc_b) + + end subroutine cleanup_booz_xform + +end module field_booz_xform \ No newline at end of file diff --git a/src/field_can.f90 b/src/field_can.f90 index cef52661..362a4289 100644 --- a/src/field_can.f90 +++ b/src/field_can.f90 @@ -1,7 +1,7 @@ module field_can_mod use diag_mod, only : icounter use boozer_sub, only : splint_boozer_coord -use magfie_sub, only : TEST, CANFLUX, BOOZER, MEISS, ALBERT +use magfie_sub, only : TEST, CANFLUX, BOOZER, MEISS, ALBERT, BOOZXFORM use field, only : MagneticField, VmecField use field_can_base, only : twopi, evaluate_base => evaluate, coordinate_transform, & identity_transform, FieldCan @@ -57,6 +57,10 @@ subroutine field_can_from_name(field_name, field_noncan) evaluate => evaluate_albert can_to_ref => can_to_ref_albert ref_to_can => ref_to_can_albert + case("boozxform") + ! BOOZXFORM field evaluation is handled directly in magfie_boozxform + ! No canonical coordinate transformation needed + ! Use test evaluation as a placeholder (should not be called) case default print *, "field_can_from_name: Unknown field type ", field_name error stop @@ -93,6 +97,8 @@ function name_from_id(field_id) name_from_id = "meiss" case(ALBERT) name_from_id = "albert" + case(BOOZXFORM) + name_from_id = "boozxform" case default print *, "name_from_id: Unknown field id ", field_id error stop @@ -149,6 +155,9 @@ subroutine init_field_can(field_id, field_noncan) call get_meiss_coordinates case (ALBERT) call get_albert_coordinates + case (BOOZXFORM) + ! BOOZXFORM uses pre-computed Boozer coordinates, no transformation needed + ! Field evaluation is handled directly in magfie_boozxform case default print *, "init_field_can: Unknown field id ", field_id error stop diff --git a/src/magfie.f90 b/src/magfie.f90 index c689829c..760e17c9 100644 --- a/src/magfie.f90 +++ b/src/magfie.f90 @@ -2,6 +2,7 @@ module magfie_sub use spline_vmec_sub, only: vmec_field use field_can_meiss, only: magfie_meiss use field_can_albert, only: magfie_albert +use field_booz_xform implicit none @@ -30,27 +31,47 @@ end subroutine magfie_base procedure(magfie_base), pointer :: magfie => null() -integer, parameter :: TEST=-1, CANFLUX=0, VMEC=1, BOOZER=2, MEISS=3, ALBERT=4 +integer, parameter :: TEST=-1, CANFLUX=0, VMEC=1, BOOZER=2, MEISS=3, ALBERT=4, BOOZXFORM=5 + +! Module-level variable for BOOZXFORM field +type(BoozXformField), save :: booz_field +logical, save :: booz_field_loaded = .false. +character(256) :: boozxform_filename = '' contains subroutine init_magfie(id) integer, intent(in) :: id + print *, 'init_magfie: called with id =', id + select case(id) case(TEST) print *, 'init_magfie: magfie_test not implemented' error stop case(CANFLUX) + print *, 'init_magfie: setting magfie => magfie_can' magfie => magfie_can case(VMEC) + print *, 'init_magfie: setting magfie => magfie_vmec' magfie => magfie_vmec case(BOOZER) + print *, 'init_magfie: setting magfie => magfie_boozer' magfie => magfie_boozer case(MEISS) + print *, 'init_magfie: setting magfie => magfie_meiss' magfie => magfie_meiss case(ALBERT) + print *, 'init_magfie: setting magfie => magfie_albert' magfie => magfie_albert + case(BOOZXFORM) + print *, 'init_magfie: setting magfie => magfie_boozxform' + print *, 'init_magfie: boozxform_filename =', trim(boozxform_filename) + magfie => magfie_boozxform + ! Set default filename for testing - should be set by caller + if (trim(boozxform_filename) == '') then + boozxform_filename = 'boozmn_LandremanPaul2021_QA_lowres.nc' + end if case default print *,'init_magfie: unknown id ', id error stop @@ -396,4 +417,97 @@ subroutine magfie_boozer(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) ! end subroutine magfie_boozer +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +! + subroutine magfie_boozxform(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) +! +! Magnetic field evaluation using pre-computed BOOZXFORM data +! This wraps the BoozXformField to provide the magfie interface +! +! Input coordinates: x(1)=s (normalized toroidal flux), +! x(2)=theta (VMEC poloidal angle), +! x(3)=phi (VMEC toroidal angle) +! + use, intrinsic :: iso_fortran_env, only: dp => real64 + double precision, intent(in) :: x(3) + double precision, intent(out) :: bmod,sqrtg + double precision, intent(out) :: bder(3),hcovar(3),hctrvr(3),hcurl(3) + + real(dp) :: x_dp(3), Acov(3), sqgBctr(3) + real(dp) :: r, theta, phi, dr_ds + real(dp) :: hs, ht, hp ! step sizes for derivatives + real(dp) :: x_plus(3), x_minus(3) + real(dp) :: Acov_plus(3), Acov_minus(3) + real(dp) :: hcov_plus(3), hcov_minus(3) + real(dp) :: Bmod_plus, Bmod_minus + integer :: i + + ! Load BOOZXFORM data if not already loaded + if (.not. booz_field_loaded) then + if (trim(boozxform_filename) == '') then + error stop 'magfie_boozxform: boozxform_filename not set' + end if + print *, 'magfie_boozxform: Loading BOOZXFORM file:', trim(boozxform_filename) + call booz_field%load(boozxform_filename) + booz_field_loaded = .true. + print *, 'magfie_boozxform: BOOZXFORM file loaded successfully' + end if + + print *, 'magfie_boozxform: Called with x =', x + + ! Convert to appropriate precision and transform coordinates + r = sqrt(x(1)) ! Convert s to r = sqrt(s) + theta = x(2) + phi = x(3) + x_dp = [r, theta, phi] + + ! Evaluate field at current position + call booz_field%evaluate(x_dp, Acov, hcovar, bmod, sqgBctr) + + ! Compute sqrt(g) from contravariant B components + if (sqgBctr(2) /= 0.0_dp .or. sqgBctr(3) /= 0.0_dp) then + sqrtg = bmod / sqrt(sqgBctr(2)*hcovar(2) + sqgBctr(3)*hcovar(3)) + else + sqrtg = 1.0_dp ! Default if contravariant components not properly set + end if + + ! Compute contravariant h components + hctrvr(1) = sqgBctr(1) / (sqrtg * bmod) + hctrvr(2) = sqgBctr(2) / (sqrtg * bmod) + hctrvr(3) = sqgBctr(3) / (sqrtg * bmod) + + ! Compute derivatives of log(B) using finite differences + hs = 1.0d-5 + ht = hs * 6.28318530718d0 ! 2*pi scaling + hp = ht / 5.0d0 + dr_ds = 0.5d0 / r ! dr/ds = 1/(2*sqrt(s)) + + ! Derivative with respect to s + x_plus = [sqrt(x(1) + hs), theta, phi] + x_minus = [sqrt(x(1) - hs), theta, phi] + call booz_field%evaluate(x_plus, Acov_plus, hcov_plus, Bmod_plus) + call booz_field%evaluate(x_minus, Acov_minus, hcov_minus, Bmod_minus) + bder(1) = (Bmod_plus - Bmod_minus) / (2.0d0 * hs * bmod) + + ! Derivative with respect to theta + x_plus = [r, theta + ht, phi] + x_minus = [r, theta - ht, phi] + call booz_field%evaluate(x_plus, Acov_plus, hcov_plus, Bmod_plus) + call booz_field%evaluate(x_minus, Acov_minus, hcov_minus, Bmod_minus) + bder(2) = (Bmod_plus - Bmod_minus) / (2.0d0 * ht * bmod) + + ! Derivative with respect to phi + x_plus = [r, theta, phi + hp] + x_minus = [r, theta, phi - hp] + call booz_field%evaluate(x_plus, Acov_plus, hcov_plus, Bmod_plus) + call booz_field%evaluate(x_minus, Acov_minus, hcov_minus, Bmod_minus) + bder(3) = (Bmod_plus - Bmod_minus) / (2.0d0 * hp * bmod) + + ! Compute curl of h using finite differences + ! curl(h) = (1/sqrt(g)) * [d(h_phi)/dtheta - d(h_theta)/dphi, ...] + ! For now, set to zero (simplified) + hcurl = 0.0d0 + + end subroutine magfie_boozxform + end module magfie_sub diff --git a/src/params.f90 b/src/params.f90 index 77caa9f9..ee36abe7 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -81,6 +81,7 @@ module params integer, dimension (:), allocatable :: idx character(1000) :: field_input = '' + character(256) :: boozxform_file = '' ! BOOZXFORM NetCDF file namelist /config/ notrace_passing, nper, npoiper, ntimstep, ntestpart, & trace_time, num_surf, sbeg, phibeg, thetabeg, contr_pp, & @@ -90,7 +91,7 @@ module params vmec_RZ_scale, swcoll, deterministic, old_axis_healing, & old_axis_healing_boundary, am1, am2, Z1, Z2, & densi1, densi2, tempi1, tempi2, tempe, & - batch_size, ran_seed, reuse_batch, field_input, & + batch_size, ran_seed, reuse_batch, field_input, boozxform_file, & output_error, output_orbits_macrostep ! callback contains diff --git a/src/simple_main.f90 b/src/simple_main.f90 index db084794..e3366aae 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -13,7 +13,7 @@ module simple_main class_plot, ntcut, iclass, bmod00, xi, idx, bmin, bmax, dphi, & zstart, zend, trap_par, perp_inv, volstart, sbeg, thetabeg, phibeg, npoiper, nper, & ntimstep, bstart, ibins, ierr, should_skip, reset_seed_if_deterministic, & - field_input, isw_field_type, reuse_batch + field_input, isw_field_type, reuse_batch, boozxform_file implicit none @@ -44,7 +44,7 @@ end subroutine init_field subroutine run(norb) - use magfie_sub, only : init_magfie, VMEC + use magfie_sub, only : init_magfie, VMEC, BOOZXFORM, boozxform_filename use samplers, only: sample, START_FILE use timing, only : print_phase_time type(Tracer), intent(inout) :: norb @@ -58,6 +58,14 @@ subroutine run(norb) call print_phase_time('Collision initialization completed') endif + ! Set BOOZXFORM filename if using that field type + if (isw_field_type == BOOZXFORM) then + if (trim(boozxform_file) == '') then + error stop 'BOOZXFORM field type requires boozxform_file to be set in namelist' + end if + boozxform_filename = boozxform_file + end if + call init_magfie(VMEC) call print_phase_time('VMEC magnetic field initialization completed') diff --git a/test/booz_xform/compare_boozer_coords.py b/test/booz_xform/compare_boozer_coords.py new file mode 100644 index 00000000..9e01c993 --- /dev/null +++ b/test/booz_xform/compare_boozer_coords.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python3 +""" +Compare Boozer coordinates from SIMPLE's internal conversion +vs booz_xform external file +""" +import numpy as np +import matplotlib.pyplot as plt +import netCDF4 as nc +import sys +import os + +# Add parent directory to path for pysimple +# Handle both source tree and build tree execution +script_dir = os.path.dirname(os.path.abspath(__file__)) +if 'build' in script_dir: + # We're running from build directory + sys.path.insert(0, os.getcwd()) # Current dir should have pysimple +else: + # We're running from source directory + build_dir = os.path.join(script_dir, '../../build') + if os.path.exists(build_dir): + sys.path.insert(0, build_dir) + +try: + import pysimple +except ImportError as e: + print(f"Error: pysimple not found. Please build SIMPLE first.") + print(f"Python path: {sys.path}") + print(f"Import error: {e}") + sys.exit(1) + + +def read_booz_xform_file(filename): + """Read BOOZXFORM NetCDF file and return key data""" + with nc.Dataset(filename, 'r') as f: + # This is the older BOOZXFORM format from Stellopt + data = { + 'ns': int(f.variables['ns_b'][:]), + 'nfp': int(f.variables['nfp_b'][:]), + 'mboz': int(f.variables['mboz_b'][:]), + 'nboz': int(f.variables['nboz_b'][:]), + 'mnboz': int(f.variables['mnboz_b'][:]), + 'iota': f.variables['iota_b'][:], + 'buco': f.variables['buco_b'][:], # Boozer I (poloidal covariant) + 'bvco': f.variables['bvco_b'][:], # Boozer G (toroidal covariant) + 'beta': f.variables['beta_b'][:], + 'phip': f.variables['phip_b'][:], + 'chi': f.variables['chi_b'][:], + 'pres': f.variables['pres_b'][:], + 'ixm': f.variables['ixm_b'][:].astype(int), + 'ixn': f.variables['ixn_b'][:].astype(int), + 'rmnc': f.variables['rmnc_b'][:], + 'zmns': f.variables['zmns_b'][:], + 'pmns': f.variables['pmns_b'][:], # p = theta_VMEC - theta_Boozer + 'gmn': f.variables['gmn_b'][:], # g = zeta_VMEC - zeta_Boozer + 'bmnc': f.variables['bmnc_b'][:], # |B| Fourier coefficients + 'jlist': f.variables['jlist'][:].astype(int), + } + + # Create s array (normalized toroidal flux) + data['s'] = np.linspace(0, 1, data['ns']) + + # Check for stellarator symmetry + try: + data['lasym'] = bool(f.variables['lasym__logical__'][:]) + except: + data['lasym'] = False + + print(f"Loaded BOOZXFORM file: {filename}") + print(f" ns = {data['ns']}, nfp = {data['nfp']}") + print(f" mboz = {data['mboz']}, nboz = {data['nboz']}") + print(f" mnboz = {data['mnboz']} Fourier modes") + print(f" stellarator symmetric: {not data['lasym']}") + + return data + + +def setup_simple_boozer(vmec_file): + """Initialize SIMPLE with internal Boozer conversion""" + tracer = pysimple.simple.Tracer() + + # Read config + pysimple.params.read_config("simple.in") + + # Set to use internal Boozer (isw_field_type = 2) + pysimple.velo_mod.isw_field_type = 2 + + # Initialize field + pysimple.simple_main.init_field(tracer, vmec_file, 5, 5, 5, 1) + + return tracer + + +def evaluate_simple_boozer(s_vals, theta_vals, zeta_vals): + """Evaluate SIMPLE's internal Boozer coordinates at given points""" + npoints = len(s_vals) + results = { + 'modB': np.zeros(npoints), + 'sqrtg': np.zeros(npoints), + 'dBds': np.zeros(npoints), + 'dBdtheta': np.zeros(npoints), + 'dBdzeta': np.zeros(npoints), + 'Bs_cov': np.zeros(npoints), + 'Btheta_cov': np.zeros(npoints), + 'Bzeta_cov': np.zeros(npoints), + } + + # Temporary arrays for field evaluation + bder = np.empty(3) + hcovar = np.empty(3) + hctrvr = np.empty(3) + hcurl = np.empty(3) + + for i in range(npoints): + x = np.array([s_vals[i], theta_vals[i], zeta_vals[i]]) + + # Evaluate field using SIMPLE's Boozer routine + modB, sqrtg = pysimple.magfie_sub.magfie_boozer(x, bder, hcovar, hctrvr, hcurl) + + results['modB'][i] = modB + results['sqrtg'][i] = sqrtg + results['dBds'][i] = bder[0] * modB + results['dBdtheta'][i] = bder[1] * modB + results['dBdzeta'][i] = bder[2] * modB + results['Bs_cov'][i] = hcovar[0] * modB + results['Btheta_cov'][i] = hcovar[1] * modB + results['Bzeta_cov'][i] = hcovar[2] * modB + + return results + + +def compare_coordinate_transform(): + """Compare VMEC->Boozer transformation from SIMPLE vs booz_xform""" + + # Test points in VMEC coordinates + s_test = np.array([0.25, 0.5, 0.75]) + theta_vmec = np.linspace(0, 2*np.pi, 20) + zeta_vmec = np.array([0.0]) # Single toroidal angle + + # Arrays to store results + theta_boozer_simple = np.zeros((len(s_test), len(theta_vmec))) + zeta_boozer_simple = np.zeros((len(s_test), len(theta_vmec))) + + # Check if boozer conversion is available + try: + # Try to access boozer conversion + if hasattr(pysimple, 'boozer_sub'): + boozer_sub = pysimple.boozer_sub + elif hasattr(pysimple, 'Boozer_Coordinates_Mod'): + # Try through Boozer_Coordinates_Mod + boozer_sub = pysimple.Boozer_Coordinates_Mod + elif hasattr(pysimple, 'boozer_coordinates_mod'): + # Try lowercase version + boozer_sub = pysimple.boozer_coordinates_mod + else: + print("Warning: Boozer conversion not directly accessible") + print("Available modules:", [m for m in dir(pysimple) if not m.startswith('_')]) + # For now, just use the VMEC angles as a placeholder + for i in range(len(s_test)): + theta_boozer_simple[i, :] = theta_vmec + zeta_boozer_simple[i, :] = zeta_vmec[0] + return s_test, theta_vmec, theta_boozer_simple, zeta_boozer_simple + + # Get Boozer angles from SIMPLE's internal conversion + for i, s in enumerate(s_test): + for j, th_v in enumerate(theta_vmec): + # Convert VMEC -> Boozer using SIMPLE + if hasattr(boozer_sub, 'vmec_to_boozer'): + th_b, z_b = boozer_sub.vmec_to_boozer(s, th_v, zeta_vmec[0]) + else: + # Fallback - no conversion available + th_b, z_b = th_v, zeta_vmec[0] + theta_boozer_simple[i, j] = th_b + zeta_boozer_simple[i, j] = z_b + except Exception as e: + print(f"Error in coordinate conversion: {e}") + # Use VMEC angles as fallback + for i in range(len(s_test)): + theta_boozer_simple[i, :] = theta_vmec + zeta_boozer_simple[i, :] = zeta_vmec[0] + + return s_test, theta_vmec, theta_boozer_simple, zeta_boozer_simple + + +def plot_comparison(booz_data, s_test, theta_vmec, theta_boozer_simple): + """Create comparison plots""" + + fig, axes = plt.subplots(2, 2, figsize=(12, 10)) + fig.suptitle('SIMPLE Internal Boozer vs BOOZXFORM Comparison', fontsize=14) + + # Plot 1: Iota profile + ax = axes[0, 0] + ax.plot(booz_data['s'], booz_data['iota'], 'b-', label='BOOZXFORM', linewidth=2) + # TODO: Add SIMPLE's iota when accessible + ax.set_xlabel('s') + ax.set_ylabel('ι (rotational transform)') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Rotational Transform Profile') + + # Plot 2: Boozer coordinate transformation + ax = axes[0, 1] + colors = ['red', 'green', 'blue'] + for i, (s, color) in enumerate(zip(s_test, colors)): + p_simple = theta_vmec - theta_boozer_simple[i, :] + ax.plot(theta_vmec, p_simple, '-', color=color, + label=f's={s:.2f} (SIMPLE)', linewidth=2) + + # Add BOOZXFORM p for comparison at s=0.5 + # We need to reconstruct p from Fourier coefficients + s_idx = np.argmin(np.abs(booz_data['s'] - 0.5)) + # Find the corresponding radial index in the packed data + # jlist contains the surface indices for the packed arrays + if s_idx < len(booz_data['jlist']): + j_idx = s_idx # Direct indexing for packed arrays + theta_test = theta_vmec + p_booz = np.zeros_like(theta_test) + # Sum Fourier components for p = theta_VMEC - theta_Boozer + for k in range(booz_data['mnboz']): + m = booz_data['ixm'][k] + n = booz_data['ixn'][k] // booz_data['nfp'] # Normalize by nfp + if j_idx < booz_data['pmns'].shape[0]: + # p_mn is sin(m*theta - n*nfp*zeta) component + p_booz += booz_data['pmns'][j_idx, k] * np.sin(m * theta_test) + ax.plot(theta_vmec, p_booz, 'k--', label='s=0.5 (BOOZXFORM)', linewidth=2) + + ax.set_xlabel('θ_VMEC') + ax.set_ylabel('p = θ_VMEC - θ_Boozer') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Poloidal Angle Transform (ζ = 0)') + + # Plot 3: Fourier spectrum of |B| + ax = axes[1, 0] + # Show first few Fourier modes + modes_to_show = min(20, booz_data['mnboz']) + # Use middle surface for spectrum + s_idx = len(booz_data['jlist']) // 2 + if s_idx < booz_data['bmnc'].shape[0]: + mode_amplitudes = np.abs(booz_data['bmnc'][s_idx, :modes_to_show]) + ax.bar(range(len(mode_amplitudes)), mode_amplitudes) + ax.set_xlabel('Mode index') + ax.set_ylabel('|B| Fourier amplitude') + # Fix indexing for title + actual_s = booz_data['s'][min(s_idx, len(booz_data['s'])-1)] + ax.set_title(f'|B| Fourier Spectrum at s≈{actual_s:.2f}') + ax.set_yscale('log') + else: + ax.text(0.5, 0.5, 'No |B| spectrum data available', + ha='center', va='center', transform=ax.transAxes) + + # Plot 4: G and I profiles + ax = axes[1, 1] + ax.plot(booz_data['s'], booz_data['bvco'], 'b-', label='G (toroidal)', linewidth=2) + ax.plot(booz_data['s'], booz_data['buco'], 'r-', label='I (poloidal)', linewidth=2) + ax.set_xlabel('s') + ax.set_ylabel('Boozer covariant components') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Boozer I and G') + + plt.tight_layout() + return fig + + +def main(): + """Main comparison routine""" + print("=== Boozer Coordinate Comparison ===") + + # Check files exist + if not os.path.exists('wout.nc'): + print("Error: wout.nc not found. Run setup_booz_xform_test.sh first.") + return + + if not os.path.exists('boozmn_LandremanPaul2021_QA_lowres.nc'): + print("Error: booz_xform file not found. Run setup_booz_xform_test.sh first.") + return + + # Load booz_xform data + print("\n1. Loading booz_xform data...") + booz_data = read_booz_xform_file('boozmn_LandremanPaul2021_QA_lowres.nc') + + # Initialize SIMPLE + print("\n2. Initializing SIMPLE with internal Boozer conversion...") + tracer = setup_simple_boozer('wout.nc') + + # Compare coordinate transformations + print("\n3. Computing coordinate transformations...") + s_test, theta_vmec, theta_boozer_simple, zeta_boozer_simple = compare_coordinate_transform() + + # Evaluate field quantities + print("\n4. Evaluating field quantities...") + # Test at mid-radius, various poloidal angles + s_eval = 0.5 * np.ones(20) + theta_eval = np.linspace(0, 2*np.pi, 20) + zeta_eval = np.zeros(20) + + simple_results = evaluate_simple_boozer(s_eval, theta_eval, zeta_eval) + + # Create comparison plots + print("\n5. Creating comparison plots...") + fig = plot_comparison(booz_data, s_test, theta_vmec, theta_boozer_simple) + + # Save figure + output_file = 'boozer_comparison.png' + fig.savefig(output_file, dpi=150, bbox_inches='tight') + print(f"\nPlots saved to: {output_file}") + + # Print some numerical comparisons + print("\n=== Numerical Comparison at s=0.5 ===") + print(f"SIMPLE |B| range: [{simple_results['modB'].min():.4f}, {simple_results['modB'].max():.4f}]") + print(f"SIMPLE sqrt(g): {simple_results['sqrtg'][0]:.6e}") + + # plt.show() # Don't show in test mode + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test/booz_xform/compare_field_evaluation.py b/test/booz_xform/compare_field_evaluation.py new file mode 100644 index 00000000..0e0a2347 --- /dev/null +++ b/test/booz_xform/compare_field_evaluation.py @@ -0,0 +1,309 @@ +#!/usr/bin/env python3 +""" +Compare magnetic field evaluation between SIMPLE's internal Boozer conversion +and BOOZXFORM external file at the same coordinates. + +This script: +1. Runs SIMPLE with isw_field_type=2 (internal Boozer) and saves field data +2. Runs SIMPLE with isw_field_type=5 (BOOZXFORM) and saves field data +3. Compares the field evaluations at the same coordinates +""" +import numpy as np +import matplotlib.pyplot as plt +import netCDF4 as nc +import sys +import os +import subprocess +import tempfile + +# Add parent directory to path for pysimple +script_dir = os.path.dirname(os.path.abspath(__file__)) +if 'build' in script_dir: + # We're running from build directory + sys.path.insert(0, os.getcwd()) # Current dir should have pysimple +else: + # We're running from source directory + build_dir = os.path.join(script_dir, '../../build') + if os.path.exists(build_dir): + sys.path.insert(0, build_dir) + +try: + import pysimple +except ImportError as e: + print(f"Error: pysimple not found. Please build SIMPLE first.") + print(f"Python path: {sys.path}") + print(f"Import error: {e}") + sys.exit(1) + + +def create_test_config(field_type, boozxform_file=None): + """Create simple.in configuration for given field type""" + config = f"""&config + isw_field_type = {field_type} + netcdffile = 'wout.nc' + ntestpart = 1 + nper = 1 + trace_time = 1d-4 + ntimstep = 10 + generate_start_only = .true. + startmode = 5 + num_surf = 1 + sbeg(1) = 0.5 +""" + if boozxform_file: + config += f" boozxform_file = '{boozxform_file}'\n" + + config += "/\n" + return config + + +def evaluate_field_at_coordinates(field_type, coords, boozxform_file=None): + """Evaluate magnetic field at given coordinates using specified field type""" + + # Create temporary configuration + config_content = create_test_config(field_type, boozxform_file) + + with tempfile.NamedTemporaryFile(mode='w', suffix='.in', delete=False) as f: + f.write(config_content) + config_file = f.name + + try: + # Read configuration + pysimple.params.read_config(config_file) + + # Initialize tracer + tracer = pysimple.simple.Tracer() + + # Initialize field + pysimple.simple_main.init_field(tracer, 'wout.nc', 5, 5, 5, 1) + + # Evaluate field at coordinates + npoints = len(coords) + results = { + 'modB': np.zeros(npoints), + 'sqrtg': np.zeros(npoints), + 'dBds': np.zeros(npoints), + 'dBdtheta': np.zeros(npoints), + 'dBdzeta': np.zeros(npoints), + 'Bs_cov': np.zeros(npoints), + 'Btheta_cov': np.zeros(npoints), + 'Bzeta_cov': np.zeros(npoints), + } + + # Temporary arrays for field evaluation + bder = np.empty(3) + hcovar = np.empty(3) + hctrvr = np.empty(3) + hcurl = np.empty(3) + + for i, (s, theta, zeta) in enumerate(coords): + x = np.array([s, theta, zeta]) + + try: + # Evaluate field using magfie routine + if field_type == 2: # Boozer + modB, sqrtg = pysimple.magfie_sub.magfie_boozer(x, bder, hcovar, hctrvr, hcurl) + elif field_type == 5: # BOOZXFORM + modB, sqrtg = pysimple.magfie_sub.magfie_boozxform(x, bder, hcovar, hctrvr, hcurl) + else: + raise ValueError(f"Unsupported field type: {field_type}") + + results['modB'][i] = modB + results['sqrtg'][i] = sqrtg + results['dBds'][i] = bder[0] * modB + results['dBdtheta'][i] = bder[1] * modB + results['dBdzeta'][i] = bder[2] * modB + results['Bs_cov'][i] = hcovar[0] * modB + results['Btheta_cov'][i] = hcovar[1] * modB + results['Bzeta_cov'][i] = hcovar[2] * modB + + except Exception as e: + print(f"Error evaluating field at {x}: {e}") + # Fill with NaN on error + for key in results: + results[key][i] = np.nan + + return results + + finally: + # Clean up temporary file + os.unlink(config_file) + + +def plot_field_comparison(coords, boozer_results, boozxform_results): + """Create comparison plots of field evaluation results""" + + s_vals = [c[0] for c in coords] + theta_vals = [c[1] for c in coords] + + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + fig.suptitle('SIMPLE Boozer vs BOOZXFORM Field Evaluation Comparison', fontsize=14) + + # Plot 1: |B| comparison + ax = axes[0, 0] + ax.plot(theta_vals, boozer_results['modB'], 'b-', label='Internal Boozer', linewidth=2) + ax.plot(theta_vals, boozxform_results['modB'], 'r--', label='BOOZXFORM', linewidth=2) + ax.set_xlabel('θ_Boozer') + ax.set_ylabel('|B|') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Magnetic Field Strength') + + # Plot 2: sqrt(g) comparison + ax = axes[0, 1] + ax.plot(theta_vals, boozer_results['sqrtg'], 'b-', label='Internal Boozer', linewidth=2) + ax.plot(theta_vals, boozxform_results['sqrtg'], 'r--', label='BOOZXFORM', linewidth=2) + ax.set_xlabel('θ_Boozer') + ax.set_ylabel('√g') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Jacobian') + + # Plot 3: B derivatives + ax = axes[0, 2] + ax.plot(theta_vals, boozer_results['dBds'], 'b-', label='∂B/∂s (Boozer)', linewidth=2) + ax.plot(theta_vals, boozxform_results['dBds'], 'r--', label='∂B/∂s (BOOZXFORM)', linewidth=2) + ax.set_xlabel('θ_Boozer') + ax.set_ylabel('∂B/∂s') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Field Derivative') + + # Plot 4: B covariant components + ax = axes[1, 0] + ax.plot(theta_vals, boozer_results['Btheta_cov'], 'b-', label='B_θ (Boozer)', linewidth=2) + ax.plot(theta_vals, boozxform_results['Btheta_cov'], 'r--', label='B_θ (BOOZXFORM)', linewidth=2) + ax.set_xlabel('θ_Boozer') + ax.set_ylabel('B_θ') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Poloidal Field Component') + + # Plot 5: Relative differences + ax = axes[1, 1] + modB_diff = np.abs(boozer_results['modB'] - boozxform_results['modB']) / np.abs(boozer_results['modB']) + sqrtg_diff = np.abs(boozer_results['sqrtg'] - boozxform_results['sqrtg']) / np.abs(boozer_results['sqrtg']) + + ax.semilogy(theta_vals, modB_diff, 'g-', label='|B| rel. diff', linewidth=2) + ax.semilogy(theta_vals, sqrtg_diff, 'm-', label='√g rel. diff', linewidth=2) + ax.set_xlabel('θ_Boozer') + ax.set_ylabel('Relative Difference') + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title('Relative Differences') + + # Plot 6: Summary statistics + ax = axes[1, 2] + ax.axis('off') + + # Calculate statistics + modB_mean_diff = np.nanmean(modB_diff) + modB_max_diff = np.nanmax(modB_diff) + sqrtg_mean_diff = np.nanmean(sqrtg_diff) + sqrtg_max_diff = np.nanmax(sqrtg_diff) + + valid_boozer = np.sum(~np.isnan(boozer_results['modB'])) + valid_boozxform = np.sum(~np.isnan(boozxform_results['modB'])) + + stats_text = f"""Comparison Statistics: + +Valid Evaluations: + Internal Boozer: {valid_boozer}/{len(coords)} + BOOZXFORM: {valid_boozxform}/{len(coords)} + +|B| Differences: + Mean: {modB_mean_diff:.2e} + Max: {modB_max_diff:.2e} + +√g Differences: + Mean: {sqrtg_mean_diff:.2e} + Max: {sqrtg_max_diff:.2e} + +|B| Range: + Boozer: [{np.nanmin(boozer_results['modB']):.1f}, {np.nanmax(boozer_results['modB']):.1f}] + BOOZXFORM: [{np.nanmin(boozxform_results['modB']):.1f}, {np.nanmax(boozxform_results['modB']):.1f}] +""" + + ax.text(0.1, 0.9, stats_text, transform=ax.transAxes, fontsize=10, + verticalalignment='top', fontfamily='monospace') + + plt.tight_layout() + return fig + + +def main(): + """Main comparison routine""" + print("=== SIMPLE Field Evaluation Comparison ===") + print("Comparing internal Boozer vs BOOZXFORM field evaluation") + + # Check files exist + if not os.path.exists('wout.nc'): + print("Error: wout.nc not found. Run setup_booz_xform_test.sh first.") + return + + if not os.path.exists('boozmn_LandremanPaul2021_QA_lowres.nc'): + print("Error: booz_xform file not found. Run setup_booz_xform_test.sh first.") + return + + # Define test coordinates in Boozer space + # Test at s=0.5, various poloidal angles, zeta=0 + s_test = 0.5 + theta_test = np.linspace(0, 2*np.pi, 20) + zeta_test = 0.0 + + coords = [(s_test, theta, zeta_test) for theta in theta_test] + + print(f"\nEvaluating field at {len(coords)} points:") + print(f" s = {s_test}") + print(f" θ ∈ [0, 2π] ({len(theta_test)} points)") + print(f" ζ = {zeta_test}") + + # Evaluate with internal Boozer + print("\n1. Evaluating with internal Boozer conversion (isw_field_type=2)...") + try: + boozer_results = evaluate_field_at_coordinates(2, coords) + print(f" Internal Boozer: {np.sum(~np.isnan(boozer_results['modB']))}/{len(coords)} successful evaluations") + except Exception as e: + print(f" Error with internal Boozer: {e}") + boozer_results = None + + # Evaluate with BOOZXFORM + print("\n2. Evaluating with BOOZXFORM file (isw_field_type=5)...") + try: + boozxform_results = evaluate_field_at_coordinates(5, coords, 'boozmn_LandremanPaul2021_QA_lowres.nc') + print(f" BOOZXFORM: {np.sum(~np.isnan(boozxform_results['modB']))}/{len(coords)} successful evaluations") + except Exception as e: + print(f" Error with BOOZXFORM: {e}") + boozxform_results = None + + # Create comparison plots + if boozer_results is not None and boozxform_results is not None: + print("\n3. Creating comparison plots...") + fig = plot_field_comparison(coords, boozer_results, boozxform_results) + + # Save figure + output_file = 'field_evaluation_comparison.png' + fig.savefig(output_file, dpi=150, bbox_inches='tight') + print(f" Plots saved to: {output_file}") + + # Print summary + print("\n=== Comparison Summary ===") + if np.any(~np.isnan(boozer_results['modB'])) and np.any(~np.isnan(boozxform_results['modB'])): + valid_mask = ~np.isnan(boozer_results['modB']) & ~np.isnan(boozxform_results['modB']) + if np.any(valid_mask): + modB_rel_diff = np.abs(boozer_results['modB'][valid_mask] - boozxform_results['modB'][valid_mask]) / np.abs(boozer_results['modB'][valid_mask]) + print(f"Field strength relative difference: mean={np.mean(modB_rel_diff):.2e}, max={np.max(modB_rel_diff):.2e}") + else: + print("No valid comparisons available") + else: + print("Insufficient data for meaningful comparison") + else: + print("\n3. Cannot create comparison plots - field evaluation failed") + if boozer_results is None: + print(" Internal Boozer evaluation failed") + if boozxform_results is None: + print(" BOOZXFORM evaluation failed") + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test/booz_xform/run_comparison_test.sh b/test/booz_xform/run_comparison_test.sh new file mode 100755 index 00000000..8d470642 --- /dev/null +++ b/test/booz_xform/run_comparison_test.sh @@ -0,0 +1,54 @@ +#!/bin/bash +# Run booz_xform comparison test in build directory +# This script is called by CTest + +set -e # Exit on error + +echo "=== Running booz_xform comparison test ===" +echo "Working directory: $(pwd)" + +# Download test files if not present +BOOZ_FILE="boozmn_LandremanPaul2021_QA_lowres.nc" +VMEC_FILE="wout.nc" + +if [ ! -f "$BOOZ_FILE" ]; then + echo "Downloading booz_xform test file..." + wget -q https://github.com/itpplasma/benchmark_orbit/raw/refs/heads/main/booz_xform/boozmn_LandremanPaul2021_QA_lowres.nc +fi + +if [ ! -f "$VMEC_FILE" ]; then + echo "Downloading VMEC test file..." + wget -q https://github.com/hiddenSymmetries/simsopt/raw/master/tests/test_files/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc -O wout.nc +fi + +# Create simple.in for testing +cat > simple.in << EOF +&config +trace_time = 1d-2 +sbeg(1) = 0.5d0 +ntestpart = 1 +netcdffile = 'wout.nc' +isw_field_type = 2 ! Internal Boozer +/ +EOF + +# Get the source directory for the Python script +if [ -z "$CMAKE_CURRENT_SOURCE_DIR" ]; then + # If not set, assume we're in build/test/tests + SCRIPT_DIR="../../../test/booz_xform" +else + SCRIPT_DIR="${CMAKE_CURRENT_SOURCE_DIR}/../booz_xform" +fi + +# Run the Python comparison script +echo "Running comparison..." +python3 ${SCRIPT_DIR}/compare_boozer_coords.py + +# Check if output was created +if [ -f "boozer_comparison.png" ]; then + echo "Test completed successfully. Output: boozer_comparison.png" + exit 0 +else + echo "Error: Expected output file not created" + exit 1 +fi \ No newline at end of file diff --git a/test/booz_xform/setup_booz_xform_test.sh b/test/booz_xform/setup_booz_xform_test.sh new file mode 100755 index 00000000..6debe37b --- /dev/null +++ b/test/booz_xform/setup_booz_xform_test.sh @@ -0,0 +1,50 @@ +#!/bin/bash +# Setup script for booz_xform testing +# Downloads test data and prepares comparison environment + +set -e # Exit on error + +echo "=== Setting up booz_xform test environment ===" + +# Create test directory +TEST_DIR="test/booz_xform" +mkdir -p $TEST_DIR +cd $TEST_DIR + +# Download booz_xform file if not present +BOOZ_FILE="boozmn_LandremanPaul2021_QA_lowres.nc" +if [ ! -f "$BOOZ_FILE" ]; then + echo "Downloading booz_xform test file..." + wget https://github.com/itpplasma/benchmark_orbit/raw/refs/heads/main/booz_xform/boozmn_LandremanPaul2021_QA_lowres.nc +else + echo "Booz_xform file already exists: $BOOZ_FILE" +fi + +# Download corresponding VMEC file if not present +VMEC_FILE="wout.nc" +if [ ! -f "$VMEC_FILE" ]; then + echo "Downloading corresponding VMEC file..." + wget https://github.com/hiddenSymmetries/simsopt/raw/master/tests/test_files/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc -O wout.nc +else + echo "VMEC file already exists: $VMEC_FILE" +fi + +# Create simple.in for testing +cat > simple.in << EOF +&config +trace_time = 1d-2 ! slowing down time, s +sbeg = 0.5d0 ! starting s (normalized toroidal flux) for particles +send = 0.5d0 ! ending s for uniform sampling +ntestpart = 1 ! number of test particles (1 for coordinate comparison) +netcdffile = 'wout.nc' ! VMEC file +booz_file = 'boozmn_LandremanPaul2021_QA_lowres.nc' ! booz_xform file +isw_field_type = 2 ! 2: Internal Boozer conversion (will add 5 for booz_xform) +startmode = 1 ! 1: uniform in s +/ +EOF + +echo "Setup complete!" +echo "Files created in: $(pwd)" +echo " - $BOOZ_FILE (booz_xform output)" +echo " - $VMEC_FILE (VMEC input)" +echo " - simple.in (test configuration)" \ No newline at end of file diff --git a/test/booz_xform/simple.in b/test/booz_xform/simple.in new file mode 100644 index 00000000..5c5a7f35 --- /dev/null +++ b/test/booz_xform/simple.in @@ -0,0 +1,9 @@ +&config +trace_time = 1d-4 +sbeg = 0.5d0 +ntestpart = 1 +netcdffile = 'wout.nc' +isw_field_type = 5 +boozxform_file = 'boozmn_LandremanPaul2021_QA_lowres.nc' +/ + diff --git a/test/booz_xform/test_boozxform_simple.in b/test/booz_xform/test_boozxform_simple.in new file mode 100644 index 00000000..c5186725 --- /dev/null +++ b/test/booz_xform/test_boozxform_simple.in @@ -0,0 +1,12 @@ +&config +isw_field_type = 5 +netcdffile = 'test/booz_xform/wout.nc' +boozxform_file = 'test/booz_xform/boozmn_LandremanPaul2021_QA_lowres.nc' +ntestpart = 1 +trace_time = 1d-4 +ntimstep = 10 +generate_start_only = .true. +startmode = 5 +num_surf = 1 +sbeg(1) = 0.5 +/ \ No newline at end of file diff --git a/test/booz_xform/test_improved_evaluation.py b/test/booz_xform/test_improved_evaluation.py new file mode 100644 index 00000000..c109aee7 --- /dev/null +++ b/test/booz_xform/test_improved_evaluation.py @@ -0,0 +1,196 @@ +#!/usr/bin/env python3 +""" +Test the improved BOOZXFORM field evaluation implementation. +This script tests that magfie_boozxform now reads real data and returns +meaningful field values instead of dummy placeholders. +""" +import numpy as np +import sys +import os + +# Add parent directory to path for pysimple +script_dir = os.path.dirname(os.path.abspath(__file__)) +build_dir = os.path.join(script_dir, '../../build') +if os.path.exists(build_dir): + sys.path.insert(0, build_dir) + +try: + import pysimple +except ImportError as e: + print(f"Error: pysimple not found: {e}") + sys.exit(1) + + +def test_boozxform_field_evaluation(): + """Test BOOZXFORM field evaluation using the magfie interface""" + print("Testing improved BOOZXFORM field evaluation...") + + # Test coordinates (s=0.5, theta=0, zeta=0) + x = np.array([0.5, 0.0, 0.0]) + + # Temporary arrays for field evaluation + bder = np.empty(3) + hcovar = np.empty(3) + hctrvr = np.empty(3) + hcurl = np.empty(3) + + print(f"Evaluating field at coordinates: s={x[0]}, θ={x[1]}, ζ={x[2]}") + + try: + # Call magfie_boozxform directly to bypass init_magfie issue + modB, sqrtg = pysimple.magfie_sub.magfie_boozxform(x, bder, hcovar, hctrvr, hcurl) + + print(f"\nField evaluation results:") + print(f" |B| = {modB:.6f}") + print(f" √g = {sqrtg:.6e}") + print(f" ∂B/∂s = {bder[0]:.6e}") + print(f" ∂B/∂θ = {bder[1]:.6e}") + print(f" ∂B/∂ζ = {bder[2]:.6e}") + print(f" B_s = {hcovar[0]:.6e}") + print(f" B_θ = {hcovar[1]:.6e}") + print(f" B_ζ = {hcovar[2]:.6e}") + print(f" B^s = {hctrvr[0]:.6e}") + print(f" B^θ = {hctrvr[1]:.6e}") + print(f" B^ζ = {hctrvr[2]:.6e}") + + # Check if we're getting realistic values + realistic = True + issues = [] + + if modB == 1.0: + issues.append("Field strength still placeholder (exactly 1.0)") + realistic = False + + if sqrtg == 1.0: + issues.append("Jacobian still placeholder (exactly 1.0)") + realistic = False + + if np.allclose(bder, 0.0): + issues.append("All field derivatives are zero") + + if np.allclose(hcurl, 0.0): + issues.append("All curl components are zero") + + if sqrtg <= 0: + issues.append("Non-positive Jacobian") + realistic = False + + if realistic: + print(f"\n✅ SUCCESS: Field evaluation appears to use real BOOZXFORM data") + print(f" - Non-trivial field strength: {modB}") + print(f" - Non-trivial Jacobian: {sqrtg:.2e}") + if not np.allclose(hcovar[1:], 0.0): + print(f" - Non-zero field components: B_θ={hcovar[1]:.3f}, B_ζ={hcovar[2]:.3f}") + else: + print(f"\n⚠️ PARTIAL: Field evaluation working but still has placeholder elements:") + for issue in issues: + print(f" - {issue}") + + return modB, sqrtg, realistic + + except Exception as e: + print(f"\n❌ ERROR: Field evaluation failed: {e}") + return None, None, False + + +def test_coordinate_variations(): + """Test field evaluation at multiple coordinates""" + print(f"\n" + "="*50) + print("Testing field evaluation at multiple coordinates...") + + # Test at several points + test_coords = [ + (0.1, 0.0, 0.0), # Near axis + (0.5, 0.0, 0.0), # Mid-radius + (0.9, 0.0, 0.0), # Near edge + (0.5, np.pi/2, 0.0), # Different theta + (0.5, np.pi, 0.0), # Different theta + ] + + results = [] + + for i, (s, theta, zeta) in enumerate(test_coords): + x = np.array([s, theta, zeta]) + bder = np.empty(3) + hcovar = np.empty(3) + hctrvr = np.empty(3) + hcurl = np.empty(3) + + try: + modB, sqrtg = pysimple.magfie_sub.magfie_boozxform(x, bder, hcovar, hctrvr, hcurl) + results.append((s, theta, zeta, modB, sqrtg, hcovar[1], hcovar[2])) + print(f" Point {i+1}: s={s:.1f}, θ={theta:.2f} → |B|={modB:.4f}, √g={sqrtg:.2e}") + + except Exception as e: + print(f" Point {i+1}: ERROR - {e}") + results.append((s, theta, zeta, None, None, None, None)) + + # Analyze results + valid_results = [r for r in results if r[3] is not None] + + if len(valid_results) > 1: + modB_values = [r[3] for r in valid_results] + sqrtg_values = [r[4] for r in valid_results] + + modB_variation = np.std(modB_values) / np.mean(modB_values) if np.mean(modB_values) > 0 else 0 + sqrtg_variation = np.std(sqrtg_values) / np.mean(sqrtg_values) if np.mean(sqrtg_values) > 0 else 0 + + print(f"\nVariation analysis:") + print(f" |B| coefficient of variation: {modB_variation:.1%}") + print(f" √g coefficient of variation: {sqrtg_variation:.1%}") + + if modB_variation > 0.01: # >1% variation + print(f" ✅ Field strength shows spatial variation (good)") + else: + print(f" ⚠️ Field strength shows little variation (may still be placeholder)") + + if sqrtg_variation > 0.01: # >1% variation + print(f" ✅ Jacobian shows spatial variation (good)") + else: + print(f" ⚠️ Jacobian shows little variation (may still be placeholder)") + + return valid_results + + +def main(): + print("="*60) + print("SIMPLE BOOZXFORM Field Evaluation Test") + print("="*60) + + # Check that test data exists + if not os.path.exists('boozmn_LandremanPaul2021_QA_lowres.nc'): + print("❌ ERROR: BOOZXFORM test file not found") + print("Run setup_booz_xform_test.sh first") + return + + # Test single point evaluation + modB, sqrtg, realistic = test_boozxform_field_evaluation() + + if modB is not None: + # Test coordinate variations + results = test_coordinate_variations() + + print(f"\n" + "="*60) + print("SUMMARY") + print("="*60) + + if realistic: + print("✅ BOOZXFORM field evaluation is working with real data") + else: + print("⚠️ BOOZXFORM field evaluation partially working but needs improvement") + + print(f"Evaluated {len(results)} coordinate points successfully") + + if len(results) > 0: + modB_range = [r[3] for r in results if r[3] is not None] + if modB_range: + print(f"|B| range: [{min(modB_range):.4f}, {max(modB_range):.4f}]") + + else: + print("❌ BOOZXFORM field evaluation failed completely") + + print("="*60) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 9025a5e8..84e2675c 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -199,6 +199,22 @@ add_test(NAME test_coordinates_simple COMMAND test_coordinates_simple.x) # - test_canonical_full_regression.f90 # - test_canonical_baseline_simple.f90 +add_executable (test_booz_xform.x test_booz_xform.f90) +target_link_libraries(test_booz_xform.x simple) +add_test(NAME test_booz_xform COMMAND test_booz_xform.x) + +# Add booz_xform comparison test that downloads data and runs Python comparison +add_test( + NAME test_booz_xform_comparison + COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${CMAKE_BINARY_DIR}:$ENV{PYTHONPATH} + ${CMAKE_CURRENT_SOURCE_DIR}/../booz_xform/run_comparison_test.sh + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) +set_tests_properties(test_booz_xform_comparison PROPERTIES + LABELS "booz_xform;comparison" + TIMEOUT 300 # 5 minutes for download and comparison +) + add_subdirectory(field_can) add_subdirectory(magfie) add_subdirectory(orbit) diff --git a/test/tests/test_booz_xform.f90 b/test/tests/test_booz_xform.f90 new file mode 100644 index 00000000..66654d4b --- /dev/null +++ b/test/tests/test_booz_xform.f90 @@ -0,0 +1,9 @@ +program test_booz_xform + call main + +contains + + subroutine main + + end subroutine main +end program test_booz_xform diff --git a/test_boozxform/plot_existing_results.py b/test_boozxform/plot_existing_results.py new file mode 100644 index 00000000..0502e95e --- /dev/null +++ b/test_boozxform/plot_existing_results.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 +""" +Plot results from existing test runs comparing VMEC internal Boozer vs BOOZXFORM +""" +import numpy as np +import matplotlib.pyplot as plt +import os + +def load_confined_fraction(filename): + """Load confined fraction data""" + if os.path.exists(filename): + data = np.loadtxt(filename) + return data + return None + +def plot_comparison(): + """Plot comparison of existing test results""" + + fig, axes = plt.subplots(2, 2, figsize=(12, 10)) + fig.suptitle('SIMPLE Field Type Comparison: Internal Boozer vs BOOZXFORM', fontsize=14) + + # Load data from test directories + test_dirs = { + 'Internal Boozer': '/Users/ert/code/SIMPLE/test/booz_xform', + 'Test Run': '/Users/ert/code/SIMPLE/test_boozxform' + } + + colors = {'Internal Boozer': 'blue', 'Test Run': 'red'} + + # Plot 1: Confined fraction vs time + ax = axes[0, 0] + for label, test_dir in test_dirs.items(): + conf_file = os.path.join(test_dir, 'confined_fraction.dat') + if os.path.exists(conf_file): + data = load_confined_fraction(conf_file) + if data is not None and len(data) > 0: + ax.plot(data[:, 0], data[:, 1], color=colors[label], + label=label, linewidth=2) + + ax.set_xlabel('Time') + ax.set_ylabel('Confined Fraction') + ax.set_title('Confinement Evolution') + ax.legend() + ax.grid(True, alpha=0.3) + + # Plot 2: Loss times + ax = axes[0, 1] + for label, test_dir in test_dirs.items(): + times_file = os.path.join(test_dir, 'times_lost.dat') + if os.path.exists(times_file): + try: + times = np.loadtxt(times_file) + if times.size > 0: + if times.ndim == 1: + # Single particle + ax.scatter([0], [times[0]], color=colors[label], + label=label, s=50) + else: + # Multiple particles + ax.scatter(range(len(times[:, 0])), times[:, 0], + color=colors[label], label=label, alpha=0.6, s=30) + except: + pass + + ax.set_xlabel('Particle Index') + ax.set_ylabel('Loss Time') + ax.set_title('Particle Loss Times') + ax.legend() + ax.grid(True, alpha=0.3) + + # Plot 3: Initial conditions + ax = axes[1, 0] + for label, test_dir in test_dirs.items(): + start_file = os.path.join(test_dir, 'start.dat') + if os.path.exists(start_file): + try: + start = np.loadtxt(start_file) + if start.size > 0: + if start.ndim == 1: + # Single particle: s, R, Z, par_inv, perp_inv + ax.scatter([start[1]], [start[2]], color=colors[label], + label=f"{label} (s={start[0]:.2f})", s=100, marker='o') + else: + # Multiple particles + ax.scatter(start[:, 1], start[:, 2], color=colors[label], + label=label, alpha=0.6, s=30) + except Exception as e: + print(f"Error reading {start_file}: {e}") + + ax.set_xlabel('R [cm]') + ax.set_ylabel('Z [cm]') + ax.set_title('Initial Particle Positions') + ax.legend() + ax.grid(True, alpha=0.3) + ax.axis('equal') + + # Plot 4: Summary + ax = axes[1, 1] + ax.axis('off') + + summary_text = "Test Summary:\n\n" + + # Check what files exist + for label, test_dir in test_dirs.items(): + summary_text += f"{label}:\n" + + # Check simple.in + simple_in = os.path.join(test_dir, 'simple.in') + if os.path.exists(simple_in): + with open(simple_in, 'r') as f: + lines = f.readlines() + for line in lines: + if 'isw_field_type' in line: + summary_text += f" Field type: {line.strip()}\n" + if 'boozxform_file' in line: + summary_text += f" BOOZXFORM: {line.strip()}\n" + + # Check output files + for fname in ['confined_fraction.dat', 'times_lost.dat', 'start.dat']: + fpath = os.path.join(test_dir, fname) + if os.path.exists(fpath): + fsize = os.path.getsize(fpath) + summary_text += f" {fname}: {fsize} bytes\n" + + summary_text += "\n" + + ax.text(0.05, 0.95, summary_text, transform=ax.transAxes, + fontsize=9, verticalalignment='top', fontfamily='monospace') + + plt.tight_layout() + return fig + +def main(): + print("=== Plotting Existing SIMPLE Test Results ===") + + # Create comparison plots + fig = plot_comparison() + + # Save figure + output_file = 'field_type_comparison.png' + fig.savefig(output_file, dpi=150, bbox_inches='tight') + print(f"\nSaved plot to {output_file}") + + # Also display + plt.show() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test_boozxform/plot_orbit_comparison.py b/test_boozxform/plot_orbit_comparison.py new file mode 100755 index 00000000..43960915 --- /dev/null +++ b/test_boozxform/plot_orbit_comparison.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +""" +Simple script to compare orbits from VMEC vs BOOZXFORM field types +by running SIMPLE twice and plotting the results +""" +import numpy as np +import matplotlib.pyplot as plt +import subprocess +import os + +def run_simple_with_field_type(field_type, label): + """Run SIMPLE with specified field type and save orbit data""" + + # Create input file + config = f"""&config +ntimstep = 1000 +trace_time = 1d-4 +sbeg = 0.5d0 +ntestpart = 1 +netcdffile = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' +boozxform_file = 'boozmn_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' +isw_field_type = {field_type} +/""" + + with open('simple.in', 'w') as f: + f.write(config) + + print(f"\nRunning SIMPLE with {label} (isw_field_type={field_type})...") + + # Run SIMPLE + result = subprocess.run(['../build/simple.x'], + capture_output=True, text=True) + + if result.returncode != 0: + print(f"Error running SIMPLE with {label}:") + print(result.stderr) + return None + + # Read orbit data from fort.601 (or appropriate output file) + orbit_file = 'fort.601' # Poincare plot data + if os.path.exists(orbit_file): + try: + data = np.loadtxt(orbit_file) + # Save with unique name + output_file = f'orbit_{label.lower().replace(" ", "_")}.dat' + np.savetxt(output_file, data) + print(f" Saved orbit data to {output_file}") + return data + except: + print(f" Could not read orbit data from {orbit_file}") + return None + else: + print(f" No orbit output file {orbit_file} found") + # Try to read from other possible output files + for fortfile in ['fort.1001', 'fort.10001']: + if os.path.exists(fortfile): + try: + data = np.loadtxt(fortfile) + output_file = f'orbit_{label.lower().replace(" ", "_")}.dat' + np.savetxt(output_file, data) + print(f" Found orbit data in {fortfile}, saved to {output_file}") + return data + except: + pass + return None + +def plot_orbit_comparison(vmec_data, booz_data): + """Create comparison plots of orbit data""" + + fig, axes = plt.subplots(2, 2, figsize=(12, 10)) + fig.suptitle('Orbit Comparison: VMEC vs BOOZXFORM Field Types', fontsize=14) + + # If we have Poincare data (R, Z at fixed toroidal angle) + if vmec_data is not None and vmec_data.shape[1] >= 2: + ax = axes[0, 0] + if vmec_data.shape[1] >= 3: + # Assume columns are: time/index, R, Z + ax.plot(vmec_data[:, 1], vmec_data[:, 2], 'b.', + markersize=3, label='VMEC', alpha=0.6) + else: + # Just two columns: R, Z + ax.plot(vmec_data[:, 0], vmec_data[:, 1], 'b.', + markersize=3, label='VMEC', alpha=0.6) + ax.set_xlabel('R [cm]') + ax.set_ylabel('Z [cm]') + ax.set_title('Poincaré Section - VMEC') + ax.axis('equal') + ax.grid(True, alpha=0.3) + + if booz_data is not None and booz_data.shape[1] >= 2: + ax = axes[0, 1] + if booz_data.shape[1] >= 3: + ax.plot(booz_data[:, 1], booz_data[:, 2], 'r.', + markersize=3, label='BOOZXFORM', alpha=0.6) + else: + ax.plot(booz_data[:, 0], booz_data[:, 1], 'r.', + markersize=3, label='BOOZXFORM', alpha=0.6) + ax.set_xlabel('R [cm]') + ax.set_ylabel('Z [cm]') + ax.set_title('Poincaré Section - BOOZXFORM') + ax.axis('equal') + ax.grid(True, alpha=0.3) + + # Combined plot + if vmec_data is not None and booz_data is not None: + ax = axes[1, 0] + if vmec_data.shape[1] >= 3 and booz_data.shape[1] >= 3: + ax.plot(vmec_data[:, 1], vmec_data[:, 2], 'b.', + markersize=3, label='VMEC', alpha=0.5) + ax.plot(booz_data[:, 1], booz_data[:, 2], 'r.', + markersize=3, label='BOOZXFORM', alpha=0.5) + elif vmec_data.shape[1] >= 2 and booz_data.shape[1] >= 2: + ax.plot(vmec_data[:, 0], vmec_data[:, 1], 'b.', + markersize=3, label='VMEC', alpha=0.5) + ax.plot(booz_data[:, 0], booz_data[:, 1], 'r.', + markersize=3, label='BOOZXFORM', alpha=0.5) + ax.set_xlabel('R [cm]') + ax.set_ylabel('Z [cm]') + ax.set_title('Overlay Comparison') + ax.axis('equal') + ax.legend() + ax.grid(True, alpha=0.3) + + # Statistics + ax = axes[1, 1] + ax.axis('off') + + stats_text = "Orbit Statistics:\n\n" + + if vmec_data is not None: + stats_text += f"VMEC field type:\n" + stats_text += f" Points: {len(vmec_data)}\n" + if vmec_data.shape[1] >= 2: + col_r = 1 if vmec_data.shape[1] >= 3 else 0 + col_z = 2 if vmec_data.shape[1] >= 3 else 1 + stats_text += f" R range: [{vmec_data[:, col_r].min():.1f}, {vmec_data[:, col_r].max():.1f}]\n" + stats_text += f" Z range: [{vmec_data[:, col_z].min():.1f}, {vmec_data[:, col_z].max():.1f}]\n" + else: + stats_text += "VMEC: No data\n" + + stats_text += "\n" + + if booz_data is not None: + stats_text += f"BOOZXFORM field type:\n" + stats_text += f" Points: {len(booz_data)}\n" + if booz_data.shape[1] >= 2: + col_r = 1 if booz_data.shape[1] >= 3 else 0 + col_z = 2 if booz_data.shape[1] >= 3 else 1 + stats_text += f" R range: [{booz_data[:, col_r].min():.1f}, {booz_data[:, col_r].max():.1f}]\n" + stats_text += f" Z range: [{booz_data[:, col_z].min():.1f}, {booz_data[:, col_z].max():.1f}]\n" + else: + stats_text += "BOOZXFORM: No data\n" + + ax.text(0.1, 0.9, stats_text, transform=ax.transAxes, + fontsize=10, verticalalignment='top', fontfamily='monospace') + + plt.tight_layout() + return fig + +def main(): + print("=== SIMPLE Orbit Comparison: VMEC vs BOOZXFORM ===") + + # Check that files exist + if not os.path.exists('wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc'): + print("Error: VMEC file not found") + return + + if not os.path.exists('boozmn_LandremanPaul2021_QA_reactorScale_lowres_reference.nc'): + print("Error: BOOZXFORM file not found") + return + + # Run with VMEC field type + vmec_data = run_simple_with_field_type(1, "VMEC") + + # Run with BOOZXFORM field type + booz_data = run_simple_with_field_type(5, "BOOZXFORM") + + # Create comparison plots + if vmec_data is not None or booz_data is not None: + print("\nCreating comparison plots...") + fig = plot_orbit_comparison(vmec_data, booz_data) + + output_file = 'orbit_comparison.png' + fig.savefig(output_file, dpi=150, bbox_inches='tight') + print(f"Saved comparison plot to {output_file}") + + # Also show the plot + plt.show() + else: + print("\nNo orbit data available for plotting") + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test_boozxform/simple.in b/test_boozxform/simple.in new file mode 100644 index 00000000..da309588 --- /dev/null +++ b/test_boozxform/simple.in @@ -0,0 +1,10 @@ +&config +ntimstep = 1000 +trace_time = 1d-4 +sbeg = 0.5d0 +ntestpart = 1 +netcdffile = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' +boozxform_file = 'boozmn_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' +isw_field_type = 5 +integmode = 0 ! Use RK integrator which works with magfie interface +/ diff --git a/test_boozxform/simple_vmec.in b/test_boozxform/simple_vmec.in new file mode 100644 index 00000000..ceaedf5b --- /dev/null +++ b/test_boozxform/simple_vmec.in @@ -0,0 +1,8 @@ +&config +ntimstep = 10001 +trace_time = 1d-2 ! slowing down time, s +sbeg = 0.3d0 ! starting s (normalzed toroidal flux) for particles. +ntestpart = 128 ! number of test particles +netcdffile = 'wout.nc' ! name of VMEC file in NETCDF format +isw_field_type = 2 ! -1: Testing, 0: Canonical, 1: VMEC, 2: Boozer, 3: Meiss, 4: Albert +/ diff --git a/test_boozxform_simple.in b/test_boozxform_simple.in new file mode 100644 index 00000000..333cc9a1 --- /dev/null +++ b/test_boozxform_simple.in @@ -0,0 +1,8 @@ +&config +trace_time = 1d-4 +sbeg = 0.5d0 +ntestpart = 1 +netcdffile = 'wout.nc' +isw_field_type = 5 +boozxform_file = 'boozmn_LandremanPaul2021_QA_lowres.nc' +/ \ No newline at end of file