Skip to content

Unify reals #33

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
krystophny merged 3 commits into from
Jan 12, 2026
Merged

Unify reals #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2532ff3
Make all reals real(dp) where dp => real64
zandivx Jan 9, 2026
78440c7
Reformat touched lines
zandivx Jan 10, 2026
2e121a6
Update test/test_reslines.f90
zandivx Jan 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions src/attenuation_factor.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,12 @@
module thetadata_mod
use iso_fortran_env, only: dp => real64

implicit none

logical :: attenuation_data_initialized = .false.
integer :: npoiarg
double precision :: argmin, argstepinv
double precision, dimension(:), allocatable :: argvals, thetavals
real(dp) :: argmin, argstepinv
real(dp), dimension(:), allocatable :: argvals, thetavals

! All variables are shared (read-only after init_attenuation_data)

Expand All @@ -15,7 +17,7 @@ subroutine init_attenuation_data()
! MUST be called BEFORE any parallel region
integer, parameter :: iunitarg = 741
integer :: i
double precision :: arglog, theta_val
real(dp) :: arglog, theta_val

if (attenuation_data_initialized) return

Expand All @@ -40,25 +42,25 @@ subroutine init_attenuation_data()
close (iunitarg)

argmin = argvals(1)
argstepinv = 1.d0/(argvals(2) - argvals(1))
argstepinv = 1.0_dp / (argvals(2) - argvals(1))

attenuation_data_initialized = .true.
end subroutine init_attenuation_data

end module thetadata_mod

subroutine attenuation_factor(D, Theta)

use iso_fortran_env, only: dp => real64
use thetadata_mod, only: npoiarg, argmin, argstepinv, argvals, thetavals
use polylag_3, only: mp, indef, plag1d

implicit none

double precision, intent(in) :: D
double precision, intent(out) :: Theta
double precision :: arglog, fun, der
real(dp), intent(in) :: D
real(dp), intent(out) :: Theta
real(dp) :: arglog, fun, der
integer, dimension(mp) :: indu
double precision, dimension(mp) :: xp, fp
real(dp), dimension(mp) :: xp, fp

arglog = log(D)

Expand Down
152 changes: 79 additions & 73 deletions src/collis_nbi.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,12 @@

module collis_alp
use iso_fortran_env, only: dp => real64

implicit none

integer, parameter :: nsorts = 3, ns = 10000
logical :: swcoll = .false.
double precision, dimension(nsorts) :: efcolf, velrat, enrat
real(dp), dimension(nsorts) :: efcolf, velrat, enrat

!$omp threadprivate (swcoll, efcolf, velrat, enrat)

Expand All @@ -29,62 +33,62 @@ subroutine coleff(p, dpp, dhh, fpeff)
implicit none

integer i
double precision, intent(in) :: p
double precision, intent(out) :: dpp, dhh, fpeff
double precision :: plim, xbeta, dp, dh, dpd
real(dp), intent(in) :: p
real(dp), intent(out) :: dpp, dhh, fpeff
real(dp) :: plim, xbeta, d_p, dh, dpd

plim = max(p, 1.d-8)
plim = max(p, 1.0e-8_dp)

dpp = 0.0d0
dhh = 0.0d0
fpeff = 0.0d0
dpp = 0.0_dp
dhh = 0.0_dp
fpeff = 0.0_dp

do i = 1, nsorts
xbeta = p*velrat(i)
xbeta = p * velrat(i)

call onseff(xbeta, dp, dh, dpd)
call onseff(xbeta, d_p, dh, dpd)

dpp = dpp + dp*efcolf(i)
dhh = dhh + dh*efcolf(i)
fpeff = fpeff + (dpd/plim - 2.0d0*dp*p*enrat(i))*efcolf(i)
dpp = dpp + d_p * efcolf(i)
dhh = dhh + dh * efcolf(i)
fpeff = fpeff + (dpd / plim - 2.0_dp * d_p * p * enrat(i)) * efcolf(i)
end do

dhh = dhh/plim**2
dhh = dhh / plim**2

return
end subroutine coleff

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)
subroutine onseff(v, d_p, dh, dpd)
! d_p - dimensionless dpp
! dh - dhh*p^2 (p - dmls)
! dpd - (1/p)(d/dp)p^2*d_p (p - dmls)

implicit none

! square root of pi
double precision, parameter :: sqp = 1.7724538d0
real(dp), parameter :: sqp = 1.7724538_dp
! cons=4./(3.*sqrt(pi))
double precision, parameter :: cons = .75225278d0
double precision, intent(in) :: v
double precision, intent(out) :: dp, dh, dpd
double precision :: v2, v3, ex, er
real(dp), parameter :: cons = 0.75225278_dp
real(dp), intent(in) :: v
real(dp), intent(out) :: d_p, dh, dpd
real(dp) :: 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
v3 = v2 * v
if (v < 0.01_dp) then
d_p = cons * (1.0_dp - 0.6_dp * v2)
dh = cons * (1.0_dp - 0.2_dp * v2)
dpd = 2.0_dp * cons * (1.0_dp - 1.2_dp * v2)
elseif (v > 6.0_dp) then
d_p = 1.0_dp / v3
dh = (1.0_dp - 0.5_dp / v2) / v
dpd = -1.0_dp / v3
else
ex = exp(-v2)/sqp
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
d_p = er / v3 - 2.0_dp * ex / v2
dh = er * (1.0_dp - 0.5_dp / v2) / v + ex / v2
dpd = 4.0_dp * ex - d_p
end if

return
Expand Down Expand Up @@ -124,43 +128,45 @@ subroutine loacol_nbi(amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2,

implicit none

double precision :: amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, ebeam, dense
double precision :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm, vti1, vti2, vte
double precision :: pi, pmass, emass, e, ev, alame, frecol_base, alami1, alami2

pi = 3.14159265358979d0
pmass = 1.6726d-24
emass = 9.1094d-28
e = 4.8032d-10
ev = 1.6022d-12

enrat(1) = ebeam/tempi1
enrat(2) = ebeam/tempi2
enrat(3) = ebeam/tempe

v0 = sqrt(2.d0*ebeam*ev/(amb*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)*Zb*Z1*(amb + am1)/(amb*tempi1 + am1*ebeam)))
alami2 = 23.d0 - log(max(epsilon(1.d0), &
sqrt(densi2*Z2**2/tempi2)*Zb*Z2*(amb + am2)/(amb*tempi2 + am2*ebeam)))
alame = 24.d0 - log(sqrt(dense)/tempe)
frecol_base = 2.d0*pi*dense*e**4*Zb**2/((amb*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
real(dp) :: amb, am1, am2, Zb, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, ebeam, dense
real(dp) :: v0, dchichi, slowrate, dchichi_norm, slowrate_norm, vti1, vti2, vte
real(dp) :: pi, pmass, emass, e, ev, alame, frecol_base, alami1, alami2

pi = 3.14159265358979_dp
pmass = 1.6726e-24_dp
emass = 9.1094e-28_dp
e = 4.8032e-10_dp
ev = 1.6022e-12_dp

enrat(1) = ebeam / tempi1
enrat(2) = ebeam / tempi2
enrat(3) = ebeam / tempe

v0 = sqrt(2.0_dp * ebeam * ev / (amb * pmass))
vti1 = sqrt(2.0_dp * tempi1 * ev / (pmass * am1))
vti2 = sqrt(2.0_dp * tempi2 * ev / (pmass * am2))
vte = sqrt(2.0_dp * tempe * ev / emass)

velrat(1) = v0 / vti1
velrat(2) = v0 / vti2
velrat(3) = v0 / vte

dense = densi1 * Z1 + densi2 * Z2
alami1 = 23.0_dp - log(max(epsilon(1.0_dp), &
sqrt(densi1 * Z1**2 / tempi1) * Zb * Z1 * (amb + am1) / (amb * &
tempi1 + am1 * ebeam)))
alami2 = 23.0_dp - log(max(epsilon(1.0_dp), &
sqrt(densi2 * Z2**2 / tempi2) * Zb * Z2 * (amb + am2) / (amb * &
tempi2 + am2 * ebeam)))
alame = 24.0_dp - log(sqrt(dense) / tempe)
frecol_base = 2.0_dp * pi * dense * e**4 * Zb**2 / ((amb * 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
end subroutine loacol_nbi

end module collis_alp
113 changes: 57 additions & 56 deletions src/datatypes.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,49 +1,50 @@
module neort_datatypes
use iso_fortran_env, only: dp => real64

implicit none

type :: magfie_params_t
real(8) :: s ! Normalised toroidal flux
real(8) :: R0 ! Major radius [cgs]
real(8) :: a ! Minor radius [cgs]
real(8) :: eps ! Inverse aspect ratio (a/R0)
real(8) :: psi_pr ! Radial derivative of poloidal flux
real(8) :: B0 ! Reference magnetic field magnitude [cgs]
real(8) :: Bthcov ! Covariant poloidal field component
real(8) :: Bphcov ! Covariant toroidal field component
real(8) :: dBthcovds ! Derivative of Bthcov with respect to s
real(8) :: dBphcovds ! Derivative of Bphcov with respect to s
real(8) :: q ! Safety factor
real(8) :: iota ! Rotational transform (1/q)
real(8) :: dVds ! Derivative of flux-surface volume with respect to s
real(8) :: M_t ! Toroidal Mach number
real(8) :: Om_tE ! Toroidal rotation frequency from electric field
real(8) :: Om_tBref ! Reference toroidal rotation frequency
real(8) :: vth ! Thermal velocity [cgs]
real(8) :: T_in_eV ! Temperature [eV]
real(8) :: m0 ! Poloidal mode number of perturbation
real(8) :: n0 ! Toroidal mode number of perturbation
real(8) :: Dp ! Thermodynamic force parameter
real(8) :: Drp_over_Dp ! Ratio of rotation force to thermodynamic force
real(8) :: etatp ! Normalized temperature gradient
real(8) :: etadt ! Normalized density gradient
real(dp) :: s ! Normalised toroidal flux
real(dp) :: R0 ! Major radius [cgs]
real(dp) :: a ! Minor radius [cgs]
real(dp) :: eps ! Inverse aspect ratio (a/R0)
real(dp) :: psi_pr ! Radial derivative of poloidal flux
real(dp) :: B0 ! Reference magnetic field magnitude [cgs]
real(dp) :: Bthcov ! Covariant poloidal field component
real(dp) :: Bphcov ! Covariant toroidal field component
real(dp) :: dBthcovds ! Derivative of Bthcov with respect to s
real(dp) :: dBphcovds ! Derivative of Bphcov with respect to s
real(dp) :: q ! Safety factor
real(dp) :: iota ! Rotational transform (1/q)
real(dp) :: dVds ! Derivative of flux-surface volume with respect to s
real(dp) :: M_t ! Toroidal Mach number
real(dp) :: Om_tE ! Toroidal rotation frequency from electric field
real(dp) :: Om_tBref ! Reference toroidal rotation frequency
real(dp) :: vth ! Thermal velocity [cgs]
real(dp) :: T_in_eV ! Temperature [eV]
real(dp) :: m0 ! Poloidal mode number of perturbation
real(dp) :: n0 ! Toroidal mode number of perturbation
real(dp) :: Dp ! Thermodynamic force parameter
real(dp) :: Drp_over_Dp ! Ratio of rotation force to thermodynamic force
real(dp) :: etatp ! Normalized temperature gradient
real(dp) :: etadt ! Normalized density gradient
logical :: pertfile ! Flag: use perturbation from file (not analytic)
logical :: nonlin ! Flag: enable nonlinear physics calculations
real(8) :: dpp ! Nonlinear diagnostic parameter
real(8) :: dhh ! Nonlinear diagnostic parameter
real(8) :: fpeff ! Effective passing fraction parameter
real(dp) :: dpp ! Nonlinear diagnostic parameter
real(dp) :: dhh ! Nonlinear diagnostic parameter
real(dp) :: fpeff ! Effective passing fraction parameter
end type magfie_params_t

type :: magfie_tensors_t
real(8) :: theta ! Boozer poloidal angle [rad]
real(8) :: bmod ! Magnetic field magnitude |B| [cgs]
real(8) :: sqrtg ! Jacobian determinant sqrt(g)
real(8) :: hder(3) ! Contravariant basis vectors
real(8) :: hcovar(3) ! Covariant basis vectors
real(8) :: hctrvr(3) ! Contravariant components of field-line curvature
real(8) :: hcurl(3) ! Components of curl of basis vectors
complex(8) :: bn ! Perturbation amplitude normalised by |B|
complex(8) :: eps_exp ! Reference analytic perturbation epsmn*exp(i*m0*theta)
real(dp) :: theta ! Boozer poloidal angle [rad]
real(dp) :: bmod ! Magnetic field magnitude |B| [cgs]
real(dp) :: sqrtg ! Jacobian determinant sqrt(g)
real(dp) :: hder(3) ! Contravariant basis vectors
real(dp) :: hcovar(3) ! Covariant basis vectors
real(dp) :: hctrvr(3) ! Contravariant components of field-line curvature
real(dp) :: hcurl(3) ! Components of curl of basis vectors
complex(dp) :: bn ! Perturbation amplitude normalised by |B|
complex(dp) :: eps_exp ! Reference analytic perturbation epsmn*exp(i*m0*theta)
end type magfie_tensors_t

type :: magfie_data_t
Expand All @@ -54,37 +55,37 @@ module neort_datatypes
end type magfie_data_t

type :: transport_summary_t
real(8) :: M_t ! Toroidal Mach number
real(8) :: Dco(2) ! Transport coefficients (D11, D12) for co-passing particles
real(8) :: Dctr(2) ! Transport coefficients (D11, D12) for counter-passing particles
real(8) :: Dt(2) ! Transport coefficients (D11, D12) for trapped particles
real(dp) :: M_t ! Toroidal Mach number
real(dp) :: Dco(2) ! Transport coefficients (D11, D12) for co-passing particles
real(dp) :: Dctr(2) ! Transport coefficients (D11, D12) for counter-passing particles
real(dp) :: Dt(2) ! Transport coefficients (D11, D12) for trapped particles
! D11: particle flux coefficient, D12: momentum flux coefficient
end type transport_summary_t

type :: torque_summary_t
logical :: has_torque ! Flag: torque calculation enabled
real(8) :: s ! Normalised toroidal flux
real(8) :: dVds ! Derivative of flux-surface volume with respect to s
real(8) :: M_t ! Toroidal Mach number
real(8) :: Tco ! Torque density from co-passing particles
real(8) :: Tctr ! Torque density from counter-passing particles
real(8) :: Tt ! Torque density from trapped particles
real(dp) :: s ! Normalised toroidal flux
real(dp) :: dVds ! Derivative of flux-surface volume with respect to s
real(dp) :: M_t ! Toroidal Mach number
real(dp) :: Tco ! Torque density from co-passing particles
real(dp) :: Tctr ! Torque density from counter-passing particles
real(dp) :: Tt ! Torque density from trapped particles
end type torque_summary_t

type :: transport_harmonic_t
integer :: mth ! Poloidal mode number associated with resonance
real(8) :: Dresco(2)
real(dp) :: Dresco(2)
! Transport coefficients (D11, D12) for co-passing particles, per harmonic
real(8) :: Dresctr(2)
real(dp) :: Dresctr(2)
! Transport coefficients (D11, D12) for counter-passing particles, per harmonic
real(8) :: Drest(2) ! Transport coefficients (D11, D12) for trapped particles, per harmonic
real(8) :: Tresco ! Torque from co-passing particles, per harmonic
real(8) :: Tresctr ! Torque from counter-passing particles, per harmonic
real(8) :: Trest ! Torque from trapped particles, per harmonic
real(8) :: vminp_over_vth ! Lower bound of passing velocity grid normalised
real(8) :: vmaxp_over_vth ! Upper bound of passing velocity grid normalised
real(8) :: vmint_over_vth ! Lower bound of trapped velocity grid normalised
real(8) :: vmaxt_over_vth ! Upper bound of trapped velocity grid normalised
real(dp) :: Drest(2) ! Transport coefficients (D11, D12) for trapped particles, per harmonic
real(dp) :: Tresco ! Torque from co-passing particles, per harmonic
real(dp) :: Tresctr ! Torque from counter-passing particles, per harmonic
real(dp) :: Trest ! Torque from trapped particles, per harmonic
real(dp) :: vminp_over_vth ! Lower bound of passing velocity grid normalised
real(dp) :: vmaxp_over_vth ! Upper bound of passing velocity grid normalised
real(dp) :: vmint_over_vth ! Lower bound of trapped velocity grid normalised
real(dp) :: vmaxt_over_vth ! Upper bound of trapped velocity grid normalised
end type transport_harmonic_t

type :: transport_data_t
Expand Down
Loading