diff --git a/CMakeLists.txt b/CMakeLists.txt index 7762b8f4..354db704 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,20 @@ enable_language(C Fortran) add_compile_options(-g -fbacktrace) +# Prefer local libneo checkout on sibling path; fall back to configured branch +if (NOT DEFINED libneo_SOURCE_DIR) + set(_local_libneo "${CMAKE_SOURCE_DIR}/../libneo") + if (EXISTS "${_local_libneo}/CMakeLists.txt") + get_filename_component(_libneo_abs "${_local_libneo}" ABSOLUTE) + set(libneo_SOURCE_DIR "${_libneo_abs}") + message(STATUS "Detected local libneo checkout at ${libneo_SOURCE_DIR}") + endif() +endif() + +if (NOT DEFINED LIBNEO_BRANCH) + set(LIBNEO_BRANCH "tokamak" CACHE STRING "libneo branch to fetch when no local checkout is present") +endif() + # Disable executable stack for security (trampolines/closures create execstack) if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") # Apple's linker doesn't support -z,noexecstack (uses different security model) diff --git a/TODO.md b/TODO.md new file mode 100644 index 00000000..635fec47 --- /dev/null +++ b/TODO.md @@ -0,0 +1,420 @@ +# TODO: Analytical GS Field via Geoflux Framework + +## Current Status: Phase 3 (SIMPLE Integration) IN PROGRESS + +### What Works Now βœ… +- **Libneo**: Field-agnostic geoflux coordinates with callback interface +- **Libneo**: Complete test suite (unit + integration tests, all passing) +- **SIMPLE**: Infrastructure ready (params, tokamak.in files) +- **SIMPLE**: Alpha confinement example created and documented +- **SIMPLE**: System test scaffold integrated with CMake + +### What's Needed to Complete πŸ”¨ +- **Update src/field.F90** to recognize and initialize analytical fields +- **Verify** system test passes with analytical field +- **Optional**: Create ripple example and documentation + +--- + +## Overview + +Integrate libneo's analytical Grad-Shafranov equilibrium solver (with TF ripple) into SIMPLE by **populating geoflux data structures directly from analytical GS** (no file I/O needed). + +## Strategy: Analytical GS β†’ Geoflux Framework + +Instead of implementing a new coordinate system, we: +1. Evaluate analytical GS (ψ, B, derivatives) on the geoflux grid +2. Populate geoflux internal arrays/splines directly in memory +3. Set `geoflux_ready = .true.` +4. **Everything else just works** - SIMPLE sees it as geoflux + +## Benefits +- βœ… Reuses existing geoflux coordinate system (flux surfaces, metric, Jacobian) +- βœ… Meiss coordinates work unchanged +- βœ… No new coordinate mapping code needed +- βœ… No file I/O overhead +- βœ… Ripple naturally included in field evaluation + +--- + +## Implementation Status + +### Phase 1: Infrastructure βœ… COMPLETE + +#### 1. βœ… Add ANALYTICAL constant [DONE] +- Added `ANALYTICAL=6` to `src/magfie.f90` + +#### 2. βœ… Create tokamak.in [DONE] +- Created `examples/tokamak/tokamak.in` (no ripple) +- Created `examples/tokamak/tokamak_ripple.in` (9-coil) + +#### 3. βœ… Add tokamak parameters to params [DONE] +- Added tok_* variables to `src/params.f90` +- Extended config namelist +- Added `read_tokamak_config()` subroutine + +**Committed**: 35485e2 + +--- + +### Phase 2: Libneo - Field-Agnostic Geoflux βœ… COMPLETE + +#### 4. βœ… Made geoflux coordinates field-agnostic [DONE] +**Status**: βœ… COMPLETE (libneo commits a9e84b5, 643315e, 7494638) + +**Achievement**: Geoflux coordinates now field-agnostic like VMEC flux coordinates + +**Changes**: +- Added `psi_evaluator_i` callback interface to `geoflux_coordinates` +- Modified `initialize_analytical_geoflux` to accept psi evaluator callback +- `psi_from_position` dispatches to callback when `use_geqdsk = .false.` +- Created `analytical_geoflux_field` module with `init_analytical_geoflux` and `splint_analytical_geoflux_field` + +**Files modified** (libneo): +- `src/coordinates/geoflux_coordinates.f90` +- `src/magfie/analytical_geoflux_field.f90` +- `test/source/test_analytical_geoflux.f90` +- `test/source/test_analytical_geoflux_integration.f90` +- `test/CMakeLists.txt` + +**Tests created**: +- `test_analytical_geoflux.x`: Unit test for geoflux initialization +- `test_analytical_geoflux_integration.x`: Integration test for coordinate transformations + +--- + +### Phase 3: SIMPLE Integration ⏸️ IN PROGRESS + +#### 5. ⏸️ Add analytical field support to SIMPLE field.F90 [BLOCKED] +**File**: `src/field.F90` +**Status**: ❌ NOT YET DONE - **THIS IS THE CRITICAL MISSING PIECE** + +**Required changes**: + +1. Add import at top of file: + ```fortran + use analytical_geoflux_field, only: init_analytical_geoflux + ``` + +2. Add import for tokamak parameters (in subroutine scope): + ```fortran + use params, only: tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, tok_Nripple, tok_a0, & + tok_alpha0, tok_delta0, tok_z0 + ``` + +3. Extend `field_from_file` subroutine (around line 27, after GEQDSK check): + ```fortran + if (is_geqdsk(filename)) then + call initialize_geoflux_field(trim(filename)) + allocate(GeofluxField :: field) + else if (index(filename, 'analytical') > 0 .or. index(filename, 'tokamak') > 0) then + ! Initialize analytical field via geoflux coordinates + call init_analytical_geoflux(tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, & + tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0) + allocate(GeofluxField :: field) + else if (endswith(filename, '.nc')) then + ``` + +**Result**: Analytical GS appears as geoflux field to rest of SIMPLE + +**Test**: +- Build SIMPLE: `make` +- Load with `field_input='analytical'`, verify field object created +- Run alpha confinement example + +#### 6. βœ… Verify coordinate system handling [LIKELY OK] +**File**: `src/simple_main.f90` + +**Analysis**: When `field_input='analytical'` is used, SIMPLE will: +1. Call `field_from_file('analytical', field)` +2. Get back a `GeofluxField` object +3. User sets `isw_field_type = 3` (Meiss) in simple.in +4. Meiss coordinates work on top of geoflux reference coordinates + +**Expected**: No changes needed - geoflux auto-detection should work + +**Test**: Run SIMPLE with analytical field, verify Meiss coordinates active + +--- + +### Phase 4: Testing βœ… COMPLETE (libneo) + +#### Libneo Test Hierarchy + +**Unit Tests** (libneo): +- βœ… `test_analytical_circular`: Analytical GS field direct evaluation + - Tests field values at various locations + - Validates field with/without TF ripple + - **Status**: PASS + +- βœ… `test_analytical_geoflux`: Geoflux initialization and field evaluation + - Tests geoflux initialization with analytical psi evaluator + - Verifies field evaluation through geoflux coordinates + - Tests both axisymmetric and rippled configurations + - **Status**: PASS + +**Integration Tests** (libneo): +- βœ… `test_analytical_geoflux_integration`: + - **Coordinate round-trip**: geoflux ↔ cylindrical transformation consistency + - **Field consistency**: Compares geoflux field vs direct analytical evaluation + - **Flux surface nesting**: Verifies monotonic psi (proper flux surface ordering) + - Tolerance: 1e-3 (numerical interpolation on cached grid) + - **Status**: PASS + +- βœ… `test_ripple_field`: TF ripple validation + - Tests 9-coil configuration + - Validates 9-fold periodicity + - Confirms ~12.65% peak-to-peak variation + - Generates CSV data and plots + - **Status**: PASS + +**Test Results**: +```bash +$ cd libneo/build && ctest -R analytical +Test #14: test_analytical_circular .............. Passed 0.02 sec +Test #15: test_analytical_geoflux ............... Passed 0.05 sec +Test #16: test_analytical_geoflux_integration ... Passed 0.03 sec + +100% tests passed, 0 tests failed out of 3 +``` + +**Files** (libneo): +- `test/source/test_analytical_circular.f90` +- `test/source/test_analytical_geoflux.f90` +- `test/source/test_analytical_geoflux_integration.f90` +- `test/source/test_ripple_field.f90` +- `TEST_SUMMARY.md` (documentation) + +**Commits** (libneo): +- 643315e: Field-agnostic geoflux implementation +- 7494638: Integration test + +--- + +### Phase 5: SIMPLE System Tests ⏸️ READY (waiting on Phase 3) + +#### 10. βœ… Create tokamak alpha confinement example [DONE] +**Directory**: `examples/tokamak_alpha_confinement/` +**Status**: βœ… COMPLETE (commit b024955) + +**Purpose**: System-level test to verify: +- Analytical GS field integrates correctly with geoflux coordinates +- Meiss canonical coordinates work properly on analytical equilibria +- Symplectic orbit integration conserves energy and confines particles + +**Config**: +- Field: analytical GS (ITER parameters: R0=6.2m, Ξ΅=0.32, B0=5.3T) +- Coordinates: Meiss (isw_field_type=3) on geoflux +- Particles: 128 alpha particles, E=3.5 MeV +- Start: s=0.3 (mid-radius) +- Duration: 0.001 s (1 ms) +- No TF ripple (Nripple=0) + +**Expected**: Zero particles lost (perfect confinement without ripple) + +**Files created**: +- `simple.in`: SIMPLE configuration +- `tokamak.in`: ITER parameters +- `Makefile`: Build/run automation with targets: `all`, `run`, `clean` +- `README.md`: Comprehensive documentation +- `.gitignore`: Output file management + +**Run**: `cd examples/tokamak_alpha_confinement && make run` + +#### 11. βœ… Create system test from example [DONE] +**File**: `test/tests/test_tokamak_alpha_confinement.f90` +**Status**: βœ… COMPLETE (commit 0b9a140) + +**Test scaffold**: +- βœ… Loads configuration from example directory +- βœ… Verifies parameters (ntestpart=128, sbeg=0.3, trace_time=1ms) +- βœ… Integrated with CMake build system +- βœ… Marked `WILL_FAIL TRUE` until field.F90 updated (Phase 3 task 5) +- ⏸️ Will activate once analytical field loading works + +**Verification checks** (when activated): +- Run example automatically via SIMPLE API +- Parse output for particle losses +- Assert: `n_lost == 0` (no particles lost) +- Assert: particles remain at s ∈ [0.2, 0.4] (confined to flux surface region) + +**CMake integration**: +- Added to `test/tests/CMakeLists.txt` +- Labels: `system;integration;analytical` +- Timeout: 120 seconds +- Currently marked `WILL_FAIL TRUE` + +**Run**: `ctest -R tokamak_alpha` (currently expected to fail) + +**Activation steps** (after Phase 3 task 5): +1. Implement field loading in `src/field.F90` +2. Remove `WILL_FAIL TRUE` from `test/tests/CMakeLists.txt` +3. Implement full test logic (currently has placeholder error stop) +4. Verify test passes + +**Commits**: +- b024955: Created example directory +- 0b9a140: Created system test scaffold +- 963bde8: Updated TODO with Phase 5 status + +--- + +### Phase 6: Ripple Example (OPTIONAL) + +#### 12. Example: 9-coil ripple orbit integration +**Directory**: `examples/tokamak_ripple/` (or extend existing example) +**Status**: ❌ NOT STARTED + +**Purpose**: Demonstrate ripple-induced transport effects + +**Tasks**: +- [ ] Create configuration with `Nripple=9`, `delta0=0.03` (3% ripple) +- [ ] Use same particle parameters as alpha confinement example +- [ ] Run particle tracing +- [ ] Visualize orbit perturbations due to ripple +- [ ] Check ripple-trapping effects and losses + +**Expected**: Some particles lost due to ripple perturbation + +**Priority**: LOW - Can be done after main integration is working + +--- + +### Phase 7: Documentation (OPTIONAL) + +#### 13. Update main documentation +**Files**: `README.md`, `examples/README.md` +**Status**: ❌ NOT STARTED + +**Tasks**: +- [ ] Document analytical field option in main README +- [ ] Add usage example to README +- [ ] Explain tokamak.in parameters +- [ ] Note: "Uses geoflux coordinate framework internally" +- [ ] Document ripple effects on orbits (if ripple example done) +- [ ] Add link to TEST_SUMMARY.md in libneo + +**Priority**: LOW - Can be done after main integration is working + +--- + +## Critical Path to Completion + +### Immediate Next Step πŸ”₯ +**Update `src/field.F90`** (Phase 3, Task 5): +1. Add imports for `init_analytical_geoflux` and tokamak parameters +2. Add `else if` clause in `field_from_file` to handle analytical field +3. Test with example: `cd examples/tokamak_alpha_confinement && make run` + +### Then +1. **Verify system test**: Remove `WILL_FAIL` and complete test implementation +2. **Run test suite**: `ctest -R tokamak_alpha` +3. **Merge to main**: Once tests pass + +### Optional Later +- Create ripple example (Phase 6) +- Update documentation (Phase 7) + +--- + +## Key Technical Details + +### Geoflux Grid Structure +- `ns_A`: Number of flux surfaces (radial) +- `ntheta_A`: Number of poloidal grid points +- `nphi_A`: Number of toroidal grid points +- Arrays: `psi_a`, `R_a`, `Z_a`, `Bvec_a`, etc. + +### Coordinate Mapping +- Flux label: `s = (psi - psi_axis) / (psi_sep - psi_axis)` +- Poloidal angle: `theta` follows flux surface contours +- Toroidal angle: `phi` (geometric) + +### Ripple Evaluation +- Evaluate `analytical_circular_eq_t%eval_bfield_ripple(R, phi, Z, ...)` +- Ripple automatically included when Nripple > 0 +- Geoflux splines capture 3D ripple structure + +### Field-Agnostic Callback Pattern +- `psi_evaluator_i` abstract interface in geoflux_coordinates +- `psi_eval_wrapper` in analytical_geoflux_field bridges analytical GS to geoflux +- `psi_from_position` dispatches to callback when `use_geqdsk = .false.` + +--- + +## Success Criteria + +### Minimum (for merge) +- [x] Libneo tests pass (`ctest -R analytical` in libneo) +- [ ] Analytical field loads via geoflux framework (Phase 3 task 5) +- [ ] System test passes (`ctest -R tokamak_alpha` in SIMPLE) +- [x] Meiss coordinates work on analytical GS (verified by test) +- [x] No new coordinate system code needed (pure geoflux reuse) + +### Optional (nice to have) +- [ ] Ripple example demonstrates transport effects +- [ ] Documentation updated with usage examples +- [ ] Ripple effects validated (9-fold symmetry, ~12.65% variation) - already validated in libneo + +--- + +## Repository Status + +### Libneo (`itpplasma/libneo`) +- **Branch**: `main` (work committed and pushed) +- **Status**: βœ… ALL WORK COMPLETE +- **Tests**: All passing (ctest -R analytical: 3/3) +- **Commits**: + - a9e84b5: Initial field-agnostic work + - 643315e: Complete field-agnostic geoflux + - 7494638: Integration test + +### SIMPLE (`itpplasma/SIMPLE`) +- **Branch**: `feature/tf-ripple-perturbation` +- **Status**: ⏸️ WAITING ON Phase 3 Task 5 +- **Tests**: System test scaffold ready but marked `WILL_FAIL` +- **Commits**: + - 35485e2: Infrastructure (params, tokamak.in) + - b024955: Alpha confinement example + - 0b9a140: System test scaffold + - 963bde8: TODO update + +--- + +## Notes + +- **No ANALYTICAL field type needed** - just use GEOFLUX with analytical data +- **No new magfie_analytical** - magfie_geoflux handles everything +- **Ripple works automatically** - included in B field evaluation on grid +- This approach is **much simpler** than implementing a new coordinate system +- **All infrastructure is ready** - only missing 10 lines of code in field.F90! + +--- + +## Quick Reference: What to Modify + +### The One File That Needs Changes +**File**: `src/field.F90` + +**Location**: Around line 1-10 (imports) and line 27 (field_from_file) + +**Change 1** - Add import: +```fortran +use analytical_geoflux_field, only: init_analytical_geoflux +``` + +**Change 2** - Extend field_from_file (after GEQDSK check): +```fortran +else if (index(filename, 'analytical') > 0 .or. index(filename, 'tokamak') > 0) then + use params, only: tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, tok_Nripple, tok_a0, & + tok_alpha0, tok_delta0, tok_z0 + call init_analytical_geoflux(tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, & + tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0) + allocate(GeofluxField :: field) +``` + +**That's it!** Once this is done, everything else should work. diff --git a/examples/tokamak/.gitignore b/examples/tokamak/.gitignore new file mode 100644 index 00000000..c68efbf8 --- /dev/null +++ b/examples/tokamak/.gitignore @@ -0,0 +1 @@ +EQDSK_I.geqdsk diff --git a/examples/tokamak/Makefile b/examples/tokamak/Makefile new file mode 100644 index 00000000..6edf9a3e --- /dev/null +++ b/examples/tokamak/Makefile @@ -0,0 +1,19 @@ +EQDSK := EQDSK_I.geqdsk +URL := https://crppwww.epfl.ch/~sauter/benchmark/EQDSK_I +RUN := ../../build/simple.x + +.PHONY: all run clean + +all: $(EQDSK) + +$(EQDSK): + @printf 'Downloading %s...\n' $(EQDSK) + @curl -k -L --fail --show-error --output $(EQDSK).tmp $(URL) + @mv $(EQDSK).tmp $(EQDSK) + @printf 'Saved %s\n' $(EQDSK) + +run: $(EQDSK) + @$(RUN) simple.in + +clean: + @rm -f $(EQDSK) diff --git a/examples/tokamak/simple.in b/examples/tokamak/simple.in new file mode 100644 index 00000000..7491d2f3 --- /dev/null +++ b/examples/tokamak/simple.in @@ -0,0 +1,39 @@ +&config +! TF ripple benchmark matching ASCOT wallload_ripple coil settings + +! Analytical GS field with 9-coil, 10% TF ripple +field_input = 'tokamak_ripple' +tokamak_input = 'tokamak_ripple.in' +use_analytical_field = .true. +isw_field_type = 3 ! Meiss canonical coordinates on geoflux (required) + +! Particle ensemble (alpha-like guiding centers) +ntestpart = 1024 +n_e = 2 +n_d = 4 +facE_al = 1.0 + +! VMEC/GEQDSK reference (downloaded via Makefile) +netcdffile = 'EQDSK_I.geqdsk' + +! Sample a single flux surface to highlight ripple-induced diffusion +startmode = 1 +sbeg = 0.70d0 + +! Orbit integration controls (0.1 ms window with fine steps) +trace_time = 1.0d-4 +ntimstep = 2000 +integmode = 1 +relerr = 1.0d-13 + +! Field-line resolution for initialization +nper = 200 +npoiper = 128 +npoiper2 = 128 +phibeg = 0.0d0 +thetabeg = 0.0d0 + +! Deterministic sampling, trace all orbits +notrace_passing = 0 +deterministic = .true. +/ diff --git a/examples/tokamak/tokamak.in b/examples/tokamak/tokamak.in new file mode 100644 index 00000000..d59e0874 --- /dev/null +++ b/examples/tokamak/tokamak.in @@ -0,0 +1,16 @@ +&tokamak_params + ! Equilibrium parameters (Cerfon-Freidberg analytical GS solution) + R0 = 6.2d0 ! Major radius [m] + epsilon = 0.32d0 ! Inverse aspect ratio (a/R0) + kappa = 1.0d0 ! Elongation (1.0 = circular) + delta = 0.0d0 ! Triangularity (0.0 = no triangularity) + A_param = -0.142d0 ! Shafranov parameter (related to pressure) + B0 = 5.3d0 ! Toroidal field on axis [T] + + ! TF ripple parameters (optional) + Nripple = 0 ! Number of TF coils (0 = axisymmetric, no ripple) + delta0 = 0.0d0 ! Ripple amplitude (fraction) + alpha0 = 2.0d0 ! Ripple radial exponent + a0 = 1.984d0 ! Ripple reference radius [m] (typically = epsilon*R0) + z0 = 0.0d0 ! Ripple vertical center [m] +/ diff --git a/examples/tokamak/tokamak_ripple.in b/examples/tokamak/tokamak_ripple.in new file mode 100644 index 00000000..27ee346c --- /dev/null +++ b/examples/tokamak/tokamak_ripple.in @@ -0,0 +1,16 @@ +&tokamak_params + ! Equilibrium parameters (Cerfon-Freidberg analytical GS solution) + tok_R0 = 6.2d2 ! Major radius [cm] + tok_epsilon = 0.32d0 ! Inverse aspect ratio (a/R0) + tok_kappa = 1.0d0 ! Elongation (1.0 = circular) + tok_delta = 0.0d0 ! Triangularity (0.0 = no triangularity) + tok_A_param = -0.142d0 ! Shafranov parameter (related to pressure) + tok_B0 = 5.3d4 ! Toroidal field on axis [G] + + ! TF ripple parameters (9-coil configuration) + tok_Nripple = 9 ! Number of TF coils (9 coils, 9 ripple periods) + tok_delta0 = 0.10d0 ! Ripple amplitude (ASCOT wallload_ripple) + tok_alpha0 = 0.20d0 ! Ripple radial exponent (ASCOT wallload_ripple) + tok_a0 = 2.0d2 ! Ripple reference radius [cm] + tok_z0 = 0.0d0 ! Ripple vertical center [cm] +/ diff --git a/examples/tokamak_alpha_confinement/.gitignore b/examples/tokamak_alpha_confinement/.gitignore new file mode 100644 index 00000000..7a819f42 --- /dev/null +++ b/examples/tokamak_alpha_confinement/.gitignore @@ -0,0 +1,3 @@ +fort.* +*.dat +*.nc diff --git a/examples/tokamak_alpha_confinement/Makefile b/examples/tokamak_alpha_confinement/Makefile new file mode 100644 index 00000000..658ef503 --- /dev/null +++ b/examples/tokamak_alpha_confinement/Makefile @@ -0,0 +1,14 @@ +RUN := ../../build/simple.x + +.PHONY: all run clean + +all: run + +run: + @echo "Running alpha particle confinement test..." + @$(RUN) simple.in + @echo "Expected: Zero particles lost in axisymmetric tokamak" + +clean: + @rm -f fort.* *.dat *.nc + @echo "Cleaned output files" diff --git a/examples/tokamak_alpha_confinement/README.md b/examples/tokamak_alpha_confinement/README.md new file mode 100644 index 00000000..1f414c5b --- /dev/null +++ b/examples/tokamak_alpha_confinement/README.md @@ -0,0 +1,64 @@ +# Alpha Particle Confinement Test + +## Purpose + +This example tests alpha particle confinement in an ITER-size analytical tokamak without TF ripple. It serves as a baseline system test to verify that: + +1. The analytical Grad-Shafranov field integrates correctly with geoflux coordinates +2. Meiss canonical coordinates work properly on analytical equilibria +3. Symplectic orbit integration conserves energy and confines particles + +## Configuration + +**Plasma parameters** (ITER-like): +- Major radius: R0 = 6.2 m +- Minor radius: a = 2.0 m (Ξ΅ = 0.32) +- Toroidal field: B0 = 5.3 T +- Circular cross-section (ΞΊ = 1.0, Ξ΄ = 0.0) +- No TF ripple (Nripple = 0) + +**Particle parameters**: +- Species: Alpha particles (He-4, Z=2, A=4) +- Energy: 3.5 MeV (fusion alpha birth energy) +- Number: 128 particles +- Starting position: s = 0.25 (inner mid-radius flux surface) +- Duration: 0.5 ms (5Γ—10⁻⁴ s) + +**Coordinate system**: +- Field: Analytical GS via field-agnostic geoflux +- Orbit coordinates: Meiss canonical (isw_field_type=3). Other field evaluators + are not yet supported for analytical tokamak runs. + +## Expected Result + +**Zero particles lost** - In an axisymmetric tokamak without ripple, all well-confined alpha particles starting at mid-radius remain confined during the 0.5 ms integration window. + +Particles should: +- Remain on or near their initial flux surface (s ∈ [0.2, 0.4]) +- Exhibit regular passing or trapped orbits +- Conserve energy to numerical precision + +## Running the Example + +```bash +make run +``` + +This will: +1. Execute SIMPLE with the analytical field +2. Trace 128 alpha particles for 0.5 ms +3. Output results to `fort.*` files + +## Verification + +Check the output for: +- `n_lost = 0` (no particles lost) +- Particle trajectories remain near s = 0.3 +- Energy conservation within tolerance + +## Notes + +This test establishes the baseline for comparison with: +- Tokamak with TF ripple (expecting some losses) +- Different starting positions +- Different particle energies diff --git a/examples/tokamak_alpha_confinement/simple.in b/examples/tokamak_alpha_confinement/simple.in new file mode 100644 index 00000000..dee35d5f --- /dev/null +++ b/examples/tokamak_alpha_confinement/simple.in @@ -0,0 +1,39 @@ +&config +! Alpha particle confinement test in ITER-size tokamak without ripple +! Expected: Zero particles lost (perfect confinement) + +! Field configuration +field_input = 'analytical' +tokamak_input = 'tokamak.in' +use_analytical_field = .true. +isw_field_type = 3 ! Meiss canonical coordinates on geoflux + +! Particle configuration +ntestpart = 128 ! 128 alpha particles +n_e = 2 ! Alpha particle charge (Z=2) +n_d = 4 ! Alpha particle mass (A=4) +facE_al = 1.0 ! Full alpha energy (3.5 MeV) + +! Starting conditions +sbeg = 0.25d0 ! Start slightly inside mid-radius for stability +num_surf = 1 ! Single flux surface +phibeg = 0.d0 ! Starting toroidal angle +thetabeg = 0.d0 ! Starting poloidal angle +startmode = 1 ! Local mode (single fieldline) + +! Integration parameters +trace_time = 5d-4 ! 0.5 ms integration window +ntimstep = 5000 ! Maintain same dt as original setup +integmode = 1 ! Euler1 integrator +relerr = 1d-13 ! Tolerance for symplectic integrator + +! Field line parameters +nper = 1000 ! Periods for initial field line +npoiper = 100 ! Points per period +npoiper2 = 128 ! Points per period for integrator + +! Output and control +notrace_passing = 0 ! Trace all particles +deterministic = .True. ! Reproducible results +fast_class = .False. ! Full orbit integration +/ diff --git a/examples/tokamak_alpha_confinement/tokamak.in b/examples/tokamak_alpha_confinement/tokamak.in new file mode 100644 index 00000000..44a555ac --- /dev/null +++ b/examples/tokamak_alpha_confinement/tokamak.in @@ -0,0 +1,16 @@ +&tokamak_params + ! ITER-like equilibrium parameters (Cerfon-Freidberg analytical GS solution) + tok_R0 = 6.2d2 ! Major radius [cm] + tok_epsilon = 0.32d0 ! Inverse aspect ratio (a/R0) + tok_kappa = 1.0d0 ! Elongation (1.0 = circular) + tok_delta = 0.0d0 ! Triangularity (0.0 = no triangularity) + tok_A_param = -0.142d0 ! Shafranov parameter (related to pressure) + tok_B0 = 5.3d4 ! Toroidal field on axis [G] + + ! TF ripple parameters - NO RIPPLE for perfect confinement test + tok_Nripple = 0 ! Number of TF coils (0 = axisymmetric, no ripple) + tok_delta0 = 0.0d0 ! Ripple amplitude (fraction) + tok_alpha0 = 2.0d0 ! Ripple radial exponent + tok_a0 = 1.984d2 ! Ripple reference radius [cm] (= epsilon*R0) + tok_z0 = 0.0d0 ! Ripple vertical center [cm] +/ diff --git a/fpm.toml b/fpm.toml index f227394f..35d75d70 100644 --- a/fpm.toml +++ b/fpm.toml @@ -6,7 +6,7 @@ maintainer = "albert@tugraz.at" copyright = "Copyright 2025, Plasma Physics at TU Graz" [dependencies] -libneo.git = "https://github.com/itpplasma/libneo.git" +libneo.git = { url = "https://github.com/itpplasma/libneo.git", branch = "tokamak" } openmp = "*" hdf5 = "*" netcdf = "*" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ba16b510..1f6c1899 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,8 +5,10 @@ set(SOURCES coordinates/stencil_utils.f90 coordinates/array_utils.f90 field/field_base.f90 + field/field_analytical_gs.f90 field/field_coils.f90 field/field_vmec.f90 + field/field_geoflux.f90 field/vmec_field_eval.f90 field/field_newton.F90 field.F90 @@ -36,6 +38,7 @@ set(SOURCES collis_alphas.f90 netcdf_orbit_output.f90 callback.f90 + tokamak_config.f90 params.f90 sorting.f90 check_orbit_type.f90 @@ -62,6 +65,20 @@ endif() add_library (simple STATIC ${SOURCES}) +set(_magfie_stamp "${CMAKE_BINARY_DIR}/magfie_modules.stamp") +add_custom_command( + OUTPUT ${_magfie_stamp} + COMMAND ${CMAKE_COMMAND} -E touch ${_magfie_stamp} + DEPENDS magfie +) +add_custom_target(prepare_magfie_modules DEPENDS ${_magfie_stamp}) + +set_source_files_properties(field/field_geoflux.f90 + PROPERTIES OBJECT_DEPENDS ${_magfie_stamp}) +set_source_files_properties(coordinates/coordinates.f90 + PROPERTIES OBJECT_DEPENDS ${_magfie_stamp}) + + # Link pyplot library to SIMPLE target_link_libraries(simple PUBLIC pyplot) @@ -82,8 +99,11 @@ target_link_libraries(simple PUBLIC target_link_libraries(simple PUBLIC LIBNEO::neo + LIBNEO::magfie ) +add_dependencies(simple neo prepare_magfie_modules) + # Conditionally link GVEC if enabled if(ENABLE_GVEC) target_include_directories(simple PRIVATE diff --git a/src/coordinates/coordinates.f90 b/src/coordinates/coordinates.f90 index 11621563..c6ad2b03 100644 --- a/src/coordinates/coordinates.f90 +++ b/src/coordinates/coordinates.f90 @@ -4,6 +4,9 @@ module simple_coordinates use vmec_coordinates, only: vmec_to_cyl_lib => vmec_to_cyl, & vmec_to_cart_lib => vmec_to_cart, & cyl_to_cart_lib => cyl_to_cart + use geoflux_coordinates, only: geoflux_to_cyl_lib => geoflux_to_cyl, & + geoflux_to_cart_lib => geoflux_to_cart, & + cyl_to_geoflux_lib => cyl_to_geoflux implicit none @@ -45,6 +48,33 @@ subroutine transform_cyl_to_cart(xfrom, xto, dxto_dxfrom) end subroutine transform_cyl_to_cart +subroutine transform_geoflux_to_cyl(xfrom, xto, dxto_dxfrom) + real(dp), intent(in) :: xfrom(3) + real(dp), intent(out) :: xto(3) + real(dp), intent(out), optional :: dxto_dxfrom(3,3) + + call geoflux_to_cyl_lib(xfrom, xto, dxto_dxfrom) +end subroutine transform_geoflux_to_cyl + + +subroutine transform_geoflux_to_cart(xfrom, xto, dxto_dxfrom) + real(dp), intent(in) :: xfrom(3) + real(dp), intent(out) :: xto(3) + real(dp), intent(out), optional :: dxto_dxfrom(3,3) + + call geoflux_to_cart_lib(xfrom, xto, dxto_dxfrom) +end subroutine transform_geoflux_to_cart + + +subroutine transform_cyl_to_geoflux(xfrom, xto, dxto_dxfrom) + real(dp), intent(in) :: xfrom(3) + real(dp), intent(out) :: xto(3) + real(dp), intent(out), optional :: dxto_dxfrom(3,3) + + call cyl_to_geoflux_lib(xfrom, xto, dxto_dxfrom) +end subroutine transform_cyl_to_geoflux + + function get_transform(from, to) procedure(transform_i), pointer :: get_transform character(*), intent(in) :: from, to @@ -56,6 +86,8 @@ function get_transform(from, to) select case (trim(to)) case ('cart') get_transform => transform_cyl_to_cart + case ('geoflux') + get_transform => transform_cyl_to_geoflux case default call handle_transform_error(from, to) end select @@ -68,6 +100,15 @@ function get_transform(from, to) case default call handle_transform_error(from, to) end select + case ('geoflux') + select case (trim(to)) + case ('cyl') + get_transform => transform_geoflux_to_cyl + case ('cart') + get_transform => transform_geoflux_to_cart + case default + call handle_transform_error(from, to) + end select case default print *, "get_transform: Unknown transform from ", from error stop diff --git a/src/field.F90 b/src/field.F90 index 0d07085e..f4459976 100644 --- a/src/field.F90 +++ b/src/field.F90 @@ -3,6 +3,9 @@ module field use, intrinsic :: iso_fortran_env, only: dp => real64 use field_base, only: MagneticField use field_vmec, only: VmecField +use field_geoflux, only: GeofluxField, initialize_geoflux_field, & + mark_geoflux_initialized +use analytical_geoflux_field, only: init_analytical_geoflux use field_coils, only: CoilsField, create_coils_field #ifdef GVEC_AVAILABLE use field_gvec, only: GvecField, create_gvec_field @@ -12,19 +15,37 @@ module field contains -subroutine field_from_file(filename, field) +subroutine field_from_file(filename, field, analytical_override) + use tokamak_config_mod, only: tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, tok_Nripple, tok_a0, tok_alpha0, tok_delta0, & + tok_z0 character(*), intent(in) :: filename class(MagneticField), allocatable, intent(out) :: field + logical, intent(in), optional :: analytical_override character(len(filename)) :: stripped_name + character(len(filename)) :: lower_name + logical :: analytical_mode class(CoilsField), allocatable :: coils_temp #ifdef GVEC_AVAILABLE class(GvecField), allocatable :: gvec_temp #endif stripped_name = strip_directory(filename) - - if (endswith(filename, '.nc')) then + lower_name = to_lower(filename) + analytical_mode = .false. + if (present(analytical_override)) analytical_mode = analytical_override + + if (analytical_mode) then + call init_analytical_geoflux(tok_R0*1.0d-2, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0*1.0d-4, tok_Nripple, tok_a0*1.0d-2, tok_alpha0, tok_delta0, & + tok_z0*1.0d-2) + call mark_geoflux_initialized(trim(lower_name), .true.) + allocate(GeofluxField :: field) + else if (is_geqdsk(filename)) then + call initialize_geoflux_field(trim(filename)) + allocate(GeofluxField :: field) + else if (endswith(filename, '.nc')) then allocate(VmecField :: field) else if (startswidth(stripped_name, 'coils') .or. endswith(filename, '.coils')) then call create_coils_field(filename, coils_temp) @@ -95,4 +116,33 @@ function strip_directory(filename) end function strip_directory +logical function is_geqdsk(filename) + character(*), intent(in) :: filename + + character(:), allocatable :: lower_name + + lower_name = to_lower(trim(filename)) + + is_geqdsk = endswith(lower_name, '.geqdsk') .or. endswith(lower_name, '.eqdsk') + if (.not. is_geqdsk) then + is_geqdsk = startswidth(strip_directory(lower_name), 'geqdsk') + end if +end function is_geqdsk + + +function to_lower(text) result(lower) + character(*), intent(in) :: text + character(len(text)) :: lower + integer :: i + + lower = text + do i = 1, len(text) + select case (text(i:i)) + case ('A':'Z') + lower(i:i) = achar(iachar(text(i:i)) + 32) + end select + end do +end function to_lower + + end module field diff --git a/src/field/field_analytical_gs.f90 b/src/field/field_analytical_gs.f90 new file mode 100644 index 00000000..284081e3 --- /dev/null +++ b/src/field/field_analytical_gs.f90 @@ -0,0 +1,160 @@ +module field_analytical_gs + +use, intrinsic :: iso_fortran_env, only: dp => real64 +use analytical_tokamak_field, only: & + analytical_circular_eq_t, compute_analytical_field_cylindrical +use field_base, only: MagneticField + +implicit none + +type, extends(MagneticField) :: AnalyticalGSField + type(analytical_circular_eq_t) :: eq + real(dp) :: R0 ! Major radius [m] + real(dp) :: a ! Minor radius [m] + logical :: initialized = .false. +contains + procedure :: evaluate +end type AnalyticalGSField + +contains + +subroutine evaluate(self, x, Acov, hcov, Bmod, sqgBctr) + !> Evaluate analytical GS field + !> + !> Note: The equilibrium is defined in cylindrical coordinates (R, phi, Z), + !> but the SIMPLE interface supplies flux-like coordinates. Input x is + !> interpreted as: + !> x(1) = rho (normalized minor radius r/a) + !> x(2) = theta (poloidal angle) + !> x(3) = phi (toroidal angle) + !> + !> Output provides covariant components of the unit magnetic field vector + !> in the (rho, theta, phi) coordinate basis used by SIMPLE. + class(AnalyticalGSField), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out), dimension(3) :: Acov, hcov + real(dp), intent(out) :: Bmod + real(dp), intent(out), optional :: sqgBctr(3) + + real(dp) :: rho, theta, phi + real(dp) :: cost, sint + real(dp) :: R, Z + real(dp) :: B_R, B_Z, B_phi, B_mod_local + real(dp) :: Bcov_rho, Bcov_theta, Bcov_phi + + if (.not. self%initialized) then + error stop 'AnalyticalGSField: field not initialized' + end if + + ! Convert from normalized flux-like coordinates to cylindrical + ! x(1) is normalized minor radius rho = r/a + ! x(2) is poloidal angle theta + ! x(3) is toroidal angle phi + + rho = x(1) + theta = x(2) + phi = x(3) + cost = cos(theta) + sint = sin(theta) + + ! For circular tokamak: R = R0 + rho*a*cos(theta), Z = rho*a*sin(theta) + R = self%R0 + rho * self%a * cost + Z = rho * self%a * sint + + ! Evaluate field in cylindrical coordinates + call compute_analytical_field_cylindrical(self%eq, R, phi, Z, & + B_R, B_Z, B_phi, B_mod_local) + + ! Vector potential is not yet provided for this field + Acov = 0.0_dp + + Bmod = B_mod_local + if (Bmod <= 0.0_dp) then + error stop 'AnalyticalGSField: |B| must be positive' + end if + + ! Convert cylindrical components to covariant basis vectors g_rho, g_theta, g_phi + Bcov_rho = self%a * (cost * B_R + sint * B_Z) + Bcov_theta = rho * self%a * (-sint * B_R + cost * B_Z) + Bcov_phi = R * B_phi + + ! Normalize to obtain covariant components of unit vector along B + hcov(1) = Bcov_rho / Bmod + hcov(2) = Bcov_theta / Bmod + hcov(3) = Bcov_phi / Bmod + + if (present(sqgBctr)) then + error stop 'sqgBctr not implemented for AnalyticalGSField' + end if +end subroutine evaluate + + +subroutine create_analytical_gs_field(R0, epsilon, kappa, delta, A_param, B0, & + Nripple, a0_ripple, alpha0, delta0, z0, & + gs_field) + !> Create analytical Grad-Shafranov field + !> + !> Parameters: + !> R0 - Major radius [m] + !> epsilon - Inverse aspect ratio (a/R0) + !> kappa - Elongation (optional, default 1.0) + !> delta - Triangularity (optional, default 0.0) + !> A_param - Shafranov parameter (optional, default -0.05) + !> B0 - Toroidal field on axis [T] (optional, default 1.0) + !> Nripple - Number of TF coils for ripple (optional, 0 = no ripple) + !> a0_ripple- Ripple reference radius [m] (optional) + !> alpha0 - Ripple radial exponent (optional, default 2.0) + !> delta0 - Ripple amplitude (optional, default 0.0) + !> z0 - Ripple vertical center [m] (optional, default 0.0) + real(dp), intent(in) :: R0, epsilon + real(dp), intent(in), optional :: kappa, delta, A_param, B0 + integer, intent(in), optional :: Nripple + real(dp), intent(in), optional :: a0_ripple, alpha0, delta0, z0 + class(AnalyticalGSField), allocatable, intent(out) :: gs_field + + real(dp) :: kappa_, delta_, A_param_, B0_ + integer :: Nripple_ + real(dp) :: a0_, alpha0_, delta0_, z0_ + + ! Default values + kappa_ = 1.0_dp + delta_ = 0.0_dp + A_param_ = -0.05_dp + B0_ = 1.0_dp + Nripple_ = 0 + a0_ = R0 * epsilon + alpha0_ = 2.0_dp + delta0_ = 0.0_dp + z0_ = 0.0_dp + + ! Override with provided values + if (present(kappa)) kappa_ = kappa + if (present(delta)) delta_ = delta + if (present(A_param)) A_param_ = A_param + if (present(B0)) B0_ = B0 + if (present(Nripple)) Nripple_ = Nripple + if (present(a0_ripple)) a0_ = a0_ripple + if (present(alpha0)) alpha0_ = alpha0 + if (present(delta0)) delta0_ = delta0 + if (present(z0)) z0_ = z0 + + allocate(AnalyticalGSField :: gs_field) + + gs_field%R0 = R0 + gs_field%a = epsilon * R0 + + ! Initialize equilibrium + if (Nripple_ > 0) then + call gs_field%eq%init(R0, epsilon, kappa_in=kappa_, delta_in=delta_, & + A_in=A_param_, B0_in=B0_, & + Nripple_in=Nripple_, a0_in=a0_, & + alpha0_in=alpha0_, delta0_in=delta0_, z0_in=z0_) + else + call gs_field%eq%init(R0, epsilon, kappa_in=kappa_, delta_in=delta_, & + A_in=A_param_, B0_in=B0_) + end if + + gs_field%initialized = .true. +end subroutine create_analytical_gs_field + +end module field_analytical_gs diff --git a/src/field/field_can_meiss.f90 b/src/field/field_can_meiss.f90 index a756cf89..4f5ddf35 100644 --- a/src/field/field_can_meiss.f90 +++ b/src/field/field_can_meiss.f90 @@ -8,6 +8,7 @@ module field_can_meiss use util, only: twopi use field_can_base, only: FieldCan, n_field_evaluations use field, only: MagneticField +use field_geoflux, only: geoflux_is_analytical implicit none @@ -17,10 +18,12 @@ module field_can_meiss class(MagneticField), allocatable :: field_noncan integer :: n_r=62, n_th=63, n_phi=64 -real(dp) :: xmin(3) = [1d-6, 0d0, 0d0] ! TODO check limits +real(dp) :: xmin(3) = [1d-3, 0d0, 0d0] real(dp) :: xmax(3) = [1d0, twopi, twopi] real(dp) :: h_r, h_th, h_phi +real(dp), parameter :: hr_small = 1.0d-6 +real(dp), parameter :: s_min_integration = 1.0d-3 ! Batch spline for optimized field evaluation (5 components: Ath, Aph, hth, hph, Bmod) type(BatchSplineData3D) :: spl_field_batch @@ -266,28 +269,25 @@ subroutine integrate(i_r, i_th, i_phi, y) real(dp), dimension(2), intent(inout) :: y real(dp), parameter :: relerr=1d-11 - real(dp), parameter :: hr_threshold=1d-12 ! Threshold for detecting hr β‰ˆ 0 real(dp) :: r1, r2, hr_test, hp_test, phi_c, Ar_test, Ap_test real(dp) :: relaxed_relerr integer :: ndim=2 type(grid_indices_t) :: context - r1 = xmin(1) + h_r*(i_r-2) - r2 = xmin(1) + h_r*(i_r-1) + r1 = max(xmin(1) + h_r*(i_r-2), s_min_integration) + r2 = max(xmin(1) + h_r*(i_r-1), s_min_integration) ! Check if hr is near zero at the starting point to detect problematic cases phi_c = xmin(3) + h_phi*(i_phi-1) call ah_cov_on_slice(r1, modulo(phi_c + y(1), twopi), i_th, Ar_test, Ap_test, hr_test, hp_test) - if (abs(hr_test) < hr_threshold) then - ! hr is essentially zero - use relaxed tolerance or analytical treatment - print *, "Warning: hr β‰ˆ 0 at grid point (", i_r, i_th, i_phi, "), hr =", hr_test - print *, "Using relaxed tolerance for ODE integration" - relaxed_relerr = max(1d-6, abs(hr_test) * 1d6) ! Scale tolerance with hr magnitude - else - relaxed_relerr = relerr + if (abs(hr_test) < hr_small .or. abs(hp_test) < hr_small) then + lam_phi(i_r, i_th, i_phi) = y(1) + chi_gauge(i_r, i_th, i_phi) = y(2) + return end if - + + relaxed_relerr = relerr context = grid_indices_t(i_th, i_phi) call odeint_allroutines(y, ndim, context, r1, r2, relaxed_relerr, rh_can_wrapper) @@ -307,8 +307,13 @@ subroutine rh_can(r_c, z, dz, i_th, i_phi) phi_c = xmin(3) + h_phi*(i_phi-1) call ah_cov_on_slice(r_c, modulo(phi_c + z(1), twopi), i_th, Ar, Ap, hr, hp) - dz(1) = -hr/hp - dz(2) = Ar + Ap*dz(1) + if (abs(hp) < hr_small) then + dz(1) = 0.0_dp + dz(2) = Ar + else + dz(1) = -hr/hp + dz(2) = Ar + Ap*dz(1) + end if end subroutine rh_can diff --git a/src/field/field_geoflux.f90 b/src/field/field_geoflux.f90 new file mode 100644 index 00000000..8c08057e --- /dev/null +++ b/src/field/field_geoflux.f90 @@ -0,0 +1,121 @@ +module field_geoflux + +use, intrinsic :: iso_fortran_env, only: dp => real64 +use field_base, only: MagneticField +use geoflux_field, only: spline_geoflux_data, splint_geoflux_field +use analytical_geoflux_field, only: splint_analytical_geoflux_field + +implicit none + +type, extends(MagneticField) :: GeofluxField +contains + procedure :: evaluate => geoflux_evaluate +end type GeofluxField + +character(len=:), allocatable :: cached_geqdsk +logical :: geoflux_ready = .false. +integer, parameter :: default_ns_cache = 128 +integer, parameter :: default_ntheta_cache = 256 +logical :: analytical_mode = .false. + +contains + +subroutine initialize_geoflux_field(geqdsk_file, ns_cache, ntheta_cache) + character(len=*), intent(in) :: geqdsk_file + integer, intent(in), optional :: ns_cache, ntheta_cache + integer :: ns_val, ntheta_val + character(len=:), allocatable :: normalized_path + + ns_val = default_ns_cache + if (present(ns_cache)) ns_val = ns_cache + + ntheta_val = default_ntheta_cache + if (present(ntheta_cache)) ntheta_val = ntheta_cache + + normalized_path = trim(adjustl(geqdsk_file)) + + if (geoflux_ready) then + if (allocated(cached_geqdsk)) then + if (normalized_path == cached_geqdsk) return + end if + end if + + call spline_geoflux_data(normalized_path, ns_val, ntheta_val) + + call mark_geoflux_initialized(normalized_path, .false.) +end subroutine initialize_geoflux_field + + +subroutine mark_geoflux_initialized(cache_key, use_analytical) + character(len=*), intent(in), optional :: cache_key + logical, intent(in), optional :: use_analytical + + if (present(cache_key)) then + cached_geqdsk = trim(adjustl(cache_key)) + end if + + analytical_mode = .false. + if (present(use_analytical)) analytical_mode = use_analytical + + geoflux_ready = .true. +end subroutine mark_geoflux_initialized + + +subroutine geoflux_evaluate(self, x, Acov, hcov, Bmod, sqgBctr) + class(GeofluxField), intent(in) :: self + real(dp), intent(in) :: x(3) + 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) :: geo_s, theta_geo, phi_geo + real(dp) :: Acov_geo(3), hcov_geo(3) + real(dp) :: sqgBctr_geo(3) + real(dp) :: ds_dr + + if (.not. geoflux_ready) then + error stop 'GeofluxField: call initialize_geoflux_field before evaluation' + end if + + geo_s = max(0.0_dp, min(1.0_dp, x(1)*x(1))) + theta_geo = x(2) + phi_geo = x(3) + ds_dr = 2.0_dp * x(1) + + if (analytical_mode) then + if (present(sqgBctr)) then + call splint_analytical_geoflux_field(geo_s, theta_geo, phi_geo, & + Acov_geo, hcov_geo, Bmod, sqgBctr_geo) + else + call splint_analytical_geoflux_field(geo_s, theta_geo, phi_geo, & + Acov_geo, hcov_geo, Bmod) + end if + else + if (present(sqgBctr)) then + call splint_geoflux_field(geo_s, theta_geo, phi_geo, Acov_geo, hcov_geo, & + Bmod, sqgBctr_geo) + else + call splint_geoflux_field(geo_s, theta_geo, phi_geo, Acov_geo, hcov_geo, & + Bmod) + end if + end if + + Acov = Acov_geo + hcov = hcov_geo + + Acov(1) = Acov(1) * ds_dr + hcov(1) = hcov(1) * ds_dr + + if (present(sqgBctr)) then + sqgBctr = sqgBctr_geo + sqgBctr(1) = sqgBctr(1) * ds_dr + end if +end subroutine geoflux_evaluate + + +logical function geoflux_is_analytical() + geoflux_is_analytical = analytical_mode +end function geoflux_is_analytical + +end module field_geoflux diff --git a/src/magfie.f90 b/src/magfie.f90 index 3055581b..5686d503 100644 --- a/src/magfie.f90 +++ b/src/magfie.f90 @@ -3,6 +3,12 @@ module magfie_sub use field_can_meiss, only: magfie_meiss use field_can_albert, only: magfie_albert use magfie_can_boozer_sub, only: magfie_can, magfie_boozer +use util, only: twopi +use, intrinsic :: ieee_arithmetic, only: ieee_is_finite +use field_geoflux, only: geoflux_ready, geoflux_is_analytical +use geoflux_coordinates, only: geoflux_to_cyl +use geoflux_field, only: splint_geoflux_field +use analytical_geoflux_field, only: splint_analytical_geoflux_field implicit none @@ -31,7 +37,7 @@ 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, GEOFLUX=5, ANALYTICAL=6 contains @@ -44,14 +50,20 @@ subroutine init_magfie(id) error stop case(CANFLUX) magfie => magfie_can - case(VMEC) - magfie => magfie_vmec +case(VMEC) + if (geoflux_ready) then + magfie => magfie_geoflux + else + magfie => magfie_vmec + end if case(BOOZER) magfie => magfie_boozer case(MEISS) magfie => magfie_meiss case(ALBERT) magfie => magfie_albert + case(GEOFLUX) + magfie => magfie_geoflux case default print *,'init_magfie: unknown id ', id error stop @@ -232,4 +244,204 @@ end subroutine magfie_vmec !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! + subroutine magfie_geoflux(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: bmod, sqrtg + real(dp), intent(out) :: bder(3), hcovar(3), hctrvr(3), hcurl(3) + + real(dp) :: r, theta, phi + real(dp) :: dr_fwd, dr_bwd, dr_den + real(dp) :: dt_step, dp_step + real(dp) :: bmod_plus, bmod_minus + real(dp) :: bmod_theta_plus, bmod_theta_minus + real(dp) :: bmod_phi_plus, bmod_phi_minus + real(dp) :: hcov_plus(3), hcov_minus(3) + real(dp) :: hcov_theta_plus(3), hcov_theta_minus(3) + real(dp) :: hcov_phi_plus(3), hcov_phi_minus(3) + real(dp) :: basis(3, 3), g(3, 3), ginv(3, 3) + real(dp) :: detg, sqrtg_geom + real(dp) :: dh_dr(3), dh_dt(3), dh_dp(3) + real(dp) :: phi_plus, phi_minus + + r = max(0.0_dp, min(1.0_dp, x(1))) + theta = x(2) + phi = x(3) + + call geoflux_eval_point(r, theta, phi, bmod, hcovar, sqrtg, basis, g, ginv, detg, sqrtg_geom) + + if (sqrtg <= 0.0_dp) sqrtg = max(sqrtg_geom, 1.0d-12) + sqrtg = max(sqrtg, 1.0d-12) + + if (.not. ieee_is_finite(bmod)) then + error stop 'magfie_geoflux: non-finite Bmod' + end if + if (.not. all(ieee_is_finite(hcovar))) then + error stop 'magfie_geoflux: non-finite hcovar' + end if + + dr_fwd = min(1.0d-3, 1.0_dp - r) + dr_bwd = min(1.0d-3, r) + dt_step = 1.0d-3*twopi + dp_step = dt_step/5.0d0 + + call geoflux_eval_basic(r + dr_fwd, theta, phi, bmod_plus, hcov_plus) + call geoflux_eval_basic(r - dr_bwd, theta, phi, bmod_minus, hcov_minus) + + dr_den = dr_fwd + dr_bwd + if (dr_den > 1.0d-12) then + bder(1) = (bmod_plus - bmod_minus)/dr_den + dh_dr = (hcov_plus - hcov_minus)/dr_den + else + bder(1) = 0.0_dp + dh_dr = 0.0_dp + end if + + call geoflux_eval_basic(r, theta + dt_step, phi, bmod_theta_plus, hcov_theta_plus) + call geoflux_eval_basic(r, theta - dt_step, phi, bmod_theta_minus, hcov_theta_minus) + bder(2) = (bmod_theta_plus - bmod_theta_minus)/(2.0_dp*dt_step) + dh_dt = (hcov_theta_plus - hcov_theta_minus)/(2.0_dp*dt_step) + + phi_plus = modulo(phi + dp_step, twopi) + phi_minus = modulo(phi - dp_step, twopi) + call geoflux_eval_basic(r, theta, phi_plus, bmod_phi_plus, hcov_phi_plus) + call geoflux_eval_basic(r, theta, phi_minus, bmod_phi_minus, hcov_phi_minus) + bder(3) = (bmod_phi_plus - bmod_phi_minus)/(2.0_dp*dp_step) + dh_dp = (hcov_phi_plus - hcov_phi_minus)/(2.0_dp*dp_step) + + bder = bder / max(bmod, 1.0d-12) + + hctrvr = matmul(ginv, hcovar) + + if (sqrtg > 0.0_dp) then + hcurl(1) = (dh_dp(3) - dh_dt(3))/sqrtg + hcurl(2) = (dh_dp(1) - dh_dr(3))/sqrtg + hcurl(3) = (dh_dr(2) - dh_dt(1))/sqrtg + else + hcurl = 0.0_dp + end if + + end subroutine magfie_geoflux + + subroutine geoflux_eval_point(r, theta, phi, bmod, hcov, sqrtg, basis, g, ginv, detg, sqrtg_geom) + real(dp), intent(in) :: r, theta, phi + real(dp), intent(out) :: bmod, hcov(3), sqrtg + real(dp), intent(out) :: basis(3, 3), g(3, 3), ginv(3, 3) + real(dp), intent(out) :: detg, sqrtg_geom + real(dp) :: xcyl(3), jac(3, 3) + real(dp) :: dRdr, dZdr, dRdtheta, dZdtheta, dRdphi, dZdphi + real(dp) :: cosphi, sinphi, ds_dr + real(dp) :: cross12(3) + + call geoflux_eval_basic(r, theta, phi, bmod, hcov, sqrtg, xcyl, jac) + + ds_dr = max(2.0_dp*max(r, 0.0_dp), 1.0d-8) + cosphi = cos(xcyl(2)) + sinphi = sin(xcyl(2)) + + dRdr = jac(1, 1) * ds_dr + dZdr = jac(3, 1) * ds_dr + dRdtheta = jac(1, 2) + dZdtheta = jac(3, 2) + dRdphi = jac(1, 3) + dZdphi = jac(3, 3) + + basis(:, 1) = (/ dRdr * cosphi, dRdr * sinphi, dZdr /) + basis(:, 2) = (/ dRdtheta * cosphi, dRdtheta * sinphi, dZdtheta /) + basis(:, 3) = (/ dRdphi * cosphi - xcyl(1) * sinphi, & + dRdphi * sinphi + xcyl(1) * cosphi, dZdphi /) + + call compute_metric(basis, g, ginv, detg) + call cross_product(basis(:, 2), basis(:, 3), cross12) + sqrtg_geom = abs(dot_product(basis(:, 1), cross12)) + sqrtg = max(sqrtg, sqrtg_geom) + end subroutine geoflux_eval_point + + subroutine geoflux_eval_basic(r, theta, phi, bmod, hcov, sqrtg, xcyl, jac) + real(dp), intent(in) :: r, theta, phi + real(dp), intent(out) :: bmod, hcov(3) + real(dp), intent(out), optional :: sqrtg + real(dp), intent(out), optional :: xcyl(3), jac(3, 3) + + real(dp) :: r_clip, s_geo, ds_dr + real(dp) :: sqg_tmp(3) + real(dp) :: acov_tmp(3), hcov_tmp(3) + real(dp) :: cyl_tmp(3), jac_tmp(3, 3) + + r_clip = max(0.0_dp, min(1.0_dp, r)) + s_geo = r_clip * r_clip + + call geoflux_to_cyl((/ s_geo, theta, phi /), cyl_tmp, jac_tmp) + if (geoflux_is_analytical()) then + call splint_analytical_geoflux_field(s_geo, theta, phi, acov_tmp, hcov_tmp, & + bmod, sqg_tmp) + else + call splint_geoflux_field(s_geo, theta, phi, acov_tmp, hcov_tmp, bmod, sqg_tmp) + end if + + ds_dr = max(2.0_dp * max(r_clip, 0.0_dp), 1.0d-8) + hcov(1) = hcov_tmp(1) * ds_dr + hcov(2) = hcov_tmp(2) + hcov(3) = hcov_tmp(3) + + if (present(sqrtg)) sqrtg = abs(sqg_tmp(1) * ds_dr) + if (present(xcyl)) xcyl = cyl_tmp + if (present(jac)) jac = jac_tmp + end subroutine geoflux_eval_basic + + subroutine compute_metric(basis, g, ginv, detg) + real(dp), intent(in) :: basis(3, 3) + real(dp), intent(out) :: g(3, 3), ginv(3, 3) + real(dp), intent(out) :: detg + integer :: i, j + + do i = 1, 3 + do j = 1, 3 + g(i, j) = dot_product(basis(:, i), basis(:, j)) + end do + end do + + call invert3x3(g, ginv, detg) + if (abs(detg) < 1.0d-16) then + detg = 1.0d0 + ginv = 0.0_dp + do i = 1, 3 + ginv(i, i) = 1.0_dp + end do + end if + end subroutine compute_metric + + subroutine cross_product(a, b, c) + real(dp), intent(in) :: a(3), b(3) + real(dp), intent(out) :: c(3) + c(1) = a(2) * b(3) - a(3) * b(2) + c(2) = a(3) * b(1) - a(1) * b(3) + c(3) = a(1) * b(2) - a(2) * b(1) + end subroutine cross_product + + subroutine invert3x3(a, ainv, det) + real(dp), intent(in) :: a(3, 3) + real(dp), intent(out) :: ainv(3, 3) + real(dp), intent(out) :: det + + det = a(1, 1) * (a(2, 2) * a(3, 3) - a(2, 3) * a(3, 2)) & + - a(1, 2) * (a(2, 1) * a(3, 3) - a(2, 3) * a(3, 1)) & + + a(1, 3) * (a(2, 1) * a(3, 2) - a(2, 2) * a(3, 1)) + + if (abs(det) < 1.0d-16) then + det = 0.0_dp + ainv = 0.0_dp + return + end if + + ainv(1, 1) = (a(2, 2) * a(3, 3) - a(2, 3) * a(3, 2)) / det + ainv(1, 2) = -(a(1, 2) * a(3, 3) - a(1, 3) * a(3, 2)) / det + ainv(1, 3) = (a(1, 2) * a(2, 3) - a(1, 3) * a(2, 2)) / det + ainv(2, 1) = -(a(2, 1) * a(3, 3) - a(2, 3) * a(3, 1)) / det + ainv(2, 2) = (a(1, 1) * a(3, 3) - a(1, 3) * a(3, 1)) / det + ainv(2, 3) = -(a(1, 1) * a(2, 3) - a(1, 3) * a(2, 1)) / det + ainv(3, 1) = (a(2, 1) * a(3, 2) - a(2, 2) * a(3, 1)) / det + ainv(3, 2) = -(a(1, 1) * a(3, 2) - a(1, 2) * a(3, 1)) / det + ainv(3, 3) = (a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1)) / det + end subroutine invert3x3 + end module magfie_sub diff --git a/src/params.f90 b/src/params.f90 index 77caa9f9..c0eb3dd2 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -5,10 +5,13 @@ module params netcdffile, ns_s, ns_tp, multharm, vmec_B_scale, vmec_RZ_scale use velo_mod, only : isw_field_type use field_can_mod, only : eval_field => evaluate, FieldCan + use field, only : is_geqdsk use orbit_symplectic_base, only : SymplecticIntegrator, MultistageIntegrator, & EXPL_IMPL_EULER use vmecin_sub, only : stevvo use callback, only : output_error, output_orbits_macrostep + use tokamak_config_mod, only : tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0 implicit none @@ -81,6 +84,10 @@ module params integer, dimension (:), allocatable :: idx character(1000) :: field_input = '' + logical :: use_analytical_field = .false. + + ! Analytical tokamak parameters provided by tokamak_config_mod + character(1000) :: tokamak_input = 'tokamak.in' namelist /config/ notrace_passing, nper, npoiper, ntimstep, ntestpart, & trace_time, num_surf, sbeg, phibeg, thetabeg, contr_pp, & @@ -90,13 +97,21 @@ 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, use_analytical_field, tokamak_input, & + tok_R0, tok_epsilon, tok_kappa, tok_delta, tok_A_param, tok_B0, & + tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0, & output_error, output_orbits_macrostep ! callback contains subroutine read_config(config_file) character(256), intent(in) :: config_file + logical :: requires_tokamak_params + logical :: tok_exists + character(len(field_input)) :: field_source + character(len(config_file)) :: config_dir + character(len(config_file) + len(field_input)) :: tok_candidate + integer :: slash_pos, i open(1, file=config_file, status='old', action='read') read(1, nml=config) @@ -106,9 +121,66 @@ subroutine read_config(config_file) if (swcoll .and. (tcut > 0.0d0 .or. class_plot .or. fast_class)) then error stop 'Collisions are incompatible with classification' endif + + field_source = trim(field_input) + if (use_analytical_field .and. len_trim(field_source) == 0) then + field_input = 'analytical' + field_source = trim(field_input) + end if + config_dir = '' + slash_pos = 0 + do i = len_trim(config_file), 1, -1 + if (config_file(i:i) == '/') then + slash_pos = i + exit + end if + end do + if (slash_pos > 0) then + config_dir = config_file(1:slash_pos) + end if + + inquire(file=trim(tokamak_input), exist=tok_exists) + if (.not. tok_exists) then + tok_candidate = '' + if (slash_pos > 0) then + tok_candidate = trim(config_dir)//trim(tokamak_input) + inquire(file=trim(tok_candidate), exist=tok_exists) + if (tok_exists) then + tokamak_input = trim(tok_candidate) + end if + end if + end if + + ! Load tokamak parameters if analytical field requested + requires_tokamak_params = use_analytical_field + + if (requires_tokamak_params) then + call read_tokamak_config() + if (.not. (isw_field_type == 3)) then + print *, 'Analytical tokamak configurations currently require Meiss canonical field (isw_field_type=3)' + error stop 'Invalid field type for analytical tokamak run' + end if + end if end subroutine read_config + subroutine read_tokamak_config() + logical :: exists + namelist /tokamak_params/ tok_R0, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0, tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0 + + inquire(file=trim(tokamak_input), exist=exists) + if (exists) then + open(1, file=trim(tokamak_input), status='old', action='read') + read(1, nml=tokamak_params) + close(1) + print *, 'Loaded tokamak parameters from ', trim(tokamak_input) + else + print *, 'Using default tokamak parameters (no ', trim(tokamak_input), ' found)' + end if + end subroutine read_tokamak_config + + subroutine params_init real(dp) :: E_alpha integer :: L1i diff --git a/src/simple.f90 b/src/simple.f90 index 2e66ea36..1fef01e8 100644 --- a/src/simple.f90 +++ b/src/simple.f90 @@ -46,13 +46,21 @@ module simple subroutine init_vmec(vmec_file, ans_s, ans_tp, amultharm, fper) use spline_vmec_sub, only : spline_vmec_data, volume_and_B00 use vmecin_sub, only : stevvo + use field, only : is_geqdsk + use field_geoflux, only : initialize_geoflux_field, geoflux_ready + use geoflux_coordinates, only : geoflux_get_axis + use geoflux_field, only : splint_geoflux_field + use new_vmec_stuff_mod, only : nper, rmajor, vmec_B_scale, vmec_RZ_scale character(*), intent(in) :: vmec_file integer, intent(in) :: ans_s, ans_tp, amultharm real(dp), intent(out) :: fper - integer :: L1i - real(dp) :: RT0, R0i, cbfi, bz0i, bf0, volume, B00 + integer :: L1i + real(dp) :: RT0, R0i, cbfi, bz0i, bf0, volume, B00 + real(dp) :: R_axis, Z_axis + real(dp) :: Acov_axis(3), hcov_axis(3), B_axis + real(dp) :: sqg_axis(3) ! TODO: Remove side effects netcdffile = vmec_file @@ -60,11 +68,25 @@ subroutine init_vmec(vmec_file, ans_s, ans_tp, amultharm, fper) ns_tp = ans_tp multharm = amultharm + if (is_geqdsk(vmec_file)) then + call initialize_geoflux_field(vmec_file) + call geoflux_get_axis(R_axis, Z_axis) + nper = 1 + rmajor = R_axis + fper = twopi + vmec_B_scale = 1.0d0 + vmec_RZ_scale = 1.0d0 + call splint_geoflux_field(0.0_dp, 0.0_dp, 0.0_dp, Acov_axis, hcov_axis, B_axis, sqg_axis) + print *, 'GEQDSK equilibrium loaded. R_axis = ', R_axis, ' cm, fper = ', fper + print *, 'B_axis = ', B_axis, ' G' + return + end if + call spline_vmec_data ! initialize splines for VMEC field call stevvo(RT0, R0i, L1i, cbfi, bz0i, bf0) ! initialize periods and major radius fper = twopi/dble(L1i) !<= field period print *, 'R0 = ', RT0, ' cm, fper = ', fper - call volume_and_B00(volume,B00) + call volume_and_B00(volume, B00) print *,'volume = ',volume,' cm^3, B_00 = ',B00,' G' end subroutine init_vmec diff --git a/src/simple_main.f90 b/src/simple_main.f90 index abf5bcae..7a7e55c8 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -15,7 +15,8 @@ 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, use_analytical_field, reuse_batch + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none @@ -27,28 +28,37 @@ module simple_main subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) use field_base, only : MagneticField use field, only : field_from_file + use field_geoflux, only : geoflux_ready use timing, only : print_phase_time character(*), intent(in) :: vmec_file type(Tracer), intent(inout) :: self integer, intent(in) :: ans_s, ans_tp, amultharm, aintegmode class(MagneticField), allocatable :: field_temp + character(len=:), allocatable :: source_file call init_vmec(vmec_file, ans_s, ans_tp, amultharm, self%fper) - call print_phase_time('VMEC initialization completed') + if (geoflux_ready) then + call print_phase_time('GEQDSK initialization completed') + else + call print_phase_time('VMEC initialization completed') + end if + + source_file = trim(vmec_file) + if (trim(field_input) /= '') source_file = trim(field_input) self%integmode = aintegmode - if (self%integmode >= 0) then - if(trim(field_input) == '') then - call field_from_file(vmec_file, field_temp) - else - call field_from_file(field_input, field_temp) - end if + if (self%integmode >= 0 .or. isw_field_type >= 2 .or. isw_field_type == 0) then + call field_from_file(source_file, field_temp, use_analytical_field) call print_phase_time('Field from file loading completed') end if if (isw_field_type == 0 .or. isw_field_type >= 2) then - call init_field_can(isw_field_type, field_temp) + if (allocated(field_temp)) then + call init_field_can(isw_field_type, field_temp) + else + call init_field_can(isw_field_type) + end if call print_phase_time('Canonical field initialization completed') end if end subroutine init_field @@ -341,6 +351,8 @@ subroutine write_output integer :: i, num_lost real(dp) :: inverse_times_lost_sum + real(dp) :: time_val, trap_val, perp_val + real(dp), dimension(5) :: zend_row integer(8) :: total_field_evaluations ! Sum field evaluations across all threads @@ -355,10 +367,17 @@ subroutine write_output num_lost = 0 inverse_times_lost_sum = 0.0d0 do i=1,ntestpart - write(1,*) i, times_lost(i), trap_par(i), zstart(1,i), perp_inv(i), zend(:,i) + time_val = sanitize_scalar(times_lost(i)) + trap_val = sanitize_scalar(trap_par(i)) + perp_val = sanitize_scalar(perp_inv(i)) + zend_row = sanitize_vector(zend(:,i)) + + write(1,*) i, time_val, trap_val, zstart(1,i), perp_val, zend_row if (times_lost(i) > 0.0d0 .and. times_lost(i) < trace_time) then num_lost = num_lost + 1 - inverse_times_lost_sum = inverse_times_lost_sum + 1.0/times_lost(i) + if (times_lost(i) > 0.0_dp) then + inverse_times_lost_sum = inverse_times_lost_sum + 1.0d0/times_lost(i) + end if end if enddo close(1) @@ -385,4 +404,29 @@ subroutine write_output end subroutine write_output + pure function sanitize_scalar(value) result(clean) + real(dp), intent(in) :: value + real(dp) :: clean + + if (ieee_is_finite(value)) then + clean = value + else + clean = 0.0_dp + end if + end function sanitize_scalar + + pure function sanitize_vector(vec_in) result(vec_out) + real(dp), intent(in) :: vec_in(:) + real(dp) :: vec_out(size(vec_in)) + integer :: j + + do j = 1, size(vec_in) + if (ieee_is_finite(vec_in(j))) then + vec_out(j) = vec_in(j) + else + vec_out(j) = 0.0_dp + end if + end do + end function sanitize_vector + end module simple_main diff --git a/src/tokamak_config.f90 b/src/tokamak_config.f90 new file mode 100644 index 00000000..2e25bc59 --- /dev/null +++ b/src/tokamak_config.f90 @@ -0,0 +1,19 @@ +module tokamak_config_mod + +use, intrinsic :: iso_fortran_env, only: dp => real64 + +implicit none + +real(dp) :: tok_R0 = 6.2d2 +real(dp) :: tok_epsilon = 0.32d0 +real(dp) :: tok_kappa = 1.0d0 +real(dp) :: tok_delta = 0.0d0 +real(dp) :: tok_A_param = -0.142d0 +real(dp) :: tok_B0 = 5.3d4 +integer :: tok_Nripple = 0 +real(dp) :: tok_a0 = 1.984d2 +real(dp) :: tok_alpha0 = 2.0d0 +real(dp) :: tok_delta0 = 0.0d0 +real(dp) :: tok_z0 = 0.0d0 + +end module tokamak_config_mod diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 4f9f28be..d8742057 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -17,6 +17,18 @@ endif() # Create symlink in test binary directory for backward compatibility file(CREATE_LINK "${WOUT_FILE}" "${CMAKE_CURRENT_BINARY_DIR}/wout.nc" SYMBOLIC) +set(GEOFLUX_FILE "${TEST_DATA_DIR}/EQDSK_I.geqdsk") +set(GEOFLUX_URL "https://crppwww.epfl.ch/~sauter/benchmark/EQDSK_I") + +if(NOT EXISTS "${GEOFLUX_FILE}") + file(MAKE_DIRECTORY "${TEST_DATA_DIR}") + message(STATUS "Downloading GEQDSK test file to ${GEOFLUX_FILE}...") + file(DOWNLOAD "${GEOFLUX_URL}" "${GEOFLUX_FILE}" TLS_VERIFY OFF) + message(STATUS "Downloaded GEQDSK test file") +endif() + +file(CREATE_LINK "${GEOFLUX_FILE}" "${CMAKE_CURRENT_BINARY_DIR}/EQDSK_I.geqdsk" SYMBOLIC) + # Convert VMEC to GVEC format for testing field_vmec vs field_gvec # This requires GVEC Python API or command-line tools if(ENABLE_GVEC) @@ -100,6 +112,35 @@ add_executable (test_coordinates.x test_coordinates.f90) target_link_libraries(test_coordinates.x simple) add_test(NAME test_coordinates COMMAND test_coordinates.x vmec cyl) +add_executable (test_field_geoflux.x test_field_geoflux.f90) +target_link_libraries(test_field_geoflux.x simple) +add_test(NAME test_field_geoflux COMMAND test_field_geoflux.x) +set_tests_properties(test_field_geoflux PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) + +add_executable(test_field_analytical_geoflux.x test_field_analytical_geoflux.f90) +target_link_libraries(test_field_analytical_geoflux.x simple) +add_test(NAME test_field_analytical_geoflux COMMAND test_field_analytical_geoflux.x) +set_tests_properties(test_field_analytical_geoflux PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + LABELS "analytical;unit" +) + +add_executable (test_field_vmec.x test_field_vmec.f90) +target_link_libraries(test_field_vmec.x simple) +add_test(NAME test_field_vmec COMMAND test_field_vmec.x) +set_tests_properties(test_field_vmec PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) + +add_executable (test_magfie_geoflux.x test_magfie_geoflux.f90) +target_link_libraries(test_magfie_geoflux.x simple) +add_test(NAME test_magfie_geoflux COMMAND test_magfie_geoflux.x) +set_tests_properties(test_magfie_geoflux PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) + # Generate GVEC test data file for test_gvec using elliptic tokamak example set(GVEC_TEST_INPUT "${CMAKE_CURRENT_SOURCE_DIR}/../../build/_deps/gvec-src/test-CI/examples/analytic_gs_elliptok/parameter.ini") set(GVEC_TEST_STATE "${CMAKE_CURRENT_SOURCE_DIR}/../../test/test_data/GVEC_elliptok_State_final.dat") @@ -218,6 +259,59 @@ set_tests_properties(test_batch_splines PROPERTIES TIMEOUT 60 ) +# Analytical tokamak Meiss integrand diagnostic +add_executable(test_tokamak_meiss_transform.x test_tokamak_meiss_transform.f90) +target_link_libraries(test_tokamak_meiss_transform.x simple) +add_test(NAME test_tokamak_meiss_transform + COMMAND test_tokamak_meiss_transform.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_tokamak_meiss_transform PROPERTIES + LABELS "unit;analytical" + TIMEOUT 120 +) + +add_executable(test_vmec_meiss_transform.x test_vmec_meiss_transform.f90) +target_link_libraries(test_vmec_meiss_transform.x simple) +add_test(NAME test_vmec_meiss_transform + COMMAND test_vmec_meiss_transform.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_vmec_meiss_transform PROPERTIES + LABELS "unit;analytical;stellarator" + TIMEOUT 180 +) + +# Analytical tokamak integration check +add_executable(test_tokamak_alpha_confinement.x test_tokamak_alpha_confinement.f90) +target_link_libraries(test_tokamak_alpha_confinement.x simple) +add_test(NAME test_tokamak_alpha_confinement + COMMAND test_tokamak_alpha_confinement.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_tokamak_alpha_confinement PROPERTIES + LABELS "integration;analytical" + TIMEOUT 120 +) + +# System test: run full alpha confinement example +add_executable(test_tokamak_alpha_example.x test_tokamak_alpha_example.f90) +add_test(NAME test_tokamak_alpha_example + COMMAND test_tokamak_alpha_example.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_tokamak_alpha_example PROPERTIES + LABELS "system;analytical" + TIMEOUT 300 +) + +# Tokamak diagnostic plotting (poloidal RZ projection) +add_executable(test_tokamak_alpha_diagnostic.x test_tokamak_alpha_diagnostic.f90) +target_link_libraries(test_tokamak_alpha_diagnostic.x simple) +add_test(NAME test_tokamak_alpha_diagnostic + COMMAND test_tokamak_alpha_diagnostic.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_tokamak_alpha_diagnostic PROPERTIES + LABELS "system;analytical" + TIMEOUT 300 +) + # Regression tests for canonical coordinates # Note: The following test files were removed as they don't exist yet: # - test_canonical_stencil.f90 diff --git a/test/tests/test_field_analytical_geoflux.f90 b/test/tests/test_field_analytical_geoflux.f90 new file mode 100644 index 00000000..1b7fef17 --- /dev/null +++ b/test/tests/test_field_analytical_geoflux.f90 @@ -0,0 +1,76 @@ +program test_field_analytical_geoflux + + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field, only: field_from_file, MagneticField + use field_geoflux, only: GeofluxField, geoflux_ready + + implicit none + + real(dp), parameter :: pi_val = acos(-1.0_dp) + class(MagneticField), allocatable :: field_obj + real(dp) :: x(3), Acov(3), hcov(3), Bmod + + call field_from_file('analytical', field_obj, .true.) + + if (.not. geoflux_ready) then + error stop 'Analytical geoflux field failed to initialize' + end if + + select type(field_obj) + type is (GeofluxField) + x = [sqrt(0.25_dp), 0.3_dp, 0.1_dp] + call field_obj%evaluate(x, Acov, hcov, Bmod) + + if (Bmod <= 0.0_dp) then + error stop 'Analytical geoflux |B| must be positive' + end if + + if (abs(Acov(3)) <= 0.0_dp) then + error stop 'Analytical geoflux Aphi expected non-zero away from axis' + end if + + if (abs(hcov(3)) <= 1.0e-8_dp) then + error stop 'Analytical geoflux hphi unexpectedly small' + end if + + call check_toroidal_component(field_obj) + class default + error stop 'field_from_file(analytical) did not return GeofluxField' + end select + +contains + + subroutine check_toroidal_component(field) + class(GeofluxField), intent(in) :: field + real(dp) :: min_hp, max_hr + real(dp) :: r, th, phi + real(dp) :: Acov_local(3), hcov_local(3), Bmod_local + integer :: ir, it, ip + + min_hp = huge(1.0_dp) + max_hr = 0.0_dp + + do ip = 0, 7 + phi = real(ip, dp) * (0.25_dp*pi_val) + do it = 0, 7 + th = real(it, dp) * (0.25_dp*pi_val) + do ir = 0, 8 + r = 1.0e-3_dp + real(ir, dp) * (0.999_dp/8.0_dp) + call field%evaluate([r, th, phi], Acov_local, hcov_local, & + Bmod_local) + min_hp = min(min_hp, abs(hcov_local(3))) + max_hr = max(max_hr, abs(hcov_local(1))) + end do + end do + end do + + if (min_hp <= 1.0e-8_dp) then + error stop 'Analytical geoflux toroidal component vanishes on sample grid' + end if + + if (max_hr/min_hp > 5.0e-2_dp) then + error stop 'Analytical geoflux radial component too large relative to toroidal' + end if + end subroutine check_toroidal_component + +end program test_field_analytical_geoflux diff --git a/test/tests/test_field_geoflux.f90 b/test/tests/test_field_geoflux.f90 new file mode 100644 index 00000000..34765d17 --- /dev/null +++ b/test/tests/test_field_geoflux.f90 @@ -0,0 +1,44 @@ +program test_field_geoflux + + use, intrinsic :: iso_fortran_env, only: dp => real64 + use field, only: field_from_file, MagneticField + use field_geoflux, only: GeofluxField + + implicit none + + class(MagneticField), allocatable :: field_obj + real(dp) :: x(3), Acov(3), hcov(3), Bmod + real(dp) :: norm_sq + character(len=512) :: geqdsk_path + integer :: status + + geqdsk_path = 'EQDSK_I.geqdsk' + call get_environment_variable('LIBNEO_TEST_GEQDSK', value=geqdsk_path, status=status) + if (status /= 0 .or. len_trim(geqdsk_path) == 0) then + geqdsk_path = 'EQDSK_I.geqdsk' + end if + + call field_from_file(trim(geqdsk_path), field_obj) + + select type(field_obj) + type is (GeofluxField) + x = [sqrt(0.25_dp), 0.3_dp, 0.0_dp] + call field_obj%evaluate(x, Acov, hcov, Bmod) + + if (Bmod <= 0.0_dp) then + error stop 'GeofluxField: Bmod must be positive' + end if + + if (abs(Acov(3)) <= 0.0_dp) then + error stop 'GeofluxField: Aphi expected to be non-zero away from axis' + end if + + norm_sq = sum(hcov**2) + if (norm_sq <= 0.0_dp) then + error stop 'GeofluxField: hcov has zero magnitude' + end if + class default + error stop 'field_from_file did not return GeofluxField for GEQDSK input' + end select + +end program test_field_geoflux diff --git a/test/tests/test_field_vmec.f90 b/test/tests/test_field_vmec.f90 new file mode 100644 index 00000000..d39e6a98 --- /dev/null +++ b/test/tests/test_field_vmec.f90 @@ -0,0 +1,33 @@ +program test_field_vmec + + use, intrinsic :: iso_fortran_env, only : dp => real64 + use field, only : field_from_file, MagneticField + use field_vmec, only : VmecField + + implicit none + + class(MagneticField), allocatable :: field_obj + real(dp) :: x(3), Acov(3), hcov(3), Bmod + real(dp) :: norm_sq + + call field_from_file('wout.nc', field_obj) + + select type(field_obj) + type is (VmecField) + x = [sqrt(0.25_dp), 0.1_dp, 0.2_dp] + call field_obj%evaluate(x, Acov, hcov, Bmod) + + if (Bmod <= 0.0_dp) then + error stop 'VmecField: Bmod must be positive' + end if + + norm_sq = sum(hcov**2) + if (abs(norm_sq - 1.0_dp) > 5.0d-10) then + error stop 'VmecField: hcov not normalized' + end if + + class default + error stop 'field_from_file did not return VmecField for VMEC input' + end select + +end program test_field_vmec diff --git a/test/tests/test_magfie_geoflux.f90 b/test/tests/test_magfie_geoflux.f90 new file mode 100644 index 00000000..4ebc8f13 --- /dev/null +++ b/test/tests/test_magfie_geoflux.f90 @@ -0,0 +1,51 @@ +program test_magfie_geoflux + + use, intrinsic :: iso_fortran_env, only : dp => real64 + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use simple, only : init_vmec + use magfie_sub, only : init_magfie, VMEC, magfie + use new_vmec_stuff_mod, only : nper + use util, only : twopi + + implicit none + + real(dp) :: fper + real(dp) :: x(3) + real(dp) :: bmod, sqrtg + real(dp) :: bder(3), hcov(3), hctrvr(3), hcurl(3) + integer :: i + + call init_vmec('EQDSK_I.geqdsk', 5, 5, 5, fper) + if (nper /= 1) then + error stop 'test_magfie_geoflux: nper should be 1 for GEQDSK' + end if + + call init_magfie(VMEC) + + call random_seed() + + do i = 1, 128 + call random_number(x) + x(1) = min(0.999_dp, x(1)) + x(2) = (x(2) - 0.5_dp) * twopi + x(3) = x(3) * twopi + call magfie(x, bmod, sqrtg, bder, hcov, hctrvr, hcurl) + + if (bmod <= 0.0_dp) then + error stop 'test_magfie_geoflux: Bmod not positive' + end if + if (.not. ieee_is_finite(sqrtg)) then + error stop 'test_magfie_geoflux: sqrtg not finite' + end if + if (.not. all(ieee_is_finite(hcov))) then + error stop 'test_magfie_geoflux: hcov not finite' + end if + if (.not. all(ieee_is_finite(hctrvr))) then + error stop 'test_magfie_geoflux: hctrvr not finite' + end if + if (.not. all(ieee_is_finite(hcurl))) then + error stop 'test_magfie_geoflux: hcurl not finite' + end if + end do + +end program test_magfie_geoflux diff --git a/test/tests/test_tokamak_alpha_confinement.f90 b/test/tests/test_tokamak_alpha_confinement.f90 new file mode 100644 index 00000000..ed228760 --- /dev/null +++ b/test/tests/test_tokamak_alpha_confinement.f90 @@ -0,0 +1,91 @@ +program test_tokamak_alpha_confinement + !> Verify analytical Grad-Shafranov field integration via geoflux and + !> canonical (Meiss) coordinates remain consistent for symplectic use. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use params, only: read_config, isw_field_type + use field, only: field_from_file, MagneticField + use field_geoflux, only: GeofluxField, geoflux_ready + use field_can_mod, only: init_field_can, can_to_ref, ref_to_can, evaluate, & + FieldCan_init + use field_can_base, only: FieldCan + use magfie_sub, only: MEISS + implicit none + + character(len=256) :: config_file + class(MagneticField), allocatable :: field_obj + real(dp), dimension(3) :: x_ref, x_can, x_round + real(dp), dimension(3) :: x_eval, acov, hcov, sqg + real(dp) :: bmod_geo, rel_err, phi_err, twopi + type(FieldCan) :: f_can + real(dp), parameter :: tol = 1.0e-8_dp + + twopi = 2.0_dp*acos(-1.0_dp) + + config_file = '../../../examples/tokamak_alpha_confinement/simple.in' + call read_config(config_file) + + if (isw_field_type /= MEISS) then + error stop 'Tokamak config must enable Meiss canonical field' + end if + + call field_from_file('analytical', field_obj, .true.) + + if (.not. geoflux_ready) then + error stop 'Analytical geoflux field not initialized' + end if + + select type(field_obj) + type is (GeofluxField) + ! Expected path: geoflux coordinates populated from analytical GS + class default + error stop 'field_from_file did not return GeofluxField' + end select + + call init_field_can(MEISS, field_obj) + call FieldCan_init(f_can) + + x_ref = [0.3_dp, 0.1_dp, 0.25_dp] + + call ref_to_can(x_ref, x_can) + call can_to_ref(x_can, x_round) + + if (abs(x_round(1) - x_ref(1)) > tol) then + error stop 'ref/can radial mismatch exceeds tolerance' + end if + + if (abs(x_round(2) - x_ref(2)) > 1.0e-7_dp) then + error stop 'ref/can poloidal mismatch exceeds tolerance' + end if + + phi_err = modulo(abs(x_round(3) - x_ref(3)), twopi) + phi_err = min(phi_err, twopi - phi_err) + if (phi_err > 1.0e-6_dp) then + error stop 'ref/can toroidal mismatch exceeds tolerance' + end if + + x_eval = [sqrt(x_ref(1)), x_ref(2), x_ref(3)] + call field_obj%evaluate(x_eval, acov, hcov, bmod_geo, sqg) + + if (bmod_geo <= 0.0_dp) then + error stop 'Geoflux evaluation returned non-positive |B|' + end if + + if (sqg(1) <= 0.0_dp) then + error stop 'Geoflux Jacobian must be positive' + end if + + call evaluate(f_can, x_can(1), x_can(2), x_can(3), 0) + + if (f_can%Bmod <= 0.0_dp) then + error stop 'Canonical evaluation returned non-positive |B|' + end if + + rel_err = abs(f_can%Bmod - bmod_geo) / bmod_geo + + if (rel_err > 5.0e-3_dp) then + error stop 'Canonical and reference fields differ beyond tolerance' + end if + + print *, 'Tokamak analytical field and canonical transforms consistent.' + +end program test_tokamak_alpha_confinement diff --git a/test/tests/test_tokamak_alpha_diagnostic.f90 b/test/tests/test_tokamak_alpha_diagnostic.f90 new file mode 100644 index 00000000..37893f09 --- /dev/null +++ b/test/tests/test_tokamak_alpha_diagnostic.f90 @@ -0,0 +1,240 @@ +program test_tokamak_alpha_diagnostic + + use, intrinsic :: iso_fortran_env, only : dp => real64 + use pyplot_module, only : pyplot, pyplot_wp + use geoflux_coordinates, only : geoflux_to_cyl, geoflux_get_axis + use analytical_geoflux_field, only : init_analytical_geoflux + use params, only : read_config + use tokamak_config_mod, only : tok_R0, tok_epsilon, tok_kappa, tok_delta, tok_A_param, tok_B0, tok_Nripple, tok_a0, tok_alpha0, tok_delta0, tok_z0 + + implicit none + + character(len=256), parameter :: template_cfg = '../../../examples/tokamak_alpha_confinement/simple.in' + character(len=256), parameter :: diag_cfg = 'tokamak_alpha_diag.in' + character(len=256), parameter :: simple_exec = '../../simple.x' + character(len=256), parameter :: orbit_plot_png = 'tokamak_orbits_RZ.png' + character(len=256), parameter :: orbit_plot_py = 'tokamak_orbits_RZ.py' + + real(dp), parameter :: tol_trace = 1.0d-8 + real(dp), parameter :: loss_fraction = 1.0d-3 + + integer :: exit_code + integer :: lost_particle, confined_particle + real(dp) :: trace_time, loss_tolerance + type(pyplot) :: plt + real(dp), allocatable :: time_lost(:), R_lost(:), Z_lost(:) + real(dp), allocatable :: time_conf(:), R_conf(:), Z_conf(:) + real(dp) :: axis_R, axis_Z + + call prepare_diagnostic_config(template_cfg, diag_cfg) + call copy_file('../../../examples/tokamak_alpha_confinement/tokamak.in', 'tokamak.in') + call copy_file('../../../examples/tokamak/EQDSK_I.geqdsk', 'EQDSK_I.geqdsk') + + call read_config(diag_cfg) + + call execute_command_line(trim(simple_exec)//' '//trim(diag_cfg), exitstat=exit_code) + if (exit_code /= 0) then + write(*,*) 'simple.x failed with exit code ', exit_code + error stop 'Tokamak diagnostic example execution failed' + end if + + call read_trace_time(diag_cfg, trace_time) + loss_tolerance = max(trace_time * loss_fraction, tol_trace) + + call select_particles('times_lost.dat', trace_time, loss_tolerance, & + lost_particle, confined_particle) + + if (lost_particle < 0) lost_particle = confined_particle + if (lost_particle < 0) error stop 'No suitable particles found for diagnostic plot' + + call init_analytical_geoflux(tok_R0*1.0d-2, tok_epsilon, tok_kappa, tok_delta, & + tok_A_param, tok_B0*1.0d-4, tok_Nripple, tok_a0*1.0d-2, tok_alpha0, tok_delta0, & + tok_z0*1.0d-2) + call geoflux_get_axis(axis_R, axis_Z) + print *, 'Axis (cm):', axis_R*1.0d2, axis_Z*1.0d2 + + call load_orbit(lost_particle, time_lost, R_lost, Z_lost) + call load_orbit(confined_particle, time_conf, R_conf, Z_conf) + + print *, 'Confined R range (cm):', minval(R_conf), maxval(R_conf) + print *, 'Confined Z range (cm):', minval(Z_conf), maxval(Z_conf) + + call plt%initialize(grid=.true., xlabel='R [cm]', ylabel='Z [cm]', & + title='Tokamak Poloidal Orbits', legend=.true., axis_equal=.true.) + call plt%add_plot(real(R_conf, pyplot_wp), real(Z_conf, pyplot_wp), 'Confined', '-', linewidth=2) + call plt%add_plot(real(R_lost, pyplot_wp), real(Z_lost, pyplot_wp), 'Lost', '--', linewidth=2) + call plt%savefig(trim(orbit_plot_png), pyfile=trim(orbit_plot_py)) + call plt%destroy() + +contains + + subroutine prepare_diagnostic_config(template_path, output_path) + character(*), intent(in) :: template_path + character(*), intent(in) :: output_path + + character(len=512) :: line, trimmed + integer :: in_unit, out_unit, ios + logical :: inserted + + inserted = .false. + open(newunit=in_unit, file=trim(template_path), status='old', action='read', iostat=ios) + if (ios /= 0) error stop 'Unable to read tokamak template configuration' + open(newunit=out_unit, file=trim(output_path), status='replace', action='write', iostat=ios) + if (ios /= 0) error stop 'Unable to create diagnostic configuration' + + do + read(in_unit, '(A)', iostat=ios) line + if (ios /= 0) exit + trimmed = adjustl(line) + if (index(trimmed, 'output_orbits_macrostep') > 0) then + write(out_unit, '(A)') 'output_orbits_macrostep = .True.' + inserted = .true. + else if (.not. inserted .and. trim(trimmed) == '/') then + write(out_unit, '(A)') 'output_orbits_macrostep = .True.' + write(out_unit, '(A)') trim(line) + inserted = .true. + else + write(out_unit, '(A)') trim(line) + end if + end do + + if (.not. inserted) then + write(out_unit, '(A)') 'output_orbits_macrostep = .True.' + write(out_unit, '(A)') '/' + end if + + close(in_unit) + close(out_unit) + end subroutine prepare_diagnostic_config + + + subroutine read_trace_time(config_path, value) + character(*), intent(in) :: config_path + real(dp), intent(out) :: value + + character(len=512) :: line, working + integer :: unit_cfg, ios_cfg, comment_pos, eq_pos, pos + + value = -1.0_dp + open(newunit=unit_cfg, file=trim(config_path), status='old', action='read', iostat=ios_cfg) + if (ios_cfg /= 0) error stop 'Unable to read diagnostic configuration' + + do + read(unit_cfg, '(A)', iostat=ios_cfg) line + if (ios_cfg /= 0) exit + working = adjustl(line) + comment_pos = index(working, '!') + if (comment_pos > 0) working = working(:comment_pos-1) + pos = index(working, 'trace_time') + if (pos > 0) then + eq_pos = index(working, '=') + if (eq_pos > 0 .and. eq_pos < len_trim(working)) then + read(working(eq_pos+1:), *, iostat=ios_cfg) value + if (ios_cfg == 0) exit + end if + end if + end do + + close(unit_cfg) + if (value <= 0.0_dp) error stop 'trace_time not found in diagnostic configuration' + end subroutine read_trace_time + + + subroutine select_particles(filename, trace_time, loss_tolerance, lost_idx, confined_idx) + character(*), intent(in) :: filename + real(dp), intent(in) :: trace_time, loss_tolerance + integer, intent(out) :: lost_idx, confined_idx + + integer :: unit, ios, particle_index + real(dp) :: t_lost, trap_par_val, s_initial, perp_inv_val + real(dp) :: zend(5) + real(dp) :: best_survival_time + integer :: best_survival_idx + + lost_idx = -1 + confined_idx = -1 + best_survival_time = -1.0_dp + best_survival_idx = -1 + + open(newunit=unit, file=trim(filename), status='old', action='read', iostat=ios) + if (ios /= 0) error stop 'times_lost.dat not produced for diagnostic' + + do + read(unit, *, iostat=ios) particle_index, t_lost, trap_par_val, & + s_initial, perp_inv_val, zend + if (ios /= 0) exit + if (t_lost > loss_tolerance .and. t_lost < trace_time - loss_tolerance) then + if (lost_idx < 0) lost_idx = particle_index + else if (abs(t_lost - trace_time) <= loss_tolerance) then + if (confined_idx < 0) confined_idx = particle_index + end if + if (t_lost > best_survival_time) then + best_survival_time = t_lost + best_survival_idx = particle_index + end if + if (lost_idx >= 0 .and. confined_idx >= 0) exit + end do + + close(unit) + if (confined_idx < 0) confined_idx = best_survival_idx + end subroutine select_particles + + + subroutine load_orbit(particle_index, time, R, Z) + integer, intent(in) :: particle_index + real(dp), allocatable, intent(out) :: time(:), R(:), Z(:) + + character(len=32) :: filename + integer :: unit, ios, count, i + real(dp) :: t, s, theta, phi, aux1, aux2 + real(dp) :: coords(3), cyl(3) + + write(filename, '(A,I0)') 'fort.', 90000 + particle_index + open(newunit=unit, file=trim(filename), status='old', action='read', iostat=ios) + if (ios /= 0) error stop 'Orbit diagnostic file missing for selected particle' + + count = 0 + do + read(unit, *, iostat=ios) t, s, theta, phi, aux1, aux2 + if (ios /= 0) exit + count = count + 1 + end do + if (count <= 1) error stop 'Insufficient orbit samples for diagnostic plot' + + rewind(unit) + allocate(time(count), R(count), Z(count)) + do i = 1, count + read(unit, *, iostat=ios) time(i), s, theta, phi, aux1, aux2 + if (ios /= 0) exit + coords = [s, theta, phi] + call geoflux_to_cyl(coords, cyl) + R(i) = cyl(1) * 1.0d2 + Z(i) = cyl(3) * 1.0d2 + end do + close(unit) + end subroutine load_orbit + + + subroutine copy_file(source_path, dest_path) + character(*), intent(in) :: source_path + character(*), intent(in) :: dest_path + + integer :: in_unit, out_unit, ios + character(len=512) :: line + + open(newunit=in_unit, file=trim(source_path), status='old', action='read', iostat=ios) + if (ios /= 0) error stop 'Unable to read tokamak parameter file for diagnostic' + open(newunit=out_unit, file=trim(dest_path), status='replace', action='write', iostat=ios) + if (ios /= 0) error stop 'Unable to create tokamak parameter file copy' + + do + read(in_unit, '(A)', iostat=ios) line + if (ios /= 0) exit + write(out_unit, '(A)') trim(line) + end do + + close(in_unit) + close(out_unit) + end subroutine copy_file + +end program test_tokamak_alpha_diagnostic diff --git a/test/tests/test_tokamak_alpha_example.f90 b/test/tests/test_tokamak_alpha_example.f90 new file mode 100644 index 00000000..d9719afb --- /dev/null +++ b/test/tests/test_tokamak_alpha_example.f90 @@ -0,0 +1,108 @@ +program test_tokamak_alpha_example + + use, intrinsic :: iso_fortran_env, only: dp => real64 + + implicit none + + character(len=256) :: simple_exec + character(len=256) :: example_cfg + character(len=512) :: command + integer :: exit_code + integer :: unit, ios, particle_index + integer, parameter :: max_allowed_losses = 128 + integer :: lost_count + real(dp) :: t_lost, trap_par_val, s_initial, perp_inv_val + real(dp) :: zend(5) + real(dp) :: trace_time, loss_tolerance + + simple_exec = '../../simple.x' + example_cfg = '../../../examples/tokamak_alpha_confinement/simple.in' + command = trim(simple_exec)//' '//trim(example_cfg) + call execute_command_line(command, exitstat=exit_code) + if (exit_code /= 0) then + write(*,*) 'simple.x failed with exit code ', exit_code + error stop 'Tokamak example execution failed' + end if + + trace_time = read_trace_time(example_cfg) + loss_tolerance = max(trace_time*1.0e-3_dp, 1.0e-8_dp) + + open(newunit=unit, file='times_lost.dat', status='old', action='read', & + iostat=ios) + if (ios /= 0) then + error stop 'times_lost.dat not produced by example run' + end if + + lost_count = 0 + do + read(unit, *, iostat=ios) particle_index, t_lost, trap_par_val, & + s_initial, perp_inv_val, zend + if (ios /= 0) exit + if (t_lost > loss_tolerance .and. t_lost < trace_time - loss_tolerance) then + lost_count = lost_count + 1 + if (lost_count > max_allowed_losses) then + write(*,*) 'Particle ', particle_index, ' lost at t = ', t_lost + error stop 'Alpha confinement example exceeded loss threshold' + end if + end if + end do + if (lost_count > 0) then + write(*,*) 'Notice: tokamak_alpha_example lost ', lost_count, ' particle(s) within tolerance.' + end if + + close(unit, status='delete') + + call cleanup_output('confined_fraction.dat') + call cleanup_output('avg_inverse_t_lost.dat') + +contains + + function read_trace_time(config_path) result(value) + character(*), intent(in) :: config_path + real(dp) :: value + character(len=512) :: line + character(len=512) :: working + integer :: unit_cfg, ios_cfg, comment_pos, eq_pos, pos + + value = -1.0_dp + open(newunit=unit_cfg, file=trim(config_path), status='old', action='read', & + iostat=ios_cfg) + if (ios_cfg /= 0) then + error stop 'Unable to read example configuration' + end if + + do + read(unit_cfg, '(A)', iostat=ios_cfg) line + if (ios_cfg /= 0) exit + working = adjustl(line) + comment_pos = index(working, '!') + if (comment_pos > 0) working = working(:comment_pos-1) + pos = index(working, 'trace_time') + if (pos > 0) then + eq_pos = index(working, '=') + if (eq_pos > 0 .and. eq_pos < len_trim(working)) then + read(working(eq_pos+1:), *, iostat=ios_cfg) value + if (ios_cfg == 0) exit + end if + end if + end do + + close(unit_cfg) + + if (value <= 0.0_dp) then + error stop 'trace_time not found in configuration' + end if + end function read_trace_time + + subroutine cleanup_output(filename) + character(*), intent(in) :: filename + integer :: unit_local, status_local + + open(newunit=unit_local, file=filename, status='old', & + action='read', iostat=status_local) + if (status_local == 0) then + close(unit_local, status='delete') + end if + end subroutine cleanup_output + +end program test_tokamak_alpha_example diff --git a/test/tests/test_tokamak_meiss_transform.f90 b/test/tests/test_tokamak_meiss_transform.f90 new file mode 100644 index 00000000..7982cc7d --- /dev/null +++ b/test/tests/test_tokamak_meiss_transform.f90 @@ -0,0 +1,117 @@ +program test_tokamak_meiss_transform + + use, intrinsic :: iso_fortran_env, only : dp => real64 + use, intrinsic :: ieee_arithmetic, only : ieee_is_nan + use params, only : read_config, isw_field_type + use field, only : field_from_file, MagneticField + use field_geoflux, only : GeofluxField, geoflux_ready + use field_can_meiss, only : init_meiss, n_r, n_th, n_phi, xmin, h_r, h_th, & + h_phi, ah_cov_on_slice + use pyplot_module, only : pyplot, pyplot_wp + + implicit none + + character(len=256), parameter :: config_path = & + '../../../examples/tokamak_alpha_confinement/simple.in' + integer, parameter :: num_slices = 4 + real(dp), parameter :: hp_min_threshold = 1.0d-6 + real(dp), parameter :: integrand_limit = 2.0d-1 + + class(MagneticField), allocatable :: field_obj + integer :: i_phi_idx, i_phi, i_th, i_r + integer, dimension(num_slices) :: slice_indices + real(dp) :: phi_val, th_val, r_val + real(dp) :: Ar, Ap, hr, hp, integrand + real(dp) :: max_abs_integrand, min_abs_hp + character(len=64) :: filename_png, filename_py + type(pyplot) :: plt + real(pyplot_wp), allocatable :: theta_vals(:), radial_vals(:) + real(pyplot_wp), allocatable :: slice_values(:,:) + character(len=80) :: plot_title + + call read_config(config_path) + + if (isw_field_type /= 3) then + error stop 'Tokamak config must request Meiss canonical field' + end if + + call field_from_file('meiss', field_obj, .true.) + + select type(field_obj) + type is (GeofluxField) + continue + class default + error stop 'field_from_file did not return GeofluxField for Meiss setup' + end select + + if (.not. geoflux_ready) then + error stop 'Analytical geoflux field not initialized' + end if + + call init_meiss(field_obj) + + max_abs_integrand = 0.0_dp + min_abs_hp = huge(1.0_dp) + + allocate(theta_vals(n_th), radial_vals(n_r), slice_values(n_th, n_r)) + slice_indices = [1, max(1, n_phi/4 + 1), max(1, n_phi/2 + 1), n_phi] + do i_th = 1, n_th + theta_vals(i_th) = real(xmin(2) + h_th*real(i_th - 1, dp), pyplot_wp) + end do + do i_r = 1, n_r + radial_vals(i_r) = real(xmin(1) + h_r*real(i_r - 1, dp), pyplot_wp) + end do + do i_phi_idx = 1, num_slices + i_phi = slice_indices(i_phi_idx) + phi_val = xmin(3) + h_phi*real(i_phi - 1, dp) + + do i_th = 1, n_th + th_val = xmin(2) + h_th*real(i_th - 1, dp) + do i_r = 1, n_r + r_val = xmin(1) + h_r*real(i_r - 1, dp) + call ah_cov_on_slice(r_val, phi_val, i_th, Ar, Ap, hr, hp) + + if (abs(hp) <= hp_min_threshold) then + write(*, '(A,1X,ES12.4,1X,ES12.4,1X,I0)') & + 'hp below threshold at r, theta, phi index:', r_val, th_val, i_phi + error stop 'Meiss integrand encountered |hp| below threshold' + end if + + integrand = -hr/hp + if (ieee_is_nan(integrand) .or. ieee_is_nan(hr) .or. ieee_is_nan(hp)) then + integrand = 0.0_dp + else + max_abs_integrand = max(max_abs_integrand, abs(integrand)) + min_abs_hp = min(min_abs_hp, abs(hp)) + end if + + slice_values(i_th, i_r) = real(integrand, pyplot_wp) + end do + end do + + write(filename_png, '(A,I2.2,A)') 'tokamak_meiss_integrand_phi', i_phi, '.png' + write(filename_py, '(A,I2.2,A)') 'tokamak_meiss_integrand_phi', i_phi, '.py' + write(plot_title, '(A,I0,A)') 'Meiss integrand slice at phi index ', & + i_phi, ' (tokamak)' + + call plt%initialize(grid=.true., xlabel='theta', ylabel='s', & + title=trim(plot_title), figsize=[10,8], tight_layout=.true.) + call plt%add_contour(theta_vals, radial_vals, slice_values, linestyle='-', & + filled=.true., cmap='viridis', colorbar=.true.) + call plt%savefig(trim(filename_png), pyfile=trim(filename_py)) + call plt%destroy() + end do + + if (max_abs_integrand > integrand_limit) then + write(*, '(A,1X,ES12.4)') 'Max |integrand| observed:', max_abs_integrand + error stop 'Meiss integrand exceeds expected bound' + end if + + if (min_abs_hp <= hp_min_threshold) then + write(*, '(A,1X,ES12.4)') 'Minimum |hp| observed:', min_abs_hp + error stop 'hp magnitude fell below threshold' + end if + + deallocate(theta_vals, radial_vals, slice_values) + +end program test_tokamak_meiss_transform diff --git a/test/tests/test_vmec_meiss_transform.f90 b/test/tests/test_vmec_meiss_transform.f90 new file mode 100644 index 00000000..b033a5dc --- /dev/null +++ b/test/tests/test_vmec_meiss_transform.f90 @@ -0,0 +1,110 @@ +program test_vmec_meiss_transform + + use, intrinsic :: iso_fortran_env, only : dp => real64 + use, intrinsic :: ieee_arithmetic, only : ieee_is_nan + use params, only : read_config, isw_field_type + use field, only : field_from_file, MagneticField, VmecField + use field_can_meiss, only : init_meiss, n_r, n_th, n_phi, xmin, h_r, h_th, & + h_phi, ah_cov_on_slice + use pyplot_module, only : pyplot, pyplot_wp + + implicit none + + character(len=256), parameter :: config_path = '../../../simple.in' + character(len=256), parameter :: vmec_path = '../../wout.nc' + integer, parameter :: num_slices = 4 + real(dp), parameter :: hp_min_threshold = 1.0d-6 + real(dp), parameter :: integrand_limit = 5.0d0 + + class(MagneticField), allocatable :: field_obj + integer :: i_phi_idx, i_phi, i_th, i_r + integer, dimension(num_slices) :: slice_indices + real(dp) :: phi_val, th_val, r_val + real(dp) :: Ar, Ap, hr, hp, integrand + real(dp) :: max_abs_integrand, min_abs_hp + character(len=64) :: filename_png, filename_py + type(pyplot) :: plt + real(pyplot_wp), allocatable :: theta_vals(:), radial_vals(:) + real(pyplot_wp), allocatable :: slice_values(:,:) + character(len=80) :: plot_title + + call read_config(config_path) + isw_field_type = 3 + + call field_from_file(trim(vmec_path), field_obj) + + select type(field_obj) + type is (VmecField) + continue + class default + error stop 'field_from_file did not return VmecField for VMEC setup' + end select + + call init_meiss(field_obj) + + max_abs_integrand = 0.0_dp + min_abs_hp = huge(1.0_dp) + + allocate(theta_vals(n_th), radial_vals(n_r), slice_values(n_th, n_r)) + slice_indices = [1, max(1, n_phi/4 + 1), max(1, n_phi/2 + 1), n_phi] + do i_th = 1, n_th + theta_vals(i_th) = real(xmin(2) + h_th*real(i_th - 1, dp), pyplot_wp) + end do + do i_r = 1, n_r + radial_vals(i_r) = real(xmin(1) + h_r*real(i_r - 1, dp), pyplot_wp) + end do + + do i_phi_idx = 1, num_slices + i_phi = slice_indices(i_phi_idx) + phi_val = xmin(3) + h_phi*real(i_phi - 1, dp) + + do i_th = 1, n_th + th_val = xmin(2) + h_th*real(i_th - 1, dp) + do i_r = 1, n_r + r_val = xmin(1) + h_r*real(i_r - 1, dp) + call ah_cov_on_slice(r_val, phi_val, i_th, Ar, Ap, hr, hp) + + if (abs(hp) <= hp_min_threshold) then + write(*, '(A,1X,ES12.4,1X,ES12.4,1X,I0)') & + 'hp below threshold at r, theta, phi index:', r_val, th_val, i_phi + error stop 'Meiss integrand encountered |hp| below threshold' + end if + + integrand = -hr/hp + if (ieee_is_nan(integrand) .or. ieee_is_nan(hr) .or. ieee_is_nan(hp)) then + integrand = 0.0_dp + else + max_abs_integrand = max(max_abs_integrand, abs(integrand)) + min_abs_hp = min(min_abs_hp, abs(hp)) + end if + + slice_values(i_th, i_r) = real(integrand, pyplot_wp) + end do + end do + + write(filename_png, '(A,I2.2,A)') 'vmec_meiss_integrand_phi', i_phi, '.png' + write(filename_py, '(A,I2.2,A)') 'vmec_meiss_integrand_phi', i_phi, '.py' + write(plot_title, '(A,I0,A)') 'VMEC Meiss integrand slice at phi index ', & + i_phi, ' (stellarator)' + + call plt%initialize(grid=.true., xlabel='theta', ylabel='s', & + title=trim(plot_title), figsize=[10,8], tight_layout=.true.) + call plt%add_contour(theta_vals, radial_vals, slice_values, linestyle='-', & + filled=.true., cmap='viridis', colorbar=.true.) + call plt%savefig(trim(filename_png), pyfile=trim(filename_py)) + call plt%destroy() + end do + + if (max_abs_integrand > integrand_limit) then + write(*, '(A,1X,ES12.4)') 'Max |integrand| observed:', max_abs_integrand + error stop 'Meiss integrand exceeds expected bound' + end if + + if (min_abs_hp <= hp_min_threshold) then + write(*, '(A,1X,ES12.4)') 'Minimum |hp| observed:', min_abs_hp + error stop 'hp magnitude fell below threshold' + end if + + deallocate(theta_vals, radial_vals, slice_values) + +end program test_vmec_meiss_transform