From feb92b1f9bb0d4b8b455afe15f22e75322833fc8 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 4 Dec 2025 16:48:10 +0100 Subject: [PATCH] Add radial plasma profiles for collision calculations Add support for radially-varying plasma profiles (density, temperature) enabling realistic collisionality calculations. New features: - profiles module with power_series and two_power profile types - Power series: f(s) = scale * sum(coef(n) * s^n) - Two-power: f(s) = scale * (1 - s^p1)^p2 (non-integer exponents) - Collision coefficients precomputed on s grid and interpolated - Backward compatible: existing simulations work unchanged Files: - src/profiles.f90: New profiles module with namelist parsing - src/collis_alphas.f90: Refactored with profile support - src/simple_main.f90: Initialize profiles before collisions - test/tests/test_profiles.f90: Unit tests for profile evaluation Closes #202 --- examples/simple_full.in | 25 ++ src/CMakeLists.txt | 1 + src/collis_alphas.f90 | 569 ++++++++++++------------- src/profiles.f90 | 134 ++++++ src/simple_main.f90 | 36 +- test/golden_record/collision/simple.in | 11 + test/tests/CMakeLists.txt | 5 + test/tests/test_profiles.f90 | 348 +++++++++++++++ 8 files changed, 825 insertions(+), 304 deletions(-) create mode 100644 src/profiles.f90 create mode 100644 test/golden_record/collision/simple.in create mode 100644 test/tests/test_profiles.f90 diff --git a/examples/simple_full.in b/examples/simple_full.in index 1ad1bda3..c08a9799 100644 --- a/examples/simple_full.in +++ b/examples/simple_full.in @@ -60,3 +60,28 @@ reuse_batch = .False. ! Reuse batch from last run. Previous indices are stored i output_orbits_macrostep = .False. ! Output orbit positions to fort.9XXXX output_error = .False. ! Output additional error diagnostics / + +&profiles +! Radial plasma profiles for collision calculations. +! profile_type can be "power_series" or "two_power". +! Default values give uniform profiles matching the old scalar parameters above. +profile_type = "two_power" + +! Scale factors (on-axis values) +Te_scale = 1.0d4 ! electron temperature in eV +Ti1_scale = 1.0d4 ! ion species 1 temperature in eV +Ti2_scale = 1.0d4 ! ion species 2 temperature in eV +ni1_scale = 0.5d20 ! ion species 1 density in m^-3 +ni2_scale = 0.5d20 ! ion species 2 density in m^-3 + +! Two-power profile: f(s) = scale * (1 - s^p1)^p2 +! p2 = 0 gives flat profile, p2 > 0 gives peaked profile +Te_p1 = 1.0d0, Te_p2 = 0.0d0 ! flat electron temperature +Ti1_p1 = 1.0d0, Ti1_p2 = 0.0d0 ! flat ion 1 temperature +Ti2_p1 = 1.0d0, Ti2_p2 = 0.0d0 ! flat ion 2 temperature +ni1_p1 = 1.0d0, ni1_p2 = 0.0d0 ! flat ion 1 density +ni2_p1 = 1.0d0, ni2_p2 = 0.0d0 ! flat ion 2 density + +! Power series profile: f(s) = scale * sum(coef(n) * s^(n-1)) +! Te_coef = 1.0d0, 0.0d0, ... ! coefficient array (up to 10 terms) +/ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 24c28c40..f88130ac 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -41,6 +41,7 @@ classification.f90 simple.f90 parse_ants.f90 + profiles.f90 collis_alphas.f90 netcdf_orbit_output.f90 netcdf_results_output.f90 diff --git a/src/collis_alphas.f90 b/src/collis_alphas.f90 index 5fe9b7a2..9e088678 100644 --- a/src/collis_alphas.f90 +++ b/src/collis_alphas.f90 @@ -1,301 +1,274 @@ - ! module collis_alp - implicit none - - ! Define real(wp) kind parameter - integer, parameter :: wp = kind(1.0d0) - integer, parameter :: nsorts=3 - real(wp), dimension(nsorts) :: efcolf,velrat,enrat - - contains - ! - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! - subroutine coleff(p,dpp,dhh,fpeff) - ! - ! Computes local values of dimensionless contravariant components - ! of collisional diffusion tensor and friction force for nonrelativistic - ! plasma. Backgound temperature is the same for all sorts. - ! - ! Input variables: - ! formal: p - dimensionless momentum module (p/(sqrt(2)*p_T) - ! common: efcolf - dmls collision frequencies - ! velrat - ratio of test species thermal velocity to - ! background species thermal velocity - ! Output variables: - ! formal: dpp - dimensionless momentum module diffusion - ! coefficient - ! dhh - dimensionless pitch angle diffusion coeff. - ! fpeff - effective dimensionless drag force (prop. to linear - ! deviation in Fokker-Planck eq.) - ! - integer :: i - real(wp) :: p,dpp,dhh,fpeff,plim,xbeta,dp,dh,dpd - ! - plim=max(p,1.d-8) - ! - dpp=0.0d0 - dhh=0.0d0 - fpeff=0.0d0 - ! - do i=1,nsorts - xbeta=p*velrat(i) - ! - call onseff(xbeta,dp,dh,dpd) - ! - dpp=dpp+dp*efcolf(i) - dhh=dhh+dh*efcolf(i) - fpeff=fpeff+(dpd/plim-2.0*dp*p*enrat(i))*efcolf(i) - enddo - ! - dhh=dhh/plim**2 - ! - return - end - ! - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! - subroutine onseff(v,dp,dh,dpd) - ! - ! dp - dimensionless dpp - ! dh - dhh*p^2 (p - dmls) - ! dpd - (1/p)(d/dp)p^2*dp (p - dmls) - ! - implicit none - ! - ! square root of pi - real(wp), parameter :: sqp=1.7724538d0 - ! cons=4./(3.*sqrt(pi)) - real(wp), parameter :: cons=.75225278d0 - real(wp) :: v,dp,dh,dpd,v2,v3,ex,er - ! - v2=v**2 - v3=v2*v - if(v.lt.0.01d0) then - dp=cons*(1.d0-0.6d0*v2) - dh=cons*(1.d0-0.2d0*v2) - dpd=2.d0*cons*(1.d0-1.2d0*v2) - elseif(v.gt.6.d0) then - dp=1.d0/v3 - dh=(1.d0-0.5d0/v2)/v - dpd=-1.d0/v3 - else - ex=exp(-v2)/sqp - er=erf(v) - dp=er/v3-2.d0*ex/v2 - dh=er*(1.d0-0.5d0/v2)/v+ex/v2 - dpd=4.d0*ex-dp - endif - ! - return - end - ! - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! FUNCTION ERF(X) - ! PARAMETER ( A1 = 0.07052 30784, A2 = 0.04228 20123, - ! , A3 = 0.00927 05272, A4 = 0.00015 10143, - ! , A5 = 0.00027 65672, A6 = 0.00004 30638 ) - ! F(T) = 1./((1.+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+T*A6))))))**4)**4 - ! W = 1. - F(ABS(X)) - ! ERF = SIGN(W,X) - ! RETURN - ! END - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! - subroutine loacol_alpha(am1,am2,Z1,Z2,densi1,densi2,tempi1,tempi2,tempe,ealpha, & - v0,dchichi,slowrate,dchichi_norm,slowrate_norm) - ! - ! Performs precomputation of the constants for Coulomb collision - ! operator for alpha-particles colliding with 2 sorts of ions and electrons - ! - ! - ! Normalisation: test particle (alpha) velocity is normalized by v0, alpha-particle - ! birth velocity, time is multiplied with v0 and has a meaning of free path of alpha - ! particle with v0. Respectively, collision frequencies have the meaning of inverse - ! mean free paths. - ! - ! Input variables: - ! formal: am1,am2 - mass numbers of the first and second background ion species - ! Z1,Z2 - charge numbers of these species - ! densi1,densi2 - densities of these species, 1/cm**3 - ! tempi1,tempi2,tempe - temperatures of two ion species and electrons, eV - ! ealpha - initial alpha particle energy, eV - ! Output variables: - ! formal: v0 - initial alpha particle velocity, cm/s - ! dchichi - pitch angle scattering frequency, $D^{\chi\chi}$, of alpha - ! particle with initial velocity, 1/s - ! slowrate - slowing down rate, $F^v / v_0$, of alpha particle with - ! initial velocity, 1/s - ! dchichi_norim - normalized pitch angle scattering frequency, 1/cm - ! slowrate_norm - normalized slowing down rate, 1/cm - ! module collis_alp: - ! efcolf - normalized collision frequencies - ! velrat - ratio of initial alpha particle velocity v0 to the - ! specific background particle thermal velocity $v_{t}=\sqrt(2T/m)$ - ! enrat - ratio of initial alpha particle energy to the background species - ! energy - ! - real(wp) :: am1,am2,Z1,Z2,densi1,densi2,tempi1,tempi2,tempe,ealpha,dense - real(wp) :: v0,dchichi,slowrate,dchichi_norm,slowrate_norm,vti1,vti2,vte - real(wp) :: pi,pmass,emass,e,ev,alame,frecol_base,alami1,alami2 - real(wp) :: p,dpp,dhh,fpeff - ! - pi=3.14159265358979d0 - pmass=1.6726d-24 - emass=9.1094d-28 - e=4.8032d-10 - ev=1.6022d-12 - ! - enrat(1)=ealpha/tempi1 - enrat(2)=ealpha/tempi2 - enrat(3)=ealpha/tempe - ! - v0=sqrt(2.d0*ealpha*ev/(4.d0*pmass)) - vti1=sqrt(2.d0*tempi1*ev/(pmass*am1)) - vti2=sqrt(2.d0*tempi2*ev/(pmass*am2)) - vte=sqrt(2.d0*tempe*ev/emass) - ! - velrat(1)=v0/vti1 - velrat(2)=v0/vti2 - velrat(3)=v0/vte - ! - dense=densi1*Z1+densi2*Z2 - alami1=23.d0-log(max(epsilon(1.d0), & - sqrt(densi1*Z1**2/tempi1)*2.d0*Z1*(4.d0+am1)/(4.d0*tempi1+am1*ealpha))) - alami2=23.d0-log(max(epsilon(1.d0), & - sqrt(densi2*Z2**2/tempi2)*2.d0*Z2*(4.d0+am2)/(4.d0*tempi2+am2*ealpha))) - alame=24.d0-log(sqrt(dense)/tempe) - frecol_base=2.d0*pi*dense*e**4*2.d0**2/((4.d0*pmass)**2*v0**3) !usual - frecol_base=frecol_base/v0 !normalized - ! - efcolf(1)=frecol_base*Z1**2*alami1*densi1/dense - efcolf(2)=frecol_base*Z2**2*alami2*densi2/dense - efcolf(3)=frecol_base*alame - ! - efcolf=efcolf*velrat - ! - p=1.d0 - ! - call coleff(p,dpp,dhh,fpeff) - ! - dchichi=dhh*v0 - slowrate=abs(fpeff)*v0 - dchichi_norm=dhh - slowrate_norm=abs(fpeff) - return - end - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! - subroutine stost(z,dtauc,iswmode,ierr) - ! - ! z(1:5) - phase space coordinates: z(1:3) - spatial position, - ! z(4) - normalized - ! z(5) - pitch parameter - ! dtauc - normalized time step: dtauc=dtau*v0 - has the dimension of length - ! iswmode - switch of the collision mode: - ! 1 - full operator (pitch-angle and energy scattering and drag) - ! 2 - energy scattering and drag only - ! 3 - drag only - ! 4 - pitch-angle scattering only - ! ierr - error code: - ! 0 - good, - ! 1 - bad argument (pitch |z(5)| > 1 ), - ! 2 - step over pitch exceeds 1 (pitch was - ! replaced by randomly distributed on [-1,1]), - ! 3 - new pitch module exceeds 1, reflection from - ! the boudary was performed, - ! 10 or >10 - new momentum module is less then - ! prescribed minimum, reflection was performed. - ! - integer :: iswmode,ierr - real(wp), parameter :: pmin=1.e-8 - real(wp) :: dtauc,p,dpp,dhh,fpeff,alam,dalam,coala - real(wp), dimension(5) :: z - real :: ur - ! - p=z(4) - ! - call coleff(p,dpp,dhh,fpeff) - ! - ierr=0 - ! - if(iswmode.eq.1.or.iswmode.eq.4) then - alam=z(5) - coala=1.d0-alam**2 - ! - if(coala.lt.0.d0) then - ierr=1 - return - endif - ! - call getran(1,ur) - ! - dalam=sqrt(2.d0*dhh*coala*dtauc)*dble(ur)-2.d0*alam*dhh*dtauc - ! - if(abs(dalam).gt.1.d0) then - ierr=2 - ! - call random_number(ur) - ! - alam=2.d0*(dble(ur)-0.5d0) - else - alam=alam+dalam - if(alam.gt.1.d0) then - ierr=3 - alam=2.d0-alam - elseif(alam.lt.-1.d0) then - ierr=3 - alam=-2.d0-alam - endif - endif - ! - z(5)=alam - if(iswmode.eq.4) return - endif - ! - if(iswmode.lt.3) then - ! - call getran(0,ur) - ! - z(4)=z(4)+sqrt(abs(2.d0*dpp*dtauc))*dble(ur)+fpeff*dtauc - else - z(4)=z(4)+fpeff*dtauc - endif - ! - if(z(4).lt.pmin) then - ierr=ierr+10 - z(4)=pmin+abs(pmin-z(4)) - endif - ! - return - end - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - subroutine getran(irand,ur) - ! - ! Produces the random number with zero deviation and unit square - ! deviation - ! - ! Input parameters: irand - 0 for continious, 1 for +1 -1, - ! Output parameters: ur - random number - ! - integer, intent(in) :: irand - real, intent(out) :: ur - ! - call random_number(ur) - ! - if(irand.eq.0) then - ! continiuos random number, constant is sqrt(12) - ur=3.464102*(ur-.5) - else - ! discrete random number - if(ur.gt..5) then - ur=1. - else - ur=-1. - endif - endif - end - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + use, intrinsic :: iso_fortran_env, only: wp => real64 + implicit none + + integer, parameter :: nsorts = 3 + integer, parameter :: N_S_GRID = 101 + + real(wp), dimension(nsorts) :: efcolf, velrat, enrat + real(wp), dimension(nsorts, N_S_GRID) :: efcolf_grid, velrat_grid, enrat_grid + +contains + + subroutine coleff(p, dpp, dhh, fpeff) + real(wp), intent(in) :: p + real(wp), intent(out) :: dpp, dhh, fpeff + + call coleff_local(p, efcolf, velrat, enrat, dpp, dhh, fpeff) + end subroutine coleff + + subroutine coleff_local(p, efcolf_loc, velrat_loc, enrat_loc, dpp, dhh, fpeff) + real(wp), intent(in) :: p + real(wp), intent(in) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) + real(wp), intent(out) :: dpp, dhh, fpeff + + integer :: i + real(wp) :: plim, xbeta, dp, dh, dpd + + plim = max(p, 1.d-8) + + dpp = 0.0d0 + dhh = 0.0d0 + fpeff = 0.0d0 + + do i = 1, nsorts + xbeta = p*velrat_loc(i) + call onseff(xbeta, dp, dh, dpd) + dpp = dpp + dp*efcolf_loc(i) + dhh = dhh + dh*efcolf_loc(i) + fpeff = fpeff + (dpd/plim - 2.0d0*dp*p*enrat_loc(i))*efcolf_loc(i) + end do + + dhh = dhh/plim**2 + end subroutine coleff_local + + subroutine onseff(v, dp, dh, dpd) + real(wp), intent(in) :: v + real(wp), intent(out) :: dp, dh, dpd + + real(wp), parameter :: sqp = 1.7724538d0 + real(wp), parameter :: cons = 0.75225278d0 + real(wp) :: v2, v3, ex, er + + v2 = v**2 + v3 = v2*v + if (v < 0.01d0) then + dp = cons*(1.d0 - 0.6d0*v2) + dh = cons*(1.d0 - 0.2d0*v2) + dpd = 2.d0*cons*(1.d0 - 1.2d0*v2) + else if (v > 6.d0) then + dp = 1.d0/v3 + dh = (1.d0 - 0.5d0/v2)/v + dpd = -1.d0/v3 + else + ex = exp(-v2)/sqp + er = erf(v) + dp = er/v3 - 2.d0*ex/v2 + dh = er*(1.d0 - 0.5d0/v2)/v + ex/v2 + dpd = 4.d0*ex - dp + end if + end subroutine onseff + + subroutine loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, & + tempe, ealpha, v0, dchichi, slowrate, & + dchichi_norm, slowrate_norm) + real(wp), intent(in) :: am1, am2, Z1, Z2, densi1, densi2 + real(wp), intent(in) :: tempi1, tempi2, tempe, ealpha + real(wp), intent(out) :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm + + call compute_collision_coeffs(am1, am2, Z1, Z2, densi1, densi2, & + tempi1, tempi2, tempe, ealpha, v0, & + efcolf, velrat, enrat) + + call compute_rates(v0, dchichi, slowrate, dchichi_norm, slowrate_norm) + end subroutine loacol_alpha + + subroutine compute_collision_coeffs(am1, am2, Z1, Z2, densi1, densi2, & + tempi1, tempi2, tempe, ealpha, v0, & + efcolf_out, velrat_out, enrat_out) + real(wp), intent(in) :: am1, am2, Z1, Z2, densi1, densi2 + real(wp), intent(in) :: tempi1, tempi2, tempe, ealpha + real(wp), intent(out) :: v0 + real(wp), intent(out) :: efcolf_out(nsorts), velrat_out(nsorts), enrat_out(nsorts) + + real(wp), parameter :: pi = 3.14159265358979d0 + real(wp), parameter :: pmass = 1.6726d-24 + real(wp), parameter :: emass = 9.1094d-28 + real(wp), parameter :: e_cgs = 4.8032d-10 + real(wp), parameter :: ev_cgs = 1.6022d-12 + + real(wp) :: dense, vti1, vti2, vte + real(wp) :: alami1, alami2, alame, frecol_base + + enrat_out(1) = ealpha/tempi1 + enrat_out(2) = ealpha/tempi2 + enrat_out(3) = ealpha/tempe + + v0 = sqrt(2.d0*ealpha*ev_cgs/(4.d0*pmass)) + vti1 = sqrt(2.d0*tempi1*ev_cgs/(pmass*am1)) + vti2 = sqrt(2.d0*tempi2*ev_cgs/(pmass*am2)) + vte = sqrt(2.d0*tempe*ev_cgs/emass) + + velrat_out(1) = v0/vti1 + velrat_out(2) = v0/vti2 + velrat_out(3) = v0/vte + + dense = densi1*Z1 + densi2*Z2 + alami1 = 23.d0 - log(max(epsilon(1.d0), & + sqrt(densi1*Z1**2/tempi1)*2.d0*Z1*(4.d0 + am1)/(4.d0*tempi1 + am1*ealpha))) + alami2 = 23.d0 - log(max(epsilon(1.d0), & + sqrt(densi2*Z2**2/tempi2)*2.d0*Z2*(4.d0 + am2)/(4.d0*tempi2 + am2*ealpha))) + alame = 24.d0 - log(sqrt(dense)/tempe) + + frecol_base = 2.d0*pi*dense*e_cgs**4*2.d0**2/((4.d0*pmass)**2*v0**3) + frecol_base = frecol_base/v0 + + efcolf_out(1) = frecol_base*Z1**2*alami1*densi1/dense + efcolf_out(2) = frecol_base*Z2**2*alami2*densi2/dense + efcolf_out(3) = frecol_base*alame + + efcolf_out = efcolf_out*velrat_out + end subroutine compute_collision_coeffs + + subroutine compute_rates(v0, dchichi, slowrate, dchichi_norm, slowrate_norm) + real(wp), intent(in) :: v0 + real(wp), intent(out) :: dchichi, slowrate, dchichi_norm, slowrate_norm + + real(wp) :: p, dpp, dhh, fpeff + + p = 1.d0 + call coleff(p, dpp, dhh, fpeff) + + dchichi = dhh*v0 + slowrate = abs(fpeff)*v0 + dchichi_norm = dhh + slowrate_norm = abs(fpeff) + end subroutine compute_rates + + subroutine init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) + use simple_profiles, only: get_plasma_params + + real(wp), intent(in) :: am1, am2, Z1, Z2, ealpha, v0 + real(wp) :: s, Te, Ti1, Ti2, ni1, ni2 + real(wp) :: densi1_loc, densi2_loc, tempi1_loc, tempi2_loc, tempe_loc + real(wp) :: v0_loc + integer :: i + + do i = 1, N_S_GRID + s = dble(i - 1)/dble(N_S_GRID - 1) + + call get_plasma_params(s, Te, Ti1, Ti2, ni1, ni2) + + tempe_loc = max(Te, 1.0d0) + tempi1_loc = max(Ti1, 1.0d0) + tempi2_loc = max(Ti2, 1.0d0) + densi1_loc = max(ni1*1.0d-6, 1.0d0) + densi2_loc = max(ni2*1.0d-6, 1.0d0) + + call compute_collision_coeffs(am1, am2, Z1, Z2, & + densi1_loc, densi2_loc, & + tempi1_loc, tempi2_loc, tempe_loc, & + ealpha, v0_loc, & + efcolf_grid(:, i), velrat_grid(:, i), enrat_grid(:, i)) + end do + end subroutine init_collision_profiles + + subroutine get_local_coeffs(s, efcolf_loc, velrat_loc, enrat_loc) + real(wp), intent(in) :: s + real(wp), intent(out) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) + + real(wp) :: s_clamped, s_norm, weight + integer :: idx_lo, idx_hi + + s_clamped = max(0.0d0, min(1.0d0, s)) + s_norm = s_clamped*dble(N_S_GRID - 1) + idx_lo = max(1, min(N_S_GRID - 1, int(s_norm) + 1)) + idx_hi = idx_lo + 1 + weight = s_norm - dble(idx_lo - 1) + + efcolf_loc = (1.0d0 - weight)*efcolf_grid(:, idx_lo) + weight*efcolf_grid(:, idx_hi) + velrat_loc = (1.0d0 - weight)*velrat_grid(:, idx_lo) + weight*velrat_grid(:, idx_hi) + enrat_loc = (1.0d0 - weight)*enrat_grid(:, idx_lo) + weight*enrat_grid(:, idx_hi) + end subroutine get_local_coeffs + + subroutine stost(z, dtauc, iswmode, ierr) + integer, intent(in) :: iswmode + real(wp), intent(in) :: dtauc + real(wp), intent(inout) :: z(5) + integer, intent(out) :: ierr + + real(wp), parameter :: pmin = 1.d-8 + real(wp) :: p, dpp, dhh, fpeff, alam, dalam, coala + real(wp) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) + real :: ur + + p = z(4) + + call get_local_coeffs(z(1), efcolf_loc, velrat_loc, enrat_loc) + call coleff_local(p, efcolf_loc, velrat_loc, enrat_loc, dpp, dhh, fpeff) + + ierr = 0 + + if (iswmode == 1 .or. iswmode == 4) then + alam = z(5) + coala = 1.d0 - alam**2 + + if (coala < 0.d0) then + ierr = 1 + return + end if + + call getran(1, ur) + + dalam = sqrt(2.d0*dhh*coala*dtauc)*dble(ur) - 2.d0*alam*dhh*dtauc + + if (abs(dalam) > 1.d0) then + ierr = 2 + call random_number(ur) + alam = 2.d0*(dble(ur) - 0.5d0) + else + alam = alam + dalam + if (alam > 1.d0) then + ierr = 3 + alam = 2.d0 - alam + else if (alam < -1.d0) then + ierr = 3 + alam = -2.d0 - alam + end if + end if + + z(5) = alam + if (iswmode == 4) return + end if + + if (iswmode < 3) then + call getran(0, ur) + z(4) = z(4) + sqrt(abs(2.d0*dpp*dtauc))*dble(ur) + fpeff*dtauc + else + z(4) = z(4) + fpeff*dtauc + end if + + if (z(4) < pmin) then + ierr = ierr + 10 + z(4) = pmin + abs(pmin - z(4)) + end if + end subroutine stost + + subroutine getran(irand, ur) + integer, intent(in) :: irand + real, intent(out) :: ur + + call random_number(ur) + + if (irand == 0) then + ur = 3.464102*(ur - 0.5) + else + if (ur > 0.5) then + ur = 1.0 + else + ur = -1.0 + end if + end if + end subroutine getran end module collis_alp diff --git a/src/profiles.f90 b/src/profiles.f90 new file mode 100644 index 00000000..deec7732 --- /dev/null +++ b/src/profiles.f90 @@ -0,0 +1,134 @@ +module simple_profiles + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + + integer, parameter :: MAX_POWER_SERIES = 10 + + integer, parameter :: POWER_SERIES = 1 + integer, parameter :: TWO_POWER = 2 + + character(len=32) :: profile_type = "power_series" + integer :: active_profile = POWER_SERIES + + real(dp) :: Te_scale = 1.0d4 + real(dp) :: Ti1_scale = 1.0d4 + real(dp) :: Ti2_scale = 1.0d4 + real(dp) :: ni1_scale = 0.5d20 + real(dp) :: ni2_scale = 0.5d20 + + real(dp) :: Te_coef(MAX_POWER_SERIES) = [1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & + 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0] + real(dp) :: Ti1_coef(MAX_POWER_SERIES) = [1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & + 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0] + real(dp) :: Ti2_coef(MAX_POWER_SERIES) = [1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & + 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0] + real(dp) :: ni1_coef(MAX_POWER_SERIES) = [1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & + 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0] + real(dp) :: ni2_coef(MAX_POWER_SERIES) = [1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & + 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0] + + real(dp) :: Te_p1 = 1.0d0, Te_p2 = 0.0d0 + real(dp) :: Ti1_p1 = 1.0d0, Ti1_p2 = 0.0d0 + real(dp) :: Ti2_p1 = 1.0d0, Ti2_p2 = 0.0d0 + real(dp) :: ni1_p1 = 1.0d0, ni1_p2 = 0.0d0 + real(dp) :: ni2_p1 = 1.0d0, ni2_p2 = 0.0d0 + + namelist /profiles/ profile_type, & + Te_scale, Ti1_scale, Ti2_scale, ni1_scale, ni2_scale, & + Te_coef, Ti1_coef, Ti2_coef, ni1_coef, ni2_coef, & + Te_p1, Te_p2, Ti1_p1, Ti1_p2, Ti2_p1, Ti2_p2, & + ni1_p1, ni1_p2, ni2_p1, ni2_p2 + +contains + + subroutine read_profiles(config_file) + character(len=*), intent(in) :: config_file + integer :: ios, unit + + open (newunit=unit, file=config_file, status='old', action='read', iostat=ios) + if (ios /= 0) return + + read (unit, nml=profiles, iostat=ios) + close (unit) + + select case (trim(profile_type)) + case ("two_power") + active_profile = TWO_POWER + case default + active_profile = POWER_SERIES + end select + end subroutine read_profiles + + subroutine get_plasma_params(s, Te, Ti1, Ti2, ni1, ni2) + real(dp), intent(in) :: s + real(dp), intent(out) :: Te, Ti1, Ti2, ni1, ni2 + + if (active_profile == TWO_POWER) then + call get_two_power(s, Te, Ti1, Ti2, ni1, ni2) + else + call get_power_series(s, Te, Ti1, Ti2, ni1, ni2) + end if + end subroutine get_plasma_params + + pure function eval_power_series(s, scale, coef) result(f) + real(dp), intent(in) :: s + real(dp), intent(in) :: scale + real(dp), intent(in) :: coef(MAX_POWER_SERIES) + real(dp) :: f + real(dp) :: s_power + integer :: n + + f = 0.0d0 + s_power = 1.0d0 + do n = 1, MAX_POWER_SERIES + f = f + coef(n)*s_power + s_power = s_power*s + end do + f = scale*f + end function eval_power_series + + pure function eval_two_power(s, scale, p1, p2) result(f) + real(dp), intent(in) :: s + real(dp), intent(in) :: scale + real(dp), intent(in) :: p1 + real(dp), intent(in) :: p2 + real(dp) :: f + real(dp) :: base + + base = 1.0d0 - s**p1 + if (base < 0.0d0) then + f = 0.0d0 + else if (abs(base) < epsilon(1.0d0)) then + if (abs(p2) < epsilon(1.0d0)) then + f = scale + else + f = 0.0d0 + end if + else + f = scale*base**p2 + end if + end function eval_two_power + + subroutine get_power_series(s, Te, Ti1, Ti2, ni1, ni2) + real(dp), intent(in) :: s + real(dp), intent(out) :: Te, Ti1, Ti2, ni1, ni2 + + Te = eval_power_series(s, Te_scale, Te_coef) + Ti1 = eval_power_series(s, Ti1_scale, Ti1_coef) + Ti2 = eval_power_series(s, Ti2_scale, Ti2_coef) + ni1 = eval_power_series(s, ni1_scale, ni1_coef) + ni2 = eval_power_series(s, ni2_scale, ni2_coef) + end subroutine get_power_series + + subroutine get_two_power(s, Te, Ti1, Ti2, ni1, ni2) + real(dp), intent(in) :: s + real(dp), intent(out) :: Te, Ti1, Ti2, ni1, ni2 + + Te = eval_two_power(s, Te_scale, Te_p1, Te_p2) + Ti1 = eval_two_power(s, Ti1_scale, Ti1_p1, Ti1_p2) + Ti2 = eval_two_power(s, Ti2_scale, Ti2_p1, Ti2_p2) + ni1 = eval_two_power(s, ni1_scale, ni1_p1, ni1_p2) + ni2 = eval_two_power(s, ni2_scale, ni2_p1, ni2_p2) + end subroutine get_two_power + +end module simple_profiles diff --git a/src/simple_main.f90 b/src/simple_main.f90 index cd7d8842..58b1f7f4 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -4,7 +4,7 @@ module simple_main use util, only: sqrt2 use simple, only: init_vmec, init_sympl, tracer_t use diag_mod, only: icounter - use collis_alp, only: loacol_alpha, stost + use collis_alp, only: loacol_alpha, stost, init_collision_profiles use samplers, only: sample use field_can_mod, only: integ_to_ref, ref_to_integ, init_field_can use callback, only: output_orbits_macrostep @@ -71,6 +71,9 @@ subroutine main call read_config(config_file) call print_phase_time('Configuration reading completed') + call read_profiles_config(config_file) + call print_phase_time('Profiles configuration reading completed') + call init_field(norb, netcdffile, ns_s, ns_tp, multharm, integmode) call print_phase_time('Field initialization completed') @@ -407,19 +410,39 @@ subroutine print_parameters print *, 'v0 = ', v0 end subroutine print_parameters + subroutine read_profiles_config(config_file) + use simple_profiles, only: read_profiles + + character(256), intent(in) :: config_file + + call read_profiles(config_file) + end subroutine read_profiles_config + subroutine init_collisions - use params, only: am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & - facE_al, dchichi, slowrate, dchichi_norm, slowrate_norm, v0 + use params, only: am1, am2, Z1, Z2, facE_al, dchichi, slowrate, & + dchichi_norm, slowrate_norm, v0 + use simple_profiles, only: Te_scale, Ti1_scale, Ti2_scale, & + ni1_scale, ni2_scale - real(dp) :: v0_coll + real(dp) :: v0_coll, ealpha + real(dp) :: densi1, densi2, tempi1, tempi2, tempe + + ealpha = 3.5d6/facE_al + densi1 = ni1_scale*1.0d-6 + densi2 = ni2_scale*1.0d-6 + tempi1 = Ti1_scale + tempi2 = Ti2_scale + tempe = Te_scale call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & - 3.5d6/facE_al, v0_coll, dchichi, slowrate, dchichi_norm, & + ealpha, v0_coll, dchichi, slowrate, dchichi_norm, & slowrate_norm) if (abs(v0_coll - v0) > 1d-6) then error stop 'simple_main.init_collisions: v0_coll != v0' end if + + call init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) end subroutine init_collisions subroutine sample_particles @@ -803,7 +826,8 @@ subroutine update_momentum(anorb, z) end subroutine update_momentum subroutine collide(z, dt) - real(dp), intent(in) :: z(5), dt + real(dp), intent(inout) :: z(5) + real(dp), intent(in) :: dt integer :: ierr_coll call stost(z, dt, 1, ierr_coll) diff --git a/test/golden_record/collision/simple.in b/test/golden_record/collision/simple.in new file mode 100644 index 00000000..03e79a2a --- /dev/null +++ b/test/golden_record/collision/simple.in @@ -0,0 +1,11 @@ +&config +multharm = 3 +contr_pp = -1e10 +trace_time = 1d-4 +sbeg = 0.3d0 +ntestpart = 32 +netcdffile = 'wout.nc' +isw_field_type = 2 +deterministic = .True. +swcoll = .True. +/ diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index d249240c..9f439df6 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -303,6 +303,11 @@ target_link_libraries(test_orbits.x simple) add_executable (test_collis.x ../../src/test_collis.f90) target_link_libraries(test_collis.x simple) +add_executable (test_profiles.x test_profiles.f90) +target_link_libraries(test_profiles.x simple) +add_test(NAME test_profiles COMMAND test_profiles.x) +set_tests_properties(test_profiles PROPERTIES LABELS "unit") + add_executable (test_samplers.x test_samplers.f90) target_link_libraries(test_samplers.x simple) diff --git a/test/tests/test_profiles.f90 b/test/tests/test_profiles.f90 new file mode 100644 index 00000000..a7241dba --- /dev/null +++ b/test/tests/test_profiles.f90 @@ -0,0 +1,348 @@ +program test_profiles + use, intrinsic :: iso_fortran_env, only: dp => real64 + use simple_profiles + use collis_alp + implicit none + + logical :: all_passed + integer :: n_failed + + all_passed = .true. + n_failed = 0 + + call test_power_series_evaluation(all_passed, n_failed) + call test_two_power_evaluation(all_passed, n_failed) + call test_flat_profile_coefficients(all_passed, n_failed) + call test_peaked_profile_collision_rates(all_passed, n_failed) + call plot_collision_frequency_comparison() + + if (all_passed) then + print *, 'All profile tests passed' + else + print *, 'FAILED: ', n_failed, ' profile tests failed' + error stop 1 + end if + +contains + + subroutine test_power_series_evaluation(passed, nfail) + logical, intent(inout) :: passed + integer, intent(inout) :: nfail + + real(dp) :: coef(MAX_POWER_SERIES), result, expected + real(dp), parameter :: tol = 1.0d-12 + + print *, 'Testing power series evaluation...' + + coef = 0.0d0 + coef(1) = 1.0d0 + coef(2) = -1.0d0 + result = eval_power_series(0.5d0, 1.0d0, coef) + expected = 0.5d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: 1-s at s=0.5 expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: 1-s at s=0.5' + end if + + coef = 0.0d0 + coef(1) = 1.0d0 + coef(2) = -2.0d0 + coef(3) = 1.0d0 + result = eval_power_series(0.5d0, 1.0d0, coef) + expected = 0.25d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: (1-s)^2 at s=0.5 expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: (1-s)^2 at s=0.5' + end if + + result = eval_power_series(0.0d0, 100.0d0, coef) + expected = 100.0d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: scale at s=0 expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: scale at s=0' + end if + end subroutine test_power_series_evaluation + + subroutine test_two_power_evaluation(passed, nfail) + logical, intent(inout) :: passed + integer, intent(inout) :: nfail + + real(dp) :: result, expected + real(dp), parameter :: tol = 1.0d-12 + + print *, 'Testing two-power evaluation...' + + result = eval_two_power(0.5d0, 1.0d0, 1.0d0, 2.0d0) + expected = 0.25d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: (1-s)^2 at s=0.5 expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: (1-s)^2 at s=0.5' + end if + + result = eval_two_power(0.5d0, 1.0d0, 1.0d0, 3.5d0) + expected = 0.5d0**3.5d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: (1-s)^3.5 at s=0.5 expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: (1-s)^3.5 at s=0.5' + end if + + result = eval_two_power(0.0d0, 10010.0d0, 1.0d0, 2.0d0) + expected = 10010.0d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: on-axis value expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: on-axis value' + end if + + result = eval_two_power(1.0d0, 10010.0d0, 1.0d0, 2.0d0) + expected = 0.0d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: edge value expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: edge value is zero' + end if + + result = eval_two_power(1.0d0, 10010.0d0, 1.0d0, 0.0d0) + expected = 10010.0d0 + if (abs(result - expected) > tol) then + print *, ' FAIL: flat profile at edge expected', expected, 'got', result + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: flat profile at edge gives scale' + end if + end subroutine test_two_power_evaluation + + subroutine test_flat_profile_coefficients(passed, nfail) + logical, intent(inout) :: passed + integer, intent(inout) :: nfail + + real(dp) :: am1, am2, Z1, Z2, ealpha, v0 + real(dp) :: densi1, densi2, tempi1, tempi2, tempe + real(dp) :: dchichi, slowrate, dchichi_norm, slowrate_norm + real(dp) :: efcolf_scalar(nsorts), velrat_scalar(nsorts), enrat_scalar(nsorts) + real(dp) :: efcolf_profile(nsorts), velrat_profile(nsorts), enrat_profile(nsorts) + real(dp), parameter :: tol = 1.0d-10 + integer :: i + + print *, 'Testing flat profile vs scalar coefficients...' + + am1 = 2.0d0 + am2 = 3.0d0 + Z1 = 1.0d0 + Z2 = 1.0d0 + densi1 = 0.5d14 + densi2 = 0.5d14 + tempi1 = 1.0d4 + tempi2 = 1.0d4 + tempe = 1.0d4 + ealpha = 3.5d6 + + call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & + ealpha, v0, dchichi, slowrate, dchichi_norm, slowrate_norm) + efcolf_scalar = efcolf + velrat_scalar = velrat + enrat_scalar = enrat + + profile_type = "two_power" + active_profile = TWO_POWER + Te_scale = tempe + Te_p1 = 1.0d0 + Te_p2 = 0.0d0 + Ti1_scale = tempi1 + Ti1_p1 = 1.0d0 + Ti1_p2 = 0.0d0 + Ti2_scale = tempi2 + Ti2_p1 = 1.0d0 + Ti2_p2 = 0.0d0 + ni1_scale = densi1*1.0d6 + ni1_p1 = 1.0d0 + ni1_p2 = 0.0d0 + ni2_scale = densi2*1.0d6 + ni2_p1 = 1.0d0 + ni2_p2 = 0.0d0 + + call init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) + + call get_local_coeffs(0.5d0, efcolf_profile, velrat_profile, enrat_profile) + + do i = 1, nsorts + if (abs(efcolf_profile(i) - efcolf_scalar(i))/abs(efcolf_scalar(i)) > tol) then + print *, ' FAIL: efcolf mismatch for species', i + print *, ' scalar:', efcolf_scalar(i), 'profile:', efcolf_profile(i) + passed = .false. + nfail = nfail + 1 + end if + if (abs(velrat_profile(i) - velrat_scalar(i))/abs(velrat_scalar(i)) > tol) then + print *, ' FAIL: velrat mismatch for species', i + passed = .false. + nfail = nfail + 1 + end if + if (abs(enrat_profile(i) - enrat_scalar(i))/abs(enrat_scalar(i)) > tol) then + print *, ' FAIL: enrat mismatch for species', i + passed = .false. + nfail = nfail + 1 + end if + end do + + print *, ' PASS: flat profile coefficients match scalar' + end subroutine test_flat_profile_coefficients + + subroutine test_peaked_profile_collision_rates(passed, nfail) + logical, intent(inout) :: passed + integer, intent(inout) :: nfail + + real(dp) :: am1, am2, Z1, Z2, ealpha, v0 + real(dp) :: densi1, densi2, tempi1, tempi2, tempe + real(dp) :: dchichi, slowrate, dchichi_norm, slowrate_norm + real(dp) :: efcolf_core(nsorts), efcolf_edge(nsorts) + real(dp) :: velrat_core(nsorts), velrat_edge(nsorts) + real(dp) :: enrat_core(nsorts), enrat_edge(nsorts) + real(dp) :: s_core, s_edge + + print *, 'Testing peaked profile produces different collision rates...' + + am1 = 2.0d0 + am2 = 3.0d0 + Z1 = 1.0d0 + Z2 = 1.0d0 + densi1 = 0.5d14 + densi2 = 0.5d14 + tempi1 = 1.0d4 + tempi2 = 1.0d4 + tempe = 1.0d4 + ealpha = 3.5d6 + + call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & + ealpha, v0, dchichi, slowrate, dchichi_norm, slowrate_norm) + + profile_type = "two_power" + active_profile = TWO_POWER + Te_scale = 10010.0d0 + Te_p1 = 1.0d0 + Te_p2 = 2.0d0 + Ti1_scale = 10010.0d0 + Ti1_p1 = 1.0d0 + Ti1_p2 = 2.0d0 + Ti2_scale = 0.0d0 + Ti2_p1 = 1.0d0 + Ti2_p2 = 0.0d0 + ni1_scale = 1.822d21 + ni1_p1 = 1.0d0 + ni1_p2 = 3.5d0 + ni2_scale = 0.0d0 + ni2_p1 = 1.0d0 + ni2_p2 = 0.0d0 + + call init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) + + s_core = 0.04d0 + s_edge = 0.64d0 + + call get_local_coeffs(s_core, efcolf_core, velrat_core, enrat_core) + call get_local_coeffs(s_edge, efcolf_edge, velrat_edge, enrat_edge) + + if (efcolf_core(1) <= efcolf_edge(1)) then + print *, ' FAIL: core collision rate should be higher than edge' + print *, ' core efcolf(1):', efcolf_core(1) + print *, ' edge efcolf(1):', efcolf_edge(1) + passed = .false. + nfail = nfail + 1 + else + print *, ' PASS: core has higher collision rate than edge' + print *, ' core/edge ratio:', efcolf_core(1)/efcolf_edge(1) + end if + end subroutine test_peaked_profile_collision_rates + + subroutine plot_collision_frequency_comparison() + use fortplot, only: figure, plot, savefig, xlabel, ylabel, title, legend + + real(dp) :: am1, am2, Z1, Z2, ealpha, v0 + real(dp) :: densi1, densi2, tempi1, tempi2, tempe + real(dp) :: dchichi, slowrate, dchichi_norm, slowrate_norm + real(dp) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts) + real(dp) :: efcolf_onaxis + integer, parameter :: npts = 101 + real(dp) :: s_arr(npts), efcolf_profile(npts), efcolf_flat(npts) + real(dp) :: s, ds + integer :: i + + print *, 'Generating collision frequency comparison plot...' + + am1 = 2.0d0 + am2 = 3.0d0 + Z1 = 1.0d0 + Z2 = 1.0d0 + ealpha = 3.5d6 + + profile_type = "two_power" + active_profile = TWO_POWER + Te_scale = 10010.0d0 + Te_p1 = 1.0d0 + Te_p2 = 2.0d0 + Ti1_scale = 10010.0d0 + Ti1_p1 = 1.0d0 + Ti1_p2 = 2.0d0 + Ti2_scale = 0.0d0 + Ti2_p1 = 1.0d0 + Ti2_p2 = 0.0d0 + ni1_scale = 1.822d21 + ni1_p1 = 1.0d0 + ni1_p2 = 3.5d0 + ni2_scale = 0.0d0 + ni2_p1 = 1.0d0 + ni2_p2 = 0.0d0 + + densi1 = ni1_scale*1.0d-6 + densi2 = ni2_scale*1.0d-6 + tempi1 = Ti1_scale + tempi2 = Ti2_scale + tempe = Te_scale + + call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & + ealpha, v0, dchichi, slowrate, dchichi_norm, slowrate_norm) + efcolf_onaxis = efcolf(1) + + call init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) + + ds = 1.0d0/dble(npts - 1) + do i = 1, npts + s = dble(i - 1)*ds + s_arr(i) = s + call get_local_coeffs(s, efcolf_loc, velrat_loc, enrat_loc) + efcolf_profile(i) = efcolf_loc(1) + efcolf_flat(i) = efcolf_onaxis + end do + + call figure() + call plot(s_arr, efcolf_profile, label="Peaked profile", linestyle="b-") + call plot(s_arr, efcolf_flat, label="Uniform (on-axis)", linestyle="r--") + call xlabel("s (normalized toroidal flux)") + call ylabel("Collision frequency [1/s]") + call title("Collision Frequency: Profile vs Uniform") + call legend() + call savefig("collision_frequency_comparison.png") + + print *, ' Saved plot to collision_frequency_comparison.png' + end subroutine plot_collision_frequency_comparison + +end program test_profiles