Write dimensional helical response output

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)

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
@@ -1507,6 +1512,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 +1549,34 @@ 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);
+    ! 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
+       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)
+    end if
+
     write (*,*) "Computing weighting coefficients..."
 
     if (allocated(C_m)) deallocate(C_m)
@@ -1615,6 +1649,25 @@ 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.
+    ! 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**(3*k + alpha) * exp(-(1+beta)*x**2) * phi_prj(m, x)
+    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**(3*k + alpha) * exp(-(exp_maxwell+beta)*x**2) &
+           * phi_prj(m, x)
+    end function am_hel_rel1
+
     function am_rel1(x)
       real(kind=dp) :: x, am_rel1
       real(kind=dp) :: exp_maxwell, gam, gam_fac
@@ -1875,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
@@ -2594,14 +2652,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
 
@@ -2622,6 +2684,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

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

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

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,12 +181,24 @@
                        ! 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.
              ! 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]

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 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")
@@ -48,3 +49,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"
+)

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

NEO-2-QL/CMakeLists.txt

@@ -31,7 +31,11 @@ add_custom_target(copy
   ${PROJECT_SOURCE_DIR}/../../TEST_RUNS/${EXE_NAME}
   )
 
-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")

NEO-2-QL/CMakeSources.in

@@ -2,10 +2,13 @@ set(NEO2_QL_SRC_FILES
   arnoldi_mod.f90
   er_rotation_mod.f90
   ntv_mod.f90
+  helical_source_mod.f90
+  helical_response_mod.f90
   propagator.f90
   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

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

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