AtHeartEngineer/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-10 07:38:00 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'kdev' into quickersquish

Browse Source
This commit is contained in:
Jeffrey Garretson
2020-10-10 19:13:24 -06:00
parent 3ae67fca3d 4472fe282e
commit d525907337
6 changed files with 75 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/base/kdefs.F90
View File

@@ -6,6 +6,7 @@ module kdefs
!Define variable precisions
integer, parameter :: sp = REAL32
integer, parameter :: dp = REAL64
integer, parameter :: qp = REAL128 !Quad precision
!integer, parameter :: ip = INT64
integer, parameter :: ip = INT32

50
src/rcm/rcm_subs.F90
View File

@@ -2270,7 +2270,7 @@ real :: v_1_1, v_1_2, v_2_1, v_2_2
V = sqrt(2*Kj/M)*100.0 !m/s->cm/s
!Timescale
tScl = cos(30.0*PI/180.0)**3.5
tScl = cos(45*PI/180.0)**3.5 !Using Smith & Bewtra 1976 scaling
Tau = tScl*1.0/(Ngeo*V*Sig)
cxrate = 1.0/Tau
@@ -2484,7 +2484,7 @@ FUNCTION FLCRat(ie,alam,vm,beq,rcurv,lossc) result(lossFLC)
real(rprec), intent(in) :: alam,vm,beq,rcurv,lossc
real(rprec) :: lossFLC
real(rprec) :: bfp,ftv,K,V,TauSS,Rgyro,eps,xSS,TauFLC
real(rprec) :: bfp,ftv,K,V,TauSS,Rgyro,eps,xSS,TauFLC,earg
bfp = beq/(sin(lossc)**2.0) !Foot point field strength, nT
ftv = (1.0/vm)**(3.0/2.0) !flux-tube volume Re/nT
@@ -2508,17 +2508,16 @@ FUNCTION FLCRat(ie,alam,vm,beq,rcurv,lossc) result(lossFLC)
eps = Rgyro/rcurv
!Chen+ 2019
!xSS = max(100.0*eps**-5.0,1.0)
!K: Mockup between Chen/Gibson, transition between eps^-5 dep. and strong scattering at kappa = sqrt(8)
xSS = max( (8.0*eps)**(-5.0), 1.0 )
!xSS = max( (8.0*eps)**(-5.0), 1.0 )
earg = eps**(-5.0)
xSS = max(100.0*earg,1.0)
!xSS = max(10.0*earg,1.0)
TauFLC = xSS*TauSS
lossFLC = 1.0/TauFLC !Rate, 1/s
! !Mixed Chen+ 2019/Gilson+ 2012
! !Above eps = 0.125, kap <= sqrt(8) use full strong-scattering
! lossFLC = (1.0/TauSS)*RampUp(eps,0.1_rp,0.15_rp)
END FUNCTION FLCRat
!=========================================================================
@@ -2561,6 +2560,9 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
!---
!Do prep work
where (eeta<0)
eeta = 0.0
endwhere
joff=jwrap-1
@@ -2580,7 +2582,6 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
fac = 1.0E-3*signbe*bir*alpha*beta*dlam*dpsi*ri**2
!---
!Get OCB
isOpen = (vm < 0)
@@ -2619,7 +2620,11 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
ie = 2 ! protons
END IF
IF (MAXVAL(eeta(:,:,kc)) < machine_tiny) CYCLE !Skip boring channels
IF (MAXVAL(eeta(:,:,kc)) < machine_tiny) then
!Skip boring channels
eeta(:,:,kc) = 0.0
CYCLE
END IF
mass_factor = SQRT (xmass(1)/xmass(ie))
!---
@@ -2663,6 +2668,10 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
!Freeze flow too close to MHD boundary
didt(iMHD-1:,:) = 0.0
djdt(iMHD+1:,:) = 0.0
!Freeze flow into the domain, only move stuff around from MHD buffer
didt(1:2,:) = 0.0
djdt(1 ,:) = 0.0
call PadClaw(didt)
call PadClaw(djdt)
@@ -2683,8 +2692,9 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
if ( L_dktime .and. (.not. isOpen(i,j)) ) then
!Do losses even in buffer region in case stuff moves in/out
r_dist = sqrt(xmin(i,j)**2+ymin(i,j)**2)
lossCX = Cexrat(ie,abs(alamc(kc))*vm(i,j),r_dist,sunspot_number, &
dktime,irdk,inrgdk,isodk,iondk)
!lossCX = Cexrat(ie,abs(alamc(kc))*vm(i,j),r_dist,sunspot_number, &
! dktime,irdk,inrgdk,isodk,iondk)
lossCX = CXKaiju(ie,abs(alamc(kc))*vm(i,j),r_dist)
!Placeholder for FLC loss, uses radcurv(i,j) [Re]
lossFLC = FLCRat(ie,alamc(kc),vm(i,j),bmin(i,j),radcurv(i,j),losscone(i,j))
endif
@@ -2703,7 +2713,7 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
!---
!Advect w/ clawpack
sumEtaBEF = sum(eeta(:,:,kc)) !Total content before clawpack
sumEtaBEF = sum(eeta(:,j1:j2,kc)) !Total content before clawpack
call rcm2claw(eeta(:,:,kc),etaC)
@@ -2730,11 +2740,11 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
call circle(eeta(:,:,kc))
!Check total content after versus before
sumEtaAFT = sum(eeta(:,:,kc))
if (sumEtaAFT>sumEtaBEF) then
!Can only increase content due to numerical shennanigans, i.e. borrowing from vacuum
eeta(:,:,kc) = (sumEtaBEF/sumEtaAFT)*eeta(:,:,kc)
endif
sumEtaAFT = sum(eeta(:,j1:j2,kc))
! if (sumEtaAFT>sumEtaBEF) then
! !Can only increase content due to numerical shennanigans, i.e. borrowing from vacuum
! eeta(:,:,kc) = (sumEtaBEF/sumEtaAFT)*eeta(:,:,kc)
! endif
if (doOCBNuke) then
!Go through and nuke any content next to open cell
@@ -2761,7 +2771,7 @@ SUBROUTINE Move_plasma_grid_MHD (dt)
!K: Added kc==1 check 8/11/20
call Kaiju_Plasmasphere_Refill(eeta(:,:,1), rmin, aloct, vm, dt)
call circle(eeta(:,:,kc)) !Probably don't need to re-circle
endif
endif
enddo !Main kc loop
@@ -3192,7 +3202,7 @@ SUBROUTINE Kaiju_Plasmasphere_Refill(eeta0,rmin,aloct,vm,idt)
dndt = 10.0**(3.48-0.331*rmin(i,j)) !cm^-3/day, Denton+ 2012 eqn 1
tau = day2s*(dppT-dpsph)/dndt
eeta0(i,j) = eeta0(i,j) + min(idt/tau,1.0)*deta !Make sure not to overfill but unlikely
eeta0(i,j) = eeta0(i,j) + min(idt/tau,1.0)*max(deta,0.0) !Make sure not to overfill but unlikely
enddo
enddo

5
src/rcm/tomhd.F90
View File

@@ -45,7 +45,7 @@
real(rprec) :: dens_plasmasphere
real(rp) :: sclmass(RCMNUMFLAV) !xmass prescaled to proton
real(rp) :: kevRCM,kevMHD
real(rp), parameter :: kRatMax = 0.5
real(rp), parameter :: kRatMax = 0.9
!AMS 04-22-2020
real(rprec) :: pressure_factor,density_factor
@@ -63,7 +63,6 @@
klow = 1
endif
!Set scaled mass by hand here to avoid precision issues
sclmass(RCMELECTRON) = mass_electron/mass_proton
sclmass(RCMPROTON) = 1.0
@@ -105,7 +104,7 @@
kevRCM = abs(alamc(kcsize))*vm(i,j)*1.0e-3 !Max keV of RCM channels here
kevMHD = DP2kT(densrcm(i,j)*rcmNScl,pressrcm(i,j)*rcmPScl) !Get keV from RCM moments
if (kevMHD > kRatMax*kevRCM) then
if (kevMHD >= kRatMax*kevRCM) then
!Effective "MHD" temperature, P=nkT_{MHD} is above kRatMax the max RCM channel energy
!This is probably bad for resolving the distribution so we do some shady cooling here

41
src/rcm/torcm.F90
View File

@@ -5,7 +5,7 @@ MODULE torcm_mod
USE rice_housekeeping_module, ONLY: use_plasmasphere,LowLatMHD,L_write_vars_debug
Use rcm_mhd_interfaces
USE rcmdefs, ONLY : RCMTOPCLOSED,RCMTOPNULL,RCMTOPOPEN,DenPP0
USE kdefs, ONLY : TINY
USE kdefs, ONLY : TINY,qp
use math, ONLY : RampDown
implicit none
@@ -211,6 +211,9 @@ MODULE torcm_mod
!-----
!Fill in grid ghosts
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j,dpp)
do j=1,jsize
do i=1,isize
if (iopen(i,j) == RCMTOPOPEN) then
@@ -340,6 +343,8 @@ MODULE torcm_mod
! compute vm and find boundaries
! (K: doing in two stages to avoid open/closed/open corner case)
vm(:,:) = big_vm
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j)
do j=jwrap,jsize
do i=isize,1,-1
if (iopen(i,j) == RCMTOPCLOSED) then
@@ -536,13 +541,13 @@ MODULE torcm_mod
real(rprec) :: dmhd,dpp,pcon,t,prcmI,prcmE,pmhdI,pmhdE,psclI,psclE
integer(iprec) :: i,j,k,kmin
real(rprec) :: xp,xm,A0,delerf,delexp
real(rprec) :: xp,xm,A0
logical :: isIon
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j,k,kmin,dmhd,dpp,pcon,prcmI,prcmE,t,pmhdI,pmhdE,psclI,psclE) &
!$OMP private(xp,xm,A0,delerf,delexp,isIon)
!$OMP private(xp,xm,A0,isIon)
do j=1,jsize
do i=1,isize
!Get corrected density from MHD
@@ -573,10 +578,9 @@ MODULE torcm_mod
xp = SQRT(ev*ABS(almmax(k))*vm(i,j)/boltz/t)
xm = SQRT(ev*ABS(almmin(k))*vm(i,j)/boltz/t)
delerf = erf(xp)-erf(xm)
delexp = (2.0/sqrt(pi)) * ( xp*exp(-xp**2.0) - xm*exp(-xm**2.0) )
eeta_new(i,j,k) = A0*( delerf - delexp )
!Use quad prec calc of erf/exp differences
eeta_new(i,j,k) = erfexpdiff(A0,xp,xm)
!Pressure contribution from this channel
pcon = pressure_factor*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
@@ -622,7 +626,26 @@ MODULE torcm_mod
END SUBROUTINE Press2eta
!Calculates difference of erfs - diff of exps, i.e. Eqn B5 from Pembroke+ 2012
function erfexpdiff(A,x,y) result(z)
real(rp), intent(in) :: A,x, y
real(rp) :: z
!QUAD precision holders
real(qp) :: xq,yq,zq,differf,diffexp
xq = x
yq = y
!Replacing erf(x)-erf(y) w/ erfc to avoid flooring to zero
differf = erfc(yq)-erfc(xq)
diffexp = xq*exp(-xq**2.0) - yq*exp(-yq**2.0)
diffexp = 2.0*diffexp/sqrt(PI)
zq = A*(differf-diffexp)
z = zq
end function erfexpdiff
!===================================================================
SUBROUTINE Set_ellipse(idim,jdim,rmin,pmin,vm,big_vm,bndloc,iopen)
@@ -693,6 +716,8 @@ MODULE torcm_mod
a = (a1 - a2)/2.
!Fill isGood array
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,ell,jm,jp)
do j=1,jdim
do i=1,idim
if (iopen(i,j) == RCMTOPOPEN) then
@@ -719,6 +744,8 @@ MODULE torcm_mod
enddo !j loop
!Now set boundary based on isGood and number of required buffer cells
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j)
do j=1,jdim
do i=idim,1,-1
!Loop from the top and find first bad cell

10
src/voltron/imagtubes.F90
View File

@@ -13,7 +13,6 @@ module imagtubes
real(rp) :: Rp_m
real(rp) :: planetM0g
logical, parameter :: doNewBeta = .false.
!Some threshold values for poisoning tubes
!TODO: Make these XML parameters
@@ -152,10 +151,10 @@ module imagtubes
ijTube%wIMAG = VaMKS/( sqrt(VaMKS**2.0 + CsMKS**2.0) + VebMKS)
!Do different calculation of beta if desired
if (doNewBeta) then
bBeta = 2.0*(CsMKS/VaMKS)**2.0 !alternative beta definition
ijTube%beta_average = bBeta
endif
! if (doNewBeta) then
! bBeta = 2.0*(CsMKS/VaMKS)**2.0 !alternative beta definition
! ijTube%beta_average = bBeta
! endif
end associate
@@ -237,7 +236,6 @@ module imagtubes
allocate(isG(Ni,Nj))
isG = .not. (RCMApp%iopen == RCMTOPOPEN)
!Smooth some tubes
call Smooth2D(RCMApp%Vol) !Flux-tube volume
call Smooth2D(RCMApp%pot) !Electrostatic potential

8
src/voltron/rcmimag.F90
View File

@@ -26,7 +26,7 @@ module rcmimag
integer, parameter, private :: MHDPad = 2 !Number of padding cells between RCM domain and MHD ingestion
logical , private :: doWolfLim = .false. !Whether to do wolf-limiting
logical , private :: doBounceDT = .true. !Whether to use Alfven bounce in dt-ingest
logical , private :: doWIMTScl = .true. !Whether to modulate ingestion timescale by wIM
logical , private :: doWIMTScl = .false. !Whether to modulate ingestion timescale by wIM
real(rp), private :: nBounce = 1.0 !Scaling factor for Alfven transit
real(rp), private :: wIM_C = 0.0 !Critical wIM for MHD ingestion inclusion
@@ -254,7 +254,7 @@ module rcmimag
call TrickyTubes(RCMApp)
!Smooth out FTV on tubes b/c RCM will take gradient
!call SmoothTubes(RCMApp)
call SmoothTubes(RCMApp)
!Advance from vApp%time to tAdv
call Tic("AdvRCM")
@@ -575,8 +575,8 @@ module rcmimag
endif
write (*, '(a,2f8.3,a)') ' @ L/MLT = ' , maxL, maxMLT, ' [deg]'
write (*, '(a, f8.3,a,f8.3,a)') ' w/ D = ' , maxD, ' (RC) / ', maxDP, ' (PSPH) [#/cc]'
!write (*, '(a,1f8.3,a)') ' w/ T = ' , maxT, ' [keV]'
write (*, '(a, f8.3,a,f8.3,a)') ' w/ T = ' , maxT, ' [keV] (RCM-Max: ', maxLam, ')'
write (*, '(a,1f8.3,a)') ' w/ T = ' , maxT, ' [keV]'
!write (*, '(a, f8.3,a,f8.3,a)') ' w/ T = ' , maxT, ' [keV] (RCM-Max: ', maxLam, ')'
write(*,'(a)',advance="no") ANSIRESET!, ''