Skip to content
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

Feature/collision refactor #326

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
Open

Feature/collision refactor #326

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
e1cafa0
LIBSTELL: Created a collision operators module for a various routines.
lazersos Dec 18, 2024
f74495d
LIBSTELL: Updated Objectlist and .dep for new collisions module.
lazersos Dec 18, 2024
d2af0b8
BEAMS3D: Migrated to the collision operator module.
lazersos Dec 18, 2024
b1b984e
BEAMS3D: Slight fix to collision operator.
lazersos Dec 19, 2024
416c90d
BEAMS3D: A bit more cleanup.
lazersos Dec 19, 2024
24ae34a
Merge commit '6f70cf2331125b5d9a33baf2af038492f692bc2e' into feature/…
kudav Jan 14, 2025
0acffbf
BEAMS3D: fixed thermalization
kudav Jan 14, 2025
21ab350
BEAMS3D: fixed fact_crit calculation
kudav Jan 14, 2025
d8162a8
BEAMS3D: removed equivalent calc. for fact_crit
kudav Jan 14, 2025
23eb49a
BEAMS3D: revert thermalization fix
kudav Jan 15, 2025
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
4 changes: 4 additions & 0 deletions BEAMS3D/BEAMS3D.dep
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,7 @@ beams3d_follow.o: \
../../LIBSTELL/$(LOCTYPE)/stel_kinds.o \
../../LIBSTELL/$(LOCTYPE)/safe_open_mod.o \
../../LIBSTELL/$(LOCTYPE)/wall_mod.o \
../../LIBSTELL/$(LOCTYPE)/collision_operators.o \
../../LIBSTELL/$(LOCTYPE)/mpi_inc.o \
beams3d_runtime.o \
beams3d_lines.o \
Expand All @@ -96,6 +97,7 @@ beams3d_follow.o: \
beams3d_follow_gc.o: \
../../LIBSTELL/$(LOCTYPE)/stel_kinds.o \
../../LIBSTELL/$(LOCTYPE)/safe_open_mod.o \
../../LIBSTELL/$(LOCTYPE)/collision_operators.o \
../../LIBSTELL/$(LOCTYPE)/mpi_inc.o \
beams3d_runtime.o \
beams3d_lines.o \
Expand All @@ -107,6 +109,7 @@ beams3d_follow_gc.o: \
beams3d_follow_fo.o: \
../../LIBSTELL/$(LOCTYPE)/stel_kinds.o \
../../LIBSTELL/$(LOCTYPE)/safe_open_mod.o \
../../LIBSTELL/$(LOCTYPE)/collision_operators.o \
../../LIBSTELL/$(LOCTYPE)/mpi_inc.o \
beams3d_runtime.o \
beams3d_lines.o \
Expand Down Expand Up @@ -276,6 +279,7 @@ beams3d_physics_mod.o: \
../../LIBSTELL/$(LOCTYPE)/mpi_params.o \
../../LIBSTELL/$(LOCTYPE)/ezspline.o \
../../LIBSTELL/$(LOCTYPE)/ezspline_obj.o \
../../LIBSTELL/$(LOCTYPE)/collision_operators.o \
beams3d_runtime.o \
beams3d_lines.o \
beams3d_grid.o
Expand Down
13 changes: 4 additions & 9 deletions BEAMS3D/Sources/beams3d_follow.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ SUBROUTINE beams3d_follow
USE beams3d_write_par
USE safe_open_mod, ONLY: safe_open
USE wall_mod, ONLY: wall_free, ihit_array, nface
USE collision_operators, ONLY: SET_CRIT_FACTOR, SET_COULOMB_FACTOR
USE mpi_inc
!-----------------------------------------------------------------------
! Local Variables
Expand Down Expand Up @@ -125,11 +126,7 @@ SUBROUTINE beams3d_follow

! Some helpers
fact_vsound = 1.5*sqrt(e_charge/plasma_mass)*therm_factor
fact_crit= SQRT(2.0*e_charge/electron_mass)*(0.75*sqrt_pi*electron_mass*plasma_Zmean/plasma_mass)**(1.0/3.0) ! Wesson pg 226 5.4.9
fact_crit_legacy = SQRT(2*e_charge/plasma_mass)*(0.75*sqrt_pi*sqrt(plasma_mass/electron_mass))**(1.0/3.0)
!fact_crit=fact_crit_pro*(plasma_Zmean/plasma_mass)**(1.0/3.0)
!fact_kick = pi2*2*SQRT(pi*1E-7*plasma_mass)*E_kick*freq_kick !old
!fact_kick = 2*freq_kick*E_kick
CALL SET_CRIT_FACTOR(plasma_Zmean,plasma_mass)

! Handle the Beam defaults
IF (lbeam) THEN
Expand Down Expand Up @@ -158,8 +155,7 @@ SUBROUTINE beams3d_follow
mybeam = Beam(i)
moment = mu_start(i)
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Save the IC of the neutral
my_end = t_end(i)
myline = i
Expand Down Expand Up @@ -193,8 +189,7 @@ SUBROUTINE beams3d_follow
myline = i
mytdex = 1
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Define neutral trajectory
myv_neut(1) = vr_start(i)*cos(phi_start(i)) - vphi_start(i)*sin(phi_start(i))
myv_neut(2) = vr_start(i)*sin(phi_start(i)) + vphi_start(i)*cos(phi_start(i))
Expand Down
10 changes: 4 additions & 6 deletions BEAMS3D/Sources/beams3d_follow_fo.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ SUBROUTINE beams3d_follow_fo
USE beams3d_write_par
USE beams3d_physics_mod, ONLY: beams3d_gc2fo, beams3d_calc_dt
USE safe_open_mod, ONLY: safe_open
USE collision_operators, ONLY: SET_COULOMB_FACTOR
USE mpi_inc
!-----------------------------------------------------------------------
! Local Variables
Expand Down Expand Up @@ -149,8 +150,7 @@ SUBROUTINE beams3d_follow_fo
lneut = .false.
! Collision parameters
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
DO ! Must do it this way becasue lbeam changes q(4) values
#if defined(NAG)
CALL D02CJF(t_nag,tf_nag,neqs_nag,q,fpart_eom,tol_nag,relab,out_beams3d_part,D02CJW,w,ier)
Expand Down Expand Up @@ -195,8 +195,7 @@ SUBROUTINE beams3d_follow_fo
lneut = .false.
! Collision parameters
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Setup DRKHVG parameters
iopt = 0
DO
Expand Down Expand Up @@ -252,8 +251,7 @@ SUBROUTINE beams3d_follow_fo
lneut = .false.
! Collision parameters
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Now handle Coordinate conversion
IF (lbeam .and. mytdex == 3) mytdex = 2 ! BEAM -> FO Run
IF (lboxsim) mytdex = 1
Expand Down
10 changes: 4 additions & 6 deletions BEAMS3D/Sources/beams3d_follow_gc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ SUBROUTINE beams3d_follow_gc
USE mpi_params ! MPI
USE beams3d_write_par
USE beams3d_physics_mod, ONLY: beams3d_calc_dt
USE collision_operators, ONLY: SET_COULOMB_FACTOR
USE safe_open_mod, ONLY: safe_open
USE mpi_inc
!-----------------------------------------------------------------------
Expand Down Expand Up @@ -142,8 +143,7 @@ SUBROUTINE beams3d_follow_gc
! Collision parameters
fact_kick = 2*E_kick*mycharge/(mymass*pi2*pi2*freq_kick*freq_kick*SQRT(pi*1E-7*plasma_mass))
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Now calc dt
CALL beams3d_calc_dt(1,q(1),q(2),q(3),dt)
tf_nag = t_nag+dt
Expand Down Expand Up @@ -192,8 +192,7 @@ SUBROUTINE beams3d_follow_gc
! Collision parameters
fact_kick = 2*E_kick*mycharge/(mymass*pi2*pi2*freq_kick*freq_kick*SQRT(pi*1E-7*plasma_mass))
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Now calc dt
CALL beams3d_calc_dt(1,q(1),q(2),q(3),dt)
tf_nag = t_nag+dt
Expand Down Expand Up @@ -265,8 +264,7 @@ SUBROUTINE beams3d_follow_gc
! Collision parameters
fact_kick = 2*E_kick*mycharge/(mymass*pi2*pi2*freq_kick*freq_kick*SQRT(pi*1E-7*plasma_mass))
fact_pa = plasma_mass/(mymass*plasma_Zmean)
fact_coul = myZ*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)

CALL SET_COULOMB_FACTOR(mymass,myZ,plasma_mass)
! Now calc dt
CALL beams3d_calc_dt(1,q(1),q(2),q(3),dt)
tf_nag = t_nag+dt
Expand Down
58 changes: 33 additions & 25 deletions BEAMS3D/Sources/beams3d_physics_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@ MODULE beams3d_physics_mod
USE EZspline_obj
USE EZspline
USE adas_mod_parallel
USE collision_operators, ONLY: COULOMB_LOG_NRL_COUNTERSTREAM,&
V_CRITICAL, V_CRITICAL_WEILAND, &
TAU_SPITZER
USE mpi_params

!-----------------------------------------------------------------
Expand Down Expand Up @@ -91,12 +94,10 @@ FUNCTION coll_op_nrl19(ne_in,te_in,vbeta_in,Zeff_in) result(slow_par)
DOUBLE PRECISION :: slow_par(3)
DOUBLE PRECISION, INTENT(in) :: ne_in, te_in, vbeta_in, Zeff_in
DOUBLE PRECISION :: ne_cm,coulomb_log
ne_cm = ne_in * 1E-6
coulomb_log = 43 - log(Zeff_in*fact_coul*sqrt(ne_cm/te_in)/(vbeta_in*vbeta_in))
!fact_crit_legacy = SQRT(2*e_charge/plasma_mass)*(0.75*sqrt_pi*sqrt(plasma_mass/electron_mass))**(1.0/3.0)
slow_par(1) = fact_crit_legacy*SQRT(te_in)
slow_par(2) = 3.777183D41*mymass*SQRT(te_in*te_in*te_in)/(ne_in*myZ*myZ*coulomb_log) ! note ne should be in m^-3 here, tau_spit
slow_par(3) =Zeff_in*fact_pa
coulomb_log = COULOMB_LOG_NRL_COUNTERSTREAM(ne_in,te_in,vbeta_in,Zeff_in)
slow_par(1) = V_CRITICAL(te_in)
slow_par(2) = TAU_SPITZER(mymass,ne_in,te_in,myZ,coulomb_log)
slow_par(3) = Zeff_in*fact_pa
END FUNCTION coll_op_nrl19

!-----------------------------------------------------------------
Expand All @@ -123,14 +124,11 @@ FUNCTION coll_op_nrl19_ie(ne_in,te_in,vbeta_in,Zeff_in) result(slow_par)
DOUBLE PRECISION :: slow_par(3)
DOUBLE PRECISION, INTENT(in) :: ne_in, te_in, vbeta_in, Zeff_in
DOUBLE PRECISION :: ne_cm, coulomb_loge, coulomb_logi
ne_cm = ne_in * 1E-6
coulomb_logi = 43 - log(Zeff_in*fact_coul*sqrt(ne_cm/te_in)/(vbeta_in*vbeta_in))
coulomb_loge=log(1.09d11 * te_in/Zeff_in/sqrt(ne_cm))
!WRITE(6,*) coulomb_loge, coulomb_logi
! Callen Ch2 pg41 eq2.135 (fact*Vtherm; Vtherm = SQRT(2*E/mass) so E in J not eV)
slow_par(1) = fact_crit*SQRT(te_in)*(coulomb_logi/coulomb_loge)**(1.0/3.0) !vcrit, the coulomb ratio is from Weiland (2018) eq.11
slow_par(2) = 3.777183D41*mymass*SQRT(te_in*te_in*te_in)/(ne_in*myZ*myZ*coulomb_loge) ! note ne should be in m^-3 here, tau_spit
slow_par(3) =Zeff_in*fact_pa
coulomb_logi = COULOMB_LOG_NRL_COUNTERSTREAM(ne_in,te_in,vbeta_in,Zeff_in)
coulomb_loge = log(1.09d11 * te_in/Zeff_in/sqrt(ne_cm))
slow_par(1) = V_CRITICAL_WEILAND(te_in,coulomb_logi,coulomb_loge)
slow_par(2) = TAU_SPITZER(mymass,ne_in,te_in,myZ,coulomb_loge)
slow_par(3) = Zeff_in*fact_pa
RETURN
END FUNCTION coll_op_nrl19_ie

Expand Down Expand Up @@ -242,7 +240,7 @@ SUBROUTINE beams3d_physics_gc(t, q)
rho_temp, omeg_temp, binv, vrot_para, vrot_perp, &
vc3_tauinv, vbeta, zeff_temp,&
sm,omega2,vrel2,bmax,bmincl,bminqu,bmin,&
zdelth,zrang
zdelth,zrang, factor_ion, factor_electron, speed3inv, speed2inv
DOUBLE PRECISION :: Ebench ! for ASCOT Benchmark
DOUBLE PRECISION :: slow_par(3), ni_temp(NION)
! For splines
Expand Down Expand Up @@ -327,6 +325,7 @@ SUBROUTINE beams3d_physics_gc(t, q)
! Apply toroidal rotation
! vrot_para: Parallel rotation velocity [m/s]
! vrot_perp: Perpendicular rotation velocity [m/s]
! speed Total particle speed [m/s]
!-----------------------------------------------------------
inv_mymass = one/mymass
vrot_para = omeg_temp*r_temp*bphi_temp*binv
Expand Down Expand Up @@ -358,37 +357,46 @@ SUBROUTINE beams3d_physics_gc(t, q)
tau_spit_inv = one/slow_par(2)
vc3_tauinv = vcrit_cube*tau_spit_inv
ELSE !Dont evaluate collisions
q(4) = vll + vrot_para
!q(4) = vll + vrot_para
RETURN
END IF

!------------------------------------------------------------
! Velocity diffusion
! https://doi.org/10.1103/PhysRev.107.1
! https://doi.org/10.1063/1.1694943
! https://doi.org/10.1016/0021-9991(81)90111-X
! https://doi.org/10.1016/j.cpc.2014.01.014
! speed3inv v**-3 [m^3/s^3]
! speed2inv v**-2 [m^3/s^3]
! ddve
!------------------------------------------------------------
ddve = zero; ddvi = zero; sigma = zero; zeta = zero
ddve = zero; ddvi = zero; sigma = zero; zeta = zero;
factor_ion = one; factor_electron = one
#if defined(B3D_VEL_DIFFUSION)
speed_cube = (speed*speed*speed)
CALL gauss_rand(1,zeta) ! A random from a standard normal (1,1)
speed3inv = 1.0/(speed*speed*speed)
speed2inv = speed*speed3inv
ddve=ABS(2*e_charge*dt*te_temp*inv_mymass*tau_spit_inv)
ddvi=ABS(2*e_charge*dt*(ti_temp*vcrit_cube*inv_mymass/speed_cube)*tau_spit_inv)
ddvi=ABS(2*e_charge*dt*ti_temp*vcrit_cube*inv_mymass*speed3inv*tau_spit_inv)
CALL gauss_rand(1,zeta) ! A random from a standard normal (1,1)
sigma = sqrt( ddve+ddvi) ! The standard deviation.
ddve=zeta*ddve/sigma
ddvi=zeta*ddvi/sigma
factor_electron = ( 1.0 - 2.0*te_temp*inv_mymass*e_charge*speed2inv)
factor_ion = ( 1.0 + ti_temp*inv_mymass*e_charge*speed2inv)
#endif
!-----------------------------------------------------------
! Viscouse Velocity Reduction
! v_s Local Sound Speed
! speed Total particle speed
! dve Speed change due to electron slowing down
! dvi Speed change due to ion slowing down
! reduction Total change in speed
! newspeed New total speed
! vfrac Ratio between new and old speed (helper)
!-----------------------------------------------------------
dve = speed*tau_spit_inv*(1-2*te_temp*inv_mymass*e_charge/speed**2.0)
dvi = vc3_tauinv/(speed*speed)*(1+ti_temp*inv_mymass*e_charge/speed**2.0)
dve = factor_electron*speed*tau_spit_inv
dvi = factor_ion*vc3_tauinv/(speed*speed)
reduction = (dve + dvi)*dt+sigma*zeta
newspeed = speed - reduction
newspeed = speed - reduction
dve=dve+ddve
dvi=dvi+ddvi
vfrac = newspeed/speed
Expand Down
3 changes: 2 additions & 1 deletion LIBSTELL/LIBSTELL.dep
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
# Microsoft Developer Studio Generated Dependency File, included by LIBSTELL.mak

collision_operators.o :

read_nescoil_mod.o : \
safe_open_mod.o \
stel_kinds.o \
Expand All @@ -8,7 +10,6 @@ read_nescoil_mod.o : \
mpi_sharmem.o \
mpi_inc.o


stellopt_globals.o : \
stel_kinds.o

Expand Down
1 change: 1 addition & 0 deletions LIBSTELL/ObjectList
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
ObjectFiles = \
collision_operators.o \
thrift_input_mod.o \
thrift_globals.o \
read_nescoil_mod.o \
Expand Down
Loading
Loading