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.

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

BEAMS3D: Added ability to specify weights in input namelist. #216

Merged
merged 2 commits into from
Jan 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
8 changes: 6 additions & 2 deletions BEAMS3D/Sources/beams3d_init.f90
Original file line number Diff line number Diff line change
Expand Up @@ -732,7 +732,8 @@ SUBROUTINE beams3d_init
vr_start = vr_start_in(1:nparticles)
vphi_start = vphi_start_in(1:nparticles)
vz_start = vz_start_in(1:nparticles)
weight = 1.0/nparticles
weight = weight_in(1:nparticles)
!weight = 1.0/nparticles
Zatom = Zatom_in(1:nparticles)
mass = mass_in(1:nparticles)
charge = charge_in(1:nparticles)
Expand All @@ -742,7 +743,10 @@ SUBROUTINE beams3d_init
nbeams = 1
charge_beams(1) = charge_in(1)
mass_beams(1) = mass_in(1)
lgc2fo_start(:) = .TRUE.
lgc2fo_start = .FALSE.
WHERE ((vr_start == 0) .and. (vphi_start == 0) .and. (vz_start == 0))
lgc2fo_start = .TRUE.
END WHERE
END IF

! Duplicate particles if requested
Expand Down
86 changes: 60 additions & 26 deletions BEAMS3D/Sources/beams3d_input_mod.f90
Original file line number Diff line number Diff line change
Expand Up @@ -77,11 +77,12 @@ MODULE beams3d_input_mod
r_start_in, phi_start_in, z_start_in, &
vll_start_in, npoinc, follow_tol, &
t_end_in, mu_start_in, charge_in, &
mass_in, Zatom_in, vc_adapt_tol, &
mass_in, Zatom_in, weight_in, &
vc_adapt_tol, &
int_type, Adist_beams, Asize_beams, &
Div_beams, E_beams, Dex_beams, &
mass_beams, charge_beams, Zatom_beams, &
r_beams, z_beams, phi_beams, s_max, TE_AUX_S, &
r_beams, z_beams, phi_beams, s_max, TE_AUX_S, &
TE_AUX_F, NE_AUX_S, NE_AUX_F, TI_AUX_S, &
TI_AUX_F, POT_AUX_S, POT_AUX_F, &
NI_AUX_S, NI_AUX_F, NI_AUX_Z, NI_AUX_M, &
Expand Down Expand Up @@ -134,6 +135,7 @@ SUBROUTINE init_beams3d_input
mass_in = -1.0
charge_in = -1.0
Zatom_in = -1.0
weight_in = 1.0

Adist_beams = 1.0_rprec
Asize_beams = -1.0_rprec
Expand Down Expand Up @@ -379,7 +381,7 @@ SUBROUTINE read_beams3d_input(filename, istat)
NI_AUX_Z(1) = 1
NI_AUX_M(1) = plasma_mass
END IF
s_max_zeff=ZEFF_AUX_S(nzeff)
s_max_zeff=ZEFF_AUX_S(nzeff)

nparticles = 0
DO WHILE ((r_start_in(nparticles+1) >= 0.0).and.(nparticles<MAXPARTICLES))
Expand Down Expand Up @@ -416,7 +418,7 @@ END SUBROUTINE read_beams3d_input
SUBROUTINE write_beams3d_namelist(iunit_out, istat)
INTEGER, INTENT(in) :: iunit_out
INTEGER, INTENT(out) :: istat
INTEGER :: ik, n
INTEGER :: ik, n, l
CHARACTER(LEN=*), PARAMETER :: outboo = "(2X,A,1X,'=',1X,L1)"
CHARACTER(LEN=*), PARAMETER :: outint = "(2X,A,1X,'=',1X,I0)"
CHARACTER(LEN=*), PARAMETER :: outflt = "(2X,A,1X,'=',1X,ES22.12E3)"
Expand All @@ -428,7 +430,7 @@ SUBROUTINE write_beams3d_namelist(iunit_out, istat)
CHARACTER(LEN=*), PARAMETER :: vecvar2 = "(2X,A,'(',I3.3,',',I3.3,')',1X,'=',1X,ES22.12E3)"
istat = 0
WRITE(iunit_out,'(A)') '&BEAMS3D_INPUT'
WRITE(iunit_out,'(A)') '!---------- General Parameters ------------'
WRITE(iunit_out,'(A)') '!---------- Background Grid Parameters ------------'
WRITE(iunit_out,outint) 'NR',nr
WRITE(iunit_out,outint) 'NZ',nz
WRITE(iunit_out,outint) 'NPHI',nphi
Expand All @@ -438,15 +440,15 @@ SUBROUTINE write_beams3d_namelist(iunit_out, istat)
WRITE(iunit_out,outflt) 'ZMAX',zmax
WRITE(iunit_out,outflt) 'PHIMIN',phimin
WRITE(iunit_out,outflt) 'PHIMAX',phimax
WRITE(iunit_out,outint) 'NPOINC',npoinc
WRITE(iunit_out,outflt) 'VC_ADAPT_TOL',vc_adapt_tol
WRITE(iunit_out,'(A)') '!---------- Marker Tracking Parameters ------------'
WRITE(iunit_out,outstr) 'INT_TYPE',TRIM(int_type)
WRITE(iunit_out,outflt) 'FOLLOW_TOL',follow_tol
WRITE(iunit_out,outflt) 'VC_ADAPT_TOL',vc_adapt_tol
WRITE(iunit_out,outint) 'NPOINC',npoinc
WRITE(iunit_out,outint) 'NPARTICLES_START',nparticles_start
WRITE(iunit_out,'(A)') '!---------- Plasma Parameters ------------'
WRITE(iunit_out,outflt) 'PLASMA_MASS',plasma_mass
WRITE(iunit_out,outflt) 'PLASMA_ZMEAN',plasma_zmean
WRITE(iunit_out,outflt) 'THERM_FACTOR',therm_factor
WRITE(iunit_out,outflt) 'LENDT_M',lendt_m
WRITE(iunit_out,outflt) 'RHO_FULLORBIT',rho_fullorbit
WRITE(iunit_out,outint) 'DUPLICATE_FACTOR',duplicate_factor
WRITE(iunit_out,'(A)') '!---------- Distribution Parameters ------------'
WRITE(iunit_out,outint) 'NRHO_DIST',ns_prof1
WRITE(iunit_out,outint) 'NTHETA_DIST',ns_prof2
Expand All @@ -461,22 +463,34 @@ SUBROUTINE write_beams3d_namelist(iunit_out, istat)
WRITE(iunit_out,outflt) 'B_KICK_MIN',B_kick_min
WRITE(iunit_out,outflt) 'B_KICK_MAX',B_kick_max
END IF
WRITE(iunit_out,'(A)') '!---------- Plasma Parameters ------------'
WRITE(iunit_out,outflt) 'PLASMA_MASS',plasma_mass
WRITE(iunit_out,outflt) 'PLASMA_ZMEAN',plasma_zmean
WRITE(iunit_out,outflt) 'THERM_FACTOR',therm_factor
WRITE(iunit_out,"(A)") '!---------- Profiles ------------'
WRITE(iunit_out,outflt) 'NE_SCALE',NE_SCALE
WRITE(iunit_out,outflt) 'TE_SCALE',TE_SCALE
WRITE(iunit_out,outflt) 'TI_SCALE',TI_SCALE
WRITE(iunit_out,outflt) 'ZEFF_SCALE',ZEFF_SCALE
WRITE(iunit_out,outflt) 'THERM_FACTOR',therm_factor
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'NE_AUX_S',(ne_aux_s(n), n=1,nne)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'NE_AUX_F',(ne_aux_f(n), n=1,nne)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TE_AUX_S',(te_aux_s(n), n=1,nte)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TE_AUX_F',(te_aux_f(n), n=1,nte)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TI_AUX_S',(ti_aux_s(n), n=1,nti)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TI_AUX_F',(ti_aux_f(n), n=1,nti)
ik = COUNT(NI_AUX_S .ge. 0)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'NI_AUX_S',(ni_aux_s(n), n=1,ik)
DO n = 1, NION
IF (ANY(NI_AUX_F(n,:)>0)) THEN
WRITE(iunit_out,"(2X,A,'(',I1.1,',:) =',4(1X,ES22.12E3))") 'TI_AUX_F',n,(ni_aux_f(n,l), l=1,ik)
END IF
END DO
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'ZEFF_AUX_S',(zeff_aux_s(n), n=1,nzeff)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'ZEFF_AUX_F',(zeff_aux_f(n), n=1,nzeff)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'POT_AUX_S',(zeff_aux_s(n), n=1,npot)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'POT_AUX_F',(zeff_aux_f(n), n=1,npot)
IF (lbeam) THEN
WRITE(iunit_out,"(A)") '!---------- Profiles ------------'
WRITE(iunit_out,outflt) 'NE_SCALE',NE_SCALE
WRITE(iunit_out,outflt) 'TE_SCALE',TE_SCALE
WRITE(iunit_out,outflt) 'TI_SCALE',TI_SCALE
WRITE(iunit_out,outflt) 'ZEFF_SCALE',ZEFF_SCALE
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'NE_AUX_S',(ne_aux_s(n), n=1,nne)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'NE_AUX_F',(ne_aux_f(n), n=1,nne)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TE_AUX_S',(te_aux_s(n), n=1,nte)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TE_AUX_F',(te_aux_f(n), n=1,nte)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TI_AUX_S',(ti_aux_s(n), n=1,nti)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'TI_AUX_F',(ti_aux_f(n), n=1,nti)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'ZEFF_AUX_S',(zeff_aux_s(n), n=1,nzeff)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'ZEFF_AUX_F',(zeff_aux_f(n), n=1,nzeff)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'POT_AUX_S',(zeff_aux_s(n), n=1,npot)
WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'POT_AUX_F',(zeff_aux_f(n), n=1,npot)
DO n = 1, nbeams
WRITE(iunit_out,"(A,I2.2)") '!---- BEAM #',n
IF (dex_beams(n)>0) &
Expand Down Expand Up @@ -508,6 +522,7 @@ SUBROUTINE write_beams3d_namelist(iunit_out, istat)
WRITE(iunit_out,"(2X,A,1X,'=',10(1X,ES22.12E3))") 'MASS_IN',(mass_in(ik), ik=1,n)
WRITE(iunit_out,"(2X,A,1X,'=',10(1X,ES22.12E3))") 'CHARGE_IN',(charge_in(ik), ik=1,n)
WRITE(iunit_out,"(2X,A,1X,'=',10(1X,ES22.12E3))") 'ZATOM_IN',(zatom_in(ik), ik=1,n)
IF (ANY(weight_in /= 1)) WRITE(iunit_out,"(2X,A,1X,'=',10(1X,ES22.12E3))") 'WEIGHT_IN',(weight_in(ik), ik=1,n)
n = COUNT(t_end_in > -1)
WRITE(iunit_out,"(2X,A,1X,'=',I0,'*',ES22.12E3)") 'T_END_IN',n,MAXVAL(t_end_in)
END IF
Expand Down Expand Up @@ -549,6 +564,9 @@ SUBROUTINE BCAST_BEAMS3D_INPUT(local_master,comm,istat)
CALL MPI_BCAST(phimax,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(vc_adapt_tol,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(plasma_mass,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(lendt_m,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(rho_fullorbit,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(duplicate_factor,1,MPI_INTEGER, local_master, comm,istat)

CALL MPI_BCAST(nte,1,MPI_INTEGER, local_master, comm,istat)
CALL MPI_BCAST(nne,1,MPI_INTEGER, local_master, comm,istat)
Expand All @@ -559,10 +577,20 @@ SUBROUTINE BCAST_BEAMS3D_INPUT(local_master,comm,istat)
CALL MPI_BCAST(NE_AUX_F,MAXPROFLEN,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(TI_AUX_S,MAXPROFLEN,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(TI_AUX_F,MAXPROFLEN,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(NI_AUX_S,MAXPROFLEN,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(NI_AUX_F,MAXPROFLEN*NION,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(te_scale,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(ne_scale,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(ti_scale,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(zeff_scale,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(fusion_scale,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(therm_factor,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max_ne,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max_te,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max_ti,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max_zeff,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(s_max_pot,1,MPI_REAL8, local_master, comm,istat)

CALL MPI_BCAST(t_end_in,MAXPARTICLES,MPI_REAL8, local_master, comm,istat)

Expand All @@ -587,10 +615,16 @@ SUBROUTINE BCAST_BEAMS3D_INPUT(local_master,comm,istat)
CALL MPI_BCAST(mass_in,nparticles,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(charge_in,nparticles,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(Zatom_in,nparticles,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(weight_in,nparticles,MPI_REAL8, local_master, comm,istat)
END IF

CALL MPI_BCAST(follow_tol,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(int_type, 256, MPI_CHARACTER, local_master, comm,istat)

CALL MPI_BCAST(E_kick,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(freq_kick,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(B_kick_min,1,MPI_REAL8, local_master, comm,istat)
CALL MPI_BCAST(B_kick_max,1,MPI_REAL8, local_master, comm,istat)
#endif
END SUBROUTINE BCAST_BEAMS3D_INPUT

Expand Down
5 changes: 3 additions & 2 deletions BEAMS3D/Sources/beams3d_runtime.f90
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
! - Box Modeling implemented
! - FIELDLINES Interface Added
! v4.05 08/25/23 - Fast Tritium only calculation added
! v4.07 01/11/24 - Added ability to specifiy weights in the input
!-----------------------------------------------------------------------
MODULE beams3d_runtime
!-----------------------------------------------------------------------
Expand Down Expand Up @@ -154,7 +155,7 @@ MODULE beams3d_runtime
REAL(rprec), DIMENSION(NION) :: NI_AUX_M
REAL(rprec), DIMENSION(MAXPARTICLES) :: r_start_in, phi_start_in, z_start_in, vll_start_in, &
& mu_start_in, charge_in, Zatom_in, mass_in, t_end_in, &
vr_start_in, vphi_start_in, vz_start_in
vr_start_in, vphi_start_in, vz_start_in, weight_in
LOGICAL, ALLOCATABLE :: lgc2fo_start(:)
REAL(rprec), ALLOCATABLE :: R_start(:), phi_start(:), Z_start(:), vll_start(:), mu_start(:), &
& mass(:), charge(:), Zatom(:), t_end(:), weight(:), vr_start(:), vphi_start(:), vz_start(:)
Expand All @@ -163,7 +164,7 @@ MODULE beams3d_runtime
CHARACTER(256) :: id_string, mgrid_string, coil_string, &
vessel_string, int_type, restart_string, bbnbi_string, eqdsk_string

REAL(rprec), PARAMETER :: BEAMS3D_VERSION = 4.05 ! this is the full orbit test version
REAL(rprec), PARAMETER :: BEAMS3D_VERSION = 4.07 ! this is the full orbit test version

!-----------------------------------------------------------------------
! Subroutines
Expand Down