From dc6146f5d2788b5cab26e802b770807cfaab678f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 10:35:35 +0200 Subject: [PATCH 01/23] Classify external-data CTests for local pipeline --- NEO-2-PAR/CMakeLists.txt | 20 ++++++++++++++++++++ NEO-2-PAR/ripple_solver.f90 | 3 +-- NEO-2-QL/CMakeLists.txt | 5 +++++ 3 files changed, 26 insertions(+), 2 deletions(-) diff --git a/NEO-2-PAR/CMakeLists.txt b/NEO-2-PAR/CMakeLists.txt index 341b6fdb..402ec02c 100644 --- a/NEO-2-PAR/CMakeLists.txt +++ b/NEO-2-PAR/CMakeLists.txt @@ -48,3 +48,23 @@ ADD_TEST(Test_par_cubic_splines-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_f ADD_TEST(Test_par_cubic_splines-Par ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "par_cubic_splines" "6") ADD_TEST(Test_par_general_bsplines-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "par_general_bsplines" "0") ADD_TEST(Test_par_general_bsplines-Par ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "par_general_bsplines" "6") + +set(EXTERNAL_DATA_TESTS + Test_Testcase1-Seq + Test_Testcase1-Par + Test_Relativistic_BeliaevBudker-Seq + Test_Relativistic_BeliaevBudker-Par + Test_Relativistic_BraamsKarney-Seq + Test_Relativistic_BraamsKarney-Par + Test_par_standard_polynomials-Seq + Test_par_standard_polynomials-Par + Test_par_quadratic_polynomials-Seq + Test_par_quadratic_polynomials-Par + Test_par_cubic_splines-Seq + Test_par_cubic_splines-Par + Test_par_general_bsplines-Seq + Test_par_general_bsplines-Par +) +set_tests_properties(${EXTERNAL_DATA_TESTS} PROPERTIES + LABELS "regression;external-data" +) diff --git a/NEO-2-PAR/ripple_solver.f90 b/NEO-2-PAR/ripple_solver.f90 index 7b32f0e2..69d60da1 100644 --- a/NEO-2-PAR/ripple_solver.f90 +++ b/NEO-2-PAR/ripple_solver.f90 @@ -33,8 +33,7 @@ SUBROUTINE ripple_solver( & lsw_save_spitf USE sparse_mod, ONLY : sparse_talk,sparse_solve_method,sparse_solve, & column_full2pointer,remap_rc,sparse_solver_test - USE mag_interface_mod, ONLY: average_bhat,average_one_over_bhat, & - surface_boozer_B00,travis_convfac, & + USE mag_interface_mod, ONLY: surface_boozer_B00,travis_convfac, & mag_magfield USE development diff --git a/NEO-2-QL/CMakeLists.txt b/NEO-2-QL/CMakeLists.txt index bebb0e56..a3759d06 100644 --- a/NEO-2-QL/CMakeLists.txt +++ b/NEO-2-QL/CMakeLists.txt @@ -34,4 +34,9 @@ add_custom_target(copy ADD_TEST(NAME Test_Build_Code COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}) ADD_TEST(Extended_QL_Multispecies ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies" "2" "QL") ADD_TEST(Extended_QL_Multispecies_Read ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies_Read" "2" "QL") +set_tests_properties( + Extended_QL_Multispecies + Extended_QL_Multispecies_Read + PROPERTIES LABELS "regression;external-data" +) #ADD_TEST(Extended_QL_Multispecies ${CMAKE_SOURCE_DIR}/../ShellScripts/test_scan.sh "Extended_QL_Scan_Kink" "QL") From b7f62d05602ce43fae70b693fc9505d2181a82fc Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 10:41:19 +0200 Subject: [PATCH 02/23] Serialize duplicate CMake build tests --- NEO-2-PAR/CMakeLists.txt | 1 + NEO-2-QL/CMakeLists.txt | 1 + 2 files changed, 2 insertions(+) diff --git a/NEO-2-PAR/CMakeLists.txt b/NEO-2-PAR/CMakeLists.txt index 402ec02c..0b548c6e 100644 --- a/NEO-2-PAR/CMakeLists.txt +++ b/NEO-2-PAR/CMakeLists.txt @@ -33,6 +33,7 @@ if (COPY_TO_ARCHIVE) endif() ADD_TEST(NAME Test_Build_Code COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}) +set_tests_properties(Test_Build_Code PROPERTIES RESOURCE_LOCK neo2_build_tree) ADD_TEST(Test_Testcase1-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Testcase1" "0") ADD_TEST(Test_Testcase1-Par ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Testcase1" "6") ADD_TEST(Test_Relativistic_BeliaevBudker-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Relativistic_BeliaevBudker" "0") diff --git a/NEO-2-QL/CMakeLists.txt b/NEO-2-QL/CMakeLists.txt index a3759d06..da30c470 100644 --- a/NEO-2-QL/CMakeLists.txt +++ b/NEO-2-QL/CMakeLists.txt @@ -32,6 +32,7 @@ add_custom_target(copy ) ADD_TEST(NAME Test_Build_Code COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}) +set_tests_properties(Test_Build_Code PROPERTIES RESOURCE_LOCK neo2_build_tree) ADD_TEST(Extended_QL_Multispecies ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies" "2" "QL") ADD_TEST(Extended_QL_Multispecies_Read ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies_Read" "2" "QL") set_tests_properties( From 1248e8bcd3f32a93d24d6bad60d0c186f7d71436 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 7 Jul 2026 20:29:46 +0200 Subject: [PATCH 03/23] Allow qflux runs with prescribed perturbation mode --- DOC/neo2.in.ql-full | 6 ++- NEO-2-QL/neo2.f90 | 21 +++++++--- NEO-2-QL/neo_magfie_perturbation.f90 | 26 +++++++++--- NEO-2-QL/ntv_mod.f90 | 18 ++++++++ TEST/CMakeLists.txt | 13 ++++++ TEST/test_ntv_perturbation_input.f90 | 63 ++++++++++++++++++++++++++++ 6 files changed, 135 insertions(+), 12 deletions(-) create mode 100644 TEST/test_ntv_perturbation_input.f90 diff --git a/DOC/neo2.in.ql-full b/DOC/neo2.in.ql-full index 36baf957..b97c481d 100644 --- a/DOC/neo2.in.ql-full +++ b/DOC/neo2.in.ql-full @@ -167,9 +167,12 @@ &ntv_input B_rho_L_loc = 0.0 ! Larmor radius times magnetic field. [0.0d0] in_file_pert = 'filename' ! File that contains perturbation of the - ! magnetic field. + ! magnetic field. Use 'none' for source-only + ! non-axisymmetric runs with prescribed m_phi. isw_mag_shear = 1 ! Turn on(=1)/off(=0) non-local computations (mag. ! shear). [0] + isw_m_phi_input = 0 ! Use m_phi_input instead of the perturbation-file mode. + ! Required when in_file_pert='none'. [0] isw_ntv_mode = 1 ! [0] isw_qflux_NA = 1 ! Turn on(=1)/off(=0) computation of non-axisymmetrix ! part. [0] @@ -178,6 +181,7 @@ ! 2 Arnoldi, 1st order ! 3 Arnoldi, 2nd order ! 4 Arnoldi, 3rd order + m_phi_input = 0 ! Toroidal mode number used when isw_m_phi_input=1. [0] MtOvR = 0.0 ! Legacy single-species Mt/R. Only used when isw_calc_Er=0; ! ignored when isw_calc_Er >= 1 (multispecies code uses Om_tE). [0.0d0] Om_tE = 0.0 ! ExB toroidal precession frequency in rad/s. diff --git a/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index e9355952..b336d01f 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -61,10 +61,11 @@ module neo2_ql USE sparse_mod, ONLY : sparse_talk,sparse_solve_method,sparse_example ! Extra input for NTV computations USE ntv_mod, ONLY : isw_ntv_mode, isw_qflux_NA, in_file_pert, & + isw_m_phi_input, m_phi_input, & MtOvR, Om_tE, B_rho_L_loc, xstart_cyl, isw_ripple_solver, & isw_calc_Er, isw_calc_MagDrift, species_tag_Vphi, & isw_Vphi_loc, Vphi, R_Vphi, Z_Vphi, boozer_theta_Vphi, & - dn_spec_ov_ds, dT_spec_ov_ds + dn_spec_ov_ds, dT_spec_ov_ds, has_perturbation_file ! derivative of iota for non-local NTV computations ! (with magnetic shear) @@ -246,7 +247,8 @@ module neo2_ql ! Extra input for NTV computation NAMELIST /ntv_input/ & isw_ntv_mode, isw_qflux_NA, in_file_pert, MtOvR, Om_tE, B_rho_L_loc, & - isw_ripple_solver, isw_mag_shear + isw_ripple_solver, isw_mag_shear, & + isw_m_phi_input, m_phi_input contains @@ -660,6 +662,8 @@ SUBROUTINE write_run_info() ! ntv_input CALL h5_define_group(h5_config_id, 'ntv_input', h5_config_group) CALL h5_add(h5_config_group, 'isw_qflux_NA', isw_qflux_NA) + CALL h5_add(h5_config_group, 'isw_m_phi_input', isw_m_phi_input) + CALL h5_add(h5_config_group, 'm_phi_input', m_phi_input) CALL h5_add(h5_config_group, 'MtOvR', MtOvR) CALL h5_add(h5_config_group, 'Om_tE', Om_tE) CALL h5_add(h5_config_group, 'B_rho_L_loc', B_rho_L_loc) @@ -891,6 +895,9 @@ subroutine set_default_values() ! ntv_input isw_ntv_mode = 0 isw_qflux_NA = 0 + in_file_pert = 'none' + isw_m_phi_input = 0 + m_phi_input = 0 MtOvR = 0.0d0 Om_tE = 0.0d0 B_rho_L_loc = 0.0d0 @@ -1158,10 +1165,12 @@ subroutine prepare_mulitspecies_scan() 'if [ ! -e ' // TRIM(ADJUSTL(in_file)) // ' ]; then ln -s ../' // & TRIM(ADJUSTL(in_file)) // ' . ; fi' call command_line_wrapper(cmd_line) - cmd_line = & - 'if [ ! -e ' // TRIM(ADJUSTL(in_file_pert)) // ' ]; then ln -s ../' // & - TRIM(ADJUSTL(in_file_pert)) // ' . ; fi' - call command_line_wrapper(cmd_line) + IF (has_perturbation_file()) THEN + cmd_line = & + 'if [ ! -e ' // TRIM(ADJUSTL(in_file_pert)) // ' ]; then ln -s ../' // & + TRIM(ADJUSTL(in_file_pert)) // ' . ; fi' + call command_line_wrapper(cmd_line) + END IF cmd_line = & 'if [ ! -e neo.in ]; then ln -s ../neo.in . ; fi' call command_line_wrapper(cmd_line) diff --git a/NEO-2-QL/neo_magfie_perturbation.f90 b/NEO-2-QL/neo_magfie_perturbation.f90 index b3500876..a9b7dc10 100644 --- a/NEO-2-QL/neo_magfie_perturbation.f90 +++ b/NEO-2-QL/neo_magfie_perturbation.f90 @@ -15,7 +15,8 @@ MODULE neo_magfie_perturbation USE partpa_mod, ONLY : bmod0 !! Modification by Andreas F. Martitsch (17.07.2014) ! Extra input for NTV computations - USE ntv_mod, ONLY : in_file_pert, m_phi + USE ntv_mod, ONLY : in_file_pert, isw_m_phi_input, m_phi_input, m_phi, & + has_perturbation_file !! End Modification by Andreas F. Martitsch (17.07.2014) ! coordinates of starting point USE mag_interface_mod, ONLY : boozer_s, boozer_theta_beg,& @@ -65,6 +66,7 @@ MODULE neo_magfie_perturbation PRIVATE r_m_pert, r_mhalf_pert, sp_index_pert REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: r_m_pert, r_mhalf_pert INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: sp_index_pert + LOGICAL :: perturbation_initialized = .FALSE. PUBLIC neo_magfie_pert PRIVATE neo_magfie_pert_a, neo_magfie_pert_b @@ -94,6 +96,18 @@ SUBROUTINE neo_read_pert() ! open Boozer file and read first quantities r_un=27052014 in_file_pert=TRIM(ADJUSTL(in_file_pert)) + + IF (.NOT. has_perturbation_file()) THEN + IF (isw_m_phi_input .EQ. 0) THEN + STOP "in_file_pert='none' requires isw_m_phi_input=1" + END IF + m_phi = m_phi_input + mnmax_pert = 0 + ns_pert = 0 + perturbation_initialized = .TRUE. + RETURN + END IF + OPEN(unit=r_un,file=in_file_pert,status='old',form='formatted') IF (inp_swi == INP_SWI_TOK_CIRC) THEN ! NEW IPP TOKAMAK @@ -206,8 +220,10 @@ SUBROUTINE neo_read_pert() PRINT *,'FATAL: There is yet no other Laboratory for the perturbed field defined!' STOP END IF + IF (isw_m_phi_input .NE. 0) m_phi = m_phi_input ! close Boozer file CLOSE (unit=r_un) + perturbation_initialized = .TRUE. RETURN END SUBROUTINE neo_read_pert @@ -306,9 +322,9 @@ SUBROUTINE neo_magfie_pert_a( x, bn_hat_pert ) ! read Boozer file and prepare spline routines (1st call) - IF (.NOT. ALLOCATED(es_pert)) THEN + IF (.NOT. perturbation_initialized) THEN CALL neo_read_pert() - CALL neo_init_spline_pert() + IF (mnmax_pert .GT. 0) CALL neo_init_spline_pert() END IF ! direct summation of Fourier components @@ -410,9 +426,9 @@ SUBROUTINE neo_magfie_pert_b( x, bn_hat_pert, dbn_hat_pert_dtheta ) ! read Boozer file and prepare spline routines (1st call) - IF (.NOT. ALLOCATED(es_pert)) THEN + IF (.NOT. perturbation_initialized) THEN CALL neo_read_pert() - CALL neo_init_spline_pert() + IF (mnmax_pert .GT. 0) CALL neo_init_spline_pert() END IF ! direct summation of Fourier components diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index f34588b0..0ed103a3 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -33,6 +33,10 @@ MODULE ntv_mod INTEGER, PUBLIC :: isw_ripple_solver !> name of perturbation file CHARACTER(len=100), PUBLIC :: in_file_pert + !> switch: use m_phi_input instead of taking m_phi from the perturbation file + INTEGER, PUBLIC :: isw_m_phi_input + !> toroidal mode number used when isw_m_phi_input is enabled + INTEGER, PUBLIC :: m_phi_input !> toroidal mach number over R_major (Mt/R). !> Only used for legacy single-species NTV output (isw_calc_Er=0). !> Ignored when isw_calc_Er >= 1; the multispecies code uses Om_tE instead. @@ -193,8 +197,22 @@ MODULE ntv_mod MODULE PROCEDURE compute_TphiNA_a END INTERFACE compute_TphiNA + PUBLIC has_perturbation_file + CONTAINS + LOGICAL FUNCTION has_perturbation_file() + CHARACTER(len=100) :: file_name + + file_name = ADJUSTL(in_file_pert) + has_perturbation_file = LEN_TRIM(file_name) .GT. 0 + IF (has_perturbation_file) THEN + IF (TRIM(file_name) .EQ. 'none') has_perturbation_file = .FALSE. + IF (TRIM(file_name) .EQ. 'None') has_perturbation_file = .FALSE. + IF (TRIM(file_name) .EQ. 'NONE') has_perturbation_file = .FALSE. + END IF + END FUNCTION has_perturbation_file + SUBROUTINE compute_Dij_norm_a(qflux_NA_in,qflux_AX_in,Dij_NA,Dij_AX) use nrtype diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index cb45b4aa..98016c75 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -50,6 +50,19 @@ set_tests_properties(er_rotation_test PROPERTIES TIMEOUT 30 ) +add_executable(test_ntv_perturbation_input test_ntv_perturbation_input.f90) +target_link_libraries(test_ntv_perturbation_input neo2_ql) + +add_test(NAME ntv_perturbation_input_test + COMMAND test_ntv_perturbation_input + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ntv_perturbation_input_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_ntv_perturbation_input.f90 b/TEST/test_ntv_perturbation_input.f90 new file mode 100644 index 00000000..9e2d6cb0 --- /dev/null +++ b/TEST/test_ntv_perturbation_input.f90 @@ -0,0 +1,63 @@ +program test_ntv_perturbation_input + use ntv_mod, only: in_file_pert, isw_m_phi_input, m_phi_input, m_phi, & + has_perturbation_file + use neo_magfie_perturbation, only: neo_read_pert + implicit none + + integer :: test_status + + test_status = 0 + + call test_none_file(test_status) + call test_file_flag(test_status) + + if (test_status == 0) then + print *, "All tests passed!" + else + print *, "Some tests failed. Status:", test_status + error stop + end if + +contains + + subroutine test_none_file(status) + integer, intent(inout) :: status + + in_file_pert = 'none' + isw_m_phi_input = 1 + m_phi_input = -1 + m_phi = 0 + + call neo_read_pert() + + if (has_perturbation_file()) then + print *, "FAIL: none should disable perturbation-file input" + status = status + 1 + else + print *, "PASS: none disables perturbation-file input" + end if + + if (m_phi /= -1) then + print *, "FAIL: prescribed m_phi was not used" + print *, " Expected:", -1 + print *, " Got:", m_phi + status = status + 1 + else + print *, "PASS: prescribed m_phi was used" + end if + end subroutine test_none_file + + subroutine test_file_flag(status) + integer, intent(inout) :: status + + in_file_pert = 'perturbation.bc' + + if (has_perturbation_file()) then + print *, "PASS: regular file name enables perturbation-file input" + else + print *, "FAIL: regular file name should enable perturbation-file input" + status = status + 1 + end if + end subroutine test_file_flag + +end program test_ntv_perturbation_input From 9cb829522427d494f23aaf17e389fd30934ecbc1 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 7 Jul 2026 22:19:08 +0200 Subject: [PATCH 04/23] QL: reject source-only NA run with no perturbation field isw_qflux_NA=1 with in_file_pert='none' assembles an identically zero non-axisymmetric source. That zero vector is the Arnoldi seed, so its normalization divides by zero (0/0 -> NaN) and zhseqr then fails to converge at every step, spinning without qflux output. Reject the configuration in check() with a clear message instead. --- NEO-2-QL/neo2.f90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index b336d01f..d8ad47db 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -934,6 +934,14 @@ subroutine check() write(*,*) 'Aborting...' stop end if + + if (isw_qflux_NA == 1 .and. .not. has_perturbation_file()) then + write(*,*) 'ERROR: non-axisymmetric solve has no source.' + write(*,*) " isw_qflux_NA=1 with in_file_pert='none' provides no" + write(*,*) ' perturbation field. Provide a perturbation file.' + write(*,*) 'Aborting...' + stop + end if end subroutine check From c1d531a4d00f87b699755f52198bc9fbc987e154 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 6 Jul 2026 21:17:07 +0200 Subject: [PATCH 05/23] QL: single-helicity misalignment drive in ripple_solver_ArnoldiO2 Add an optional right-hand side to the non-axisymmetric solve: the drive of a single-helicity perturbation given by a radial corrugation field (delta B^s/B_0^phi)_mn and an electrostatic potential harmonic Phi_mn. The vpar corrugation piece enters the A1/A2 force channels with velocity weights a3m and x^2 (sigma-odd), the ExB piece with weights x^-1 and a1m (sigma-even); the two new velocity moments are computed in compute_sources. Input via new &ntv_input variables (isw_hel_drive, m_theta_hel, hel_brad_re/im, hel_phim_re/im), default off. Band profiles, moments and dimensional prefactors follow a CAS-checked derivation (see issue #118). Off by default: golden-record QL summary.h5 is reproduced to 1e-15 against a main build. --- COMMON/collop_compute.f90 | 44 +++++ COMMON/rkstep_mod.f90 | 4 + DOC/neo2.in.ql-full | 11 ++ NEO-2-QL/neo2.f90 | 20 +- NEO-2-QL/ntv_mod.f90 | 12 ++ NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 | 186 +++++++++++++++++- 6 files changed, 273 insertions(+), 4 deletions(-) diff --git a/COMMON/collop_compute.f90 b/COMMON/collop_compute.f90 index 7a42d8b3..71b936ed 100644 --- a/COMMON/collop_compute.f90 +++ b/COMMON/collop_compute.f90 @@ -1507,6 +1507,7 @@ end function phim_phimp_rel1 end subroutine compute_Minv subroutine compute_sources(asource_s, weightlag_s, weightden_s, weightparflow_s, weightenerg_s) + use rkstep_mod, only : asource_hel real(kind=dp), dimension(:,:), intent(out) :: asource_s, weightlag_s real(kind=dp), dimension(:), intent(out) :: weightden_s, weightparflow_s, weightenerg_s real(kind=dp) :: res_int @@ -1543,6 +1544,31 @@ subroutine compute_sources(asource_s, weightlag_s, weightden_s, weightparflow_s, end if call chop(asource_s) + + ! Moments of the single-helicity misalignment drive (weights x^-1, x^2); + ! always computed by quadrature, no precomputed-matrix-element path. + if (allocated(asource_hel)) deallocate(asource_hel) + allocate(asource_hel(0:lagmax, 2)) + if (isw_relativistic .eq. 0) then + do k = 1, 2 + do m = 0, lagmax + asource_hel(m, k) = integrate(am_hel, 0d0) + end do + end do + else + do k = 1, 2 + do m = 0, lagmax + asource_hel(m, k) = norm_maxwell * integrate(am_hel_rel1, 0d0) + end do + end do + end if + if (make_ortho) then + do k = 1, 2 + asource_hel(:,k) = matmul(M_transform_inv, asource_hel(:,k)) + end do + end if + call chop(asource_hel) + write (*,*) "Computing weighting coefficients..." if (allocated(C_m)) deallocate(C_m) @@ -1615,6 +1641,24 @@ function am(x) am = pi**(-3d0/2d0) * x**(4+alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x) * x**(2*k - 1 - 5*kdelta(3,k)) end function am + ! Misalignment-drive weights x^(3k-4): k=1 -> x^-1, k=2 -> x^2 + function am_hel(x) + real(kind=dp) :: x, am_hel + + am_hel = pi**(-3d0/2d0) * x**(4+alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x) * x**(3*k - 4) + end function am_hel + + function am_hel_rel1(x) + real(kind=dp) :: x, am_hel_rel1 + real(kind=dp) :: exp_maxwell, gam + + exp_maxwell = 2d0/(sqrt(1+2*x**2/rmu) + 1) + gam = sqrt(1 + 2*x**2/rmu) + + am_hel_rel1 = 1d0/gam * pi**(-3d0/2d0) * x**(4+alpha) * exp(-(exp_maxwell+beta)*x**2) & + * phi_prj(m, x) * x**(3*k - 4) + end function am_hel_rel1 + function am_rel1(x) real(kind=dp) :: x, am_rel1 real(kind=dp) :: exp_maxwell, gam, gam_fac diff --git a/COMMON/rkstep_mod.f90 b/COMMON/rkstep_mod.f90 index 7c094ad5..c6317196 100644 --- a/COMMON/rkstep_mod.f90 +++ b/COMMON/rkstep_mod.f90 @@ -5,6 +5,10 @@ MODULE rkstep_mod CHARACTER(len=3), DIMENSION(:), ALLOCATABLE :: species_tags DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: Amm DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: asource + ! Velocity moments of the single-helicity misalignment drive + ! (column 1: weight x^-1, ExB piece / A1 channel; + ! column 2: weight x^2, vpar corrugation piece / A2 channel) + DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: asource_hel DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: anumm_a DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: denmm_a DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: anumm_aa diff --git a/DOC/neo2.in.ql-full b/DOC/neo2.in.ql-full index b97c481d..2b666a0f 100644 --- a/DOC/neo2.in.ql-full +++ b/DOC/neo2.in.ql-full @@ -188,6 +188,17 @@ ! isw_calc_Er=1: computed from ambipolarity (this value ignored). ! isw_calc_Er=2: read from here (single-surface) or multispec HDF5 ! (profile runs with isw_multispecies_init > 0). [0.0d0] + isw_hel_drive = 0 ! Turn on(=1)/off(=0) the single-helicity misalignment + ! drive (additional RHS of the non-axisymmetric equation + ! set: radial corrugation + electrostatic harmonic). [0] + m_theta_hel = 0 ! Poloidal mode number m of the misalignment-drive + ! harmonic; toroidal mode number n is m_phi from the file + ! or m_phi_input. [0] + hel_brad_re = 0.0 ! Re of the corrugation amplitude (delta B^s/B_0^phi)_mn, + ! dimensionless flux-function amplitude. [0.0d0] + hel_brad_im = 0.0 ! Im of the corrugation amplitude. [0.0d0] + hel_phim_re = 0.0 ! Re of the electrostatic harmonic Phi_mn in statvolt. [0.0d0] + hel_phim_im = 0.0 ! Im of the electrostatic harmonic Phi_mn in statvolt. [0.0d0] / &plotting plot_gauss = 0 ! Plotting of gauss function in flint. [0] diff --git a/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index d8ad47db..a4ce52d7 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -65,7 +65,9 @@ module neo2_ql MtOvR, Om_tE, B_rho_L_loc, xstart_cyl, isw_ripple_solver, & isw_calc_Er, isw_calc_MagDrift, species_tag_Vphi, & isw_Vphi_loc, Vphi, R_Vphi, Z_Vphi, boozer_theta_Vphi, & - dn_spec_ov_ds, dT_spec_ov_ds, has_perturbation_file + dn_spec_ov_ds, dT_spec_ov_ds, & + isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im, has_perturbation_file ! derivative of iota for non-local NTV computations ! (with magnetic shear) @@ -248,7 +250,9 @@ module neo2_ql NAMELIST /ntv_input/ & isw_ntv_mode, isw_qflux_NA, in_file_pert, MtOvR, Om_tE, B_rho_L_loc, & isw_ripple_solver, isw_mag_shear, & - isw_m_phi_input, m_phi_input + isw_m_phi_input, m_phi_input, & + isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im contains @@ -669,6 +673,12 @@ SUBROUTINE write_run_info() CALL h5_add(h5_config_group, 'B_rho_L_loc', B_rho_L_loc) CALL h5_add(h5_config_group, 'isw_ripple_solver', isw_ripple_solver) CALL h5_add(h5_config_group, 'isw_mag_shear', isw_mag_shear) + CALL h5_add(h5_config_group, 'isw_hel_drive', isw_hel_drive) + CALL h5_add(h5_config_group, 'm_theta_hel', m_theta_hel) + CALL h5_add(h5_config_group, 'hel_brad_re', hel_brad_re) + CALL h5_add(h5_config_group, 'hel_brad_im', hel_brad_im) + CALL h5_add(h5_config_group, 'hel_phim_re', hel_phim_re) + CALL h5_add(h5_config_group, 'hel_phim_im', hel_phim_im) CALL h5_close_group(h5_config_group) CALL h5_close(h5_config_id) @@ -903,6 +913,12 @@ subroutine set_default_values() B_rho_L_loc = 0.0d0 isw_ripple_solver = 1 isw_mag_shear = 0 + isw_hel_drive = 0 + m_theta_hel = 0 + hel_brad_re = 0.0d0 + hel_brad_im = 0.0d0 + hel_phim_re = 0.0d0 + hel_phim_im = 0.0d0 !! End Modification by Andreas F. Martitsch (14.07.2015) use_fpol = .false. ! to switch off some features, which are untested with neo-2 diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index 0ed103a3..2ad32565 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -47,6 +47,18 @@ MODULE ntv_mod REAL(kind=dp), PUBLIC :: Om_tE !> Larmor radius associated with $B_{00}^{Booz}$ (rho_L_loc) times B REAL(kind=dp), PUBLIC :: B_rho_L_loc + !> switch: turn on(=1)/off(=0) the single-helicity misalignment drive + !> (additional RHS of the non-axisymmetric equation set) + INTEGER, PUBLIC :: isw_hel_drive + !> poloidal mode number m of the misalignment-drive harmonic + !> (toroidal mode number n is m_phi from the file or m_phi_input) + INTEGER, PUBLIC :: m_theta_hel + !> complex amplitude (delta B^s/B_0^phi)_mn of the radial corrugation, + !> dimensionless flux-function amplitude + REAL(kind=dp), PUBLIC :: hel_brad_re, hel_brad_im + !> complex amplitude Phi_mn of the electrostatic perturbation harmonic + !> in statvolt + REAL(kind=dp), PUBLIC :: hel_phim_re, hel_phim_im ! ADDITIONAL INPUT FOR MULTI-SPECIES COMPUTATIONS (neo2.in) diff --git a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 index d85cb5b9..396a666e 100644 --- a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 +++ b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 @@ -17,7 +17,7 @@ subroutine ripple_solver_ArnoldiO2( & use collisionality_mod, only: collpar, conl_over_mfp, isw_lorentz, & isw_energy, isw_integral, isw_axisymm, & !<-in Winny isw_momentum, nvel, num_spec, lsw_multispecies, & - T_spec, m_spec + T_spec, m_spec, z_spec ! collpar - the same as $\kappa$ - inverse mean-free path times 4 !! Modifications by Andreas F. Martitsch (01.04.2015) ! @@ -77,7 +77,9 @@ subroutine ripple_solver_ArnoldiO2( & use ntv_mod, only: isw_qflux_NA, MtOvR, Om_tE, B_rho_L_loc, & m_phi, qflux_symm, eps_M_2_val, av_gphph_val, av_inv_bhat_val, & qflux_symm_allspec, qflux_ntv_allspec, & - MtOvR_spec, isw_calc_Er, B_rho_L_loc_spec, isw_calc_MagDrift + MtOvR_spec, isw_calc_Er, B_rho_L_loc_spec, isw_calc_MagDrift, & + isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im, clight_hel => c, echarge_hel => e use ntv_mod, only: get_Er, get_B_rho_L_loc use er_rotation_mod, only: Om_tE_to_MtOvR_spec !! End Modification by Andreas F. Martitsch (14.07.2015) @@ -218,6 +220,12 @@ subroutine ripple_solver_ArnoldiO2( & complex(dp), dimension(:, :), allocatable :: q_rip_1 complex(dp), dimension(:, :), allocatable :: q_rip_incompress complex(dp), dimension(:, :), allocatable :: q_rip_parflow + ! Single-helicity misalignment drive: band profiles along the field line + ! (q_hel_b: vpar corrugation piece, q_hel_e: ExB piece) + complex(dp), dimension(:, :), allocatable :: q_hel_b, q_hel_e + complex(dp) :: amp_hel_b, amp_hel_e, hel_phase, hel_phase_fac + real(dp) :: vt_hel, rho_hel + logical :: hel_drive_active real(dp), dimension(:, :, :), allocatable :: rhs_mat_energ real(dp), dimension(:, :, :), allocatable :: rhs_mat_energ2 !NTV real(dp), dimension(:, :, :), allocatable :: pleg_bra, pleg_ket @@ -284,6 +292,8 @@ subroutine ripple_solver_ArnoldiO2( & logical :: lsw_debug_distfun = .FALSE. logical :: addboucol = .false. + hel_drive_active = .FALSE. + ! multi-species part - MPI rank determines species ispec = mpro%getRank() print *, "Species: ", ispec @@ -923,6 +933,10 @@ subroutine ripple_solver_ArnoldiO2( & allocate (q_rip_1(npart + 2, ibeg:iend)) allocate (q_rip_incompress(npart + 2, ibeg:iend)) allocate (q_rip_parflow(npart + 2, ibeg:iend)) + allocate (q_hel_b(npart + 2, ibeg:iend)) + allocate (q_hel_e(npart + 2, ibeg:iend)) + q_hel_b = (0.d0, 0.d0) + q_hel_e = (0.d0, 0.d0) allocate (convol_flux(npart + 1, ibeg:iend), convol_curr(npart + 1, ibeg:iend)) allocate (convol_flux_0(npart + 1, ibeg:iend)) allocate (ind_start(ibeg:iend)) @@ -1008,6 +1022,7 @@ subroutine ripple_solver_ArnoldiO2( & deallocate (fun_coef, ttmp_mat) deallocate (rhs_mat_energ2) !NTV deallocate (q_rip, q_rip_1, q_rip_incompress, q_rip_parflow) + deallocate (q_hel_b, q_hel_e) deallocate (convol_flux, convol_curr, convol_flux_0) deallocate (pleg_bra, pleg_ket, scalprod_pleg) write (*, *) 'ERROR: maxval(phi_divide) > 1, returning.' @@ -3180,6 +3195,45 @@ subroutine ripple_solver_ArnoldiO2( & q_rip(:, istep, 2) = -2.d0*q_rip(:, istep, 2)*bnoverb0(istep) end do + ! Single-helicity misalignment drive: fill the band profiles along the + ! field line. Amplitudes and band integrals follow the CAS-checked + ! derivation referenced in the pull request: + ! vpar corrugation piece: band source cB (eta_k - eta_{k-1})/rho_alpha, + ! sigma-odd, velocity weights a3m (A1 channel) and x^2 (A2 channel); + ! ExB piece: band source 2 (lambda_{k-1} - lambda_k) vE^s_m / + ! (Bhat h^phi vT rho_alpha), sigma-even, velocity weights x^-1 + ! (A1 channel) and a1m (A2 channel); + ! vE^s_m z e bmod0/(2 T c) reduces to the bmod0-free hat form used + ! for amp_hel_e below. + if (isw_hel_drive .ne. 0) then + hel_drive_active = .TRUE. + hel_phase_fac = imun*(DBLE(m_theta_hel)*aiota + DBLE(m_phi)) + vt_hel = SQRT(2.d0*T_spec(ispec)/m_spec(ispec)) + rho_hel = vt_hel*m_spec(ispec)*clight_hel & + /(z_spec(ispec)*echarge_hel*(bmod0*1.0d4)) + amp_hel_b = CMPLX(hel_brad_re, hel_brad_im, kind=dp)/rho_hel + amp_hel_e = imun*CMPLX(hel_phim_re, hel_phim_im, kind=dp) & + *z_spec(ispec)*echarge_hel/(2.d0*T_spec(ispec)) & + *(DBLE(m_phi)*bcovar_theta_hat - DBLE(m_theta_hel)*bcovar_phi_hat) & + /(boozer_psi_pr_hat*denomjac) + do istep = ibeg, iend + npassing = npl(istep) + hel_phase = EXP(hel_phase_fac*phi_mfl(istep)) + q_hel_b(:, istep) = (0.d0, 0.d0) + q_hel_e(:, istep) = (0.d0, 0.d0) + q_hel_b(1:npassing, istep) = amp_hel_b*hel_phase & + *(eta(1:npassing) - eta(0:npassing - 1)) + q_hel_b(npassing + 1, istep) = amp_hel_b*hel_phase & + *(1.d0/bhat_mfl(istep) - eta(npassing)) + q_hel_e(1:npassing, istep) = amp_hel_e*hel_phase & + *2.d0*(alambd(0:npassing - 1, istep) - alambd(1:npassing, istep)) & + /(bhat_mfl(istep)*h_phi_mfl(istep)) + q_hel_e(npassing + 1, istep) = amp_hel_e*hel_phase & + *2.d0*alambd(npassing, istep) & + /(bhat_mfl(istep)*h_phi_mfl(istep)) + end do + end if + if (nobounceaver) then !! Modifications by Andreas F. Martitsch (18.08.2014) ! derivative along the periodic Boozer angle theta has @@ -3312,6 +3366,7 @@ subroutine ripple_solver_ArnoldiO2( & deallocate (enu_coef2, dellampow2, rhs_mat_energ2) deallocate (npl, rhs_mat_fzero, rhs_mat_lorentz, rhs_mat_energ, q_rip, q_rip_1) deallocate (q_rip_incompress, q_rip_parflow) + deallocate (q_hel_b, q_hel_e) deallocate (fun_coef, ttmp_mat) deallocate (convol_flux, convol_curr, ind_start, convol_flux_0) deallocate (delt_pos, delt_neg, fact_pos_b, fact_neg_b, fact_pos_e, fact_neg_e) @@ -4189,8 +4244,135 @@ subroutine source_flux end do end do + if (hel_drive_active .and. iplot .ne. 1) then + call add_hel_source(q_hel_b, asource(0:lag, 2), 1, -1) + call add_hel_source(q_hel_b, asource_hel(0:lag, 2), 3, -1) + call add_hel_source(q_hel_e, asource_hel(0:lag, 1), 1, 1) + call add_hel_source(q_hel_e, asource(0:lag, 1), 3, 1) + end if + end subroutine source_flux +!------------------------------------------------------------------------ + subroutine add_hel_source(q_hel, amom, lcol, isig) + +! Adds one misalignment-drive piece (band profile q_hel, velocity moments +! amom) to source_vector column lcol, following the integration-scheme +! structure of the standard drives assembled above. isig = +1 for +! sigma-even pieces (sign pattern of columns 1 and 3), isig = -1 for +! sigma-odd pieces (sign pattern of column 2). + + complex(dp), dimension(:, ibeg:), intent(in) :: q_hel + real(dp), dimension(0:), intent(in) :: amom + integer, intent(in) :: lcol, isig + + integer :: is, mh, kh, kh_prev, np, np_prev, np_next + real(dp) :: sg + + sg = DBLE(isig) + + do is = ibeg, iend + np = npl(is) + do mh = 0, lag + kh = ind_start(is) + 2*(np + 1)*mh + + if (is .GT. ibeg) then + np_prev = npl(is - 1) + kh_prev = ind_start(is - 1) + 2*(np_prev + 1)*mh + if (MOD(is - ibeg, 2) .EQ. 1) then + np_next = npl(is + 1) + source_vector(kh + 1:kh + np + 1, lcol) & + = source_vector(kh + 1:kh + np + 1, lcol) & + + amom(mh)/1.5d0*q_hel(1:np + 1, is)*fact_pos_e(is) + source_vector(kh + 1:kh + np_prev + 1, lcol) & + = source_vector(kh + 1:kh + np_prev + 1, lcol) & + + amom(mh)/2.4d0*q_hel(1:np_prev + 1, is - 1)*fact_pos_b(is - 1) + source_vector(kh + 1:kh + np_next + 1, lcol) & + = source_vector(kh + 1:kh + np_next + 1, lcol) & + - amom(mh)/12d0*q_hel(1:np_next + 1, is + 1)*fact_pos_e(is + 1) + else + np_next = npl(is - 2) + source_vector(kh + 1:kh + np + 1, lcol) & + = source_vector(kh + 1:kh + np + 1, lcol) & + + amom(mh)/2.4d0*q_hel(1:np + 1, is)*fact_pos_e(is) + if (np_prev .LE. np) then + source_vector(kh + 1:kh + np_prev + 1, lcol) & + = source_vector(kh + 1:kh + np_prev + 1, lcol) & + + amom(mh)/1.5d0*q_hel(1:np_prev + 1, is - 1)*fact_pos_b(is - 1) + else + source_vector(kh + 1:kh + np + 1, lcol) & + = source_vector(kh + 1:kh + np + 1, lcol) & + + amom(mh)/1.5d0*q_hel(1:np + 1, is - 1)*fact_pos_b(is - 1) + source_vector(kh_prev + np_prev + 2, lcol) & + = source_vector(kh_prev + np_prev + 2, lcol) & + + amom(mh)/1.5d0*q_hel(np_prev + 1, is - 1)*fact_pos_b(is - 1) + end if + if (np_next .LE. np) then + source_vector(kh + 1:kh + np_next + 1, lcol) & + = source_vector(kh + 1:kh + np_next + 1, lcol) & + - amom(mh)/12d0*q_hel(1:np_next + 1, is - 2)*fact_pos_b(is - 2) + else + source_vector(kh + 1:kh + np + 1, lcol) & + = source_vector(kh + 1:kh + np + 1, lcol) & + - amom(mh)/12d0*q_hel(1:np + 1, is - 2)*fact_pos_b(is - 2) + source_vector(kh_prev + np_prev + 2, lcol) & + = source_vector(kh_prev + np_prev + 2, lcol) & + - amom(mh)/12d0*q_hel(np_next + 1, is - 2)*fact_pos_b(is - 2) + end if + end if + end if + + if (is .LT. iend) then + np_prev = npl(is + 1) + kh_prev = ind_start(is + 1) + 2*(np_prev + 1)*mh + if (MOD(is - ibeg, 2) .EQ. 1) then + np_next = npl(is - 1) + source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + + sg*amom(mh)/1.5d0*q_hel(np + 1:1:-1, is)*fact_neg_e(is) + source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & + = source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & + + sg*amom(mh)/2.4d0*q_hel(np_prev + 1:1:-1, is + 1)*fact_neg_b(is + 1) + source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & + = source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & + - sg*amom(mh)/12d0*q_hel(np_next + 1:1:-1, is - 1)*fact_neg_e(is - 1) + else + np_next = npl(is + 2) + source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + + sg*amom(mh)/2.4d0*q_hel(np + 1:1:-1, is)*fact_neg_e(is) + if (np_prev .LE. np) then + source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & + = source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & + + sg*amom(mh)/1.5d0*q_hel(np_prev + 1:1:-1, is + 1)*fact_neg_b(is + 1) + else + source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + + sg*amom(mh)/1.5d0*q_hel(np + 1:1:-1, is + 1)*fact_neg_b(is + 1) + source_vector(kh_prev + np_prev + 1, lcol) & + = source_vector(kh_prev + np_prev + 1, lcol) & + + sg*amom(mh)/1.5d0*q_hel(np_prev + 1, is + 1)*fact_neg_b(is + 1) + end if + if (np_next .LE. np) then + source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & + = source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & + - sg*amom(mh)/12d0*q_hel(np_next + 1:1:-1, is + 2)*fact_neg_b(is + 2) + else + source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & + - sg*amom(mh)/12d0*q_hel(np + 1:1:-1, is + 2)*fact_neg_b(is + 2) + source_vector(kh_prev + np_prev + 1, lcol) & + = source_vector(kh_prev + np_prev + 1, lcol) & + - sg*amom(mh)/12d0*q_hel(np_next + 1, is + 2)*fact_neg_b(is + 2) + end if + end if + end if + + end do + end do + + end subroutine add_hel_source + !------------------------------------------------------------------------ subroutine add_f01_source From a3fb8741a3de8bfdf7264b5d56f0e6a2ddf74d0c Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 6 Jul 2026 21:39:17 +0200 Subject: [PATCH 06/23] Harden misalignment drive against invalid configurations Review findings on the drive commit: fold the x^(3k-4) weight into a single power so the x=0 quadrature node stays finite; skip the moment quadrature on the precomputed-matrix-element path and reject isw_hel_drive at run time when the moments are unavailable; reject the unvalidated relativistic branch, reconstruction mode and unset species data; deallocate asource_hel in collop_unload. --- COMMON/collision_operator_mems.f90 | 1 + COMMON/collop_compute.f90 | 50 ++++++++++--------- NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 | 16 ++++++ 3 files changed, 44 insertions(+), 23 deletions(-) diff --git a/COMMON/collision_operator_mems.f90 b/COMMON/collision_operator_mems.f90 index d1f061dd..fb1bf362 100644 --- a/COMMON/collision_operator_mems.f90 +++ b/COMMON/collision_operator_mems.f90 @@ -502,6 +502,7 @@ subroutine collop_unload() if (allocated(denmm_a)) deallocate(denmm_a) if (allocated(Inbi_lmmp_a)) deallocate(Inbi_lmmp_a) if (allocated(asource)) deallocate(asource) + if (allocated(asource_hel)) deallocate(asource_hel) if (allocated(weightlag)) deallocate(weightlag) if (allocated(x1mm)) deallocate(x1mm) if (allocated(x2mm)) deallocate(x2mm) diff --git a/COMMON/collop_compute.f90 b/COMMON/collop_compute.f90 index 71b936ed..eab77366 100644 --- a/COMMON/collop_compute.f90 +++ b/COMMON/collop_compute.f90 @@ -1546,28 +1546,31 @@ subroutine compute_sources(asource_s, weightlag_s, weightden_s, weightparflow_s, call chop(asource_s) ! Moments of the single-helicity misalignment drive (weights x^-1, x^2); - ! always computed by quadrature, no precomputed-matrix-element path. - if (allocated(asource_hel)) deallocate(asource_hel) - allocate(asource_hel(0:lagmax, 2)) - if (isw_relativistic .eq. 0) then - do k = 1, 2 - do m = 0, lagmax - asource_hel(m, k) = integrate(am_hel, 0d0) + ! computed by quadrature only (left unallocated on the precomputed + ! matrix-element path; the solver rejects the drive in that case). + if (.not. precomp) then + if (allocated(asource_hel)) deallocate(asource_hel) + allocate(asource_hel(0:lagmax, 2)) + if (isw_relativistic .eq. 0) then + do k = 1, 2 + do m = 0, lagmax + asource_hel(m, k) = integrate(am_hel, 0d0) + end do end do - end do - else - do k = 1, 2 - do m = 0, lagmax - asource_hel(m, k) = norm_maxwell * integrate(am_hel_rel1, 0d0) + else + do k = 1, 2 + do m = 0, lagmax + asource_hel(m, k) = norm_maxwell * integrate(am_hel_rel1, 0d0) + end do end do - end do - end if - if (make_ortho) then - do k = 1, 2 - asource_hel(:,k) = matmul(M_transform_inv, asource_hel(:,k)) - end do + end if + if (make_ortho) then + do k = 1, 2 + asource_hel(:,k) = matmul(M_transform_inv, asource_hel(:,k)) + end do + end if + call chop(asource_hel) end if - call chop(asource_hel) write (*,*) "Computing weighting coefficients..." @@ -1641,11 +1644,12 @@ function am(x) am = pi**(-3d0/2d0) * x**(4+alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x) * x**(2*k - 1 - 5*kdelta(3,k)) end function am - ! Misalignment-drive weights x^(3k-4): k=1 -> x^-1, k=2 -> x^2 + ! Misalignment-drive weights x^(3k-4): k=1 -> x^-1, k=2 -> x^2. + ! Powers folded into one factor so the x=0 quadrature node stays finite. function am_hel(x) real(kind=dp) :: x, am_hel - am_hel = pi**(-3d0/2d0) * x**(4+alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x) * x**(3*k - 4) + am_hel = pi**(-3d0/2d0) * x**(3*k + alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x) end function am_hel function am_hel_rel1(x) @@ -1655,8 +1659,8 @@ function am_hel_rel1(x) exp_maxwell = 2d0/(sqrt(1+2*x**2/rmu) + 1) gam = sqrt(1 + 2*x**2/rmu) - am_hel_rel1 = 1d0/gam * pi**(-3d0/2d0) * x**(4+alpha) * exp(-(exp_maxwell+beta)*x**2) & - * phi_prj(m, x) * x**(3*k - 4) + am_hel_rel1 = 1d0/gam * pi**(-3d0/2d0) * x**(3*k + alpha) * exp(-(exp_maxwell+beta)*x**2) & + * phi_prj(m, x) end function am_hel_rel1 function am_rel1(x) diff --git a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 index 396a666e..7ba6358e 100644 --- a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 +++ b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 @@ -3206,6 +3206,22 @@ subroutine ripple_solver_ArnoldiO2( & ! vE^s_m z e bmod0/(2 T c) reduces to the bmod0-free hat form used ! for amp_hel_e below. if (isw_hel_drive .ne. 0) then + if (.not. allocated(asource_hel)) then + stop 'isw_hel_drive requires quadrature collision moments '// & + '(incompatible with precomputed matrix elements)' + end if + if (isw_relativistic .ne. 0) then + stop 'isw_hel_drive: relativistic drive moments not validated' + end if + if (iplot .eq. 1) then + stop 'isw_hel_drive: reconstruction mode not supported' + end if + if (T_spec(ispec) .le. 0.d0 .or. m_spec(ispec) .le. 0.d0) then + stop 'isw_hel_drive: invalid species temperature or mass' + end if + if (z_spec(ispec) .eq. 0.d0) then + stop 'isw_hel_drive: invalid species charge number' + end if hel_drive_active = .TRUE. hel_phase_fac = imun*(DBLE(m_theta_hel)*aiota + DBLE(m_phi)) vt_hel = SQRT(2.d0*T_spec(ispec)/m_spec(ispec)) From 7a868ae263444d772a17f69871900717f1777455 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 10:43:37 +0200 Subject: [PATCH 07/23] Test helical-drive source activation --- NEO-2-QL/ntv_mod.f90 | 12 ++++++- TEST/CMakeLists.txt | 13 ++++++++ TEST/test_ntv_hel_drive_input.f90 | 52 +++++++++++++++++++++++++++++++ 3 files changed, 76 insertions(+), 1 deletion(-) create mode 100644 TEST/test_ntv_hel_drive_input.f90 diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index 2ad32565..6ff3e4b5 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -209,10 +209,20 @@ MODULE ntv_mod MODULE PROCEDURE compute_TphiNA_a END INTERFACE compute_TphiNA - PUBLIC has_perturbation_file + PUBLIC has_perturbation_file, has_helical_drive_source CONTAINS + LOGICAL FUNCTION has_helical_drive_source() + has_helical_drive_source = .FALSE. + IF (isw_hel_drive .EQ. 0) RETURN + + IF (hel_brad_re .NE. 0.0_dp) has_helical_drive_source = .TRUE. + IF (hel_brad_im .NE. 0.0_dp) has_helical_drive_source = .TRUE. + IF (hel_phim_re .NE. 0.0_dp) has_helical_drive_source = .TRUE. + IF (hel_phim_im .NE. 0.0_dp) has_helical_drive_source = .TRUE. + END FUNCTION has_helical_drive_source + LOGICAL FUNCTION has_perturbation_file() CHARACTER(len=100) :: file_name diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 98016c75..b36a2a9f 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -63,6 +63,19 @@ set_tests_properties(ntv_perturbation_input_test PROPERTIES TIMEOUT 30 ) +add_executable(test_ntv_hel_drive_input test_ntv_hel_drive_input.f90) +target_link_libraries(test_ntv_hel_drive_input neo2_ql) + +add_test(NAME ntv_hel_drive_input_test + COMMAND test_ntv_hel_drive_input + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ntv_hel_drive_input_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_ntv_hel_drive_input.f90 b/TEST/test_ntv_hel_drive_input.f90 new file mode 100644 index 00000000..d605d2c9 --- /dev/null +++ b/TEST/test_ntv_hel_drive_input.f90 @@ -0,0 +1,52 @@ +program test_ntv_hel_drive_input + use nrtype, only: dp + use ntv_mod, only: isw_hel_drive, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im, has_helical_drive_source + implicit none + + integer :: failures + + failures = 0 + + call set_drive(0, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp) + call expect_source(.false., failures) + + call set_drive(1, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp) + call expect_source(.false., failures) + + call set_drive(1, 1.0e-3_dp, 0.0_dp, 0.0_dp, 0.0_dp) + call expect_source(.true., failures) + + call set_drive(1, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp) + call expect_source(.true., failures) + + if (failures == 0) then + print *, 'All tests passed!' + else + error stop 'helical-drive source detection failed' + end if + +contains + + subroutine set_drive(switch, brad_re, brad_im, phim_re, phim_im) + integer, intent(in) :: switch + real(dp), intent(in) :: brad_re, brad_im, phim_re, phim_im + + isw_hel_drive = switch + hel_brad_re = brad_re + hel_brad_im = brad_im + hel_phim_re = phim_re + hel_phim_im = phim_im + end subroutine set_drive + + subroutine expect_source(expected, status) + logical, intent(in) :: expected + integer, intent(inout) :: status + + if (has_helical_drive_source() .neqv. expected) then + print *, 'FAIL: unexpected helical-drive source state' + status = status + 1 + end if + end subroutine expect_source + +end program test_ntv_hel_drive_input From b00a3ca685b2e35e41a5394c7dd0711374e4f525 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 7 Jul 2026 22:20:59 +0200 Subject: [PATCH 08/23] QL: allow source-only NA run when the helical drive is the source The zero-source guard added upstream rejects isw_qflux_NA=1 with in_file_pert='none'. With the misalignment drive, that configuration is valid: the drive supplies the non-axisymmetric source. Relax the guard so it fires only when the drive is off or all its amplitudes are zero. --- NEO-2-QL/neo2.f90 | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index a4ce52d7..0838bb62 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -952,11 +952,17 @@ subroutine check() end if if (isw_qflux_NA == 1 .and. .not. has_perturbation_file()) then - write(*,*) 'ERROR: non-axisymmetric solve has no source.' - write(*,*) " isw_qflux_NA=1 with in_file_pert='none' provides no" - write(*,*) ' perturbation field. Provide a perturbation file.' - write(*,*) 'Aborting...' - stop + if (isw_hel_drive == 0 .or. & + (hel_phim_re == 0.0d0 .and. hel_phim_im == 0.0d0 .and. & + hel_brad_re == 0.0d0 .and. hel_brad_im == 0.0d0)) then + write(*,*) 'ERROR: non-axisymmetric solve has no source.' + write(*,*) " isw_qflux_NA=1 with in_file_pert='none' provides no" + write(*,*) ' perturbation field, and the helical drive is off or zero.' + write(*,*) ' Provide a perturbation file, or set isw_hel_drive with a' + write(*,*) ' non-zero hel_phim/hel_brad amplitude.' + write(*,*) 'Aborting...' + stop + end if end if end subroutine check From 8ea3eec23c2f59d27831c8abb2bf3cf2a519f7d7 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 10:49:51 +0200 Subject: [PATCH 09/23] Test source-only non-axisymmetric guard --- NEO-2-QL/neo2.f90 | 7 +++---- TEST/CMakeLists.txt | 23 +++++++++++++++++++++++ TEST/test_ntv_source_guard.f90 | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 58 insertions(+), 4 deletions(-) create mode 100644 TEST/test_ntv_source_guard.f90 diff --git a/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index 0838bb62..a3506b51 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -67,7 +67,8 @@ module neo2_ql isw_Vphi_loc, Vphi, R_Vphi, Z_Vphi, boozer_theta_Vphi, & dn_spec_ov_ds, dT_spec_ov_ds, & isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & - hel_phim_re, hel_phim_im, has_perturbation_file + hel_phim_re, hel_phim_im, has_perturbation_file, & + has_helical_drive_source ! derivative of iota for non-local NTV computations ! (with magnetic shear) @@ -952,9 +953,7 @@ subroutine check() end if if (isw_qflux_NA == 1 .and. .not. has_perturbation_file()) then - if (isw_hel_drive == 0 .or. & - (hel_phim_re == 0.0d0 .and. hel_phim_im == 0.0d0 .and. & - hel_brad_re == 0.0d0 .and. hel_brad_im == 0.0d0)) then + if (.not. has_helical_drive_source()) then write(*,*) 'ERROR: non-axisymmetric solve has no source.' write(*,*) " isw_qflux_NA=1 with in_file_pert='none' provides no" write(*,*) ' perturbation field, and the helical drive is off or zero.' diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index b36a2a9f..4b6fc233 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -76,6 +76,29 @@ set_tests_properties(ntv_hel_drive_input_test PROPERTIES TIMEOUT 30 ) +add_executable(test_ntv_source_guard test_ntv_source_guard.f90) +target_link_libraries(test_ntv_source_guard neo2_ql) + +add_test(NAME ntv_zero_source_guard_test + COMMAND test_ntv_source_guard zero + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ntv_zero_source_guard_test PROPERTIES + PASS_REGULAR_EXPRESSION "ERROR: non-axisymmetric solve has no source." + FAIL_REGULAR_EXPRESSION "FAIL;NaN;zhseqr;returned from rejected source" + TIMEOUT 5 +) + +add_test(NAME ntv_driven_source_guard_test + COMMAND test_ntv_source_guard driven + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ntv_driven_source_guard_test PROPERTIES + PASS_REGULAR_EXPRESSION "Driven source accepted!" + FAIL_REGULAR_EXPRESSION "FAIL;NaN;zhseqr;non-axisymmetric solve has no source" + TIMEOUT 5 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_ntv_source_guard.f90 b/TEST/test_ntv_source_guard.f90 new file mode 100644 index 00000000..0be73384 --- /dev/null +++ b/TEST/test_ntv_source_guard.f90 @@ -0,0 +1,32 @@ +program test_ntv_source_guard + use neo2_ql, only: init, check + use ntv_mod, only: isw_qflux_NA, in_file_pert, isw_hel_drive, & + hel_brad_re, hel_brad_im, hel_phim_re, hel_phim_im + implicit none + + character(len=16) :: test_case + + call get_command_argument(1, test_case) + call init() + + isw_qflux_NA = 1 + in_file_pert = 'none' + isw_hel_drive = 0 + hel_brad_re = 0.0d0 + hel_brad_im = 0.0d0 + hel_phim_re = 0.0d0 + hel_phim_im = 0.0d0 + + select case (trim(test_case)) + case ('zero') + call check() + error stop 'returned from rejected source' + case ('driven') + isw_hel_drive = 1 + hel_phim_re = 1.0d0 + call check() + print *, 'Driven source accepted!' + case default + error stop 'FAIL: unknown test case' + end select +end program test_ntv_source_guard From 4b00b887b6a9c37928b6575130b9ad494d2cdc2b Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 11:00:20 +0200 Subject: [PATCH 10/23] Preserve full-step endpoints in mode-3 refinement --- COMMON/flint_mod.f90 | 6 ++- TEST/CMakeLists.txt | 13 ++++++ TEST/test_phi_divider_parity.f90 | 73 ++++++++++++++++++++++++++++++++ 3 files changed, 90 insertions(+), 2 deletions(-) create mode 100644 TEST/test_phi_divider_parity.f90 diff --git a/COMMON/flint_mod.f90 b/COMMON/flint_mod.f90 index 25b86dc0..6b7640de 100644 --- a/COMMON/flint_mod.f90 +++ b/COMMON/flint_mod.f90 @@ -899,7 +899,7 @@ SUBROUTINE phi_divider(u_eta,phi_eta_ind) !<-in Winny INTEGER, INTENT(inout) :: phi_eta_ind(0:u_eta,2) INTEGER :: ub - INTEGER :: iphi,idiv,cphi,cdiv,ieta + INTEGER :: iphi,idiv,cphi,cdiv,ieta,refinement_scale INTEGER, ALLOCATABLE :: phi_eta_ind_loc(:,:) REAL(kind=dp) :: hphi REAL(kind=dp), ALLOCATABLE :: phi(:) @@ -910,6 +910,8 @@ SUBROUTINE phi_divider(u_eta,phi_eta_ind) !<-in Winny phi_eta_ind_loc = phi_eta_ind ub = UBOUND(fieldpropagator%coords%x2,1) + refinement_scale = 1 + IF (MOD(SUM(phi_divide(1:ub)),2) .NE. 0) refinement_scale = 2 ALLOCATE(phi(0:ub)) phi = fieldpropagator%coords%x2 IF (ALLOCATED(phiarr)) DEALLOCATE(phiarr) @@ -918,7 +920,7 @@ SUBROUTINE phi_divider(u_eta,phi_eta_ind) !<-in Winny CALL set_new(phi_p,phi(0)) cphi = 0 ! counter philoop: DO iphi = 1,ub - cdiv = phi_divide(iphi) + cdiv = refinement_scale*phi_divide(iphi) hphi = phi(iphi) - phi(iphi-1) divloop: DO idiv = 1,cdiv cphi = cphi + 1 diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 4b6fc233..a79e5720 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -99,6 +99,19 @@ set_tests_properties(ntv_driven_source_guard_test PROPERTIES TIMEOUT 5 ) +add_executable(test_phi_divider_parity test_phi_divider_parity.f90) +target_link_libraries(test_phi_divider_parity common) + +add_test(NAME phi_divider_parity_test + COMMAND test_phi_divider_parity + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(phi_divider_parity_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_phi_divider_parity.f90 b/TEST/test_phi_divider_parity.f90 new file mode 100644 index 00000000..7facb71c --- /dev/null +++ b/TEST/test_phi_divider_parity.f90 @@ -0,0 +1,73 @@ +program test_phi_divider_parity + use device_mod, only: fieldpropagator + use flint_mod, only: phi_divide, phi_divider, phiarr + implicit none + + integer, parameter :: ninterval = 4 + integer :: forward_map(0:ninterval), reverse_map(0:ninterval) + real(kind(1.0d0)) :: forward_nodes(0:ninterval), reverse_nodes(0:ninterval) + integer :: i, failures + + failures = 0 + call refined_node_map([2, 1, 2, 2], forward_map, forward_nodes) + call refined_node_map([2, 2, 1, 2], reverse_map, reverse_nodes) + + do i = 0, ninterval + if (forward_nodes(i) /= real(i, kind(1.0d0)) .or. & + reverse_nodes(i) /= real(i, kind(1.0d0))) then + print *, 'FAIL: refinement dropped an original grid point' + failures = failures + 1 + end if + if (reverse_map(i) /= forward_map(ninterval) - forward_map(ninterval - i)) then + print *, 'FAIL: reversed refinement changed the mapped grid' + failures = failures + 1 + end if + end do + + if (any(forward_map(1:) <= forward_map(:ninterval - 1)) .or. & + any(reverse_map(1:) <= reverse_map(:ninterval - 1))) then + print *, 'FAIL: refinement did not preserve grid ordering' + failures = failures + 1 + end if + + if (mod(forward_map(ninterval), 2) /= 0) then + print *, 'FAIL: refinement made the final Runge-Kutta point a half step' + failures = failures + 1 + end if + + if (failures == 0) then + print *, 'All tests passed!' + else + error stop 'phi-divider parity contract failed' + end if + +contains + + subroutine refined_node_map(divisions, node_map, node_values) + integer, intent(in) :: divisions(ninterval) + integer, intent(out) :: node_map(0:ninterval) + real(kind(1.0d0)), intent(out) :: node_values(0:ninterval) + + integer :: i, phi_eta_ind(0:ninterval, 2) + + allocate(fieldpropagator) + allocate(fieldpropagator%coords) + allocate(fieldpropagator%coords%x2(0:ninterval)) + fieldpropagator%coords%x2 = [(real(i, kind(1.0d0)), i = 0, ninterval)] + + allocate(phi_divide(ninterval)) + phi_divide = divisions + do i = 0, ninterval + phi_eta_ind(i, :) = i + end do + + call phi_divider(ninterval, phi_eta_ind) + node_map = phi_eta_ind(:, 1) + node_values = phiarr(node_map) + + deallocate(phiarr, phi_divide) + deallocate(fieldpropagator%coords%x2) + deallocate(fieldpropagator%coords) + deallocate(fieldpropagator) + end subroutine refined_node_map +end program test_phi_divider_parity From 691a5f55a070f92cc5fa23f394c609372251a0d8 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 11:25:58 +0200 Subject: [PATCH 11/23] Report propagator join solver failures --- NEO-2-QL/CMakeSources.in | 1 + NEO-2-QL/join_diagnostics_mod.f90 | 27 +++++++++++++++++++++++++++ NEO-2-QL/join_ripples.f90 | 16 +++++++++++++--- TEST/CMakeLists.txt | 13 +++++++++++++ TEST/test_join_failure_diagnostic.f90 | 17 +++++++++++++++++ 5 files changed, 71 insertions(+), 3 deletions(-) create mode 100644 NEO-2-QL/join_diagnostics_mod.f90 create mode 100644 TEST/test_join_failure_diagnostic.f90 diff --git a/NEO-2-QL/CMakeSources.in b/NEO-2-QL/CMakeSources.in index 9717a122..23382ac4 100644 --- a/NEO-2-QL/CMakeSources.in +++ b/NEO-2-QL/CMakeSources.in @@ -6,6 +6,7 @@ set(NEO2_QL_SRC_FILES flint.f90 ripple_solver_axi_test.f90 ripple_solver_ArnoldiOrder2_test.f90 + join_diagnostics_mod.f90 join_ripples_int.f90 join_ripples.f90 join_ends.f90 diff --git a/NEO-2-QL/join_diagnostics_mod.f90 b/NEO-2-QL/join_diagnostics_mod.f90 new file mode 100644 index 00000000..87edcfb2 --- /dev/null +++ b/NEO-2-QL/join_diagnostics_mod.f90 @@ -0,0 +1,27 @@ +module join_diagnostics_mod + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + implicit none + private + + public :: report_join_failure + +contains + + subroutine report_join_failure(stage, info, old_tags, new_tags, ndim, ndim1, & + matrix, rhs, ierr) + integer, intent(in) :: stage, info, old_tags(2), new_tags(2), ndim, ndim1 + real(kind(1.0d0)), intent(in) :: matrix(:, :), rhs(:, :) + integer, intent(out) :: ierr + + write(*, '(a,i0)') 'join_ripples: DGBSV stage=', stage + write(*, '(a,i0)') 'join_ripples: DGBSV info=', info + write(*, '(a,2(1x,i0))') 'join_ripples: old propagator tags=', old_tags + write(*, '(a,2(1x,i0))') 'join_ripples: new propagator tags=', new_tags + write(*, '(a,2(1x,i0))') 'join_ripples: dimensions=', ndim, ndim1 + write(*, '(a,i0,a,i0)') 'join_ripples: finite matrix entries=', & + count(ieee_is_finite(matrix)), '/', size(matrix) + write(*, '(a,i0,a,i0)') 'join_ripples: finite rhs entries=', & + count(ieee_is_finite(rhs)), '/', size(rhs) + ierr = stage + end subroutine report_join_failure +end module join_diagnostics_mod diff --git a/NEO-2-QL/join_ripples.f90 b/NEO-2-QL/join_ripples.f90 index f8028148..2897b462 100644 --- a/NEO-2-QL/join_ripples.f90 +++ b/NEO-2-QL/join_ripples.f90 @@ -19,6 +19,7 @@ SUBROUTINE join_ripples(ierr) USE propagator_mod USE lapack_band USE collisionality_mod, ONLY : isw_lorentz + USE join_diagnostics_mod, ONLY : report_join_failure USE development IMPLICIT NONE @@ -138,7 +139,10 @@ SUBROUTINE join_ripples(ierr) CALL gbsv(ndim,ndim,amat,ipivot,bvec_lapack,info) IF(info.NE.0) THEN - ierr=2 + CALL report_join_failure(2,info, & + [o%fieldpropagator_tag_s,o%fieldpropagator_tag_e], & + [n%fieldpropagator_tag_s,n%fieldpropagator_tag_e], & + ndim,ndim1,amat,bvec_lapack,ierr) RETURN ENDIF @@ -179,7 +183,10 @@ SUBROUTINE join_ripples(ierr) ! solution bvec_lapack -> $\ifour{\Fmat}{+}{}{o+1}{l,l}$ CALL gbsv(ndim,ndim,amat,ipivot,bvec_lapack,info) IF(info.NE.0) THEN - ierr=3 + CALL report_join_failure(3,info, & + [o%fieldpropagator_tag_s,o%fieldpropagator_tag_e], & + [n%fieldpropagator_tag_s,n%fieldpropagator_tag_e], & + ndim,ndim1,amat,bvec_lapack,ierr) RETURN ENDIF @@ -232,7 +239,10 @@ SUBROUTINE join_ripples(ierr) ! solution bvec_lapack -> $\ifour{\Fmat}{+}{}{o+1}{r,l}$ CALL gbsv(ndim,ndim,amat,ipivot,bvec_lapack,info) IF(info.NE.0) THEN - ierr=4 + CALL report_join_failure(4,info, & + [o%fieldpropagator_tag_s,o%fieldpropagator_tag_e], & + [n%fieldpropagator_tag_s,n%fieldpropagator_tag_e], & + ndim,ndim1,amat,bvec_lapack,ierr) RETURN ENDIF diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index a79e5720..7e694737 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -112,6 +112,19 @@ set_tests_properties(phi_divider_parity_test PROPERTIES TIMEOUT 30 ) +add_executable(test_join_failure_diagnostic test_join_failure_diagnostic.f90) +target_link_libraries(test_join_failure_diagnostic neo2_ql) + +add_test(NAME join_failure_diagnostic_test + COMMAND test_join_failure_diagnostic + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(join_failure_diagnostic_test PROPERTIES + PASS_REGULAR_EXPRESSION "join_ripples: DGBSV info=2" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_join_failure_diagnostic.f90 b/TEST/test_join_failure_diagnostic.f90 new file mode 100644 index 00000000..c373e620 --- /dev/null +++ b/TEST/test_join_failure_diagnostic.f90 @@ -0,0 +1,17 @@ +program test_join_failure_diagnostic + use join_diagnostics_mod, only: report_join_failure + use lapack_band, only: gbsv + implicit none + + real(kind(1.0d0)) :: matrix(2, 2), rhs(2, 1) + integer :: info, ierr, pivot(2) + + matrix = reshape([1.0d0, 2.0d0, 2.0d0, 4.0d0], shape(matrix)) + rhs(:, 1) = [1.0d0, 2.0d0] + call gbsv(2, 2, matrix, pivot, rhs, info) + + if (info /= 2) error stop 'FAIL: singular solve did not return pivot 2' + call report_join_failure(2, info, [10, 20], [21, 30], 2, 2, & + matrix, rhs, ierr) + if (ierr /= 2) error stop 'FAIL: join stage was not retained' +end program test_join_failure_diagnostic From d614da521e96690aeb51a7067be91fbe94531c57 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 11:33:23 +0200 Subject: [PATCH 12/23] Use final geometry for NTV output --- NEO-2-QL/ntv_mod.f90 | 23 +++++++++++++++++++---- NEO-2-QL/propagator.f90 | 2 +- TEST/CMakeLists.txt | 13 +++++++++++++ TEST/test_ntv_output_geometry.f90 | 19 +++++++++++++++++++ 4 files changed, 52 insertions(+), 5 deletions(-) create mode 100644 TEST/test_ntv_output_geometry.f90 diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index 6ff3e4b5..b9910869 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -128,6 +128,7 @@ MODULE ntv_mod INTEGER, PUBLIC :: m_phi ! Local copy of y-vector (see definition in rhs_kin.f90) REAL(kind=dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: y_ntv_mod + PUBLIC :: resolve_ntv_output_geometry PUBLIC write_ntv_output @@ -213,6 +214,21 @@ MODULE ntv_mod CONTAINS + SUBROUTINE resolve_ntv_output_geometry(y, final_y) + REAL(kind=dp), DIMENSION(:), ALLOCATABLE, INTENT(out) :: y + REAL(kind=dp), DIMENSION(:), OPTIONAL, INTENT(in) :: final_y + + IF (PRESENT(final_y)) THEN + ALLOCATE(y(SIZE(final_y))) + y = final_y + ELSE + IF (.NOT. ALLOCATED(y_ntv_mod)) & + ERROR STOP 'NTV output geometry is not available' + ALLOCATE(y(SIZE(y_ntv_mod))) + y = y_ntv_mod + END IF + END SUBROUTINE resolve_ntv_output_geometry + LOGICAL FUNCTION has_helical_drive_source() has_helical_drive_source = .FALSE. IF (isw_hel_drive .EQ. 0) RETURN @@ -756,7 +772,7 @@ SUBROUTINE compute_Dijab_a(qflux_ab_NA_in, qflux_ab_AX_in, & END SUBROUTINE compute_Dijab_a - SUBROUTINE write_ntv_output_a(qflux_in) + SUBROUTINE write_ntv_output_a(qflux_in, final_y) use nrtype USE neo_control, ONLY: lab_swi @@ -772,6 +788,7 @@ SUBROUTINE write_ntv_output_a(qflux_in) ! input: ! ---------------------------------------------------------------! REAL(kind=dp), DIMENSION(:,:), INTENT(in) :: qflux_in + REAL(kind=dp), DIMENSION(:), OPTIONAL, INTENT(in) :: final_y ! ---------------------------------------------------------------! ! local definitions: ! ---------------------------------------------------------------! @@ -809,9 +826,7 @@ SUBROUTINE write_ntv_output_a(qflux_in) REAL(kind=dp) :: bcovar_phi, bcovar_tht, avbhat2, avb2, avbhat CHARACTER(len=30) :: file_name - ! copy y-vector (see definition in rhs_kin.f90) - ALLOCATE(y(SIZE(y_ntv_mod,1))) - y = y_ntv_mod + CALL resolve_ntv_output_geometry(y, final_y) ! R0, aiota, av_nabla_stor : rt0 = device%r0 diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index 7de989a6..bbeabfe4 100644 --- a/NEO-2-QL/propagator.f90 +++ b/NEO-2-QL/propagator.f90 @@ -953,7 +953,7 @@ SUBROUTINE diag_propagator_res(iend) IF(lsw_multispecies) THEN ! multi-species output CALL write_multispec_output() ELSE ! single-species output - CALL write_ntv_output(prop_a%p%qflux) + CALL write_ntv_output(prop_a%p%qflux, y) END IF END IF !! End Modification by Andreas F. Martitsch (14.07.2015) diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 7e694737..9e0aee10 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -125,6 +125,19 @@ set_tests_properties(join_failure_diagnostic_test PROPERTIES TIMEOUT 30 ) +add_executable(test_ntv_output_geometry test_ntv_output_geometry.f90) +target_link_libraries(test_ntv_output_geometry neo2_ql) + +add_test(NAME ntv_output_geometry_test + COMMAND test_ntv_output_geometry + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ntv_output_geometry_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_ntv_output_geometry.f90 b/TEST/test_ntv_output_geometry.f90 new file mode 100644 index 00000000..5517ddb6 --- /dev/null +++ b/TEST/test_ntv_output_geometry.f90 @@ -0,0 +1,19 @@ +program test_ntv_output_geometry + use ntv_mod, only: y_ntv_mod, resolve_ntv_output_geometry + implicit none + + real(kind(1.0d0)), allocatable :: output_y(:) + + if (allocated(y_ntv_mod)) deallocate(y_ntv_mod) + call resolve_ntv_output_geometry(output_y, [1.0d0, 2.0d0, 3.0d0]) + if (any(output_y /= [1.0d0, 2.0d0, 3.0d0])) & + error stop 'FAIL: final geometry was not selected' + + allocate(y_ntv_mod(2)) + y_ntv_mod = [4.0d0, 5.0d0] + call resolve_ntv_output_geometry(output_y) + if (any(output_y /= y_ntv_mod)) & + error stop 'FAIL: legacy exchange geometry was not selected' + + print *, 'All tests passed!' +end program test_ntv_output_geometry From 377226441848b64f39e97462937797e10342bab5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 11:43:10 +0200 Subject: [PATCH 13/23] Synchronize final NTV output geometry --- NEO-2-QL/ntv_mod.f90 | 14 +++++++------- TEST/test_ntv_output_geometry.f90 | 5 +++++ 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index b9910869..d36d69e3 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -219,14 +219,14 @@ SUBROUTINE resolve_ntv_output_geometry(y, final_y) REAL(kind=dp), DIMENSION(:), OPTIONAL, INTENT(in) :: final_y IF (PRESENT(final_y)) THEN - ALLOCATE(y(SIZE(final_y))) - y = final_y - ELSE - IF (.NOT. ALLOCATED(y_ntv_mod)) & - ERROR STOP 'NTV output geometry is not available' - ALLOCATE(y(SIZE(y_ntv_mod))) - y = y_ntv_mod + IF (ALLOCATED(y_ntv_mod)) DEALLOCATE(y_ntv_mod) + ALLOCATE(y_ntv_mod(SIZE(final_y))) + y_ntv_mod = final_y END IF + IF (.NOT. ALLOCATED(y_ntv_mod)) & + ERROR STOP 'NTV output geometry is not available' + ALLOCATE(y(SIZE(y_ntv_mod))) + y = y_ntv_mod END SUBROUTINE resolve_ntv_output_geometry LOGICAL FUNCTION has_helical_drive_source() diff --git a/TEST/test_ntv_output_geometry.f90 b/TEST/test_ntv_output_geometry.f90 index 5517ddb6..c17c2848 100644 --- a/TEST/test_ntv_output_geometry.f90 +++ b/TEST/test_ntv_output_geometry.f90 @@ -8,7 +8,12 @@ program test_ntv_output_geometry call resolve_ntv_output_geometry(output_y, [1.0d0, 2.0d0, 3.0d0]) if (any(output_y /= [1.0d0, 2.0d0, 3.0d0])) & error stop 'FAIL: final geometry was not selected' + if (.not. allocated(y_ntv_mod)) & + error stop 'FAIL: final geometry did not initialize output state' + if (any(y_ntv_mod /= [1.0d0, 2.0d0, 3.0d0])) & + error stop 'FAIL: final geometry did not synchronize output state' + deallocate(y_ntv_mod) allocate(y_ntv_mod(2)) y_ntv_mod = [4.0d0, 5.0d0] call resolve_ntv_output_geometry(output_y) From 69fbb86dbce74526093b95936b7c052e4a7fccaf Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 12:11:29 +0200 Subject: [PATCH 14/23] Handle collision-kernel zero endpoints --- COMMON/collop_compute.f90 | 16 ++++++++++++---- TEST/CMakeLists.txt | 13 +++++++++++++ TEST/test_collop_zero_endpoint.f90 | 19 +++++++++++++++++++ 3 files changed, 44 insertions(+), 4 deletions(-) create mode 100644 TEST/test_collop_zero_endpoint.f90 diff --git a/COMMON/collop_compute.f90 b/COMMON/collop_compute.f90 index eab77366..b56aa4de 100644 --- a/COMMON/collop_compute.f90 +++ b/COMMON/collop_compute.f90 @@ -2642,14 +2642,18 @@ subroutine precompute_inner_kernel() do x_loop = number_points_inner_kernel-1,0,-1 ! For numerical reasons the negative l terms are treated separately, to avoid ! infinity/nans as result. - if (l < 0) then + if (l <= 0) then recurence_factor = 1.0d0 else recurence_factor = ((x_loop*delta_x)/((x_loop+1)*delta_x))**l end if - table_inner_kernel_x_to_infty(x_loop, l, m) = & - & recurence_factor*table_inner_kernel_x_to_infty(x_loop+1, l, m) & - & + integrate_param(kernel_x_to_infty, x_loop*delta_x, x_loop*delta_x, (x_loop+1)*delta_x) + if ((x_loop == 0) .and. (l > 0)) then + table_inner_kernel_x_to_infty(x_loop, l, m) = 0.0d0 + else + table_inner_kernel_x_to_infty(x_loop, l, m) = & + & recurence_factor*table_inner_kernel_x_to_infty(x_loop+1, l, m) & + & + integrate_param(kernel_x_to_infty, x_loop*delta_x, x_loop*delta_x, (x_loop+1)*delta_x) + end if end do end do @@ -2670,6 +2674,10 @@ function kernel_x_to_infty(xp, xval) ! infinity/nans as result. if (l < 0) then kernel_x_to_infty = xp**(abs(l)+1) * exp(-xp**2) * phi_exp(m, xp) + else if (l == 0) then + kernel_x_to_infty = xp * exp(-xp**2) * phi_exp(m, xp) + else if (xval == 0.0_dp) then + kernel_x_to_infty = 0.0_dp else kernel_x_to_infty = xp * (xval/xp)**(l) * exp(-xp**2) * phi_exp(m, xp) end if diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 9e0aee10..28cc27cf 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -138,6 +138,19 @@ set_tests_properties(ntv_output_geometry_test PROPERTIES TIMEOUT 30 ) +add_executable(test_collop_zero_endpoint test_collop_zero_endpoint.f90) +target_link_libraries(test_collop_zero_endpoint common) + +add_test(NAME collop_zero_endpoint_test + COMMAND test_collop_zero_endpoint + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(collop_zero_endpoint_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_collop_zero_endpoint.f90 b/TEST/test_collop_zero_endpoint.f90 new file mode 100644 index 00000000..76fca932 --- /dev/null +++ b/TEST/test_collop_zero_endpoint.f90 @@ -0,0 +1,19 @@ +program test_collop_zero_endpoint + use, intrinsic :: ieee_arithmetic, only: ieee_get_flag, ieee_invalid, & + ieee_set_flag, ieee_set_halting_mode + use collop_compute, only: init_collop + use rkstep_mod, only: lag, leg + implicit none + + logical :: invalid_raised + + lag = 0 + leg = 1 + call ieee_set_flag(ieee_invalid, .false.) + call ieee_set_halting_mode(ieee_invalid, .true.) + call init_collop(0, 0, 0.0d0, 0.0d0) + call ieee_get_flag(ieee_invalid, invalid_raised) + + if (invalid_raised) error stop 'FAIL: zero endpoint raised IEEE invalid' + print *, 'All tests passed!' +end program test_collop_zero_endpoint From 03f5b9dacdbf7c882fac01c37972e8e266575351 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 12:32:02 +0200 Subject: [PATCH 15/23] Skip exact-zero energy scattering integral --- COMMON/collop_compute.f90 | 12 +++++++++++- TEST/test_collop_zero_endpoint.f90 | 19 ++++++++++++++++++- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/COMMON/collop_compute.f90 b/COMMON/collop_compute.f90 index b56aa4de..7a3531d4 100644 --- a/COMMON/collop_compute.f90 +++ b/COMMON/collop_compute.f90 @@ -119,6 +119,7 @@ module collop_compute integer(HID_T) :: h5id_matelem logical :: precomp=.false. logical :: make_ortho=.true. + logical :: expansion_has_constant_mode=.false. interface chop module procedure chop_0 @@ -305,20 +306,24 @@ subroutine init_collop(collop_base_prj, collop_base_exp, scalprod_alpha, scalpro stop end select + expansion_has_constant_mode = .false. select case (collop_base_exp) case (0) + expansion_has_constant_mode = .true. write (*,*) "Using Laguerre polynomials as collision operator expansion base." init_phi_exp => init_phi_laguerre phi_exp => phi_laguerre d_phi_exp => d_phi_laguerre dd_phi_exp => dd_phi_laguerre case (1) + expansion_has_constant_mode = .true. write (*,*) "Using standard polynomials as collision operator expansion base." init_phi_exp => init_phi_polynomial phi_exp => phi_polynomial d_phi_exp => d_phi_polynomial dd_phi_exp => dd_phi_polynomial case (2) + expansion_has_constant_mode = .true. write (*,*) "Using squared polynomials as collision operator expansion base." init_phi_exp => init_phi_polynomial_2 phi_exp => phi_polynomial_2 @@ -1923,7 +1928,12 @@ subroutine compute_energyscattering(denmm_s) if (isw_relativistic .eq. 0) then do m = 0, lagmax do mp = 0, lagmax - denmm_s(m+1, mp+1) = 3d0/(4d0 * pi) * integrate(integrand, 0d0) + if ((mp == 0) .and. expansion_has_constant_mode .and. & + (T_a == T_b)) then + denmm_s(m+1, mp+1) = 0d0 + else + denmm_s(m+1, mp+1) = 3d0/(4d0 * pi) * integrate(integrand, 0d0) + end if end do end do elseif (isw_relativistic .ge. 1) then diff --git a/TEST/test_collop_zero_endpoint.f90 b/TEST/test_collop_zero_endpoint.f90 index 76fca932..c01ab53b 100644 --- a/TEST/test_collop_zero_endpoint.f90 +++ b/TEST/test_collop_zero_endpoint.f90 @@ -1,19 +1,36 @@ program test_collop_zero_endpoint use, intrinsic :: ieee_arithmetic, only: ieee_get_flag, ieee_invalid, & ieee_set_flag, ieee_set_halting_mode - use collop_compute, only: init_collop + use collop_compute, only: compute_energyscattering, gamma_ab, & + init_collop, isw_relativistic, make_ortho, T_a, T_b use rkstep_mod, only: lag, leg implicit none logical :: invalid_raised + real(kind(1.0d0)) :: denmm(1, 1) lag = 0 leg = 1 call ieee_set_flag(ieee_invalid, .false.) call ieee_set_halting_mode(ieee_invalid, .true.) call init_collop(0, 0, 0.0d0, 0.0d0) + gamma_ab = 1.0d0 + isw_relativistic = 0 + make_ortho = .false. + T_a = 1.0d0 + T_b = 1.0d0 + call compute_energyscattering(denmm) call ieee_get_flag(ieee_invalid, invalid_raised) if (invalid_raised) error stop 'FAIL: zero endpoint raised IEEE invalid' + if (denmm(1, 1) /= 0.0d0) error stop 'FAIL: constant basis energy scattering is not zero' + + gamma_ab = sqrt(2.0d0) + T_a = 2.0d0 + call compute_energyscattering(denmm) + call ieee_get_flag(ieee_invalid, invalid_raised) + + if (invalid_raised) error stop 'FAIL: unequal-temperature energy scattering raised IEEE invalid' + if (denmm(1, 1) == 0.0d0) error stop 'FAIL: unequal-temperature energy scattering was skipped' print *, 'All tests passed!' end program test_collop_zero_endpoint From 6e30459047a026f52da453eb9af59889107647a6 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 12:40:40 +0200 Subject: [PATCH 16/23] Reject zero Lorentz join columns --- NEO-2-QL/join_diagnostics_mod.f90 | 23 ++++++++++++++++++++++- NEO-2-QL/join_ripples.f90 | 12 +++++++++++- TEST/test_join_failure_diagnostic.f90 | 11 ++++++++++- 3 files changed, 43 insertions(+), 3 deletions(-) diff --git a/NEO-2-QL/join_diagnostics_mod.f90 b/NEO-2-QL/join_diagnostics_mod.f90 index 87edcfb2..41bd162e 100644 --- a/NEO-2-QL/join_diagnostics_mod.f90 +++ b/NEO-2-QL/join_diagnostics_mod.f90 @@ -3,7 +3,7 @@ module join_diagnostics_mod implicit none private - public :: report_join_failure + public :: report_join_failure, validate_join_normalization contains @@ -24,4 +24,25 @@ subroutine report_join_failure(stage, info, old_tags, new_tags, ndim, ndim1, & count(ieee_is_finite(rhs)), '/', size(rhs) ierr = stage end subroutine report_join_failure + + subroutine validate_join_normalization(facnorm, direction, column, old_tags, & + new_tags, npass_l, npass_r, ierr) + real(kind(1.0d0)), intent(in) :: facnorm + character(len=1), intent(in) :: direction + integer, intent(in) :: column, old_tags(2), new_tags(2), npass_l, npass_r + integer, intent(out) :: ierr + + if (ieee_is_finite(facnorm) .and. facnorm /= 0.0d0) then + ierr = 0 + return + end if + + write(*, '(a,a)') 'join_ripples: normalization direction=', direction + write(*, '(a,i0)') 'join_ripples: normalization column=', column + write(*, '(a,es23.15)') 'join_ripples: normalization factor=', facnorm + write(*, '(a,2(1x,i0))') 'join_ripples: old propagator tags=', old_tags + write(*, '(a,2(1x,i0))') 'join_ripples: new propagator tags=', new_tags + write(*, '(a,2(1x,i0))') 'join_ripples: pass dimensions=', npass_l, npass_r + ierr = 5 + end subroutine validate_join_normalization end module join_diagnostics_mod diff --git a/NEO-2-QL/join_ripples.f90 b/NEO-2-QL/join_ripples.f90 index 2897b462..689c48a2 100644 --- a/NEO-2-QL/join_ripples.f90 +++ b/NEO-2-QL/join_ripples.f90 @@ -19,7 +19,7 @@ SUBROUTINE join_ripples(ierr) USE propagator_mod USE lapack_band USE collisionality_mod, ONLY : isw_lorentz - USE join_diagnostics_mod, ONLY : report_join_failure + USE join_diagnostics_mod, ONLY : report_join_failure, validate_join_normalization USE development IMPLICIT NONE @@ -291,11 +291,21 @@ SUBROUTINE join_ripples(ierr) IF(isw_lorentz.EQ.1) THEN DO i=1,o%p%npass_l facnorm=SUM(o%p%amat_p_p(:,i))+SUM(o%p%amat_p_m(:,i)) + CALL validate_join_normalization(facnorm,'p',i, & + [o%fieldpropagator_tag_s,o%fieldpropagator_tag_e], & + [n%fieldpropagator_tag_s,n%fieldpropagator_tag_e], & + o%p%npass_l,o%p%npass_r,ierr) + IF(ierr.NE.0) RETURN o%p%amat_p_p(:,i)=o%p%amat_p_p(:,i)/facnorm o%p%amat_p_m(:,i)=o%p%amat_p_m(:,i)/facnorm ENDDO DO i=1,o%p%npass_r facnorm=SUM(o%p%amat_m_p(:,i))+SUM(o%p%amat_m_m(:,i)) + CALL validate_join_normalization(facnorm,'m',i, & + [o%fieldpropagator_tag_s,o%fieldpropagator_tag_e], & + [n%fieldpropagator_tag_s,n%fieldpropagator_tag_e], & + o%p%npass_l,o%p%npass_r,ierr) + IF(ierr.NE.0) RETURN o%p%amat_m_p(:,i)=o%p%amat_m_p(:,i)/facnorm o%p%amat_m_m(:,i)=o%p%amat_m_m(:,i)/facnorm ENDDO diff --git a/TEST/test_join_failure_diagnostic.f90 b/TEST/test_join_failure_diagnostic.f90 index c373e620..0b9a9bee 100644 --- a/TEST/test_join_failure_diagnostic.f90 +++ b/TEST/test_join_failure_diagnostic.f90 @@ -1,5 +1,6 @@ program test_join_failure_diagnostic - use join_diagnostics_mod, only: report_join_failure + use join_diagnostics_mod, only: report_join_failure, & + validate_join_normalization use lapack_band, only: gbsv implicit none @@ -14,4 +15,12 @@ program test_join_failure_diagnostic call report_join_failure(2, info, [10, 20], [21, 30], 2, 2, & matrix, rhs, ierr) if (ierr /= 2) error stop 'FAIL: join stage was not retained' + + call validate_join_normalization(0.0d0, 'p', 2, [12, 13], [13, 13], & + 24, 24, ierr) + if (ierr /= 5) error stop 'FAIL: zero normalization was not rejected' + + call validate_join_normalization(1.0d0, 'p', 2, [12, 13], [13, 13], & + 24, 24, ierr) + if (ierr /= 0) error stop 'FAIL: finite normalization was rejected' end program test_join_failure_diagnostic From 37d7076c13f26a08eb4147e0b95aab4f461a50f0 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 12:57:10 +0200 Subject: [PATCH 17/23] Restrict Arnoldi solver to global route --- NEO-2-QL/flint.f90 | 29 ++++++++++++++++++++++++++++- TEST/CMakeLists.txt | 13 +++++++++++++ TEST/test_ripple_solver_route.f90 | 15 +++++++++++++++ 3 files changed, 56 insertions(+), 1 deletion(-) create mode 100644 TEST/test_ripple_solver_route.f90 diff --git a/NEO-2-QL/flint.f90 b/NEO-2-QL/flint.f90 index d2fb44ab..3654ed03 100644 --- a/NEO-2-QL/flint.f90 +++ b/NEO-2-QL/flint.f90 @@ -1,3 +1,22 @@ +MODULE ripple_solver_route_mod + IMPLICIT NONE + PRIVATE + + PUBLIC :: uses_single_propagator_route + +CONTAINS + + PURE LOGICAL FUNCTION uses_single_propagator_route(magnetic_device, & + isw_axisymm, mag_magfield) + INTEGER, INTENT(in) :: magnetic_device, isw_axisymm, mag_magfield + + uses_single_propagator_route = & + (magnetic_device .EQ. 0 .AND. isw_axisymm .EQ. 1) .OR. & + mag_magfield .EQ. 0 + END FUNCTION uses_single_propagator_route + +END MODULE ripple_solver_route_mod + !> this routine now walks through the specified fieldpropagators !> (see main routine about how to set proptag_start and proptag_end) SUBROUTINE flint(eta_part_globalfac,eta_part_globalfac_p,eta_part_globalfac_t, & @@ -21,6 +40,8 @@ SUBROUTINE flint(eta_part_globalfac,eta_part_globalfac_p,eta_part_globalfac_t, & prop_reconstruct, prop_reconstruct_levels USE collisionality_mod, ONLY : isw_lorentz,isw_axisymm,y_axi_averages, & isw_momentum + USE ntv_mod, ONLY : isw_ripple_solver + USE ripple_solver_route_mod, ONLY : uses_single_propagator_route USE rkstep_mod, ONLY : lag USE collisionality_mod, ONLY : collpar,collpar_min,collpar_max @@ -53,6 +74,7 @@ SUBROUTINE flint(eta_part_globalfac,eta_part_globalfac_p,eta_part_globalfac_t, & INTEGER :: x2_ub,i_min_sav,clear_old_ripple INTEGER :: proptag_first,proptag_last INTEGER :: start_at_begin + LOGICAL :: single_propagator_route INTEGER :: i_construct,i_construct_s,i_construct_e,i_construct_d INTEGER :: eta_min_loc(1) @@ -1909,7 +1931,12 @@ SUBROUTINE flint(eta_part_globalfac,eta_part_globalfac_p,eta_part_globalfac_t, & prop_ibegperiod = 1 prop_count_call = 0 PRINT *, 'Do the real computation - ripple_solver' - IF ( (magnetic_device .EQ. 0 .AND. isw_axisymm .EQ. 1) .OR. mag_magfield .EQ. 0 ) THEN + single_propagator_route = uses_single_propagator_route( & + magnetic_device,isw_axisymm,mag_magfield) + IF (isw_ripple_solver .EQ. 3 .AND. .NOT. single_propagator_route) THEN + ERROR STOP 'isw_ripple_solver=3 requires the single-propagator global route' + END IF + IF (single_propagator_route) THEN IF (mag_magfield .EQ. 0) THEN PRINT *, 'Use only one propagator for homogeneous field' ELSE diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 28cc27cf..0642a80e 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -125,6 +125,19 @@ set_tests_properties(join_failure_diagnostic_test PROPERTIES TIMEOUT 30 ) +add_executable(test_ripple_solver_route test_ripple_solver_route.f90) +target_link_libraries(test_ripple_solver_route neo2_ql) + +add_test(NAME ripple_solver_route_test + COMMAND test_ripple_solver_route + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(ripple_solver_route_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_ntv_output_geometry test_ntv_output_geometry.f90) target_link_libraries(test_ntv_output_geometry neo2_ql) diff --git a/TEST/test_ripple_solver_route.f90 b/TEST/test_ripple_solver_route.f90 new file mode 100644 index 00000000..4e1d3665 --- /dev/null +++ b/TEST/test_ripple_solver_route.f90 @@ -0,0 +1,15 @@ +program test_ripple_solver_route + use ripple_solver_route_mod, only: uses_single_propagator_route + implicit none + + if (.not. uses_single_propagator_route(0, 1, 1)) & + error stop 'FAIL: axisymmetric tokamak route was rejected' + if (.not. uses_single_propagator_route(1, 0, 0)) & + error stop 'FAIL: homogeneous-field route was rejected' + if (uses_single_propagator_route(0, 0, 1)) & + error stop 'FAIL: non-axisymmetric tokamak route was accepted' + if (uses_single_propagator_route(1, 0, 1)) & + error stop 'FAIL: general multi-propagator route was accepted' + + print *, 'All tests passed!' +end program test_ripple_solver_route From 12aebf2e8dfc588b0d3e1ac3e74b790a66c356ba Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 13:22:30 +0200 Subject: [PATCH 18/23] Remove duplicate top-level build test --- NEO-2-PAR/CMakeLists.txt | 2 +- NEO-2-QL/CMakeLists.txt | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/NEO-2-PAR/CMakeLists.txt b/NEO-2-PAR/CMakeLists.txt index 0b548c6e..9bca09f8 100644 --- a/NEO-2-PAR/CMakeLists.txt +++ b/NEO-2-PAR/CMakeLists.txt @@ -33,7 +33,7 @@ if (COPY_TO_ARCHIVE) endif() ADD_TEST(NAME Test_Build_Code COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}) -set_tests_properties(Test_Build_Code PROPERTIES RESOURCE_LOCK neo2_build_tree) +set_tests_properties(Test_Build_Code PROPERTIES RUN_SERIAL TRUE) ADD_TEST(Test_Testcase1-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Testcase1" "0") ADD_TEST(Test_Testcase1-Par ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Testcase1" "6") ADD_TEST(Test_Relativistic_BeliaevBudker-Seq ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Relativistic_BeliaevBudker" "0") diff --git a/NEO-2-QL/CMakeLists.txt b/NEO-2-QL/CMakeLists.txt index da30c470..ed7c2359 100644 --- a/NEO-2-QL/CMakeLists.txt +++ b/NEO-2-QL/CMakeLists.txt @@ -31,8 +31,6 @@ add_custom_target(copy ${PROJECT_SOURCE_DIR}/../../TEST_RUNS/${EXE_NAME} ) -ADD_TEST(NAME Test_Build_Code COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}) -set_tests_properties(Test_Build_Code PROPERTIES RESOURCE_LOCK neo2_build_tree) ADD_TEST(Extended_QL_Multispecies ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies" "2" "QL") ADD_TEST(Extended_QL_Multispecies_Read ${PROJECT_SOURCE_DIR}/../ShellScripts/test_full_reconstruct.sh "Extended_QL_Multispecies_Read" "2" "QL") set_tests_properties( From 47343f3b639408475b45a48335376de3d1d1c0ae Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 13:23:12 +0200 Subject: [PATCH 19/23] Run helical source through local solver --- NEO-2-QL/CMakeSources.in | 1 + NEO-2-QL/helical_source_mod.f90 | 145 ++++++++++++++++++ NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 | 135 ++-------------- NEO-2-QL/ripple_solver_axi_test.f90 | 118 ++++++++++---- TEST/CMakeLists.txt | 13 ++ TEST/test_helical_source.f90 | 57 +++++++ 6 files changed, 320 insertions(+), 149 deletions(-) create mode 100644 NEO-2-QL/helical_source_mod.f90 create mode 100644 TEST/test_helical_source.f90 diff --git a/NEO-2-QL/CMakeSources.in b/NEO-2-QL/CMakeSources.in index 23382ac4..eace6793 100644 --- a/NEO-2-QL/CMakeSources.in +++ b/NEO-2-QL/CMakeSources.in @@ -2,6 +2,7 @@ set(NEO2_QL_SRC_FILES arnoldi_mod.f90 er_rotation_mod.f90 ntv_mod.f90 + helical_source_mod.f90 propagator.f90 flint.f90 ripple_solver_axi_test.f90 diff --git a/NEO-2-QL/helical_source_mod.f90 b/NEO-2-QL/helical_source_mod.f90 new file mode 100644 index 00000000..6000aab3 --- /dev/null +++ b/NEO-2-QL/helical_source_mod.f90 @@ -0,0 +1,145 @@ +module helical_source_mod + use nrtype, only: dp + implicit none + private + + public :: add_helical_source + +contains + + subroutine add_helical_source(source, profile, moments, column, sigma_sign, & + ibeg, iend, lag, npl, ind_start, fact_pos_b, fact_pos_e, & + fact_neg_b, fact_neg_e) + complex(dp), intent(inout) :: source(:, :) + complex(dp), intent(in) :: profile(:, ibeg:) + real(dp), intent(in) :: moments(0:) + integer, intent(in) :: column, sigma_sign, ibeg, iend, lag + integer, intent(in) :: npl(ibeg:iend), ind_start(ibeg:iend) + real(dp), intent(in) :: fact_pos_b(ibeg:iend), fact_pos_e(ibeg:iend) + real(dp), intent(in) :: fact_neg_b(ibeg:iend), fact_neg_e(ibeg:iend) + + call add_forward_source(source, profile, moments, column, ibeg, iend, & + lag, npl, ind_start, fact_pos_b, fact_pos_e) + call add_backward_source(source, profile, moments, column, sigma_sign, & + ibeg, iend, lag, npl, ind_start, fact_neg_b, fact_neg_e) + end subroutine add_helical_source + + subroutine add_forward_source(source, profile, moments, column, ibeg, iend, & + lag, npl, ind_start, fact_b, fact_e) + complex(dp), intent(inout) :: source(:, :) + complex(dp), intent(in) :: profile(:, ibeg:) + real(dp), intent(in) :: moments(0:), fact_b(ibeg:iend), fact_e(ibeg:iend) + integer, intent(in) :: column, ibeg, iend, lag + integer, intent(in) :: npl(ibeg:iend), ind_start(ibeg:iend) + integer :: i, m, k, k_prev, np, np_prev, np_next + + do i = ibeg, iend + np = npl(i) + do m = 0, lag + k = ind_start(i) + 2*(np + 1)*m + if (i <= ibeg) cycle + np_prev = npl(i - 1) + k_prev = ind_start(i - 1) + 2*(np_prev + 1)*m + if (mod(i - ibeg, 2) == 1) then + np_next = npl(i + 1) + source(k + 1:k + np + 1, column) = source(k + 1:k + np + 1, column) & + + moments(m)/1.5d0*profile(1:np + 1, i)*fact_e(i) + source(k + 1:k + np_prev + 1, column) = source(k + 1:k + np_prev + 1, column) & + + moments(m)/2.4d0*profile(1:np_prev + 1, i - 1)*fact_b(i - 1) + source(k + 1:k + np_next + 1, column) = source(k + 1:k + np_next + 1, column) & + - moments(m)/12d0*profile(1:np_next + 1, i + 1)*fact_e(i + 1) + else + np_next = npl(i - 2) + source(k + 1:k + np + 1, column) = source(k + 1:k + np + 1, column) & + + moments(m)/2.4d0*profile(1:np + 1, i)*fact_e(i) + call add_forward_neighbor(source, profile, moments(m), column, k, k_prev, & + np, np_prev, i - 1, ibeg, 1.5d0, fact_b(i - 1)) + call add_forward_neighbor(source, profile, -moments(m), column, k, k_prev, & + np, np_next, i - 2, ibeg, 12d0, fact_b(i - 2)) + end if + end do + end do + end subroutine add_forward_source + + subroutine add_forward_neighbor(source, profile, moment, column, k, k_prev, & + np, np_neighbor, i_neighbor, ibeg, divisor, factor) + complex(dp), intent(inout) :: source(:, :) + complex(dp), intent(in) :: profile(:, ibeg:) + real(dp), intent(in) :: moment, divisor, factor + integer, intent(in) :: column, k, k_prev, np, np_neighbor, i_neighbor, ibeg + + if (np_neighbor <= np) then + source(k + 1:k + np_neighbor + 1, column) = & + source(k + 1:k + np_neighbor + 1, column) & + + moment/divisor*profile(1:np_neighbor + 1, i_neighbor)*factor + else + source(k + 1:k + np + 1, column) = source(k + 1:k + np + 1, column) & + + moment/divisor*profile(1:np + 1, i_neighbor)*factor + source(k_prev + np_neighbor + 2, column) = source(k_prev + np_neighbor + 2, column) & + + moment/divisor*profile(np_neighbor + 1, i_neighbor)*factor + end if + end subroutine add_forward_neighbor + + subroutine add_backward_source(source, profile, moments, column, sigma_sign, & + ibeg, iend, lag, npl, ind_start, fact_b, fact_e) + complex(dp), intent(inout) :: source(:, :) + complex(dp), intent(in) :: profile(:, ibeg:) + real(dp), intent(in) :: moments(0:), fact_b(ibeg:iend), fact_e(ibeg:iend) + integer, intent(in) :: column, sigma_sign, ibeg, iend, lag + integer, intent(in) :: npl(ibeg:iend), ind_start(ibeg:iend) + integer :: i, m, k, k_prev, np, np_prev, np_next + real(dp) :: signed_moment + + do i = ibeg, iend + np = npl(i) + do m = 0, lag + k = ind_start(i) + 2*(np + 1)*m + if (i >= iend) cycle + signed_moment = real(sigma_sign, dp)*moments(m) + np_prev = npl(i + 1) + k_prev = ind_start(i + 1) + 2*(np_prev + 1)*m + if (mod(i - ibeg, 2) == 1) then + np_next = npl(i - 1) + source(k + np + 2:k + 2*np + 2, column) = & + source(k + np + 2:k + 2*np + 2, column) & + + signed_moment/1.5d0*profile(np + 1:1:-1, i)*fact_e(i) + source(k + 2*np + 2 - np_prev:k + 2*np + 2, column) = & + source(k + 2*np + 2 - np_prev:k + 2*np + 2, column) & + + signed_moment/2.4d0*profile(np_prev + 1:1:-1, i + 1)*fact_b(i + 1) + source(k + 2*np + 2 - np_next:k + 2*np + 2, column) = & + source(k + 2*np + 2 - np_next:k + 2*np + 2, column) & + - signed_moment/12d0*profile(np_next + 1:1:-1, i - 1)*fact_e(i - 1) + else + np_next = npl(i + 2) + source(k + np + 2:k + 2*np + 2, column) = & + source(k + np + 2:k + 2*np + 2, column) & + + signed_moment/2.4d0*profile(np + 1:1:-1, i)*fact_e(i) + call add_backward_neighbor(source, profile, signed_moment, column, k, k_prev, & + np, np_prev, i + 1, ibeg, 1.5d0, fact_b(i + 1)) + call add_backward_neighbor(source, profile, -signed_moment, column, k, k_prev, & + np, np_next, i + 2, ibeg, 12d0, fact_b(i + 2)) + end if + end do + end do + end subroutine add_backward_source + + subroutine add_backward_neighbor(source, profile, moment, column, k, k_prev, & + np, np_neighbor, i_neighbor, ibeg, divisor, factor) + complex(dp), intent(inout) :: source(:, :) + complex(dp), intent(in) :: profile(:, ibeg:) + real(dp), intent(in) :: moment, divisor, factor + integer, intent(in) :: column, k, k_prev, np, np_neighbor, i_neighbor, ibeg + + if (np_neighbor <= np) then + source(k + 2*np + 2 - np_neighbor:k + 2*np + 2, column) = & + source(k + 2*np + 2 - np_neighbor:k + 2*np + 2, column) & + + moment/divisor*profile(np_neighbor + 1:1:-1, i_neighbor)*factor + else + source(k + np + 2:k + 2*np + 2, column) = & + source(k + np + 2:k + 2*np + 2, column) & + + moment/divisor*profile(np + 1:1:-1, i_neighbor)*factor + source(k_prev + np_neighbor + 1, column) = source(k_prev + np_neighbor + 1, column) & + + moment/divisor*profile(np_neighbor + 1, i_neighbor)*factor + end if + end subroutine add_backward_neighbor +end module helical_source_mod diff --git a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 index 7ba6358e..26bb7111 100644 --- a/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 +++ b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 @@ -92,6 +92,7 @@ subroutine ripple_solver_ArnoldiO2( & use collop use arnoldi_mod, only: iterator, f_init_arnoldi, & & lsw_write_flux_surface_distribution, write_flux_surface_distribution + use helical_source_mod, only: add_helical_source implicit none complex(dp), parameter :: imun = (0.d0, 1.d0) @@ -4261,133 +4262,23 @@ subroutine source_flux end do if (hel_drive_active .and. iplot .ne. 1) then - call add_hel_source(q_hel_b, asource(0:lag, 2), 1, -1) - call add_hel_source(q_hel_b, asource_hel(0:lag, 2), 3, -1) - call add_hel_source(q_hel_e, asource_hel(0:lag, 1), 1, 1) - call add_hel_source(q_hel_e, asource(0:lag, 1), 3, 1) + call add_helical_source(source_vector, q_hel_b, asource(0:lag, 2), & + 1, -1, ibeg, iend, lag, npl, ind_start, fact_pos_b, & + fact_pos_e, fact_neg_b, fact_neg_e) + call add_helical_source(source_vector, q_hel_b, asource_hel(0:lag, 2), & + 3, -1, ibeg, iend, lag, npl, ind_start, fact_pos_b, & + fact_pos_e, fact_neg_b, fact_neg_e) + call add_helical_source(source_vector, q_hel_e, asource_hel(0:lag, 1), & + 1, 1, ibeg, iend, lag, npl, ind_start, fact_pos_b, & + fact_pos_e, fact_neg_b, fact_neg_e) + call add_helical_source(source_vector, q_hel_e, asource(0:lag, 1), & + 3, 1, ibeg, iend, lag, npl, ind_start, fact_pos_b, & + fact_pos_e, fact_neg_b, fact_neg_e) end if end subroutine source_flux !------------------------------------------------------------------------ - subroutine add_hel_source(q_hel, amom, lcol, isig) - -! Adds one misalignment-drive piece (band profile q_hel, velocity moments -! amom) to source_vector column lcol, following the integration-scheme -! structure of the standard drives assembled above. isig = +1 for -! sigma-even pieces (sign pattern of columns 1 and 3), isig = -1 for -! sigma-odd pieces (sign pattern of column 2). - - complex(dp), dimension(:, ibeg:), intent(in) :: q_hel - real(dp), dimension(0:), intent(in) :: amom - integer, intent(in) :: lcol, isig - - integer :: is, mh, kh, kh_prev, np, np_prev, np_next - real(dp) :: sg - - sg = DBLE(isig) - - do is = ibeg, iend - np = npl(is) - do mh = 0, lag - kh = ind_start(is) + 2*(np + 1)*mh - - if (is .GT. ibeg) then - np_prev = npl(is - 1) - kh_prev = ind_start(is - 1) + 2*(np_prev + 1)*mh - if (MOD(is - ibeg, 2) .EQ. 1) then - np_next = npl(is + 1) - source_vector(kh + 1:kh + np + 1, lcol) & - = source_vector(kh + 1:kh + np + 1, lcol) & - + amom(mh)/1.5d0*q_hel(1:np + 1, is)*fact_pos_e(is) - source_vector(kh + 1:kh + np_prev + 1, lcol) & - = source_vector(kh + 1:kh + np_prev + 1, lcol) & - + amom(mh)/2.4d0*q_hel(1:np_prev + 1, is - 1)*fact_pos_b(is - 1) - source_vector(kh + 1:kh + np_next + 1, lcol) & - = source_vector(kh + 1:kh + np_next + 1, lcol) & - - amom(mh)/12d0*q_hel(1:np_next + 1, is + 1)*fact_pos_e(is + 1) - else - np_next = npl(is - 2) - source_vector(kh + 1:kh + np + 1, lcol) & - = source_vector(kh + 1:kh + np + 1, lcol) & - + amom(mh)/2.4d0*q_hel(1:np + 1, is)*fact_pos_e(is) - if (np_prev .LE. np) then - source_vector(kh + 1:kh + np_prev + 1, lcol) & - = source_vector(kh + 1:kh + np_prev + 1, lcol) & - + amom(mh)/1.5d0*q_hel(1:np_prev + 1, is - 1)*fact_pos_b(is - 1) - else - source_vector(kh + 1:kh + np + 1, lcol) & - = source_vector(kh + 1:kh + np + 1, lcol) & - + amom(mh)/1.5d0*q_hel(1:np + 1, is - 1)*fact_pos_b(is - 1) - source_vector(kh_prev + np_prev + 2, lcol) & - = source_vector(kh_prev + np_prev + 2, lcol) & - + amom(mh)/1.5d0*q_hel(np_prev + 1, is - 1)*fact_pos_b(is - 1) - end if - if (np_next .LE. np) then - source_vector(kh + 1:kh + np_next + 1, lcol) & - = source_vector(kh + 1:kh + np_next + 1, lcol) & - - amom(mh)/12d0*q_hel(1:np_next + 1, is - 2)*fact_pos_b(is - 2) - else - source_vector(kh + 1:kh + np + 1, lcol) & - = source_vector(kh + 1:kh + np + 1, lcol) & - - amom(mh)/12d0*q_hel(1:np + 1, is - 2)*fact_pos_b(is - 2) - source_vector(kh_prev + np_prev + 2, lcol) & - = source_vector(kh_prev + np_prev + 2, lcol) & - - amom(mh)/12d0*q_hel(np_next + 1, is - 2)*fact_pos_b(is - 2) - end if - end if - end if - - if (is .LT. iend) then - np_prev = npl(is + 1) - kh_prev = ind_start(is + 1) + 2*(np_prev + 1)*mh - if (MOD(is - ibeg, 2) .EQ. 1) then - np_next = npl(is - 1) - source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - + sg*amom(mh)/1.5d0*q_hel(np + 1:1:-1, is)*fact_neg_e(is) - source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & - = source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & - + sg*amom(mh)/2.4d0*q_hel(np_prev + 1:1:-1, is + 1)*fact_neg_b(is + 1) - source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & - = source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & - - sg*amom(mh)/12d0*q_hel(np_next + 1:1:-1, is - 1)*fact_neg_e(is - 1) - else - np_next = npl(is + 2) - source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - + sg*amom(mh)/2.4d0*q_hel(np + 1:1:-1, is)*fact_neg_e(is) - if (np_prev .LE. np) then - source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & - = source_vector(kh + 2*np + 2 - np_prev:kh + 2*np + 2, lcol) & - + sg*amom(mh)/1.5d0*q_hel(np_prev + 1:1:-1, is + 1)*fact_neg_b(is + 1) - else - source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - + sg*amom(mh)/1.5d0*q_hel(np + 1:1:-1, is + 1)*fact_neg_b(is + 1) - source_vector(kh_prev + np_prev + 1, lcol) & - = source_vector(kh_prev + np_prev + 1, lcol) & - + sg*amom(mh)/1.5d0*q_hel(np_prev + 1, is + 1)*fact_neg_b(is + 1) - end if - if (np_next .LE. np) then - source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & - = source_vector(kh + 2*np + 2 - np_next:kh + 2*np + 2, lcol) & - - sg*amom(mh)/12d0*q_hel(np_next + 1:1:-1, is + 2)*fact_neg_b(is + 2) - else - source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - = source_vector(kh + np + 2:kh + 2*np + 2, lcol) & - - sg*amom(mh)/12d0*q_hel(np + 1:1:-1, is + 2)*fact_neg_b(is + 2) - source_vector(kh_prev + np_prev + 1, lcol) & - = source_vector(kh_prev + np_prev + 1, lcol) & - - sg*amom(mh)/12d0*q_hel(np_next + 1, is + 2)*fact_neg_b(is + 2) - end if - end if - end if - - end do - end do - - end subroutine add_hel_source !------------------------------------------------------------------------ subroutine add_f01_source diff --git a/NEO-2-QL/ripple_solver_axi_test.f90 b/NEO-2-QL/ripple_solver_axi_test.f90 index c9e749f5..f0153838 100644 --- a/NEO-2-QL/ripple_solver_axi_test.f90 +++ b/NEO-2-QL/ripple_solver_axi_test.f90 @@ -54,7 +54,7 @@ SUBROUTINE ripple_solver( & USE collisionality_mod, ONLY : collpar,conl_over_mfp,isw_lorentz, & isw_energy,isw_integral,isw_axisymm, & !<-in Winny isw_momentum,vel_distri_swi,vel_num,vel_max, & - nvel,vel_array,num_spec + nvel,vel_array,num_spec,isw_relativistic,T_spec,m_spec,z_spec ! collpar - the same as $\kappa$ - inverse mean-free path times 4 USE lapack_band @@ -87,6 +87,13 @@ SUBROUTINE ripple_solver( & deallocate_ntv_eqmat use mpiprovider_module, only : mpro USE collop + USE helical_source_mod, ONLY : add_helical_source + USE neo_magfie, ONLY : boozer_iota,boozer_curr_tor_hat, & + boozer_curr_pol_hat,boozer_psi_pr_hat + USE ntv_mod, ONLY : isw_hel_drive,isw_m_phi_input,m_phi_input, & + m_theta_hel,hel_brad_re,hel_brad_im,hel_phim_re,hel_phim_im, & + clight_hel => c,echarge_hel => e + USE partpa_mod, ONLY : bmod0 !! End Modification by Andreas F. Martitsch (28.07.2015) IMPLICIT NONE @@ -149,6 +156,7 @@ SUBROUTINE ripple_solver( & DOUBLE PRECISION :: diflam,diflampow,coefdir DOUBLE PRECISION :: coefenu,coefenu_averb !!!term[1] DOUBLE PRECISION :: alambd_save1,alambd_save2,alambd_save3 + DOUBLE PRECISION :: denomjac,vt_hel,rho_hel DOUBLE PRECISION :: diflam_flux,coef_cf DOUBLE PRECISION :: hxeta,fun_bound_new,amin2ovb DOUBLE PRECISION :: a11,a12,a21,a22,determ,deleta_b,dum_phi,dum_a3 @@ -218,6 +226,11 @@ SUBROUTINE ripple_solver( & DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bvec_sp,bvec_iter,bvec_lor DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: bvec_prev DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: flux_vector,source_vector + COMPLEX(kind=dp), DIMENSION(:,:), ALLOCATABLE :: helical_source + COMPLEX(kind=dp), DIMENSION(:,:), ALLOCATABLE :: q_hel_b,q_hel_e + COMPLEX(kind=dp) :: amp_hel_b,amp_hel_e,hel_phase,hel_phase_fac + LOGICAL :: hel_drive_active + INTEGER :: n_hel INTEGER :: isw_lor,isw_ene,isw_intp INTEGER, DIMENSION(:), ALLOCATABLE :: npl DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: rhs_mat_fzero @@ -860,6 +873,10 @@ SUBROUTINE ripple_solver( & ALLOCATE(rhs_mat_lorentz(5,npart+1,ibeg:iend)) ALLOCATE(rhs_mat_energ(4,npart+1,ibeg:iend)) ALLOCATE(q_rip(npart+2,ibeg:iend,0:2)) + ALLOCATE(q_hel_b(npart+2,ibeg:iend),q_hel_e(npart+2,ibeg:iend)) + q_hel_b=(0.d0,0.d0) + q_hel_e=(0.d0,0.d0) + hel_drive_active=.FALSE. ALLOCATE(convol_flux(npart+1,ibeg:iend),convol_curr(npart+1,ibeg:iend)) ALLOCATE(ind_start(ibeg:iend)) if (iplot.EQ.1) then @@ -895,7 +912,7 @@ SUBROUTINE ripple_solver( & ierr=3 DEALLOCATE(deriv_coef,npl) DEALLOCATE(rhs_mat_lorentz,rhs_mat_energ) - DEALLOCATE(q_rip) + DEALLOCATE(q_rip,q_hel_b,q_hel_e) DEALLOCATE(convol_flux,convol_curr) DEALLOCATE(pleg_bra,pleg_ket,scalprod_pleg) return @@ -1213,6 +1230,49 @@ SUBROUTINE ripple_solver( & ENDDO + IF(isw_hel_drive.NE.0) THEN + IF(isw_m_phi_input.NE.1) & + ERROR STOP 'mode-1 helical drive requires explicit m_phi_input' + IF(.NOT.ALLOCATED(asource_hel)) & + ERROR STOP 'helical drive requires quadrature collision moments' + IF(isw_relativistic.NE.0) & + ERROR STOP 'helical drive relativistic moments are not validated' + IF(iplot.EQ.1) & + ERROR STOP 'mode-1 helical drive does not support reconstruction mode' + IF(T_spec(ispec).LE.0.d0.OR.m_spec(ispec).LE.0.d0.OR.z_spec(ispec).EQ.0.d0) & + ERROR STOP 'helical drive species data are invalid' + n_hel=m_phi_input + aiota=boozer_iota + denomjac=aiota*boozer_curr_tor_hat+boozer_curr_pol_hat + IF(denomjac.EQ.0.d0.OR.boozer_psi_pr_hat.EQ.0.d0) & + ERROR STOP 'helical drive geometry denominator is zero' + hel_drive_active=.TRUE. + hel_phase_fac=(0.d0,1.d0)*(DBLE(m_theta_hel)*aiota+DBLE(n_hel)) + vt_hel=SQRT(2.d0*T_spec(ispec)/m_spec(ispec)) + rho_hel=vt_hel*m_spec(ispec)*clight_hel & + /(z_spec(ispec)*echarge_hel*(bmod0*1.d4)) + amp_hel_b=CMPLX(hel_brad_re,hel_brad_im,kind=dp)/rho_hel + amp_hel_e=(0.d0,1.d0)*CMPLX(hel_phim_re,hel_phim_im,kind=dp) & + *z_spec(ispec)*echarge_hel/(2.d0*T_spec(ispec)) & + *(DBLE(n_hel)*boozer_curr_tor_hat-DBLE(m_theta_hel)*boozer_curr_pol_hat) & + /(boozer_psi_pr_hat*denomjac) + DO istep=ibeg,iend + npassing=npl(istep) + hel_phase=EXP(hel_phase_fac*phi_mfl(istep)) + q_hel_b(:,istep)=(0.d0,0.d0) + q_hel_e(:,istep)=(0.d0,0.d0) + q_hel_b(1:npassing,istep)=amp_hel_b*hel_phase & + *(eta(1:npassing)-eta(0:npassing-1)) + q_hel_b(npassing+1,istep)=amp_hel_b*hel_phase & + *(1.d0/bhat_mfl(istep)-eta(npassing)) + q_hel_e(1:npassing,istep)=amp_hel_e*hel_phase*2.d0 & + *(alambd(0:npassing-1,istep)-alambd(1:npassing,istep)) & + /(bhat_mfl(istep)*h_phi_mfl(istep)) + q_hel_e(npassing+1,istep)=amp_hel_e*hel_phase*2.d0 & + *alambd(npassing,istep)/(bhat_mfl(istep)*h_phi_mfl(istep)) + ENDDO + ENDIF + DEALLOCATE(amat,bvec_lapack,ipivot) !! Modification by Andreas F. Martitsch (16.09.2015) @@ -1580,6 +1640,25 @@ SUBROUTINE ripple_solver( & enddo enddo + IF(hel_drive_active.AND.iplot.NE.1) THEN + ALLOCATE(helical_source(n_2d_size,3)) + helical_source=(0.d0,0.d0) + CALL add_helical_source(helical_source,q_hel_b,asource(0:lag,2), & + 1,-1,ibeg,iend,lag,npl,ind_start,fact_pos_b,fact_pos_e, & + fact_neg_b,fact_neg_e) + CALL add_helical_source(helical_source,q_hel_b,asource_hel(0:lag,2), & + 3,-1,ibeg,iend,lag,npl,ind_start,fact_pos_b,fact_pos_e, & + fact_neg_b,fact_neg_e) + CALL add_helical_source(helical_source,q_hel_e,asource_hel(0:lag,1), & + 1,1,ibeg,iend,lag,npl,ind_start,fact_pos_b,fact_pos_e, & + fact_neg_b,fact_neg_e) + CALL add_helical_source(helical_source,q_hel_e,asource(0:lag,1), & + 3,1,ibeg,iend,lag,npl,ind_start,fact_pos_b,fact_pos_e, & + fact_neg_b,fact_neg_e) + source_vector=source_vector+REAL(helical_source,dp) + DEALLOCATE(helical_source) + ENDIF + ! Determine the size of arrays (number of non-zero elements): @@ -2741,36 +2820,21 @@ SUBROUTINE ripple_solver( & qflux_allspec=qflux_allspec_tmp IF(ALLOCATED(qflux_allspec_tmp)) DEALLOCATE(qflux_allspec_tmp) IF(mpro%getrank() .EQ. 0) THEN - ! D11 - PRINT *,'qflux(1,1,0,0):' - PRINT *,qflux_allspec(1,1,0,0) - PRINT *,'qflux(1,1,1,0):' - PRINT *,qflux_allspec(1,1,1,0) - PRINT *,'qflux(1,1,0,1):' - PRINT *,qflux_allspec(1,1,0,1) - PRINT *,'qflux(1,1,1,1):' - PRINT *,qflux_allspec(1,1,1,1) - ! D12 - PRINT *,'qflux(1,3,0,0):' - PRINT *,qflux_allspec(1,3,0,0) - PRINT *,'qflux(1,3,1,0):' - PRINT *,qflux_allspec(1,3,1,0) - PRINT *,'qflux(1,3,0,1):' - PRINT *,qflux_allspec(1,3,0,1) - PRINT *,'qflux(1,3,1,1):' - PRINT *,qflux_allspec(1,3,1,1) + PRINT *,'qflux(1,1,:,:):' + PRINT *,qflux_allspec(1,1,:,:) + PRINT *,'qflux(1,3,:,:):' + PRINT *,qflux_allspec(1,3,:,:) OPEN(070915,file='qflux_symm_allspec.dat') WRITE(070915,*) boozer_s, collpar, & - qflux_allspec(1,1,0,0), qflux_allspec(1,1,1,0), & - qflux_allspec(1,1,0,1), qflux_allspec(1,1,1,1), & - qflux_allspec(1,3,0,0), qflux_allspec(1,3,1,0), & - qflux_allspec(1,3,0,1), qflux_allspec(1,3,1,1) + ((qflux_allspec(1,1,ispecp,ispecpp),ispecp=0,num_spec-1), & + ispecpp=0,num_spec-1), & + ((qflux_allspec(1,3,ispecp,ispecpp),ispecp=0,num_spec-1), & + ispecpp=0,num_spec-1) CLOSE(070915) END IF !! End Modification by Andreas F. Martitsch (23.08.2015) call cpu_time(time2) -RETURN do m=0,lag do kk=1,3 @@ -2796,7 +2860,7 @@ SUBROUTINE ripple_solver( & if(isw_intp.eq.1) deallocate(bvec_iter,bvec_lor,bvec_prev) DEALLOCATE(deriv_coef,enu_coef,alambd,Vg_vp_over_B,scalprod_pleg) DEALLOCATE(alampow,vrecurr,dellampow,convol_polpow,pleg_bra,pleg_ket) - DEALLOCATE(npl,rhs_mat_fzero,rhs_mat_lorentz,rhs_mat_energ,q_rip) + DEALLOCATE(npl,rhs_mat_fzero,rhs_mat_lorentz,rhs_mat_energ,q_rip,q_hel_b,q_hel_e) DEALLOCATE(convol_flux,convol_curr,ind_start) DEALLOCATE(phi_mfl,bhat_mfl,geodcu_mfl,h_phi_mfl,dlogbdphi_mfl,eta) DEALLOCATE(delt_pos,delt_neg,fact_pos_b,fact_neg_b,fact_pos_e,fact_neg_e) @@ -3024,7 +3088,7 @@ SUBROUTINE ripple_solver( & PRINT *,' ' DEALLOCATE(deriv_coef,enu_coef,alambd,Vg_vp_over_B,scalprod_pleg) DEALLOCATE(alampow,vrecurr,dellampow,convol_polpow,pleg_bra,pleg_ket) - DEALLOCATE(npl,rhs_mat_fzero,rhs_mat_lorentz,rhs_mat_energ,q_rip) + DEALLOCATE(npl,rhs_mat_fzero,rhs_mat_lorentz,rhs_mat_energ,q_rip,q_hel_b,q_hel_e) DEALLOCATE(convol_flux,convol_curr,ind_start) DEALLOCATE(delt_pos,delt_neg,fact_pos_b,fact_neg_b,fact_pos_e,fact_neg_e) diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 0642a80e..6dd331f9 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -138,6 +138,19 @@ set_tests_properties(ripple_solver_route_test PROPERTIES TIMEOUT 30 ) +add_executable(test_helical_source test_helical_source.f90) +target_link_libraries(test_helical_source neo2_ql) + +add_test(NAME helical_source_test + COMMAND test_helical_source + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(helical_source_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_ntv_output_geometry test_ntv_output_geometry.f90) target_link_libraries(test_ntv_output_geometry neo2_ql) diff --git a/TEST/test_helical_source.f90 b/TEST/test_helical_source.f90 new file mode 100644 index 00000000..bacb77d3 --- /dev/null +++ b/TEST/test_helical_source.f90 @@ -0,0 +1,57 @@ +program test_helical_source + use helical_source_mod, only: add_helical_source + implicit none + + integer, parameter :: dp = kind(1.0d0), ibeg = 0, iend = 4 + complex(dp) :: even_source(24, 3), odd_source(24, 3), profile(3, ibeg:iend) + complex(dp) :: forward_even(24, 3), forward_odd(24, 3) + complex(dp) :: backward_even(24, 3), backward_odd(24, 3) + complex(dp) :: weighted_sum + real(dp) :: factors(ibeg:iend), moments(0:0), zero_factors(ibeg:iend) + integer :: i, ind_start(ibeg:iend), npl(ibeg:iend) + + even_source = (0.0d0, 0.0d0) + odd_source = (0.0d0, 0.0d0) + forward_even = (0.0d0, 0.0d0) + forward_odd = (0.0d0, 0.0d0) + backward_even = (0.0d0, 0.0d0) + backward_odd = (0.0d0, 0.0d0) + zero_factors = 0.0d0 + factors = [1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0] + moments = 2.0d0 + npl = [1, 2, 1, 2, 1] + ind_start = [0, 4, 10, 14, 20] + do i = ibeg, iend + profile(:, i) = cmplx([real(i + 1, dp), real(2*i + 1, dp), real(3*i + 2, dp)], & + [real(3*i + 1, dp), real(i + 2, dp), real(2*i + 4, dp)], kind=dp) + end do + + call add_helical_source(even_source, profile, moments, 1, 1, ibeg, iend, & + 0, npl, ind_start, factors, factors, factors, factors) + call add_helical_source(odd_source, profile, moments, 1, -1, ibeg, iend, & + 0, npl, ind_start, factors, factors, factors, factors) + call add_helical_source(forward_even, profile, moments, 1, 1, ibeg, iend, & + 0, npl, ind_start, factors, factors, zero_factors, zero_factors) + call add_helical_source(forward_odd, profile, moments, 1, -1, ibeg, iend, & + 0, npl, ind_start, factors, factors, zero_factors, zero_factors) + call add_helical_source(backward_even, profile, moments, 1, 1, ibeg, iend, & + 0, npl, ind_start, zero_factors, zero_factors, factors, factors) + call add_helical_source(backward_odd, profile, moments, 1, -1, ibeg, iend, & + 0, npl, ind_start, zero_factors, zero_factors, factors, factors) + + if (any(even_source(:, 2:3) /= (0.0d0, 0.0d0))) & + error stop 'FAIL: helical source changed an unselected force column' + if (maxval(abs(forward_even - forward_odd)) > 1.0d-14) & + error stop 'FAIL: sigma sign changed the forward source' + if (maxval(abs(backward_even + backward_odd)) > 1.0d-14) & + error stop 'FAIL: sigma-odd backward source has the wrong sign' + if (maxval(abs(aimag(even_source))) == 0.0d0) & + error stop 'FAIL: complex source quadrature was discarded' + weighted_sum = sum([(even_source(i, 1)*real(i, dp), i = 1, 24)]) + if (abs(weighted_sum - cmplx(11874.666666666666d0, 14272.0d0, dp)) > 1.0d-10) & + error stop 'FAIL: shared source stencil changed' + if (abs(sum(abs(even_source(:, 1))**2) - 111813.33333333333d0) > 1.0d-9) & + error stop 'FAIL: shared source stencil norm changed' + + print *, 'All tests passed!' +end program test_helical_source From 36d2cde68f8c77e72e96a605144903da3584ce4b Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 14:11:10 +0200 Subject: [PATCH 20/23] Write dimensional helical response output --- NEO-2-QL/CMakeSources.in | 1 + NEO-2-QL/helical_response_mod.f90 | 33 ++++++++++++ NEO-2-QL/ntv_mod.f90 | 20 +++---- NEO-2-QL/propagator.f90 | 87 ++++++++++++++++++++++++++++++- TEST/CMakeLists.txt | 13 +++++ TEST/test_helical_response.f90 | 32 ++++++++++++ 6 files changed, 175 insertions(+), 11 deletions(-) create mode 100644 NEO-2-QL/helical_response_mod.f90 create mode 100644 TEST/test_helical_response.f90 diff --git a/NEO-2-QL/CMakeSources.in b/NEO-2-QL/CMakeSources.in index eace6793..de5add8c 100644 --- a/NEO-2-QL/CMakeSources.in +++ b/NEO-2-QL/CMakeSources.in @@ -3,6 +3,7 @@ set(NEO2_QL_SRC_FILES er_rotation_mod.f90 ntv_mod.f90 helical_source_mod.f90 + helical_response_mod.f90 propagator.f90 flint.f90 ripple_solver_axi_test.f90 diff --git a/NEO-2-QL/helical_response_mod.f90 b/NEO-2-QL/helical_response_mod.f90 new file mode 100644 index 00000000..f9b27ad5 --- /dev/null +++ b/NEO-2-QL/helical_response_mod.f90 @@ -0,0 +1,33 @@ +module helical_response_mod + use nrtype, only: dp + use ntv_mod, only: c, e + implicit none + private + + character(len=*), parameter, public :: D_CURRENT_UNIT = 'cm^2 G/s' + public :: dimensional_current_coefficients + public :: reconstruct_complex_response + +contains + + subroutine dimensional_current_coefficients(gamma, temperature, charge_number, & + gamma_current, d_current) + real(dp), intent(in) :: gamma(3, 3), temperature, charge_number + real(dp), intent(out) :: gamma_current(2), d_current(2) + real(dp) :: scale + + if (temperature <= 0.0_dp) error stop 'helical response temperature must be positive' + if (charge_number == 0.0_dp) error stop 'helical response charge must be nonzero' + + gamma_current = gamma(3, 1:2) + scale = 2.0_dp*temperature*c/(charge_number*e) + d_current = gamma_current*scale + end subroutine dimensional_current_coefficients + + pure subroutine reconstruct_complex_response(real_run, imaginary_run, baseline, response) + real(dp), intent(in) :: real_run(:), imaginary_run(:), baseline(:) + complex(dp), intent(out) :: response(:) + + response = cmplx(real_run - baseline, -(imaginary_run - baseline), kind=dp) + end subroutine reconstruct_complex_response +end module helical_response_mod diff --git a/NEO-2-QL/ntv_mod.f90 b/NEO-2-QL/ntv_mod.f90 index d36d69e3..6724cc8a 100644 --- a/NEO-2-QL/ntv_mod.f90 +++ b/NEO-2-QL/ntv_mod.f90 @@ -1733,42 +1733,42 @@ SUBROUTINE write_multispec_output_a() CALL h5_add(h5id_multispec, 'D12_AX', D12_AX, LBOUND(D12_AX), UBOUND(D12_AX), & & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D13_AX', D13_AX, LBOUND(D13_AX), UBOUND(D13_AX), & - & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D21_AX', D21_AX, LBOUND(D21_AX), UBOUND(D21_AX), & & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D22_AX', D22_AX, LBOUND(D22_AX), UBOUND(D22_AX), & & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D23_AX', D23_AX, LBOUND(D23_AX), UBOUND(D23_AX), & - & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D31_AX', D31_AX, LBOUND(D31_AX), UBOUND(D31_AX), & - & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D32_AX', D32_AX, LBOUND(D32_AX), UBOUND(D32_AX), & - & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D33_AX', D33_AX, LBOUND(D33_AX), UBOUND(D33_AX), & - & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for axisymmetric solution', unit='cm^2 G^2/s') !-> non-axisymmetric solution CALL h5_add(h5id_multispec, 'D11_NA', D11_NA, LBOUND(D11_NA), UBOUND(D11_NA), & & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D12_NA', D12_NA, LBOUND(D12_NA), UBOUND(D12_NA), & & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D13_NA', D13_NA, LBOUND(D13_NA), UBOUND(D13_NA), & - & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D21_NA', D21_NA, LBOUND(D21_NA), UBOUND(D21_NA), & & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D22_NA', D22_NA, LBOUND(D22_NA), UBOUND(D22_NA), & & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') CALL h5_add(h5id_multispec, 'D23_NA', D23_NA, LBOUND(D23_NA), UBOUND(D23_NA), & - & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D31_NA', D31_NA, LBOUND(D31_NA), UBOUND(D31_NA), & - & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D32_NA', D32_NA, LBOUND(D32_NA), UBOUND(D32_NA), & - & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2 G/s') CALL h5_add(h5id_multispec, 'D33_NA', D33_NA, LBOUND(D33_NA), UBOUND(D33_NA), & - & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2/s') + & comment='dimensional diffusion coefficient for non-axisymmetric solution', unit='cm^2 G^2/s') ! add ExB rotation frequency and species Mach numbers IF (isw_calc_Er .GE. 1) THEN diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index bbeabfe4..add0af59 100644 --- a/NEO-2-QL/propagator.f90 +++ b/NEO-2-QL/propagator.f90 @@ -72,7 +72,7 @@ MODULE propagator_mod ! --------------------------------------------------------------------------- ! --------------------------------------------------------------------------- - use hdf5_tools, only : HID_T + use hdf5_tools, only : HID_T, h5_add, h5_close, h5_create USE binarysplit_mod @@ -899,6 +899,8 @@ SUBROUTINE diag_propagator_res(iend) call h5_close(h5id) + call write_helical_current_response(gamma_out) + OPEN(uw,file='fulltransp.dat',status='replace') WRITE (uw,'(6(1x,i4),1000(1x,e18.5))') & full_version, & @@ -963,6 +965,89 @@ SUBROUTINE diag_propagator_res(iend) DEALLOCATE(y) END SUBROUTINE diag_propagator_res + SUBROUTINE write_helical_current_response(gamma) + USE collisionality_mod, ONLY : num_spec, z_spec, m_spec, T_spec + USE helical_response_mod, ONLY : dimensional_current_coefficients + USE ntv_mod, ONLY : c, e, has_perturbation_file, isw_hel_drive, & + isw_qflux_NA, isw_ripple_solver, m_phi_input, m_theta_hel, & + hel_brad_re, hel_brad_im, hel_phim_re, hel_phim_im + USE partpa_mod, ONLY : bmod0 + + REAL(kind=dp), INTENT(in) :: gamma(3,3) + REAL(kind=dp) :: gamma_current(2), d_current(2), bref, vt, rho + INTEGER :: baseline_exact, requires_off_reference + CHARACTER(len=16) :: quadrature + + IF(isw_hel_drive.EQ.0.OR.isw_ripple_solver.NE.1.OR.num_spec.NE.1) RETURN + + CALL dimensional_current_coefficients(gamma,T_spec(0),z_spec(0), & + gamma_current,d_current) + bref=bmod0*1.0e4_dp + vt=SQRT(2.0_dp*T_spec(0)/m_spec(0)) + rho=vt*m_spec(0)*c/(z_spec(0)*e*bref) + baseline_exact=0 + IF(isw_qflux_NA.EQ.1.AND..NOT.has_perturbation_file()) baseline_exact=1 + requires_off_reference=1-baseline_exact + quadrature='mixed' + IF(hel_brad_im.EQ.0.0_dp.AND.hel_phim_im.EQ.0.0_dp) quadrature='real' + IF(hel_brad_re.EQ.0.0_dp.AND.hel_phim_re.EQ.0.0_dp) quadrature='imaginary' + + CALL h5_create('helical_current_response.h5',h5id,1) + CALL write_helical_response_fields(gamma_current,d_current,bref,vt,rho, & + baseline_exact,requires_off_reference,quadrature) + CALL h5_close(h5id) + END SUBROUTINE write_helical_current_response + + SUBROUTINE write_helical_response_fields(gamma_current,d_current,bref,vt, & + rho,baseline_exact,requires_off_reference,quadrature) + USE collisionality_mod, ONLY : num_spec, species_tag, z_spec, m_spec, & + n_spec, T_spec + USE device_mod, ONLY : surface + USE helical_response_mod, ONLY : D_CURRENT_UNIT + USE mag_interface_mod, ONLY : boozer_s + USE ntv_mod, ONLY : isw_hel_drive, isw_m_phi_input, isw_qflux_NA, & + isw_ripple_solver, m_phi_input, m_theta_hel, hel_brad_re, & + hel_brad_im, hel_phim_re, hel_phim_im + + REAL(kind=dp), INTENT(in) :: gamma_current(2),d_current(2),bref,vt,rho + INTEGER, INTENT(in) :: baseline_exact,requires_off_reference + CHARACTER(len=*), INTENT(in) :: quadrature + + CALL h5_add(h5id,'schema_version',1) + CALL h5_add(h5id,'response_kind','total_source_quadrature') + CALL h5_add(h5id,'input_amplitude_quadrature',TRIM(quadrature)) + CALL h5_add(h5id,'phase_convention','exp(+i(m_theta*theta+m_phi*phi))') + CALL h5_add(h5id,'isw_ripple_solver',isw_ripple_solver) + CALL h5_add(h5id,'isw_hel_drive',isw_hel_drive) + CALL h5_add(h5id,'isw_qflux_NA',isw_qflux_NA) + CALL h5_add(h5id,'isw_m_phi_input',isw_m_phi_input) + CALL h5_add(h5id,'num_spec',num_spec) + CALL h5_add(h5id,'m_theta',m_theta_hel) + CALL h5_add(h5id,'m_phi',m_phi_input) + CALL h5_add(h5id,'boozer_s',boozer_s,unit='1') + CALL h5_add(h5id,'aiota',surface%aiota,unit='1') + CALL h5_add(h5id,'lambda_parallel', & + DBLE(m_theta_hel)*surface%aiota+DBLE(m_phi_input),unit='1') + CALL h5_add(h5id,'species_tag',species_tag(0)) + CALL h5_add(h5id,'charge_number',z_spec(0),unit='1') + CALL h5_add(h5id,'mass',m_spec(0),unit='g') + CALL h5_add(h5id,'density',n_spec(0),unit='1/cm^3') + CALL h5_add(h5id,'temperature',T_spec(0),unit='erg') + CALL h5_add(h5id,'Bref',bref,unit='G') + CALL h5_add(h5id,'thermal_velocity',vt,unit='cm/s') + CALL h5_add(h5id,'signed_larmor_radius',rho,unit='cm') + CALL h5_add(h5id,'hel_brad_re',hel_brad_re,unit='1') + CALL h5_add(h5id,'hel_brad_im',hel_brad_im,unit='1') + CALL h5_add(h5id,'hel_phim_re',hel_phim_re,unit='statV') + CALL h5_add(h5id,'hel_phim_im',hel_phim_im,unit='statV') + CALL h5_add(h5id,'gamma_current_A1',gamma_current(1),unit='1') + CALL h5_add(h5id,'gamma_current_A2',gamma_current(2),unit='1') + CALL h5_add(h5id,'D31',d_current(1),unit=D_CURRENT_UNIT) + CALL h5_add(h5id,'D32',d_current(2),unit=D_CURRENT_UNIT) + CALL h5_add(h5id,'zero_source_baseline_exact',baseline_exact) + CALL h5_add(h5id,'requires_off_reference',requires_off_reference) + END SUBROUTINE write_helical_response_fields + ! --------------------------------------------------------------------------- SUBROUTINE diag_propagator_dis diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 6dd331f9..e572e0ac 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -151,6 +151,19 @@ set_tests_properties(helical_source_test PROPERTIES TIMEOUT 30 ) +add_executable(test_helical_response test_helical_response.f90) +target_link_libraries(test_helical_response neo2_ql) + +add_test(NAME helical_response_test + COMMAND test_helical_response + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(helical_response_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_ntv_output_geometry test_ntv_output_geometry.f90) target_link_libraries(test_ntv_output_geometry neo2_ql) diff --git a/TEST/test_helical_response.f90 b/TEST/test_helical_response.f90 new file mode 100644 index 00000000..8d74f75c --- /dev/null +++ b/TEST/test_helical_response.f90 @@ -0,0 +1,32 @@ +program test_helical_response + use helical_response_mod, only: D_CURRENT_UNIT, & + dimensional_current_coefficients, reconstruct_complex_response + implicit none + + integer, parameter :: dp = kind(1.0d0) + real(dp) :: gamma(3, 3), gamma_current(2), d_current(2) + real(dp) :: real_run(2), imaginary_run(2), baseline(2) + complex(dp) :: response(2) + + gamma = 0.0_dp + gamma(3, 1:2) = [2.0_dp, -3.0_dp] + call dimensional_current_coefficients(gamma, 4.0_dp, -2.0_dp, & + gamma_current, d_current) + + if (any(gamma_current /= [2.0_dp, -3.0_dp])) & + error stop 'FAIL: current channels do not map to gamma rows 3 and columns 1:2' + if (maxval(abs(d_current - [-4.9931712191872093e20_dp, & + 7.489756828780813e20_dp])) > 1.0e6_dp) & + error stop 'FAIL: dimensional D31/D32 normalization changed' + if (D_CURRENT_UNIT /= 'cm^2 G/s') & + error stop 'FAIL: D31/D32 unit metadata is wrong' + real_run = [5.0_dp, 8.0_dp] + imaginary_run = [1.0_dp, 7.0_dp] + baseline = [2.0_dp, 3.0_dp] + call reconstruct_complex_response(real_run, imaginary_run, baseline, response) + if (maxval(abs(response - [cmplx(3.0_dp, 1.0_dp, dp), & + cmplx(5.0_dp, -4.0_dp, dp)])) > 1.0e-14_dp) & + error stop 'FAIL: real-imaginary quadrature reconstruction has the wrong sign' + + print *, 'All tests passed!' +end program test_helical_response From 723b5ccea66af0e2f53ae744fb58ac0dd0997010 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 15:30:11 +0200 Subject: [PATCH 21/23] Guard spline weight scans --- COMMON/spline_cof.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index 8897a2d1..d16daec0 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -916,7 +916,8 @@ END FUNCTION f ! find first regular point shift = 1 - DO WHILE ( ( shift <= SIZE(w) ) .AND. ( w(shift) == 0 ) ) + DO WHILE ( shift <= SIZE(w) ) + IF ( w(shift) /= 0 ) EXIT shift = shift + 1 END DO @@ -1424,7 +1425,8 @@ end function f ! find first regular point shift = 1 - do while ( ( shift <= size(w) ) .AND. ( w(shift) == 0 ) ) + do while ( shift <= size(w) ) + if ( w(shift) /= 0 ) exit shift = shift + 1 end do From 3350b7309a857b52a299bfb38439d539f7b34f7f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 15:30:23 +0200 Subject: [PATCH 22/23] Allocate reconstructed qflux input --- NEO-2-QL/propagator.f90 | 2 ++ TEST/CMakeLists.txt | 13 +++++++++++ TEST/test_propagator_io.f90 | 43 +++++++++++++++++++++++++++++++++++++ 3 files changed, 58 insertions(+) create mode 100644 TEST/test_propagator_io.f90 diff --git a/NEO-2-QL/propagator.f90 b/NEO-2-QL/propagator.f90 index add0af59..d9aa5f45 100644 --- a/NEO-2-QL/propagator.f90 +++ b/NEO-2-QL/propagator.f90 @@ -2331,6 +2331,8 @@ SUBROUTINE read_propagator_cont(o,prop_type,prop_start,prop_end,prop_showall_in) ! qflux IF (prop_showall .EQ. 1) THEN + IF (ALLOCATED(o%p%qflux)) DEALLOCATE(o%p%qflux) + ALLOCATE(o%p%qflux(3,3)) READ(prop_unit,*) o%p%qflux END IF diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index e572e0ac..5f9b4f1a 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -164,6 +164,19 @@ set_tests_properties(helical_response_test PROPERTIES TIMEOUT 30 ) +add_executable(test_propagator_io test_propagator_io.f90) +target_link_libraries(test_propagator_io neo2_ql) + +add_test(NAME propagator_io_test + COMMAND test_propagator_io + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(propagator_io_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_ntv_output_geometry test_ntv_output_geometry.f90) target_link_libraries(test_ntv_output_geometry neo2_ql) diff --git a/TEST/test_propagator_io.f90 b/TEST/test_propagator_io.f90 new file mode 100644 index 00000000..16d697ec --- /dev/null +++ b/TEST/test_propagator_io.f90 @@ -0,0 +1,43 @@ +program test_propagator_io + use, intrinsic :: iso_fortran_env, only: dp => real64 + use propagator_mod, only: propagator, read_propagator_content, & + write_propagator_content + implicit none + + type(propagator), pointer :: loaded, written + integer :: i, unit + real(dp) :: expected(3, 3) + + allocate (loaded, written) + written%nr_joined = 0 + written%fieldpropagator_tag_s = 41 + written%fieldpropagator_tag_e = 41 + written%fieldperiod_tag_s = 7 + written%fieldperiod_tag_e = 7 + written%phi_l = 0.25_dp + written%phi_r = 0.5_dp + written%p%npart = 1 + written%p%npass_l = 1 + written%p%npass_r = 1 + written%p%nvelocity = 0 + written%p%eta_boundary_l = 1.0_dp + written%p%eta_boundary_r = 1.0_dp + expected = reshape([(real(i, dp), i = 1, 9)], shape(expected)) + allocate (written%p%qflux(3, 3)) + written%p%qflux = expected + + call write_propagator_content(written, 3) + call read_propagator_content(loaded, 3, 41, 41, 1) + + if (.not. allocated(loaded%p%qflux)) & + error stop 'FAIL: qflux was not allocated while reading propagator content' + if (any(shape(loaded%p%qflux) /= [3, 3])) & + error stop 'FAIL: reconstructed qflux has the wrong shape' + if (any(loaded%p%qflux /= expected)) & + error stop 'FAIL: reconstructed qflux values changed during serialization' + + open (newunit=unit, file='propagator_41_41.prop', status='old') + close (unit, status='delete') + deallocate (loaded, written) + print *, 'All tests passed!' +end program test_propagator_io From 4c257e9b462defb61c614b9b1b15681acbcaeef2 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 11 Jul 2026 15:30:37 +0200 Subject: [PATCH 23/23] Preserve helical source during replay --- NEO-2-QL/helical_source_mod.f90 | 30 +++++++++++++++++++++++ NEO-2-QL/ripple_solver_axi_test.f90 | 12 ++++++---- TEST/test_helical_source.f90 | 37 ++++++++++++++++++++++++++++- 3 files changed, 74 insertions(+), 5 deletions(-) diff --git a/NEO-2-QL/helical_source_mod.f90 b/NEO-2-QL/helical_source_mod.f90 index 6000aab3..c6b9f979 100644 --- a/NEO-2-QL/helical_source_mod.f90 +++ b/NEO-2-QL/helical_source_mod.f90 @@ -4,9 +4,39 @@ module helical_source_mod private public :: add_helical_source + public :: apply_reconstructed_incoming_rows contains + subroutine apply_reconstructed_incoming_rows(source, flux_left, flux_right, & + ibeg, iend, lag, npl, ind_start) + real(dp), intent(inout) :: source(:, :) + real(dp), intent(in) :: flux_left(:, :), flux_right(:, :) + integer, intent(in) :: ibeg, iend, lag + integer, intent(in) :: npl(ibeg:iend), ind_start(ibeg:iend) + integer :: k_left, k_right, m, nleft, nright + + nleft = npl(ibeg) + 1 + nright = npl(iend) + 1 + if (size(flux_left, 1) /= nleft*(lag + 1)) & + error stop 'reconstructed left boundary has the wrong size' + if (size(flux_right, 1) /= nright*(lag + 1)) & + error stop 'reconstructed right boundary has the wrong size' + if (size(flux_left, 2) /= size(source, 2)) & + error stop 'reconstructed left boundary has the wrong columns' + if (size(flux_right, 2) /= size(source, 2)) & + error stop 'reconstructed right boundary has the wrong columns' + + do m = 0, lag + k_left = ind_start(ibeg) + 2*nleft*m + k_right = ind_start(iend) + 2*nright*m + source(k_left + 1:k_left + nleft, :) = & + flux_left(nleft*m + 1:nleft*(m + 1), :) + source(k_right + nright + 1:k_right + 2*nright, :) = & + flux_right(nright*(m + 1):nright*m + 1:-1, :) + end do + end subroutine apply_reconstructed_incoming_rows + subroutine add_helical_source(source, profile, moments, column, sigma_sign, & ibeg, iend, lag, npl, ind_start, fact_pos_b, fact_pos_e, & fact_neg_b, fact_neg_e) diff --git a/NEO-2-QL/ripple_solver_axi_test.f90 b/NEO-2-QL/ripple_solver_axi_test.f90 index f0153838..66bb002d 100644 --- a/NEO-2-QL/ripple_solver_axi_test.f90 +++ b/NEO-2-QL/ripple_solver_axi_test.f90 @@ -87,7 +87,8 @@ SUBROUTINE ripple_solver( & deallocate_ntv_eqmat use mpiprovider_module, only : mpro USE collop - USE helical_source_mod, ONLY : add_helical_source + USE helical_source_mod, ONLY : add_helical_source, & + apply_reconstructed_incoming_rows USE neo_magfie, ONLY : boozer_iota,boozer_curr_tor_hat, & boozer_curr_pol_hat,boozer_psi_pr_hat USE ntv_mod, ONLY : isw_hel_drive,isw_m_phi_input,m_phi_input, & @@ -868,6 +869,8 @@ SUBROUTINE ripple_solver( & ALLOCATE(convol_polpow(0:legmax,1:npart+3)) ! terms[2,3] ALLOCATE(pleg_bra(0:legmax,1:npart+1,ibeg:iend)) ! terms[2,3] ALLOCATE(pleg_ket(0:legmax,1:npart+1,ibeg:iend)) ! terms[2,3] + pleg_bra=0.d0 + pleg_ket=0.d0 ALLOCATE(npl(ibeg:iend)) ALLOCATE(rhs_mat_fzero(4,ibeg:iend,0:1)) ALLOCATE(rhs_mat_lorentz(5,npart+1,ibeg:iend)) @@ -1237,8 +1240,6 @@ SUBROUTINE ripple_solver( & ERROR STOP 'helical drive requires quadrature collision moments' IF(isw_relativistic.NE.0) & ERROR STOP 'helical drive relativistic moments are not validated' - IF(iplot.EQ.1) & - ERROR STOP 'mode-1 helical drive does not support reconstruction mode' IF(T_spec(ispec).LE.0.d0.OR.m_spec(ispec).LE.0.d0.OR.z_spec(ispec).EQ.0.d0) & ERROR STOP 'helical drive species data are invalid' n_hel=m_phi_input @@ -1640,7 +1641,7 @@ SUBROUTINE ripple_solver( & enddo enddo - IF(hel_drive_active.AND.iplot.NE.1) THEN + IF(hel_drive_active) THEN ALLOCATE(helical_source(n_2d_size,3)) helical_source=(0.d0,0.d0) CALL add_helical_source(helical_source,q_hel_b,asource(0:lag,2), & @@ -1657,6 +1658,9 @@ SUBROUTINE ripple_solver( & fact_neg_b,fact_neg_e) source_vector=source_vector+REAL(helical_source,dp) DEALLOCATE(helical_source) + IF(iplot.EQ.1.AND.isw_axisymm.NE.1) & + CALL apply_reconstructed_incoming_rows(source_vector,flux_pl,flux_mr, & + ibeg,iend,lag,npl,ind_start) ENDIF diff --git a/TEST/test_helical_source.f90 b/TEST/test_helical_source.f90 index bacb77d3..ad919f79 100644 --- a/TEST/test_helical_source.f90 +++ b/TEST/test_helical_source.f90 @@ -1,14 +1,21 @@ program test_helical_source - use helical_source_mod, only: add_helical_source + use helical_source_mod, only: add_helical_source, & + apply_reconstructed_incoming_rows implicit none integer, parameter :: dp = kind(1.0d0), ibeg = 0, iend = 4 + integer, parameter :: replay_ibeg = 0, replay_iend = 2, replay_lag = 1 complex(dp) :: even_source(24, 3), odd_source(24, 3), profile(3, ibeg:iend) complex(dp) :: forward_even(24, 3), forward_odd(24, 3) complex(dp) :: backward_even(24, 3), backward_odd(24, 3) complex(dp) :: weighted_sum real(dp) :: factors(ibeg:iend), moments(0:0), zero_factors(ibeg:iend) integer :: i, ind_start(ibeg:iend), npl(ibeg:iend) + real(dp) :: replay_source(24, 3), original_source(24, 3) + real(dp) :: flux_left(4, 3), flux_right(4, 3) + integer :: replay_ind_start(replay_ibeg:replay_iend) + integer :: replay_npl(replay_ibeg:replay_iend) + logical :: incoming_row(24) even_source = (0.0d0, 0.0d0) odd_source = (0.0d0, 0.0d0) @@ -53,5 +60,33 @@ program test_helical_source if (abs(sum(abs(even_source(:, 1))**2) - 111813.33333333333d0) > 1.0d-9) & error stop 'FAIL: shared source stencil norm changed' + replay_source = reshape([(real(i, dp), i = 1, size(replay_source))], & + shape(replay_source)) + original_source = replay_source + flux_left = reshape([(100.0_dp + real(i, dp), i = 1, size(flux_left))], & + shape(flux_left)) + flux_right = reshape([(200.0_dp + real(i, dp), i = 1, size(flux_right))], & + shape(flux_right)) + replay_npl = 1 + replay_ind_start = [0, 8, 16] + call apply_reconstructed_incoming_rows(replay_source, flux_left, flux_right, & + replay_ibeg, replay_iend, replay_lag, replay_npl, replay_ind_start) + + if (any(replay_source(1:2, :) /= flux_left(1:2, :))) & + error stop 'FAIL: first left incoming Laguerre block was not restored' + if (any(replay_source(5:6, :) /= flux_left(3:4, :))) & + error stop 'FAIL: second left incoming Laguerre block was not restored' + if (any(replay_source(19:20, :) /= flux_right(2:1:-1, :))) & + error stop 'FAIL: first right incoming Laguerre block was not restored' + if (any(replay_source(23:24, :) /= flux_right(4:3:-1, :))) & + error stop 'FAIL: second right incoming Laguerre block was not restored' + incoming_row = .false. + incoming_row([1, 2, 5, 6, 19, 20, 23, 24]) = .true. + do i = 1, size(replay_source, 1) + if (incoming_row(i)) cycle + if (any(replay_source(i, :) /= original_source(i, :))) & + error stop 'FAIL: reconstruction changed an interior or outgoing row' + end do + print *, 'All tests passed!' end program test_helical_source