JHUAPL/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-09 18:57:53 -05:00
Code Issues Packages Projects Releases 1 Wiki Activity

Merge branch 'development' into asd

Browse Source
This commit is contained in:
Jeffrey Garretson
2023-08-10 10:22:23 -06:00
parent 35dac4b115 d2a23a5270
commit 095a570d89
17 changed files with 424 additions and 156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
kaipy/kdefs.py
View File

@@ -48,14 +48,18 @@ RNeptuneXE = 3.860 # Rx = X*Re
#------
#Helio
#------
Rsolar = 6.956E5 #[km] Solar radius
kbltz = 1.38e-16 #Boltzmann constant [erg/K]
mp = 1.67e-24 #Proton mass in grams
Tsolar = 25.38 #Siderial solar rotation period
JD2MJD = 2400000.5 #Conversion from JD to MJD: MJD = JD - 2400000.5
Rsolar = 6.957e10 # [cm] Solar radius
kbltz = 1.38e-16 # [erg/K] Boltzmann constant
mp = 1.67e-24 # [g] Proton mass
Tsolar = 25.38 # [days] Siderial solar rotation period
JD2MJD = 2400000.5 # Conversion from JD to MJD: MJD = JD - 2400000.5
Day2s = 86400. # [s] Conversion days => seconds
vc_cgs = 2.99792458e10 # [cm/s] speed of light
#------
#Output defaults
#------
barLen = 30
barLab = 30
barLab = 30

95
scripts/preproc/wsa2gamera.py
View File

@@ -19,19 +19,13 @@ parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
# Read params from config file
prm = params.params(args.ConfigFileName)
Ng = prm.Nghost
gamma = prm.gamma
# Normalization parameters
# remember to use the same units in gamera
B0 = prm.B0
n0 = prm.n0
V0 = B0/np.sqrt(4*np.pi*mp*n0)
TCS = prm.TCS #Temperature in the current sheet for pressure balance calculation
nCS = prm.nCS #Density in the current sheet for pressure balance calculation
TCS = prm.TCS # Temperature in the current sheet for pressure balance calculation
nCS = prm.nCS # Density in the current sheet for pressure balance calculation
# Grid parameters
tMin = prm.tMin
@@ -42,6 +36,12 @@ Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
#conversions from wsa to gamera units
cms2kms = 1.e-5 # cm/s => km/s
Gs2nT = 1.e5 # Gs => nT
# Conversion for E field 1 statV/cm = 3.e7 mV/m
eScl = 3.e7
ffits = prm.wsaFile
# Generate spherical helio grid
@@ -53,15 +53,16 @@ gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
if os.path.exists(prm.gameraGridFile):
print("Grid file heliogrid.h5 is ready!")
# Read and normalize WSA
# Read WSA
jd_c,phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa,T_wsa = wsa.read(ffits,prm.densTempInfile,prm.normalized)
bi_wsa /= B0
n_wsa /= (mp*n0)
v_wsa /= V0
#convert julian date in the center of the WSA map into modified julian date
mjd_c = jd_c - JD2MJD
# Units of WSA input
# bi_wsa in [Gs]
# v_wsa in [cm/s]
# n_wsa in [g cm-3]
# T_wsa in [K]
# convert julian date in the center of the WSA map into modified julian date
mjd_c = jd_c - JD2MJD
# Get GAMERA grid for further interpolation
with h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile),'r') as f:
@@ -79,12 +80,20 @@ zc = 0.125*(z[:-1,:-1,:-1]+z[:-1,1:,:-1]+z[:-1,:-1,1:]+z[:-1,1:,1:]
# radius of the inner boundary
R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2)
r = np.sqrt(x**2+y**2+z**2)
# Calculate phi and theta in physical domain (excluding ghost cells)
P = np.arctan2(y[Ng:-Ng-1,Ng:-Ng-1,:],x[Ng:-Ng-1,Ng:-Ng-1,:])
P = np.arctan2(y[Ng:-Ng,Ng:-Ng,:],x[Ng:-Ng,Ng:-Ng,:])
P[P<0]=P[P<0]+2*np.pi
P = P % (2*np.pi) # sometimes the very first point may be a very
#P = P % (2*np.pi) # sometimes the very first point may be a very
# small negative number, which the above call sets
# to 2*pi. This takes care of it.
T = np.arccos(z[Ng:-Ng,Ng:-Ng,:]/r[Ng:-Ng,Ng:-Ng,:])
#grid for inner i-ghost region; output to innerbc.h5
P_out = P[:,:,0:Ng+1]
T_out = T[:,:,0:Ng+1]
R_out = r[Ng:-Ng,Ng:-Ng,0:Ng+1]
# Calculate r, phi and theta coordinates of cell centers in physical domain (excluding ghost cells)
Rc = np.sqrt(xc[Ng:-Ng,Ng:-Ng,:]**2+yc[Ng:-Ng,Ng:-Ng,:]**2+zc[Ng:-Ng,Ng:-Ng,:]**2)
@@ -115,8 +124,8 @@ rho = f(Pc[:,0,0],Tc[0,:,0])
# Not interpolating temperature, but calculating from the total pressure balance
# AFTER interpolating br and rho to the gamera grid
# n_CS*k*T_CS = n*k*T + Br^2/8pi
temp = (nCS*kbltz*TCS - (br*B0)**2/8./np.pi)/(rho*n0)/kbltz
# note, keep temperature in K (pressure is normalized in wsa.F90)
temp = (nCS*kbltz*TCS - br**2/8./np.pi)*Mp_cgs/rho/kbltz
#temperature in [K]
#check
#print ("Max and min of temperature in MK")
@@ -130,8 +139,12 @@ temp = (nCS*kbltz*TCS - (br*B0)**2/8./np.pi)/(rho*n0)/kbltz
fbi = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],br,kx=1,ky=1)
fv = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],vr,kx=1,ky=1)
br_kface = fbi(P[:,0,0],Tc[0,:,0])
vr_kface = fv (P[:,0,0],Tc[0,:,0])
br_kface = fbi(P[:-1,0,0],Tc[0,:,0]) #(Nk,Nj)
vr_kface = fv (P[:-1,0,0],Tc[0,:,0]) #(Nk,Nj)
# before applying scaling inside ghost region
# get br values to the left of an edge for E_theta calculation
br_kedge = np.roll(br,1, axis=1)
# Scale inside ghost region
(vr,vr_kface,rho,temp,br,br_kface) = [np.dstack(Ng*[var]) for var in (vr,vr_kface,rho,temp,br,br_kface)]
@@ -139,19 +152,33 @@ rho*=(R0/Rc[0,0,:Ng])**2
br*=(R0/Rc[0,0,:Ng])**2
br_kface*=(R0/Rc[0,0,:Ng])**2
# Calculating electric field component on k_edges
# E_theta = B_phi*Vr = - Omega*R*sin(theta)/Vr*Br * Vr = - Omega*R*sin(theta)*Br
omega=2*np.pi/Tsolar
et_kedge = - omega*R0*np.sin(Tc[:,:,Ng-1])*br_kface[:,:,-1]
# Calculating E-field component on k_edges in [mV/m]
# E_theta = B_phi*Vr/c = - Omega*R*sin(theta)/Vr*Br * Vr/c = - Omega*R*sin(theta)*Br/c
omega = 2*np.pi/(Tsolar*Day2s) # [1/s]
# Theta at centers of k-faces (== theta at kedges)
Tcf = 0.25*(T[:,:-1,:-1] + T[:,1:,1:] + T[:,:-1,1:] + T[:,1:,:-1])
et_kedge = - omega*R0*Rsolar*np.sin(Tcf[:-1,:,Ng-1])*br_kedge/vc_cgs #[statV/cm]
# Unit conversion agreement. Input to GAMERA innerbc.h5 has units V[km/s], Rho[cm-3], T[K], B[nT], E[mV/m]
vr *= cms2kms
vr_kface *= cms2kms
rho /= Mp_cgs
br *= Gs2nT
br_kface *= Gs2nT
et_kedge *= eScl
# v, rho, br are normalized, temp is in [K]
with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
hf.attrs["MJD"] = mjd_c
hf.create_dataset("vr",data=vr)
hf.create_dataset("vr_kface",data=vr_kface)
hf.create_dataset("rho",data=rho)
hf.create_dataset("temp",data=temp)
hf.create_dataset("br",data=br)
hf.create_dataset("br_kface",data=br_kface)
#hf.create_dataset("et_kedge",data=et_kedge)
hf.create_dataset("X", data=P_out)
hf.create_dataset("Y", data=T_out)
hf.create_dataset("Z", data=R_out)
grname = "Step#0"
grp = hf.create_group(grname)
grp.attrs.create("MJD", mjd_c)
grp.create_dataset("vr",data=vr)
grp.create_dataset("vr_kface",data=vr_kface) # size (Nk,Nj,Ng)
grp.create_dataset("rho",data=rho)
grp.create_dataset("temp",data=temp)
grp.create_dataset("br",data=br)
grp.create_dataset("br_kface",data=br_kface) # size (Nk,Nj,Ng)
grp.create_dataset("et_kedge",data=et_kedge) # size (Nk, Nj)
hf.close()

5
scripts/quicklook/mixpic.py
View File

@@ -152,7 +152,10 @@ if __name__ == "__main__":
nStp = sorted(sIds)[-1]
if debug:
print("nStp = %s" % nStp)
foundT = T[nStp]
idxStp = np.where(sIds==nStp)[0][0]
if debug:
print("idxStp = %s" % idxStp)
foundT = T[idxStp]
if debug:
print("foundT = %s" % foundT)
print('Found time:', Time(foundT, format='mjd').iso)

1
src/base/defs/chmpdefs.F90
View File

@@ -48,6 +48,7 @@ module chmpdefs
logical :: isRestart=.false.
logical :: isMAGE = .false.
logical :: doDip = .false.
!Output info
type (IOClock_T) :: IO

7
src/base/kdefs.F90
View File

@@ -66,7 +66,11 @@ module kdefs
real(rp), parameter :: REarth = Re_cgs*1.0e-2 !m
real(rp), parameter :: RionE = 6.5 ! Earth Ionosphere radius in 1000 km
real(rp), parameter :: EarthPsi0 = 92.4 ! Corotation potential [kV]
! Earth corotation potential
!real(rp), parameter :: EarthPsi0 = 92.4 ! !!OUTDATED, based on LFM magntic field strength
!real(rp), parameter :: EarthPsi0 = 87.62 ! Calculated using B = 0.2961737 Gauss
real(rp), parameter :: EarthPsi0 = REarth**2 * (2.0*PI/86400.0) * (EarthM0g*G2T) * 1.0e-3 ! [kV]
! Everyone should ideally get Psi0 from planet object, but things specific to Earth should still be able to rely on EarthPsi0
!Saturn
real(rp), parameter :: SaturnM0g = 0.21 !Gauss
@@ -84,6 +88,7 @@ module kdefs
!Helio constants
real(rp), parameter :: Rsolar = 6.956D5 ! [km] Solar radius
real(rp), parameter :: Msolar = 1.98847D30 ! [kg] Solar mass
real(rp), parameter :: Tsolar_synodic = 27.28 ![days] synodic Tsolar
!Numbered accessors
!Directions

52
src/base/planethelper.F90
View File

@@ -17,6 +17,7 @@ module planethelper
character(len=strLen) :: pID !Planet ID string
logical :: doCorot
real(rp) :: RIon
real(rp) :: period ! Rotation period in seconds
if (present(pStrO)) then
pID = pStrO
@@ -33,7 +34,9 @@ module planethelper
planet%ri_m = RionE*1.e6 ! Defined in kdefs in 10000km
planet%grav = 9.807
call xmlInp%Set_val(planet%magMoment, "/Kaiju/Gamera/prob/M0", EarthM0g)
planet%psiCorot = EarthPsi0
call xmlInp%Set_val(period, "/Kaiju/Gamera/prob/period", 86400.0)
planet%psiCorot = CorotPotential(planet%rp_m, period, planet%magMoment)
write(*,*)" Calculated corotation potential [kV]: ", planet%psiCorot
planet%doGrav = .true.
case("Saturn","saturn","SATURN")
planet%rp_m = RSaturnXE*REarth
@@ -127,6 +130,47 @@ module planethelper
write(*,*) "-------------"
end subroutine
subroutine writePlanetParams(planet, doIOp, h5File, gStrO)
!! Write planet parameters as attributes to a "Planet" group
!! Include a gStrO if you want it to be nested within some other group for some reason
type(planet_T), intent(in) :: planet
logical, intent(in) :: doIOp
character(len=strLen) :: h5File
character(len=strLen), optional :: gStrO
type(IOVAR_T), dimension(8) :: IOVars
call clearIO(IOVars)
! All attributes
call AddOutVar(IOVars, "Name", planet%name)
call AddOutVar(IOVars, "Rad_surface" , planet%rp_m, uStr="m")
call AddOutVar(IOVars, "Rad_ionosphere", planet%ri_m, uStr="m")
call AddOutVar(IOVars, "Gravity", planet%grav, uStr="m/s^2")
call AddOutVar(IOVars, "Mag Moment", planet%magMoment, uStr="Gauss")
call AddOutVar(IOVars, "Psi Corot", planet%psiCorot, uStr="kV")
if (planet%doGrav) then
call AddOutVar(IOVars, "doGrav", "True")
else
call AddOutVar(IOVars, "doGrav", "False")
endif
if(present(gStrO)) then
call WriteVars(IOVars, doIOp, h5File, gStrO, "Planet")
else
call WriteVars(IOVars, doIOp, h5File, "Planet")
endif
end subroutine writePlanetParams
function CorotPotential(Rp_m, period, Bmag) result(cPot)
!! Calculates corotation potential, assuming dipole and rotational axes are aligned
real(rp), intent(in) :: Rp_m ! Planetary radius in meters
real(rp), intent(in) :: period ! Rotation period in seconds
real(rp), intent(in) :: Bmag ! Magnetic moment in Gauss
real(rp) :: cPot ! Corotation potential in kV
cPot = Rp_m**2.0 * (2.0*PI/period) * (Bmag*G2T) * 1.0e-3 ! [kV]
end function CorotPotential
! Helpful thermo conversion functions
@@ -255,6 +299,12 @@ module planethelper
end function MirrorRatio
function DipColat2L(colat) result(L)
real(rp), intent(in) :: colat
real(rp) :: L
L = 1.0/sin(colat)**2
end function DipColat2L
!Calculate FTV of dipole, Rx/nT
!M0g is optional mag moment in Gauss, otherwise use Earth
function DipFTV_L(L,M0gO) result(V)

163
src/gamera/ICs/helio/wsa.F90
View File

@@ -14,16 +14,22 @@ module usergamic
enum, bind(C)
! variables passed via innerbc file
! Br, Vr, Rho, Temperature, Br @ kface, Vr @ kface
enumerator :: BRIN=1,VRIN,RHOIN,TIN,BRKFIN,VRKFIN
! Br, Vr, Rho, Temperature, Br @ kface, Vr @ kface, et @ kedge
enumerator :: BRIN=1,VRIN,RHOIN,TIN,BRKFIN,VRKFIN,ETKEDGE
endenum
integer, private, parameter :: NVARSIN=6 ! SHOULD be the same as the number of vars in the above enumerator
integer, private, parameter :: NVARSIN=7 ! SHOULD be the same as the number of vars in the above enumerator
real(rp), dimension(:,:,:,:), allocatable :: ibcVars
real(rp), dimension(:,:), allocatable :: ibcEt
!Various global would go here
real (rp) :: Rho0, P0, Vslow,Vfast, wScl, Cs0, B0, MJD_c
!Scaling from innerbc to CGS
!Elena: May be move these conversion factors to kdefs?
real (rp) :: km2cm = 1.e5
real (rp) :: nT2Gs = 1.e-5
! things we keep reusing
real(rp), dimension(NDIM) :: xyz,xyz0,rHat,phiHat
@@ -50,7 +56,7 @@ module usergamic
character(len=strLen) :: wsaFile
! use this to fix the Efield at the inner boundary
! real(rp), allocatable :: inEijk(:,:,:,:)
real(rp), allocatable :: inEijk(:,:,:,:)
! type for solar wind BC
type, extends(innerIBC_T) :: SWInnerBC_T
@@ -67,6 +73,8 @@ module usergamic
contains
!> Initialize the run
!> Set initial conditions
subroutine initUser(Model,Grid,State,inpXML)
type(Model_T), intent(inout) :: Model
type(Grid_T), intent(inout) :: Grid
@@ -80,13 +88,14 @@ module usergamic
integer :: n1, n2
integer :: kg, ke, kb, jg
real(rp) :: a
real(rp) :: Tsolar_synodic
real(rp) :: Br_left
real(rp) :: ibcVarsStaticBRKFN !For brkfn
real(rp) :: ibcVarsStatic !For br only
! if (.not.allocated(inEijk)) allocate(inEijk(1,Grid%jsg:Grid%jeg+1,Grid%ksg:Grid%keg+1,1:NDIM))
if (.not.allocated(inEijk)) allocate(inEijk(1,Grid%jsg:Grid%jeg+1,Grid%ksg:Grid%keg+1,1:NDIM))
inEijk = 0.0
! set units and other thins, like Tsolar
! set units and other things, like Tsolar
call setHeliosphere(Model,inpXML,Tsolar)
! grab inner
@@ -136,21 +145,18 @@ module usergamic
Grid%ksDT = Grid%ks
Grid%keDT = Grid%ke
! Add gravity
! eHack => EFix
! Model%HackE => eHack
! tsHack => PerStep
! Model%HackStep => tsHack
!Fixing Efield at the inner boundary
eHack => EFix
Model%HackE => eHack
!Write MJD_center_of_WSA_map to the root of H5 output
Model%HackIO_0 => writeMJDcH5Root
! Write MJD_center_of_WSA_map to the root of H5 output
Model%HackIO_0 => writeMJDcH5Root
! everybody reads WSA data
call readIBC(wsaFile)
call readIBC(wsaFile,Model)
!Set MJD0 which corresponds to sim time 0
!MJD0 is MJD_center_of_WSA_map - Tsolar_synodic/2; Tsolar_synodic = 27.28
![EP] TO DO: add synodic Tsolar to constants
Tsolar_synodic = 27.28
Model%MJD0 = MJD_c - Tsolar_synodic/2.
!if not restart set State%time according the tSpin
@@ -167,13 +173,13 @@ module usergamic
do k=Grid%ks,Grid%ke
kg = k+Grid%ijkShift(KDIR)
! map rotating to static grid
!initial spin of a map by Tsolar/2
call mapKinit(kg,ke,kb,a)
do j=Grid%js,Grid%je
jg = j+Grid%ijkShift(JDIR)
do i=Grid%isg,Grid%ieg+1
do i=Grid%isg,Grid%ieg+1
ibcVarsStatic = 1._rp
if ( (jg>=Grid%js).and.(jg<=size(ibcVars,2)) ) then
@@ -189,15 +195,13 @@ module usergamic
!
! note also that ibcVars(Model%Ng,:,:,:) corresponds to the center of the first ghost cell (just below the boundary)
!State%magFlux(i,j,k,IDIR) = ibcVars(Model%Ng,j+Grid%ijkShift(JDIR),k+Grid%ijkShift(KDIR),BRIN)*Rfactor**2*Grid%face(Grid%is,j,k,IDIR)
State%magFlux(i,j,k,IDIR) = ibcVarsStatic*Rfactor**2*Grid%face(Grid%is,j,k,IDIR)
State%magFlux(i,j,k,IDIR) = ibcVarsStatic*Rfactor**2*Grid%face(Grid%is,j,k,IDIR)
enddo
enddo
enddo
!Local functions
!NOTE: Don't put BCs here as they won't be visible after the initialization call
contains
!> Set initial conditions
subroutine GasIC(x,y,z,D,Vx,Vy,Vz,P)
real (rp), intent(in) :: x,y,z
real (rp), intent(out) :: D,Vx,Vy,Vz,P
@@ -219,6 +223,8 @@ module usergamic
Vz = r_unit(ZDIR)*Vslow
end subroutine GasIC
!> Rotates initial WSA map so that at simulation time 0
!> the start date of the map faces Earth (-X axis)
subroutine mapKinit(k,ke,kb,a)
! find the lower and upper k-index of the rotating cell on the static grid
integer, intent(in) :: k
@@ -242,22 +248,27 @@ module usergamic
end subroutine initUser
!Inner-I BC for WSA-Gamera
!> Setting inner-I boundary conditions in four ghost i-layers of cells
!> Linear interpolation from rotating to inertial grid
!> Setting cell-centered variables and magnetic fluxes at cell faces
!> Calculate user defined Efield at the inner boundary
subroutine wsaBC(bc,Model,Grid,State)
class(SWInnerBC_T), intent(inout) :: bc
type(Model_T), intent(in) :: Model
type(Grid_T), intent(in) :: Grid
type(State_T), intent(inout) :: State
procedure(HackE_T), pointer :: eHack
! local variables
integer :: i,j,k, kb, ke
integer :: kg, jg, ig ! global indices (in preparation for MPI)
integer :: kg, jg, ig ! global indices
integer :: var ! ibcVar variable number
real(rp) :: a
real(rp) :: ibcVarsStatic(NVARSIN)
real(rp) :: R, Theta, Phi
real(rp) :: Theta_kf, R_kf ! kface
real(rp), dimension(NVAR) :: conVar, pVar
real(rp) :: Br_left
!i-boundaries (IN)
!$OMP PARALLEL DO default(shared) &
@@ -286,6 +297,9 @@ module usergamic
do var=1,NVARSIN
ibcVarsStatic(var) = a*ibcVars(ig,jg,kb,var)+(1-a)*ibcVars(ig,jg,ke,var)
end do
! use Br value from the left (kb) cell of the rotating grid to calculate Ej
! Ensures constant Ej between Br updates and linearly changing Br in time.
Br_left = ibcVars(ig,jg,kb,BRIN)
end if
! do cell centered things for cell-centers only
@@ -299,7 +313,6 @@ module usergamic
phiHat = [-sin(phi),cos(phi),0._rp]
! NOTE, WSA data were already scaled appropriately in the python code
! TODO: save them in the innerbc hdf file and set in helioutils appropriately
!Set primitives
pVar(VELX:VELZ) = rHat*ibcVarsStatic(VRIN)
@@ -311,28 +324,38 @@ module usergamic
call CellP2C(Model,pVar,conVar)
State%Gas(i,j,k,:,BLK) = conVar
! note, don't need cc Bxyz because we run flux2field through ghosts
! note, no need to set cc Bxyz because we run flux2field through ghosts
end if
! also need theta at k-face for k-flux
! although we're assuming theta and R don't change from cell to k-face center,
! also need to calculate theta at k-face for k-flux
! although we're assuming theta and R don't change from cell center to k-face center,
! xyzcc used under if statement above is only defined for cell centers so need to define it here
xyz0 = Grid%xfc(i,j,k,:,KDIR) !just reusing a temp var here.
R_kf = norm2(xyz0)
Theta_kf = acos(xyz0(ZDIR)/R_kf)
! Set magnetic fluxes
! note scaling for Bi. See note above in InitUser
Rfactor = Rbc/norm2(Grid%xfc(Grid%is,j,k,:,IDIR))
State%magFlux(i,j,k,IDIR) = ibcVarsStatic(BRIN)*Rfactor**2*Grid%face(Grid%is,j,k,IDIR)
State%magFlux(i,j,k,JDIR) = 0.0
State%magFlux(i,j,k,KDIR) = - 2*PI/Tsolar*R_kf*sin(Theta_kf)/ibcVarsStatic(VRKFIN)*ibcVarsStatic(BRKFIN)*Grid%face(i,j,k,KDIR)
State%magFlux(i,j,k,KDIR) = -2*PI/Tsolar*R_kf*sin(Theta_kf)/ibcVarsStatic(VRKFIN)*ibcVarsStatic(BRKFIN)*Grid%face(i,j,k,KDIR)
! Calculate user`s electric fields at the inner boundary
if (ig == Model%Ng) then
! calculate E_theta at edges using Br_left value defined above
inEijk(Grid%is,j,k,JDIR) = -2*PI/Tsolar*Rbc*sin(Theta_kf)*Br_left*Rfactor**2
end if
end do
end do
end do
contains
contains
!> Find the lower and upper k-index for interpolation of the rotating cell on the static grid
subroutine mapK(k,ke,kb,a)
! find the lower and upper k-index of the rotating cell on the static grid
integer, intent(in) :: k
integer, intent(out) :: ke,kb ! upper (ke) and lower (kb) indices
real(rp), intent(out):: a ! interp coefficient
@@ -354,25 +377,29 @@ module usergamic
end subroutine wsaBC
subroutine eFix(Model,Gr,State)
!> Set user defined tangential Efield to enforce rotating WSA map
subroutine EFix(Model,Gr,State)
type(Model_T), intent(in) :: Model
type(Grid_T), intent(in) :: Gr
type(Grid_T), intent(inout) :: Gr
type(State_T), intent(inout) :: State
! see example of how to do this in voltron/ICs/earthcmi.F90
!Grid%externalBCs(1)%p
! State%Efld(Gr%is ,:,:,JDIR:KDIR) = inEijk(1,:,:,JDIR:KDIR)*Gr%edge(Gr%is ,:,:,JDIR:KDIR)
if (Gr%hasLowerBC(IDIR)) then
State%Efld(Gr%is,:,:,JDIR) = inEijk(1,:,:,JDIR)*Gr%edge(Gr%is ,:,:,JDIR)
State%Efld(Gr%is,:,:,KDIR) = 0.
endif
end subroutine EFix
end subroutine eFix
subroutine readIBC(ibcH5)
!> Read in innerbc.h5 file with inner boundary conditions
subroutine readIBC(ibcH5,Model)
character(len=*), intent(in) :: ibcH5
type(Model_T), intent(in) :: Model
logical :: fExist
integer :: i,nvar,dims(3)
integer :: i,nvar,dims(3), dimsE(2)
integer, parameter :: MAXIOVAR = 50
type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
character(len=10) :: Grp = "Step#0"
!Reset IO chain
call ClearIO(IOVars)
@@ -385,20 +412,27 @@ module usergamic
endif
!Setup input chain
call AddInVar(IOVars,"vr")
call AddInVar(IOVars,"vr_kface")
call AddInVar(IOVars,"rho")
call AddInVar(IOVars,"temp")
call AddInVar(IOVars,"br")
call AddInVar(IOVars,"br_kface")
call AddInVar(IOVars,"vr", vSclO=km2cm/Model%Units%gv0)
call AddInVar(IOVars,"vr_kface", vSclO=km2cm/Model%Units%gv0)
call AddInVar(IOVars,"rho", vSclO=1/Model%Units%gD0)
call AddInVar(IOVars,"temp") !Temp in K
call AddInVar(IOVars,"br", vSclO=nT2Gs/Model%Units%gB0)
call AddInVar(IOVars,"br_kface", vSclO=nT2Gs/Model%Units%gB0)
call AddInVar(IOVars,"et_kedge", vSclO=1.e-3/(Model%Units%gv0*1.e-2*Model%Units%gB0*1.e-4))
call AddInVar(IOVars,"MJD")
call ReadVars(IOVars,.false.,ibcH5) !Don't use io precision
call ReadVars(IOVars,.false.,ibcH5,Grp) !Don't use io precision
! NOTE, assuming they all have the same dimesnions here (NO2,NJ,NK)
! see wsa2gamera
dims=IOVars(1)%dims(1:3) ! i,j,k
dimsE(:) = dims(2:3) ! (Nj,Nk) for E_theta
if (.not.allocated(ibcVars)) allocate(ibcVars(dims(1),dims(2),dims(3),NVARSIN))
if (.not.allocated(ibcEt))allocate(ibcEt(dimsE(1),dimsE(2)))
! Zero out
ibcVars = 0.
ibcEt = 0.
do i=1,NVARSIN
select case (i)
@@ -414,22 +448,31 @@ module usergamic
nvar= FindIO(IOVars,"br_kface")
case (VRKFIN)
nvar= FindIO(IOVars,"vr_kface")
case (ETKEDGE)
nvar= FindIO(IOVars,"et_kedge")
end select
ibcVars(:,:,:,i) = reshape(IOVars(nvar)%data,dims)
if (i <= NVARSIN-1) then
ibcVars(:,:,:,i) = reshape(IOVars(nvar)%data*IOVars(nvar)%scale,dims)
else if (i == ETKEDGE) then
ibcEt = reshape(IOVars(nvar)%data*IOVars(nvar)%scale,dimsE)
end if
end do
!reading modified julian date from innerbc
MJD_c = GetIOReal(IOVars,"MJD")
!reading modified julian date from innerbc
MJD_c = GetIOReal(IOVars,"MJD")
end subroutine readIBC
subroutine writeMJDcH5Root(Model,Grid,IOVars)
type(Model_T), intent(in) :: Model
type(Grid_T) , intent(in) :: Grid
type(IOVAR_T), dimension(:), intent(inout) :: IOVars
!> Writing Modified Julian Date at the center of WSA map to the root of h5 file
!> MJDc is required to set the GAMERA Helio coordinate system
subroutine writeMJDcH5Root(Model,Grid,IOVars)
type(Model_T), intent(in) :: Model
type(Grid_T) , intent(in) :: Grid
type(IOVAR_T), dimension(:), intent(inout) :: IOVars
call AddOutVar(IOVars,"MJDc", MJD_c)
call AddOutVar(IOVars,"MJDc", MJD_c)
end subroutine writeMJDcH5Root
end subroutine writeMJDcH5Root
end module usergamic

21
src/gamera/apps/helioutils.F90
View File

@@ -12,7 +12,7 @@ module helioutils
implicit none
! normalization
real(rp), private :: gD0, gB0, gx0, gT0, gv0, gP0
real(rp), private :: gD0, gB0, gx0, gT0, gv0, gP0, gE0
type, extends(innerJBC_T) :: helioInnerJBC_T
contains
@@ -41,9 +41,9 @@ module helioutils
logical :: doSpin
! normalization
gD0=200. ! [/cc]
gB0=1.e-3 ! [Gs], 100 nT
gx0=Rsolar*1.e5 ! [cm], solar radius
gD0 = 200. ! [/cc]
gB0 = 1.e-3 ! [Gs], 100 nT
gx0 = Rsolar*1.e5 ! [cm], solar radius
! for Ohelio case
!gD0=10. ! [/cc]
!gB0=5.e-5 ! [Gs], 5 nT
@@ -54,6 +54,9 @@ module helioutils
gT0 = gx0/gv0 ! [s] ~ 1.25 hour for above values
gP0 = gB0**2/(4*pi) ! [erg/cm3]
!normalization of E-field
gE0 = gB0*1.e-4*gv0*1.e-2 ! B0[T]*V0[m/s]
! Use gamma=1.5 for SW calculations (set in xml, but defaults to 1.5 here)
call inpXML%Set_Val(Model%gamma,"physics/gamma",1.5_rp)
call inpXML%Set_Val(Tsolar,"prob/Tsolar",25.38_rp) ! siderial solar rotation in days
@@ -64,12 +67,10 @@ module helioutils
if (doSpin) then
call inpXML%Set_Val(tSpin,"spinup/tSpin",200.0_rp) ! set in [hr] in xml
!Rewind Gamera helio time to negative tSpin
if (.not. Model%isRestart) then
Model%t = -tSpin*3600./gT0
call inpXML%Set_Val(tIO,"spinup/tIO",0.0_rp) !Time of first restart and output
Model%IO%tRes = tIO
Model%IO%tOut = tIO
endif
Model%t = -tSpin*3600./gT0
call inpXML%Set_Val(tIO,"spinup/tIO",0.0_rp) !Time of first restart and output
Model%IO%tRes = tIO
Model%IO%tOut = tIO
endif
endif

8
src/gamera/fields.F90
View File

@@ -312,12 +312,16 @@ module fields
bR = 0.0
VolB = 0.0
!Get stencils for this sweep
!dT1 faces in dT2 direction
!Get stencils for this sweep
!dT1 faces in dT2 direction
call recLoadBlockI(Model,Gr,AreaB ,Gr%Face(:,:,:,dT1),iB,j,k,iMax,dT2)
call recLoadBlockI(Model,Gr, MagB(:,:,1),bFluxn (:,:,:,dT1),iB,j,k,iMax,dT2)
call recLoadBlockI(Model,Gr,MagB(:,:,1),bFluxHf (:,:,:,dT1),iB,j,k,iMax,dT2)
call limLoadBlockI(Model,Gr,nAreaB ,Gr%Face(:,:,:,dT1),iB,j,k,iMax,dT2)
call limLoadBlockI(Model,Gr,nMagB(:,:,1),bFluxHf(:,:,:,dT1),iB,j,k,iMax,dT2)
call limLoadBlockI(Model,Gr,nMagB(:,:,1),bFluxn (:,:,:,dT1),iB,j,k,iMax,dT2)
!Split into L/Rs
if (doBdA) then

6
src/gamera/gioH5.F90
View File

@@ -736,7 +736,7 @@ module gioH5
Model%ts = 0
Model%t = tReset
else
Model%IO%nOut = GetIOInt(IOVars,"nOut")
Model%IO%nOut = GetIOInt(IOVars,"nOut") + 1
Model%IO%nRes = GetIOInt(IOVars,"nRes") + 1
Model%ts = GetIOInt(IOVars,"ts")
Model%t = GetIOReal(IOVars,"t")
@@ -788,8 +788,8 @@ module gioH5
endif !Gas0
!Do touchup to data structures
Model%IO%tOut = floor(Model%t/Model%IO%dtOut)*Model%IO%dtOut
Model%IO%tRes = Model%t + Model%IO%dtRes
Model%IO%tOut = floor(Model%t/Model%IO%dtOut)*Model%IO%dtOut + Model%IO%dtOut
Model%IO%tRes = floor(Model%t/Model%IO%dtRes)*Model%IO%dtRes + Model%IO%dtRes
end subroutine readH5Restart

58
src/remix/mixconductance.F90
View File

@@ -369,37 +369,35 @@ module mixconductance
real(rp) :: mad,Ttmp
real(rp) :: temp(G%Np,G%Nt)
MaxIter = 5
MaxIter = 15
! Iterative diffusion algorithm to smooth out IM_GTYPE: 0/1 are boundary cells. Evolve with nine-cell mean.
! Otherwise it only has three values, 0, 0.5, and 1.0,
! which causes discontinuities in the merged precipitation.
gtype_RCM = St%Vars(:,:,IM_GTYPE) ! supposed to be between 0 and 1.
call FixPole(G,gtype_RCM)
do it=1,MaxIter
mad = 0.D0 ! max abs difference from last iteration.
temp = gtype_RCM
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,jm1,jp1,im1,ip1,Ttmp) &
!$OMP reduction(max:mad)
do j=1,G%Nt ! use open BC for lat.
if(j==1) then
jm1 = 1
else
jm1 = j-1
endif
if(j==G%Nt) then
jp1 = G%Nt
else
jp1 = j+1
endif
do i=1,G%Np ! use periodic BC for lon.
if(i==1) then
im1 = G%Np
else
im1 = i-1
endif
if(i==G%Np) then
ip1 = 1
else
ip1 = i+1
endif
jm1 = j-1
jp1 = j+1
if (j == 1) jm1 = 1
if (j == G%Nt) jp1 = G%Nt
im1 = i-1
ip1 = i+1
if (i == 1) im1 = G%Np
if (i == G%Np) ip1 = 1
if(St%Vars(i,j,IM_GTYPE)>0.01 .and. St%Vars(i,j,IM_GTYPE)<0.99) then
! Use 0.01/0.99 as the boundary for this numerical diffusion because
! the interpolation from RCM to REMIX would slightly modify the 0/1 boundary.
@@ -411,12 +409,26 @@ module mixconductance
endif
enddo
enddo
if(mad<0.05) exit
call FixPole(G,gtype_RCM)
if(mad<0.025) exit
enddo
gtype_RCM = min(gtype_RCM,1.0)
gtype_RCM = min(gtype_RCM,1.0)
gtype_RCM = max(gtype_RCM,0.0)
end subroutine conductance_IM_GTYPE
!Enforce pole condition that all values at the same point are equal
subroutine FixPole(Gr,Q)
type(mixGrid_T), intent(in) :: Gr
real(rp) , intent(inout) :: Q(Gr%Np,Gr%Nt)
real(rp) :: Qavg
Qavg = sum(Q(:,2))/Gr%Np
Q(:,1) = Qavg
end subroutine FixPole
subroutine conductance_rcmhd(conductance,G,St)
type(mixConductance_T), intent(inout) :: conductance
type(mixGrid_T), intent(in) :: G

34
src/voltron/ICs/earthcmi.F90
View File

@@ -459,9 +459,7 @@ module uservoltic
!dA = Grid%face(ig,j,k,d)/Grid%face(Grid%is,jp,kp,d)
dA = 1.0 !Using dA=1 for smoother magflux stencil
if((d == IDIR .and. j .le. Grid%jsg .and. k .le. Grid%ksg) .or. &
(d == JDIR .and. ig .le. Grid%isg .and. k .le. Grid%ksg) .or. &
(d == KDIR .and. ig .le. Grid%isg .and. j .le. Grid%jsg)) then
if (isGoodFaceIJK(Model,Grid,ig,j,k,d)) then
if ( isLowLat(Grid%xfc(ig,j,k,:,d),llBC) ) then
!State%magFlux(ig,j,k,d) = 0.0
State%magFlux(ig,j,k,d) = dApm(d)*dA*State%magFlux(Grid%is,jp,kp,d)
@@ -469,14 +467,40 @@ module uservoltic
State%magFlux(ig,j,k,d) = dApm(d)*dA*State%magFlux(Grid%is,jp,kp,d)
endif
endif
enddo
enddo !n loop (ig)
enddo !j loop
enddo !k loop
contains
function isGoodFaceIJK(Model,Grid,i,j,k,d) result(isGood)
type(Model_T), intent(in) :: Model
type(Grid_T) , intent(in) :: Grid
integer , intent(in) :: i,j,k,d
logical :: isGood
select case(d)
case(IDIR)
isGood = (i >= Grid%isg) .and. (i <= Grid%ieg+1) &
.and. (j >= Grid%jsg) .and. (j <= Grid%jeg ) &
.and. (k >= Grid%ksg) .and. (k <= Grid%keg )
case(JDIR)
isGood = (i >= Grid%isg) .and. (i <= Grid%ieg ) &
.and. (j >= Grid%jsg) .and. (j <= Grid%jeg+1) &
.and. (k >= Grid%ksg) .and. (k <= Grid%keg )
case(KDIR)
isGood = (i >= Grid%isg) .and. (i <= Grid%ieg ) &
.and. (j >= Grid%jsg) .and. (j <= Grid%jeg ) &
.and. (k >= Grid%ksg) .and. (k <= Grid%keg+1)
end select
end function isGoodFaceIJK
function isLowLat(xyz,llBC)
real(rp), dimension(NDIM), intent(in) :: xyz
real(rp), intent(in) :: llBC

9
src/voltron/ebsquish.F90
View File

@@ -25,7 +25,6 @@ module ebsquish
real(rp), parameter, private :: startEps = 0.05
real(rp), parameter, private :: rEps = 0.125
real(rp), parameter, private :: ShueScl = 1.25 !Safety factor for Shue MP
logical , parameter, private :: doDipTest = .false.
contains
!Adjust the starting indices of the squish blocks according to the passed in array values
@@ -203,7 +202,7 @@ module ebsquish
call GetSquishBds(vApp,ksB,keB)
call Tic("SQ-Project")
if (doDipTest) write(*,*) "Using fake projection for testing!"
if (ebModel%doDip) write(*,*) "Using fake projection for testing!"
!$OMP PARALLEL DO default(shared) collapse(2) &
!$OMP schedule(dynamic) &
@@ -470,10 +469,10 @@ module ebsquish
isGood = inShueMP_SM(xyz,ShueScl)
if (.not. isGood) return
if (doDipTest) then
if (ebApp%ebModel%doDip) then
xyz0 = DipoleShift(xyz,norm2(xyz)+startEps)
x1 = InvLatitude(xE)
x2 = atan2(xE(YDIR),xE(XDIR))
x1 = InvLatitude(xyz0)
x2 = atan2(xyz0(YDIR),xyz0(XDIR))
if (x2 < 0) x2 = x2 + 2*PI
return
endif

5
src/voltron/voltapp.F90
View File

@@ -151,8 +151,8 @@ module voltapp
if(gApp%Model%isRestart) then
call readVoltronRestart(vApp, xmlInp)
vApp%IO%tOut = floor(vApp%time/vApp%IO%dtOut)*vApp%IO%dtOut
vApp%IO%tRes = vApp%time + vApp%IO%dtRes
vApp%IO%tOut = floor(vApp%time/vApp%IO%dtOut)*vApp%IO%dtOut + vApp%IO%dtOut
vApp%IO%tRes = floor(vApp%time/vApp%IO%dtRes)*vApp%IO%dtRes + vApp%IO%dtRes
vApp%IO%tsNext = vApp%ts
if(vApp%isSeparate) then
gApp%Model%ts = vApp%ts
@@ -653,6 +653,7 @@ module voltapp
Model%doMHD = .true.
call inpXML%Set_Val(Model%epsds,'tracer/epsds',1.0e-2)
call setBackground(Model,inpXML)
call inpXML%Set_Val(Model%doDip,'tracer/doDip',.false.)
!Initialize ebState
!CHIMP grid is initialized from Gamera's active corners

2
src/voltron/voltio.F90
View File

@@ -276,7 +276,7 @@ module voltio
stop
endif
vApp%IO%nOut = GetIOInt(IOVars,"nOut")
vApp%IO%nOut = GetIOInt(IOVars,"nOut") + 1
vApp%IO%nRes = GetIOInt(IOVars,"nRes") + 1
vApp%ts = GetIOInt(IOVars,"ts")
vApp%MJD = GetIOReal(IOVars,"MJD")

13
tests/voltron/testplanetunits.pf
View File

@@ -5,6 +5,7 @@ module testplanetunits
use chmpunits
use planethelper
use uservoltic
use xml_input
use ioH5
use kdefs
use msphutils, only : RadIonosphere
@@ -126,11 +127,21 @@ contains
type(voltApp_T) :: vApp
type(gamApp_T) :: gApp
character(len=strLen) :: xmlName = 'cmriD_NewPlanet.xml'
type(XML_Input_T) :: inpXML
real(rp) :: corotXML
call initGamera(gApp, initUser, xmlName)
call initVoltron(vApp, gApp, xmlName)
! For this test, we should get corot potential right from xml cause we don't have a real pre-set value for a fake planet.
! This does check if Voltron and Gamera are getting the right potential from xml
! Currently, we know RCM does not, and that kai2geo doesn't need to cause its only for Earth
! so, RCM and kai2geo should only be used for Earth Earth
inpXML = New_XML_Input(trim(xmlName),"Kaiju/GAMERA",.true.)
call inpXML%Set_Val(corotXML , "prob/Psi0", EarthM0g)
! Run tests
call checkVoltron(vApp, 1.3*REarth, 1.01*1.3*REarth, 12.0_rp, 0.3*EarthM0g, 0.1*EarthPsi0, .true.)
call checkVoltron(vApp, 1.3*REarth, 1.01*1.3*REarth, 12.0_rp, 0.3*EarthM0g, corotXML, .true.)
call checkGamera(vApp, gApp)
call checkRCM(vApp)
!Does REMIX need to be checked? Only relies on msphutil's RadIonosphere() (2021/10/12)

85
tests/voltron_mpi/testhelpers.pf
View File

@@ -168,7 +168,90 @@ contains
! now calculate local-only squish
voltAppMpi%ebTrcApp%ebSquish%myNumBlocks = -1 ! I do all blocks
voltAppMpi%ebTrcApp%ebSquish%myFirstBlock = 0 ! ensure I start with the first block
voltAppMpi%ebTrcApp%ebSquish%myFirstBlock = 1 ! ensure I start with the first block
call DeepUpdate(voltAppMpi, voltAppMpi%gAppLocal)
sqErr = 0
maxSqErr = 0
do i=voltAppMpi%ebTrcApp%ebState%ebGr%is,voltAppMpi%iDeep+1
do j=voltAppMpi%ebTrcApp%ebState%ebGr%js,voltAppMpi%ebTrcApp%ebState%ebGr%je+1
do k=voltAppMpi%ebTrcApp%ebState%ebGr%ks,voltAppMpi%ebTrcApp%ebState%ebGr%ke+1
posErr = abs(distSquish(i,j,k,:) - voltAppMpi%chmp2mhd%xyzSquish(i,j,k,:))
if(posErr(2) > PI) posErr(2) = 2*PI - posErr(2)
sqErr = sqErr + NORM2(posErr)
maxSqErr = MAX(maxSqErr, NORM2(posErr))
enddo
enddo
enddo
! debugging...
!call ClearIO(IOVars)
!call AddOutVar(IOVars,"distSquish",distSquish)
!call AddOutVar(IOVars,"localSquish",voltAppMpi%chmp2mhd%xyzSquish)
!call WriteVars(IOVars,.false.,"SquishTestData.h5")
! ...debugging
deallocate(distSquish)
! Floating point error rounding
@assertLessThan(sqErr, 1.0e-15_rp, trim("Total Distributed Squish Error Too Large"))
@assertLessThan(maxSqErr, 1.0e-15_rp, trim("Per Cell Distributed Squish Error Too Large"))
else
! helpers
helperQuit = .false.
do while(.not. helperQuit)
call helpVoltron(voltAppMpi, helperQuit)
end do
endif
endif
end subroutine
@test(npes=[8])
subroutine testHelpSquishDip(this)
class (MpiTestMethod), intent(inout) :: this
real(rp), dimension(:,:,:,:), allocatable :: distSquish
real(rp), dimension(2) :: posErr
real(rp) :: sqErr, maxSqErr
integer :: i,j,k
logical :: helperQuit
! debugging
!type(IOVAR_T), dimension(12) :: IOVars
write (*,'(a,I0)') 'Testing HelpSquish ',this%getNumProcesses()
call manualSetup(this, 'testHelpersSquish_4.xml')
if(isGamera) then
call runApplication()
else
voltAppMpi%ebTrcApp%ebModel%doDip = .true.
call runApplication()
endif
! Now compare distributed squish results to local-only squish results
if(.not. isGamera) then
if(.not. voltAppMpi%amHelper) then
! call doImag once to update rTrc and nTrc before doing any squishing
call PreDeep(voltAppMpi, voltAppMpi%gAppLocal)
call DoImag(voltAppMpi)
! calculate distributed squish from most recent results
call startDeep(voltAppMpi)
do while(deepInProgress(voltAppMpi))
call doDeepBlock(voltAppMpi)
enddo
call vhReqHelperQuit(voltAppMpi) ! tell the helpers they're done
allocate(distSquish, MOLD=voltAppMpi%chmp2mhd%xyzSquish)
distSquish = voltAppMpi%chmp2mhd%xyzSquish
! now calculate local-only squish
voltAppMpi%ebTrcApp%ebSquish%myNumBlocks = -1 ! I do all blocks
voltAppMpi%ebTrcApp%ebSquish%myFirstBlock = 1 ! ensure I start with the first block
call DeepUpdate(voltAppMpi, voltAppMpi%gAppLocal)
sqErr = 0