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 7a42d8b3..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 @@ -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 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/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/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 diff --git a/DOC/neo2.in.ql-full b/DOC/neo2.in.ql-full index 36baf957..2b666a0f 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,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] diff --git a/NEO-2-PAR/CMakeLists.txt b/NEO-2-PAR/CMakeLists.txt index 341b6fdb..9bca09f8 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 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" +) 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..ed7c2359 100644 --- a/NEO-2-QL/CMakeLists.txt +++ b/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") diff --git a/NEO-2-QL/CMakeSources.in b/NEO-2-QL/CMakeSources.in index 9717a122..de5add8c 100644 --- a/NEO-2-QL/CMakeSources.in +++ b/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 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/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/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/join_diagnostics_mod.f90 b/NEO-2-QL/join_diagnostics_mod.f90 new file mode 100644 index 00000000..41bd162e --- /dev/null +++ b/NEO-2-QL/join_diagnostics_mod.f90 @@ -0,0 +1,48 @@ +module join_diagnostics_mod + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + implicit none + private + + public :: report_join_failure, validate_join_normalization + +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 + + 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 f8028148..689c48a2 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, validate_join_normalization 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 @@ -281,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/NEO-2-QL/neo2.f90 b/NEO-2-QL/neo2.f90 index e9355952..a3506b51 100644 --- a/NEO-2-QL/neo2.f90 +++ b/NEO-2-QL/neo2.f90 @@ -61,10 +61,14 @@ 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, & + isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im, has_perturbation_file, & + has_helical_drive_source ! derivative of iota for non-local NTV computations ! (with magnetic shear) @@ -246,7 +250,10 @@ 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, & + isw_hel_drive, m_theta_hel, hel_brad_re, hel_brad_im, & + hel_phim_re, hel_phim_im contains @@ -660,11 +667,19 @@ 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) 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) @@ -891,11 +906,20 @@ 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 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 @@ -927,6 +951,18 @@ subroutine check() write(*,*) 'Aborting...' stop end if + + if (isw_qflux_NA == 1 .and. .not. has_perturbation_file()) 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.' + 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 @@ -1158,10 +1194,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..6724cc8a 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. @@ -43,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) @@ -112,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 @@ -193,8 +210,47 @@ MODULE ntv_mod MODULE PROCEDURE compute_TphiNA_a END INTERFACE compute_TphiNA + PUBLIC has_perturbation_file, has_helical_drive_source + 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 + 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() + 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 + + 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 @@ -716,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 @@ -732,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: ! ---------------------------------------------------------------! @@ -769,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 @@ -1678,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 7de989a6..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, & @@ -953,7 +955,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) @@ -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/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 b/NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90 index d85cb5b9..26bb7111 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) @@ -90,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) @@ -218,6 +221,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 +293,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 +934,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 +1023,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 +3196,61 @@ 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 + 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)) + 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 +3383,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 +4261,25 @@ subroutine source_flux end do end do + if (hel_drive_active .and. iplot .ne. 1) then + 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_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 cb45b4aa..e572e0ac 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -50,6 +50,146 @@ 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_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_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_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_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_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_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_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) + +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_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..c01ab53b --- /dev/null +++ b/TEST/test_collop_zero_endpoint.f90 @@ -0,0 +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: 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 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 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 diff --git a/TEST/test_join_failure_diagnostic.f90 b/TEST/test_join_failure_diagnostic.f90 new file mode 100644 index 00000000..0b9a9bee --- /dev/null +++ b/TEST/test_join_failure_diagnostic.f90 @@ -0,0 +1,26 @@ +program test_join_failure_diagnostic + use join_diagnostics_mod, only: report_join_failure, & + validate_join_normalization + 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' + + 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 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 diff --git a/TEST/test_ntv_output_geometry.f90 b/TEST/test_ntv_output_geometry.f90 new file mode 100644 index 00000000..c17c2848 --- /dev/null +++ b/TEST/test_ntv_output_geometry.f90 @@ -0,0 +1,24 @@ +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' + 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) + 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 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 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 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 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