diff --git a/kaiju.sublime-project b/kaiju.sublime-project
index ac201879..b25bd630 100644
--- a/kaiju.sublime-project
+++ b/kaiju.sublime-project
@@ -26,4 +26,13 @@
"trim_automatic_white_space": true,
},
+ "build_systems":
+ [
+ {
+ "file_regex": "^[ ]*File \"(...*?)\", line ([0-9]*)",
+ "name": "Anaconda Python Builder",
+ "selector": "source.python",
+ "shell_cmd": "\"python\" -u \"$file\""
+ }
+ ],
}
diff --git a/kaipy/gamhelio/ConfigScripts/startup.config b/kaipy/gamhelio/ConfigScripts/startup.config
index 05e0b937..23a6948b 100644
--- a/kaipy/gamhelio/ConfigScripts/startup.config
+++ b/kaipy/gamhelio/ConfigScripts/startup.config
@@ -1,4 +1,9 @@
-#Modify if needed paths to a grid file, output file and WSA fits file
+;Comments and definitions:
+;If needed, modify the paths to the grid file, output innerbc file and WSA fits file
+;tMin and tMax set the range for theta [tMin, tMax]*pi
+;Rin and Rout are inner and outer boundaries in the radial direction
+;Ni, Nj, Nk set the number of cells in r, theta, phi directions
+;Nghost is the number of ghost cells
[Gamera]
gameraGridFile = heliogrid.h5
@@ -9,28 +14,29 @@ IbcDir = ./
[Grid]
tMin = 0.1
tMax = 0.9
-Rin = 21.5
-Rout = 215.
-Ni = 128
-Nj = 64
-Nk = 128
+Rin = 21.5
+Rout = 220.
+Ni = 128
+Nj = 64
+Nk = 128
[WSA]
;wsafile is the path to the WSA fits file relative to $KAIJUHOME
-;Helio test uses Carrington Rotation 2193
+;Helio test uses WSA file for Carrington Rotation 2193, by default
wsafile = examples/helio/vel_201708132000R002_ahmi.fits
density_temperature_infile = no
-gauss_smooth_width = 0 ; 8
-plots = yes
-normalized = no
+gauss_smooth_width = 0 ; 8
+normalized = no
[Constants]
-gamma = 1.5 ;1.05
-NO2 = 4
+gamma = 1.5
+Nghost = 4
+Tsolar = 25.38
[Normalization]
-B0 = 1.e-3 ; 100 nT
-n0 = 200. ; 200/cc
+B0 = 1.e-3 ; in [Gs] equals to 100 [nT]
+n0 = 200. ; in [cm-3]
+T0 = 1.e6 ; in [K]
diff --git a/kaipy/gamhelio/helioViz.py b/kaipy/gamhelio/helioViz.py
index 6ba8215d..db718880 100644
--- a/kaipy/gamhelio/helioViz.py
+++ b/kaipy/gamhelio/helioViz.py
@@ -9,8 +9,7 @@ import kaipy.kaiViz as kv
import kaipy.gamhelio.heliosphere as hsph
import os
-#Tsolar = 25.38
-Tsolar = 1.e6
+Tsolar = 25.38
VMax = 800.
VMin = 300.
@@ -35,8 +34,8 @@ TMin = 0.2
TMax = 2.
TCM = "copper"
-T0Min = 0.05
-T0Max = 0.15
+T0Min = 0.01
+T0Max = 0.25
BMax = 150.
BMin = -150.
@@ -46,6 +45,8 @@ BCM = "coolwarm"
B0Min = -4.
B0Max = 4.
+
+colorProf = "tab:orange"
#Function to Add different size options to argument
#not used for helio right now
def AddSizeArgs(parser):
@@ -58,13 +59,15 @@ def AddSizeArgs(parser):
def GetSizeBds(pic):
if (pic == "pic1" or pic == "pic2"):
#for inner helio
- xyBds = [-216.,216.,-216.,216.]
+ xyBds = [-220.,220.,-220.,220.]
#for 1-10 au helio
#xyBds = [-10.,10.,-10.,10.]
elif (pic == "pic3"):
xyBds = [0.,360.,-75.,75.]
elif (pic == "pic4"):
xyBds = [0.,360.,-90.,90.]
+ elif (pic == "pic5"):
+ xyBds = [20.,220.,1.,2000.]
else:
print ("No pic type specified.")
return xyBds
@@ -416,6 +419,53 @@ def PlotiSlTemp(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
Ax.set_ylabel('Latitude')
return Temp
+#Plot Density as a function of distance
+def PlotDensityProf(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
+ if (doClear):
+ Ax.clear()
+
+ D = gsph.RadProfDen(nStp)
+ rad = gsph.RadialProfileGrid()
+
+ Ax.plot(rad,D,colorProf)
+
+ if (doDeco):
+ Ax.set_xlabel('Radial distance [R_sun]')
+ Ax.set_ylabel('Density [cm-3]')
+ Ax.set_ylim(250.,450.)
+ Ax.set_xlim(20.,220.)
+ #Ax.yaxis.tick_right()
+ #Ax.yaxis.set_label_position('right')
+ return D
+
+#Plot speed as a function of distance
+def PlotSpeedProf(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
+ if (doClear):
+ Ax.clear()
+ V = gsph.RadProfSpeed(nStp)
+ rad = gsph.RadialProfileGrid()
+ Ax.plot(rad,V,colorProf)
+
+ if (doDeco):
+ Ax.set_xlabel('Radial distance [R_sun]')
+ Ax.set_ylabel('Speed [km/s]')
+ Ax.set_ylim(600.,750.)
+ Ax.set_xlim(20.,220.)
+ return V
+
+def PlotFluxProf(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
+ if (doClear):
+ Ax.clear()
+ F = gsph.RadProfFlux(nStp)
+ rad = gsph.RadialProfileGrid()
+ Ax.plot(rad,F,colorProf)
+
+ if (doDeco):
+ Ax.set_xlabel('Radial distance [R_sun]')
+ Ax.set_ylabel('RhoVr^2')
+ Ax.set_ylim(180000.,280000.)
+ Ax.set_xlim(20.,220.)
+ return F
#Adds MPI contours
#this function is from magnetosphere Viz script. PlotMPI is not used for helio as of now
diff --git a/kaipy/gamhelio/heliosphere.py b/kaipy/gamhelio/heliosphere.py
index 550cefd9..ff2d6ec8 100644
--- a/kaipy/gamhelio/heliosphere.py
+++ b/kaipy/gamhelio/heliosphere.py
@@ -9,12 +9,13 @@ import timeit
#Object to pull from MPI/Serial heliosphere runs (H5 data), extends base
-ffam = "monospace"
-dLabC = "black" #Default label color
+ffam = "monospace"
+dLabC = "black" #Default label color
dLabFS = "medium" #Default label size
-dBoxC = "lightgrey" #Default box color
-TINY = 1.0e-8
-rmStr = "mixtest"
+dBoxC = "lightgrey" #Default box color
+TINY = 1.0e-8
+rmStr = "mixtest"
+MK = 1.e6 #MegaKelvin
#Adapted to helio grid
class GamsphPipe(GameraPipe):
@@ -28,9 +29,9 @@ class GamsphPipe(GameraPipe):
self.vScl = 150. #-> km/s
self.tScl = 4637. #->seconds
self.dScl = 200. #cm-3
- self.TScl = 1.e-6/4/np.pi/200./1.38e-16/1.e6 #in MK
+ self.TScl = 1.e-6/4/np.pi/200./kbltz/MK #in MK
- # [OHelio]
+ # units for OHelio
#self.bScl = 5. #->nT
#self.vScl = 34.5 #-> km/s
#self.tScl = 1.4e8/34.5
@@ -46,23 +47,25 @@ class GamsphPipe(GameraPipe):
#inner boundary distance
self.R0 = self.xxc[0,0]
-
+
+ #j and k for radial profile
+ self.jRad = self.Nj//2
+ self.kRad = self.Nk//4
+
def OpenPipe(self,doVerbose=True):
GameraPipe.OpenPipe(self,doVerbose)
if (self.UnitsID != "CODE"):
- self.bScl = 1.0 #->nT
- self.vScl = 1.0 #-> km/s
- self.tScl = 1.0 #->Seconds
- # [EP] added
- self.dScl = 1.0
- self.TScl = 1.0
+ self.bScl = 1.0 #->nT
+ self.vScl = 1.0 #-> km/s
+ self.tScl = 1.0 #-> Seconds
+ self.dScl = 1.0 #-> cm-3
+ self.TScl = 1.0/kbltz/MK #-> MKelvin
#Rescale time
self.T = self.tScl*self.T
Neq_a = self.Nj//2 #cell above eq plane
- print (Neq_a)
Nr = self.Ni
Np = self.Nk
@@ -109,6 +112,30 @@ class GamsphPipe(GameraPipe):
Qj[:,:] = 0.5*( Q[:,ja,:] + Q[:,jb,:] )
return Qj
+
+ def RadialProfileGrid(self):
+ self.GetGrid(doVerbose=True)
+ #set j and k for a radial direction
+ #Nk2 = self.Nk//2
+ #Nj2 = self.Nj//2
+ #cell corners
+ x = self.X [:,:,:]
+ y = self.Y [:,:,:]
+ z = self.Z [:,:,:]
+ #cell centers
+ x_c = 0.125*(x[:-1,:-1,:-1]+x[:-1,:-1,1:]+x[:-1,1:,:-1]+x[:-1,1:,1:]+
+ x[1:,:-1,:-1]+x[1:,:-1,1:]+x[1:,1:,:-1]+x[1:,1:,1:])
+ y_c = 0.125*(y[:-1,:-1,:-1]+y[:-1,:-1,1:]+y[:-1,1:,:-1]+y[:-1,1:,1:]+
+ y[1:,:-1,:-1]+y[1:,:-1,1:]+y[1:,1:,:-1]+y[1:,1:,1:])
+ z_c = 0.125*(z[:-1,:-1,:-1]+z[:-1,:-1,1:]+z[:-1,1:,:-1]+z[:-1,1:,1:]+
+ z[1:,:-1,:-1]+z[1:,:-1,1:]+z[1:,1:,:-1]+z[1:,1:,1:])
+ #radius of cell centers
+ jR = self.jRad
+ kR = self.kRad
+ r = np.sqrt(x_c[:,jR,kR]**2.0 + y_c[:,jR,kR]**2.0 + z_c[:,jR,kR]**2.)
+
+ return r
+
#NOT USED merid plane Y=0
def MeridGrid(self):
#Get Grid
@@ -212,6 +239,83 @@ class GamsphPipe(GameraPipe):
#print ('jd_c = ', jd_c)
return Qi
+ #Var along 1D radial line
+ def RadialProfileVar(self,vID,sID=None,vScl=None,doVerb=True):
+ #Get full 3D variable first
+ Q = self.GetVar(vID,sID,vScl,doVerb)
+
+ #set j and k for a radial profile
+ jR = self.jRad
+ kR = self.kRad
+ Nr = self.Ni
+
+ Qi = np.zeros(Nr)
+ #variable in a cell center
+ Qi[:] = Q[:,jR,kR]
+
+ return Qi
+
+ #Radial Profile: Normalized Density
+ def RadProfDen(self,s0=0):
+ D = self.RadialProfileVar("D", s0)
+
+ self.GetGrid(doVerbose=True)
+ jR = self.jRad
+ kR = self.kRad
+ x = self.X [:,:,:]
+ y = self.Y [:,:,:]
+ z = self.Z [:,:,:]
+
+ x_c = 0.125*(x[:-1,:-1,:-1]+x[:-1,:-1,1:]+x[:-1,1:,:-1]+x[:-1,1:,1:]+
+ x[1:,:-1,:-1]+x[1:,:-1,1:]+x[1:,1:,:-1]+x[1:,1:,1:])
+ y_c = 0.125*(y[:-1,:-1,:-1]+y[:-1,:-1,1:]+y[:-1,1:,:-1]+y[:-1,1:,1:]+
+ y[1:,:-1,:-1]+y[1:,:-1,1:]+y[1:,1:,:-1]+y[1:,1:,1:])
+ z_c = 0.125*(z[:-1,:-1,:-1]+z[:-1,:-1,1:]+z[:-1,1:,:-1]+z[:-1,1:,1:]+
+ z[1:,:-1,:-1]+z[1:,:-1,1:]+z[1:,1:,:-1]+z[1:,1:,1:])
+
+ r = np.sqrt(x_c[:,jR,kR]**2.0 + y_c[:,jR,kR]**2.0+z_c[:,jR,kR]**2.)
+ Norm = r**2./r[0]/r[0]
+
+ D = D*Norm*self.dScl
+ return D
+
+ #Radial Profile: Speed
+ def RadProfSpeed(self,s0=0):
+ Vx = self.RadialProfileVar("Vx", s0)
+ Vy = self.RadialProfileVar("Vy", s0)
+ Vz = self.RadialProfileVar("Vz", s0)
+
+ MagV = self.vScl*np.sqrt(Vx**2.0+Vy**2.0+Vz**2.0)
+ return MagV
+
+ #Radial Profile: Normalized Flux rho*V*r^2
+ def RadProfFlux(self,s0=0):
+ D = self.RadialProfileVar("D", s0)
+ Vx = self.RadialProfileVar("Vx", s0)
+ Vy = self.RadialProfileVar("Vy", s0)
+ Vz = self.RadialProfileVar("Vz", s0)
+
+ self.GetGrid(doVerbose=True)
+ jR = self.jRad
+ kR = self.kRad
+ x = self.X [:,:,:]
+ y = self.Y [:,:,:]
+ z = self.Z [:,:,:]
+
+ x_c = 0.125*(x[:-1,:-1,:-1]+x[:-1,:-1,1:]+x[:-1,1:,:-1]+x[:-1,1:,1:]+
+ x[1:,:-1,:-1]+x[1:,:-1,1:]+x[1:,1:,:-1]+x[1:,1:,1:])
+ y_c = 0.125*(y[:-1,:-1,:-1]+y[:-1,:-1,1:]+y[:-1,1:,:-1]+y[:-1,1:,1:]+
+ y[1:,:-1,:-1]+y[1:,:-1,1:]+y[1:,1:,:-1]+y[1:,1:,1:])
+ z_c = 0.125*(z[:-1,:-1,:-1]+z[:-1,:-1,1:]+z[:-1,1:,:-1]+z[:-1,1:,1:]+
+ z[1:,:-1,:-1]+z[1:,:-1,1:]+z[1:,1:,:-1]+z[1:,1:,1:])
+
+ r = np.sqrt(x_c[:,jR,kR]**2.0 + y_c[:,jR,kR]**2.0+z_c[:,jR,kR]**2.)
+
+ Norm = r[:]**2./r[0]/r[0]
+
+ Flux = D*Norm*self.dScl*self.vScl*np.sqrt(Vx**2.0+Vy**2.0+Vz**2.0)
+ return Flux
+
#Speed at 1 AU
def iSliceMagV(self,s0=0):
Vx = self.iSliceVar("Vx",s0) #Unscaled
diff --git a/kaipy/gamhelio/wsa2gamera/params.py b/kaipy/gamhelio/wsa2gamera/params.py
index 6d997488..ef39f23d 100644
--- a/kaipy/gamhelio/wsa2gamera/params.py
+++ b/kaipy/gamhelio/wsa2gamera/params.py
@@ -12,15 +12,17 @@ class params():
self.wsaFile = config['WSA']['wsafile']
self.gaussSmoothWidth = config.getint('WSA','gauss_smooth_width')
- self.plots = config.getboolean('WSA','plots')
+ #self.plots = config.getboolean('WSA','plots')
self.densTempInfile = config.getboolean('WSA','density_temperature_infile')
self.normalized = config.getboolean('WSA','normalized')
self.gamma = config.getfloat('Constants','gamma')
- self.NO2 = config.getint('Constants','NO2')
+ self.Nghost = config.getint('Constants','Nghost')
+ self.Tsolar = config.getfloat('Constants','Tsolar')
self.B0 = config.getfloat('Normalization','B0')
self.n0 = config.getfloat('Normalization','n0')
+ self.T0 = config.getfloat('Normalization','T0')
self.tMin = config.getfloat('Grid','tMin')
self.tMax = config.getfloat('Grid','tMax')
diff --git a/kaipy/kdefs.py b/kaipy/kdefs.py
index 2a5e0bd3..89af4ec6 100644
--- a/kaipy/kdefs.py
+++ b/kaipy/kdefs.py
@@ -2,3 +2,4 @@
import numpy as np
EarthM0g = 0.31 #Gauss
+kbltz = 1.38e-16 #Boltzmann constant [erg/K]
diff --git a/kaipy/satcomp/scutils.py b/kaipy/satcomp/scutils.py
index 72a7ce0c..5aa9ab5f 100644
--- a/kaipy/satcomp/scutils.py
+++ b/kaipy/satcomp/scutils.py
@@ -495,13 +495,21 @@ def createHelioInputFiles(data,scDic,scId,mjd0,sec0,fdir,ftag,numSegments):
else:
print('Coordinate system transformation failed')
return
- elapsedSecs = [(tt - t[0]).seconds for tt in t]
+ #
+ from math import sqrt, pi
+ L0 = 6.955e10
+ Mp = 1.67e-24
+ in2cms = 1e-3/sqrt(4*pi*200*Mp)
+ in2s = L0/in2cms
+ elapsed = [(tt - t[0]).seconds/in2s for tt in t]
+ #
scTrackName = os.path.join(fdir,scId+".sc.h5")
with h5py.File(scTrackName,'w') as hf:
- hf.create_dataset("T" ,data=elapsedSecs)
- hf.create_dataset("X" ,data=x)
- hf.create_dataset("Y" ,data=y)
- hf.create_dataset("Z" ,data=z)
+ hf.create_dataset("T" ,data=elapsed)
+ # Reverse x for gamhelio frame.
+ hf.create_dataset("X" ,data=-c.cartesian.x)
+ hf.create_dataset("Y" ,data=c.cartesian.y)
+ hf.create_dataset("Z" ,data=c.cartesian.z)
if scId in ["ACE"]:
chimpxml = genHelioSCXML(fdir,ftag,
scid=scId,h5traj=os.path.basename(scTrackName),numSegments=0)
diff --git a/scripts/datamodel/helioSatComp.py b/scripts/datamodel/helioSatComp.py
index bbfd43dc..eb1724e8 100755
--- a/scripts/datamodel/helioSatComp.py
+++ b/scripts/datamodel/helioSatComp.py
@@ -160,16 +160,16 @@ if __name__ == "__main__":
# Use the first (non zero?) time as the inital MJD.
# N.B. THIS SKIPS THE FIRST SIMULATION STEP SINCE IT TYPICALLY HAS
# gamT[0] = 0.
- # loc = np.argwhere(gamT > 0.0)[0][0]
- # t0 = gamUT[loc] # First non-0 time
- t0 = gamUT[0]
+ loc = np.argwhere(gamT > 0.0)[0][0]
+ t0 = gamUT[loc] # First non-0 time
+ # t0 = gamUT[0]
t1 = gamUT[-1]
if debug:
print("t0 = %s" % t0)
print("t1 = %s" % t1)
# Compute the time interval (in integer code units) between step 1 and step 2.
- # deltaT = np.round(gamT[loc + 1] - gamT[loc])
+ deltaT = np.round(gamT[loc + 1] - gamT[loc])
# if debug:
# print("deltaT = %s" % deltaT)
#
@@ -177,14 +177,14 @@ if __name__ == "__main__":
#
# Save the (float) MJD of the first step.
- # mjdFileStart = gamMJD[loc]
- mjdFileStart = gamMJD[0]
+ mjdFileStart = gamMJD[loc]
+ # mjdFileStart = gamMJD[0]
if debug:
print("mjdFileStart = %s" % mjdFileStart)
# Save the elapsed simulation time in code units of the first step.
- # secFileStart = gamT[loc]
- secFileStart = gamT[0]
+ secFileStart = gamT[loc]
+ # secFileStart = gamT[0]
if debug:
print("secFileStart = %s" % secFileStart)
diff --git a/scripts/datamodel/sunpy.ipynb b/scripts/datamodel/sunpy.ipynb
new file mode 100644
index 00000000..a4954cea
--- /dev/null
+++ b/scripts/datamodel/sunpy.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a442a301447dc28f4f41d6f88accb9bd3a5ac6f3340dcf4298832117ec05ec8b
+size 7775
diff --git a/scripts/preproc/wsa2TDgamera.py b/scripts/preproc/wsa2TDgamera.py
index b59fb1c1..5f79414f 100755
--- a/scripts/preproc/wsa2TDgamera.py
+++ b/scripts/preproc/wsa2TDgamera.py
@@ -270,6 +270,7 @@ with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
#not interpolating temperature, calculating sound speed cs
#assuming uniform total pressure Rho_max*k*T0 = p+Br^2/8pi
+ #TODO: Check Temp calculation
T0 = 0.9e6
Rho0 = 1100.*mp #density in the HCS
#cs = np.sqrt(prm.gamma/rho*(rho.max()*1.38e-16*T0/1.67e-24-br**2/8/np.pi))
diff --git a/scripts/preproc/wsa2gamera.py b/scripts/preproc/wsa2gamera.py
index 4a3ddb12..ca9c7aa4 100755
--- a/scripts/preproc/wsa2gamera.py
+++ b/scripts/preproc/wsa2gamera.py
@@ -12,22 +12,33 @@ import kaipy.gamhelio.lib.wsa as wsa
import kaipy.gamera.gamGrids as gg
-#----------- PARSE ARGUMENTS ---------#
+# Parse arguments
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
-#----------- PARSE ARGUMENTS ---------#
+
+# constants
+# TODO: move to kaipy/kdefs.py
+mp = 1.67e-24
+kblts = 1.38e-16
+nCS = 1100. #for T calculation from pressure balance
# Read params from config file
prm = params.params(args.ConfigFileName)
-Ng=prm.NO2
+Ng=prm.Nghost
gamma = prm.gamma
+Tsolar = prm.Tsolar
+
+# Normalization parameters
+# remember to use the same units in gamera
B0 = prm.B0
n0 = prm.n0
-T0 = 3.44e6 #2.88e6
+V0 = B0/np.sqrt(4*np.pi*mp*n0)
+#Temperature in the Current Sheet (for T calculation from pressure balance)
+TCS = prm.T0
-#grid parameters
+# Grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
@@ -38,51 +49,31 @@ Nk = prm.Nk
ffits = os.path.join(os.getenv('KAIJUHOME'),prm.wsaFile)
-# constants
-mp = 1.67e-24
-
-#----------GENERATE HELIO GRID------
-
-print("Generating gamera-helio grid ...")
+# Generate spherical helio grid
+print("Generating gamera-helio grid Ni = %d, Nj = %d, Nk = %d " % (Ni, Nj, Nk))
X3,Y3,Z3 = gg.GenKSph(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
-
-#to generate non-uniform grid for GL cme (more fine in region 0.1-0.3 AU)
-#X3,Y3,Z3 = gg.GenKSphNonUGL(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
-print("Gamera-helio grid ready!")
-
-#----------GENERATE HELIO GRID------
+if os.path.exists(prm.gameraGridFile):
+ print("Grid file heliogrid.h5 is ready!")
-############### WSA STUFF #####################
+# Read and normalize 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)
-
-# convert the units; remember to use the same units in gamera
-# TODO: probably store units in the h5 file?
-# B0 = 1.e-3 Gs
-# n0 = 200./cc
-
-V0 = B0/np.sqrt(4*np.pi*mp*n0)
-
bi_wsa /= B0
n_wsa /= (mp*n0)
v_wsa /= V0
-#convert julian date from wsa fits into modified julian date
+#convert julian date in the center of the WSA map into modified julian date
mjd_c = jd_c - 2400000.5
-# keep temperature in K
-############### WSA STUFF #####################
-############### GAMERA STUFF #####################
-
-# GAMERA GRID
+# Get GAMERA grid for further interpolation
with h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile),'r') as f:
x=f['X'][:]
y=f['Y'][:]
z=f['Z'][:]
-
+# Cell centers, note order of indexes [k,j,i]
xc = 0.125*(x[:-1,:-1,:-1]+x[:-1,1:,:-1]+x[:-1,:-1,1:]+x[:-1,1:,1:]
+x[1:,:-1,:-1]+x[1:,1:,:-1]+x[1:,:-1,1:]+x[1:,1:,1:])
yc = 0.125*(y[:-1,:-1,:-1]+y[:-1,1:,:-1]+y[:-1,:-1,1:]+y[:-1,1:,1:]
@@ -90,15 +81,17 @@ yc = 0.125*(y[:-1,:-1,:-1]+y[:-1,1:,:-1]+y[:-1,:-1,1:]+y[:-1,1:,1:]
zc = 0.125*(z[:-1,:-1,:-1]+z[:-1,1:,:-1]+z[:-1,:-1,1:]+z[:-1,1:,1:]
+z[1:,:-1,:-1]+z[1:,1:,:-1]+z[1:,:-1,1:]+z[1:,1:,1:])
-# remove the ghosts from angular dimensions
-R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2) # radius of the inner boundary
+# radius of the inner boundary
+R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**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[P<0]=P[P<0]+2*np.pi
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.
+# 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)
Pc = np.arctan2(yc[Ng:-Ng,Ng:-Ng,:],xc[Ng:-Ng,Ng:-Ng,:])
Pc[Pc<0]=Pc[Pc<0]+2*np.pi
@@ -108,7 +101,7 @@ Tc = np.arccos(zc[Ng:-Ng,Ng:-Ng,:]/Rc)
fbi = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,bi_wsa.T,kx=1,ky=1)
br = fbi(Pc[:,0,0],Tc[0,:,0])
-############### SMOOTHING #####################
+# Smoothing
if not prm.gaussSmoothWidth==0:
import astropy
from astropy.convolution import convolve,Gaussian2DKernel
@@ -117,21 +110,22 @@ if not prm.gaussSmoothWidth==0:
br =astropy.convolution.convolve(br,gauss,boundary='extend')
-############### INTERPOLATE AND DUMP #####################
+# Interpolate to Gamera grid
fv = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,v_wsa.T,kx=1,ky=1)
vr = fv(Pc[:,0,0],Tc[0,:,0])
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa.T,kx=1,ky=1)
rho = f(Pc[:,0,0],Tc[0,:,0])
-#f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa.T,kx=1,ky=1)
-#temp = f(Pc[:,0,0],Tc[0,:,0])
-temp = 1.*T0/rho + (1.**2-(br)**2)*V0**2 / 2e8/1.38 * 1.67/rho # *****
-temp_T = temp.T
+# 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*kblts*TCS - (br*B0)**2/8./np.pi)/(rho*n0)/kblts
+# note, keep temperature in K (pressure is normalized in wsa.F90)
-pressure = ((br)**2)*V0**2 /2.*mp*n0 *0.1 + (n0*rho* temp)*1.38e-16 *0.1
-pressure_therm = (n0*rho* temp)*1.38e-16 * 0.1
-pressure_B = (br)**2 *V0**2 / 2.*mp*n0 *0.1
+#check
+#print ("Max and min of temperature in MK")
+#print (np.amax(temp)*1.e-6, np.amin(temp)*1.e-6)
# note, redefining interpolation functions we could also
# interpolate from bi_wsa as above, but then we would have to
@@ -145,11 +139,22 @@ br_kface = fbi(P[:,0,0],Tc[0,:,0])
vr_kface = fv (P[:,0,0],Tc[0,:,0])
# Scale inside ghost region
-(vr,vr_kface,rho,temp,br,br_kface) = [np.dstack(prm.NO2*[var]) for var in (vr,vr_kface,rho,temp,br,br_kface)]
+(vr,vr_kface,rho,temp,br,br_kface) = [np.dstack(Ng*[var]) for var in (vr,vr_kface,rho,temp,br,br_kface)]
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]
+
+#check
+#print ("Br kface ", br_kface.shape)
+#print ("E theta ", et_kedge.shape)
+
+# Save to file
+# 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)
@@ -158,4 +163,5 @@ with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
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.close()
diff --git a/scripts/quicklook/heliopic.py b/scripts/quicklook/heliopic.py
index 6735280b..1a97cd20 100755
--- a/scripts/quicklook/heliopic.py
+++ b/scripts/quicklook/heliopic.py
@@ -63,8 +63,10 @@ if __name__ == "__main__":
figSz = (10,12.5)
elif (pic == "pic3"):
figSz = (10,6.5)
- else:
+ elif (pic == "pic4"):
figSz = (10,6.)
+ elif (pic == "pic5"):
+ figSz = (12.,12.)
#======
#Init data
gsph = hsph.GamsphPipe(fdir,ftag,doFast=doFast)
@@ -77,7 +79,7 @@ if __name__ == "__main__":
#Setup figure
fig = plt.figure(figsize=figSz)
- if (pic != "pic4"):
+ if (pic == "pic1" or pic == "pic2" or pic == "pic3"):
gs = gridspec.GridSpec(4,6,height_ratios=[20,1,20,1])
#plots. Two rows of two plots
AxL0 = fig.add_subplot(gs[0,0:3])
@@ -92,11 +94,15 @@ if __name__ == "__main__":
AxC1_1 = fig.add_subplot(gs[3,0:3])
AxC2_1 = fig.add_subplot(gs[3,3:])
- else:
+ elif (pic == "pic4"):
gs = gridspec.GridSpec(2,1,height_ratios=[20,1])
Ax = fig.add_subplot(gs[0,0])
AxC = fig.add_subplot(gs[1,0])
-
+ elif (pic == "pic5"):
+ gs = gridspec.GridSpec(2,2)
+ Ax = fig.add_subplot(gs[0,0])
+ AxC = fig.add_subplot(gs[0,1])
+ AxC1 = fig.add_subplot(gs[1,0])
if (pic == "pic1"):
hviz.PlotEqMagV(gsph,nStp,xyBds,AxL0,AxC1_0)
@@ -118,6 +124,10 @@ if __name__ == "__main__":
hviz.PlotiSlBr(gsph,nStp,xyBds,AxR1,AxC2_1)
elif (pic == "pic4"):
hviz.PlotiSlBrRotatingFrame(gsph,nStp,xyBds,Ax,AxC)
+ elif (pic == "pic5"):
+ hviz.PlotDensityProf(gsph,nStp,xyBds,Ax)
+ hviz.PlotSpeedProf(gsph,nStp,xyBds,AxC)
+ hviz.PlotFluxProf(gsph,nStp,xyBds,AxC1)
else:
print ("Pic is empty. Choose pic1 or pic2 or pic3")
@@ -126,7 +136,7 @@ if __name__ == "__main__":
gsph.AddTime(nStp,AxL0,xy=[0.025,0.875],fs="x-large")
elif (pic == "pic3"):
gsph.AddTime(nStp,AxL0,xy=[0.015,0.82],fs="small")
- elif (pic == "pic4"):
+ elif (pic == "pic4" or pic == "pic5"):
gsph.AddTime(nStp,Ax,xy=[0.015,0.92],fs="small")
else:
print ("Pic is empty. Choose pic1 or pic2 or pic3")
diff --git a/src/base/kdefs.F90 b/src/base/kdefs.F90
index 11141620..719e4087 100644
--- a/src/base/kdefs.F90
+++ b/src/base/kdefs.F90
@@ -34,9 +34,11 @@ module kdefs
real(rp), parameter :: Re_cgs = 6.3781D8 ![cm] Earth's radius
real(rp), parameter :: Me_cgs = 9.1093837015D-28 ![g] Electron mass
real(rp), parameter :: Mp_cgs = 1.67262192369D-24 ![g] Proton mass
+ real(rp), parameter :: G_cgs = 6.6726D-8 ![cm^3/g/s^2], Gravitational constant (per NRL plasma formulary'21)
!MKS Constants
real(rp), parameter :: vc_mks = vc_cgs*(1.0e-2) ![m/s], Speed of light
+ real(rp), parameter :: G_mks = 6.6726D-11 ![m^3/kg/s^2], Gravitational constant (per NRL plasma formulary'21)
!Helper conversions
real(rp), parameter :: G2nT = 1.0E+5 !Gauss->nT
@@ -80,7 +82,8 @@ module kdefs
real(rp), parameter :: RNeptuneXE = 3.860 !Rx = X*Re
!Helio constants
- real(rp), parameter :: Rsolar = 6.956D5 ! [km] Solar radius
+ real(rp), parameter :: Rsolar = 6.956D5 ! [km] Solar radius
+ real(rp), parameter :: Msolar = 1.98847D30 ! [kg] Solar mass
!Numbered accessors
!Directions
diff --git a/src/chimp/chmpunits.F90 b/src/chimp/chmpunits.F90
index aed50d78..9d33eac4 100644
--- a/src/chimp/chmpunits.F90
+++ b/src/chimp/chmpunits.F90
@@ -212,11 +212,11 @@ module chmpunits
!Field: 100 nT
!Gamera units for heliosphere runs
L0 = 6.955e+10 !Rs in cm
- in2cms = 1.0e-3/sqrt(4*PI*200*Mp_cgs) !150e+5 cm/s
- in2G = 1.0e-3 !in [G]
- in2s = L0/in2cms ! time in s
- M0g = 0.0
- inPScl = 1.0e-6*1.0e+8/4/pi !Pressure unit B[G]^2/4pi *1.e8 in [nPa]
+ in2cms = 1.0e+5 ! km/s -> cm/s
+ in2G = 1.0e-5 ! nT -> Gs
+ in2s = 1.0 ! already in s
+ M0g = 0.0
+ inPScl = 1.0e+8 !erg/cm3 -> [nPa]
rClosed = 21.5 !Radius of inner boundary in units of grid length
case("LFM")
L0 = Re_cgs !Using scaled grid
diff --git a/src/gamera/ICs/helio/wsa.F90 b/src/gamera/ICs/helio/wsa.F90
index 7082efe9..a047f7cf 100644
--- a/src/gamera/ICs/helio/wsa.F90
+++ b/src/gamera/ICs/helio/wsa.F90
@@ -125,8 +125,6 @@ module usergamic
Grid%keDT = Grid%ke
! Add gravity
-! tsHack => PerStep
-! Model%HackStep => tsHack
! eHack => EFix
! Model%HackE => eHack
@@ -207,9 +205,6 @@ module usergamic
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,k,jg,kg,ke,kb,a,var,xyz,R,Theta,Phi,rHat,phiHat) &
!$OMP private(ibcVarsStatic,pVar,conVar,xyz0,R_kf,Theta_kf)
-
-
-
do k=Grid%ksg,Grid%keg+1 ! note, going all the way to last face for mag fluxes
kg = k+Grid%ijkShift(KDIR)
! map rotating to static grid
@@ -301,51 +296,6 @@ module usergamic
end subroutine wsaBC
-
- ! !Do per timestep, includes lazy gravitational force term
- ! subroutine PerStep(Model,Gr,State)
- ! type(Model_T), intent(in) :: Model
- ! type(Grid_T), intent(inout) :: Gr
- ! type(State_T), intent(inout) :: State
-
- ! integer :: i,j,k
-
- ! real(rp), dimension(NDIM) :: xyz, Vxyz, rHat
- ! real(rp), dimension(NVAR) :: pW,pCon
- ! real(rp) :: D,IntE,r
- ! real(rp) :: GM0
-
- ! !Scaling for gravitational force
- ! GM0 = UN/(UL**3*UB**2)*6.67408*1.99*4*pi*1.67/6.955/10 ! 2.74e4cm/s^2
-
- ! !Add grav force
- ! !$OMP PARALLEL DO default(shared) &
- ! !$OMP private(i,j,k,xyz,rHat,Vxyz,pW,pCon,r,D,IntE)
- ! do k=Gr%ksg,Gr%keg
- ! do j=Gr%jsg,Gr%jeg
- ! do i=Gr%isg,Gr%ieg
- ! xyz = Gr%xyzcc(i,j,k,:)
- ! r = norm2(xyz)
- ! rHat = xyz/r
-
- ! pCon = State%Gas(i,j,k,:,BLK)
- ! call CellC2P(Model,pCon,pW)
-
- ! D = pW(DEN)
- ! IntE = pW(PRESSURE)/(Model%gamma-1)
- ! Vxyz = pW(VELX:VELZ)
- ! Vxyz = Vxyz - Model%dt*GM0*rHat/(r*r)
-
- ! !Reset conserved State
- ! pCon(DEN) = D
- ! pCon(MOMX:MOMZ) = D*Vxyz
- ! pCon(ENERGY) = IntE + 0.5*D*dot_product(Vxyz,Vxyz)
- ! State%Gas(i,j,k,:,BLK) = pCon
- ! enddo
- ! enddo
- ! enddo
- ! end subroutine PerStep
-
subroutine eFix(Model,Gr,State)
type(Model_T), intent(in) :: Model
type(Grid_T), intent(in) :: Gr
diff --git a/src/gamera/ICs/wsa.F90 b/src/gamera/ICs/wsa.F90
new file mode 100644
index 00000000..7d7d621b
--- /dev/null
+++ b/src/gamera/ICs/wsa.F90
@@ -0,0 +1,417 @@
+module usergamic
+ use kdefs
+ use gamtypes
+ use gambctypes
+ use gamutils
+ use math
+ use gridutils
+ use xml_input
+ use bcs
+ use ioH5
+ use helioutils
+
+
+ implicit none
+
+ enum, bind(C)
+ ! variables passed via innerbc file
+ ! Br, Vr, Rho, Temperature, Br @ kface, Vr @ kface
+ enumerator :: BRIN=1,VRIN,RHOIN,TIN,BRKFIN,VRKFIN
+ endenum
+
+
+ integer, private, parameter :: NVARSIN=6 ! SHOULD be the same as the number of vars in the above enumerator
+ real(rp), dimension(:,:,:,:), allocatable :: ibcVars
+
+ !Various global would go here
+ real (rp) :: Rho0, P0, Vslow,Vfast, wScl, Cs0, B0, MJD_c
+
+ ! things we keep reusing
+ real(rp), dimension(NDIM) :: xyz,xyz0,rHat,phiHat
+ real(rp) :: Rfactor
+
+ ! global grid
+ integer :: gNkp
+
+ ! FIXME
+ ! setting it here temporarily. eventually need to read from HDF or something
+ ! note also that this is slightly incorrect, since Rbc below is used as the radius
+ ! of the center of the first ghost cell.
+
+ !for inner heliosphere
+ real(rp) :: Rbc = 21.5
+
+ real(rp) :: Tsolar ! Solar rotation period, defined in apps/helioutils.F90
+
+ character(len=strLen) :: wsaFile
+
+ ! use this to fix the Efield at the inner boundary
+! real(rp), allocatable :: inEijk(:,:,:,:)
+
+ ! type for solar wind BC
+ type, extends(innerIBC_T) :: SWInnerBC_T
+
+ !Main electric field structures
+ real(rp), allocatable, dimension(:,:,:,:) :: inEijk,inExyz
+
+ contains
+
+! procedure :: doInit => InitIonInner
+ ! TODO: this shoudl be made generic (wsa, mas, etc.) How
+ procedure :: doBC => wsaBC
+ end type SWInnerBC_T
+
+ contains
+
+ subroutine initUser(Model,Grid,State,inpXML)
+ type(Model_T), intent(inout) :: Model
+ type(Grid_T), intent(inout) :: Grid
+ type(State_T), intent(inout) :: State
+ type(XML_Input_T), intent(in) :: inpXML
+ procedure(GasIC_T), pointer :: Wxyz
+ procedure(HackStep_T), pointer :: tsHack
+ procedure(HackE_T), pointer :: eHack
+
+ integer :: i,j,k,nvar,nr,d
+ integer :: n1, n2
+
+! if (.not.allocated(inEijk)) allocate(inEijk(1,Grid%jsg:Grid%jeg+1,Grid%ksg:Grid%keg+1,1:NDIM))
+
+ ! set units and other thins, like Tsolar
+ call setHeliosphere(Model,inpXML,Tsolar)
+
+ ! grab inner
+ call inpXML%Set_Val(wsaFile,"prob/wsaFile","innerbc.h5" )
+
+
+ ! compute global Nkp
+ gNkp = Grid%Nkp*Grid%NumRk
+
+ ! initial conditions
+ ! TODO: change using norm. units in Model, set in helioutils
+ Cs0 = 0.267 ! 40 km/s
+ Vslow = 1.33 ! 200 km/s
+ Vfast = 5.33 ! 800 km/s
+ !for inner helio
+ B0 = 2.0 ! 200 nT
+ Rho0 = 1.0 ! 200/cc
+ P0 = 1.0e-4*Rho0*Cs0**2.0/Model%gamma
+
+ ![OHelio] for 1-10 au helio
+ !Cs0 = 0.78 !27 km/s
+ !B0 = 1. ! 5nT
+ !Rho0 = 1. ! 10/cc
+ !Vslow = 8.5 !300 km/s
+
+ ! deallocate default BCs
+ ! required because defaults are triply-periodic
+ ! need to wipe
+ call WipeBCs(Model,Grid)
+
+ !Set BCs
+ allocate(SWInnerBC_T :: Grid%externalBCs(1)%p)
+ allocate(helioOuterIBC_T :: Grid%externalBCs(2)%p)
+ allocate(helioInnerJBC_T :: Grid%externalBCs(3)%p)
+ allocate(helioOuterJBC_T :: Grid%externalBCs(4)%p)
+ allocate(periodicInnerKBC_T :: Grid%externalBCs(5)%p)
+ allocate(periodicOuterKBC_T :: Grid%externalBCs(6)%p)
+
+ !Set DT bounds
+ Grid%isDT = Grid%is
+ Grid%ieDT = Grid%ie
+ Grid%jsDT = Grid%js
+ Grid%jeDT = Grid%je
+ Grid%ksDT = Grid%ks
+ Grid%keDT = Grid%ke
+
+ ! Add gravity
+! eHack => EFix
+! Model%HackE => eHack
+! tsHack => PerStep
+! Model%HackStep => tsHack
+
+ ! everybody reads WSA data
+ call readIBC(wsaFile)
+
+ Model%MJD0 = MJD_c
+
+ !Map IC to grid
+ Wxyz => GasIC
+ call GasIC2State(Model,Grid,State,Wxyz)
+
+ ! NOTE, not filling ghost cells here
+ ! relying on J, K boundary conditions
+ ! first zero everything out
+ State%magFlux = 0.0
+
+ do k=Grid%ks,Grid%ke
+ do j=Grid%js,Grid%je
+ do i=Grid%isg,Grid%ieg+1
+ ! FIXME: calc Rbc appropriately above!
+ Rfactor = Rbc/norm2(Grid%xfc(Grid%is,j,k,:,IDIR))
+ ! note scaling br to the first active face
+ ! this is opposite to what we did on the python side
+ ! not elegant, but otherwise we'd have to store both br and br_iface in teh innerbc.h5 file
+ !
+ ! 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)
+ enddo
+ enddo
+ enddo
+
+ !Local functions
+ !NOTE: Don't put BCs here as they won't be visible after the initialization call
+
+ contains
+ 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
+ real (rp) :: r, theta, phi
+ real (rp) :: r_unit(NDIM)
+
+ D = Rho0
+ P = P0
+
+ !Calculate radial hat vector
+ r = sqrt(x**2.0 + y**2.0 + z**2.0)
+ theta = acos(z/R)
+ phi = atan2(y,x)
+ r_unit = [x,y,z]/r
+
+ !Set primitives (already have D/P)
+ Vx = r_unit(XDIR)*Vslow
+ Vy = r_unit(YDIR)*Vslow
+ Vz = r_unit(ZDIR)*Vslow
+ end subroutine GasIC
+
+ end subroutine initUser
+
+ !Inner-I BC for WSA-Gamera
+ 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
+
+ ! local variables
+ integer :: i,j,k, kb, ke
+ integer :: kg, jg, ig ! global indices (in preparation for MPI)
+ 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
+
+ !i-boundaries (IN)
+ !$OMP PARALLEL DO default(shared) &
+ !$OMP private(i,j,k,jg,kg,ke,kb,a,var,xyz,R,Theta,Phi,rHat,phiHat) &
+ !$OMP private(ibcVarsStatic,pVar,conVar,xyz0,R_kf,Theta_kf)
+
+
+
+ do k=Grid%ksg,Grid%keg+1 ! note, going all the way to last face for mag fluxes
+ kg = k+Grid%ijkShift(KDIR)
+ ! map rotating to static grid
+ call mapK(kg,ke,kb,a)
+
+ do j=Grid%js,Grid%je+1 ! note, going all the way to last face for mag fluxes
+ jg = j+Grid%ijkShift(JDIR)
+
+ do ig=1,Model%Ng
+ i=ig-Model%Ng
+ ! note, ibcVars are not defined in the global j-corners or k-corners
+ ! access to k-corners is fine because mapK function will always map into active domain (FIXME: check for k=nk+1 face!)
+ ! access to j-corners gets into unallocated space for ibcVars. Use this trick instead:
+ ! set everything arbitrarily to 1. in the global corners (on low and high j boundaries)
+ ! we then apply the j-boundary after i boundary anyway, so the corners will be overwritten
+ ibcVarsStatic = 1._rp
+
+ ! otherwise
+ ! interpolate linearly from rotating to inertial frame
+ if ( (jg>=Grid%js).and.(jg<=size(ibcVars,2)) ) then
+ do var=1,NVARSIN
+ ibcVarsStatic(var) = a*ibcVars(ig,jg,kb,var)+(1-a)*ibcVars(ig,jg,ke,var)
+ end do
+ end if
+
+ ! do cell centered things for cell-centers only
+ if ( (j/=Grid%jeg+1).and.(k/=Grid%keg+1) ) then
+ ! various geometrical quantities for the cell center
+ xyz = Grid%xyzcc(i,j,k,:)
+ R = norm2(xyz)
+ Theta = acos(xyz(3)/R)
+ Phi = atan2(xyz(2),xyz(1))
+ rHat = xyz/R
+ 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)
+ ! note conversion to my units with B0^2/4pi in the denominator
+ pVar(PRESSURE) = ibcVarsStatic(RHOIN)*Model%Units%gD0*Kbltz*ibcVarsStatic(TIN)/(Model%Units%gP0)
+ pVar(DEN) = ibcVarsStatic(RHOIN)
+
+ !Swap prim->con in ghost variables
+ call CellP2C(Model,pVar,conVar)
+ State%Gas(i,j,k,:,BLK) = conVar
+
+ ! note, don't need 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,
+ ! 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)
+
+ ! 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)
+ end do
+ end do
+ end do
+
+ contains
+ 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
+
+ real(rp) :: kprime
+
+ kprime = modulo(k-gNkp*State%time/Tsolar,real(gNkp,kind=rp))
+
+ if ((kprime.ge.0).and.(kprime.lt.1)) then
+ kb=gNkp
+ ke=1
+ else
+ kb=floor(kprime)
+ ke=ceiling(kprime)
+ endif
+ a = ke-kprime
+
+ end subroutine mapK
+
+ end subroutine wsaBC
+
+
+ ! !Do per timestep, includes lazy gravitational force term
+ ! subroutine PerStep(Model,Gr,State)
+ ! type(Model_T), intent(in) :: Model
+ ! type(Grid_T), intent(inout) :: Gr
+ ! type(State_T), intent(inout) :: State
+
+ ! integer :: i,j,k
+
+ ! real(rp), dimension(NDIM) :: xyz, Vxyz, rHat
+ ! real(rp), dimension(NVAR) :: pW,pCon
+ ! real(rp) :: D,IntE,r
+ ! real(rp) :: GM0
+
+ ! !Scaling for gravitational force
+ ! GM0 = UN/(UL**3*UB**2)*6.67408*1.99*4*pi*1.67/6.955/10 ! 2.74e4cm/s^2
+
+ ! !Add grav force
+ ! !$OMP PARALLEL DO default(shared) &
+ ! !$OMP private(i,j,k,xyz,rHat,Vxyz,pW,pCon,r,D,IntE)
+ ! do k=Gr%ksg,Gr%keg
+ ! do j=Gr%jsg,Gr%jeg
+ ! do i=Gr%isg,Gr%ieg
+ ! xyz = Gr%xyzcc(i,j,k,:)
+ ! r = norm2(xyz)
+ ! rHat = xyz/r
+
+ ! pCon = State%Gas(i,j,k,:,BLK)
+ ! call CellC2P(Model,pCon,pW)
+
+ ! D = pW(DEN)
+ ! IntE = pW(PRESSURE)/(Model%gamma-1)
+ ! Vxyz = pW(VELX:VELZ)
+ ! Vxyz = Vxyz - Model%dt*GM0*rHat/(r*r)
+
+ ! !Reset conserved State
+ ! pCon(DEN) = D
+ ! pCon(MOMX:MOMZ) = D*Vxyz
+ ! pCon(ENERGY) = IntE + 0.5*D*dot_product(Vxyz,Vxyz)
+ ! State%Gas(i,j,k,:,BLK) = pCon
+ ! enddo
+ ! enddo
+ ! enddo
+ ! end subroutine PerStep
+
+ subroutine eFix(Model,Gr,State)
+ type(Model_T), intent(in) :: Model
+ type(Grid_T), intent(in) :: 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)
+
+
+ end subroutine eFix
+
+ subroutine readIBC(ibcH5)
+ character(len=*), intent(in) :: ibcH5
+ logical :: fExist
+ integer :: i,nvar,dims(3)
+ integer, parameter :: MAXIOVAR = 50
+ type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
+
+
+ !Reset IO chain
+ call ClearIO(IOVars)
+
+ inquire(file=ibcH5,exist=fExist)
+ if (.not. fExist) then
+ !Error out and leave
+ write(*,*) 'Unable to open innerbc file, exiting'
+ stop
+ 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,"MJD")
+
+ call ReadVars(IOVars,.false.,ibcH5) !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
+ if (.not.allocated(ibcVars)) allocate(ibcVars(dims(1),dims(2),dims(3),NVARSIN))
+
+ do i=1,NVARSIN
+ select case (i)
+ case (BRIN)
+ nvar= FindIO(IOVars,"br")
+ case (VRIN)
+ nvar= FindIO(IOVars,"vr")
+ case (RHOIN)
+ nvar= FindIO(IOVars,"rho")
+ case (TIN)
+ nvar= FindIO(IOVars,"temp")
+ case (BRKFIN)
+ nvar= FindIO(IOVars,"br_kface")
+ case (VRKFIN)
+ nvar= FindIO(IOVars,"vr_kface")
+ end select
+
+ ibcVars(:,:,:,i) = reshape(IOVars(nvar)%data,dims)
+ end do
+ !reading modified julian date from innerbc
+ MJD_c = GetIOReal(IOVars,"MJD")
+
+ end subroutine readIBC
+end module usergamic
diff --git a/src/gamera/apps/helioutils.F90 b/src/gamera/apps/helioutils.F90
index 02d77cc8..25636fc7 100644
--- a/src/gamera/apps/helioutils.F90
+++ b/src/gamera/apps/helioutils.F90
@@ -7,6 +7,7 @@ module helioutils
use math
use gridutils
use output
+ use ioclock
implicit none
@@ -35,18 +36,21 @@ module helioutils
type(XML_Input_T), intent(in) :: inpXML
real(rp),intent(out) :: Tsolar ! Solar rotation period
+ type(IOClock_T) :: clockScl
+
! normalization
gD0=200. ! [/cc]
- !gD0=10. ! [/cc] Ohelio
gB0=1.e-3 ! [Gs], 100 nT
- !gB0=5.e-5 ! [Gs], 5 nT Ohelio
gx0=Rsolar*1.e5 ! [cm], solar radius
- !gx0 = 1.496e13 ! 1 AU in cm Ohelio
+ ! for Ohelio case
+ !gD0=10. ! [/cc]
+ !gB0=5.e-5 ! [Gs], 5 nT
+ !gx0 = 1.496e13 ! [cm], 1 AU
! get the necessary units
gv0 = gB0/sqrt(4*pi*gD0*mp_cgs) ! [cm/s] ~ 154km/s for gD0=200. and gB0 = 1.e-3
gT0 = gx0/gv0 ! [s] ~ 1.25 hour for above values
- gP0 = gB0**2/(4*pi) ! [erg/cm3]
+ gP0 = gB0**2/(4*pi) ! [erg/cm3]
! 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)
@@ -56,23 +60,23 @@ module helioutils
Tsolar = Tsolar*24.*3600./gt0
!Add gravity if required
- ! TODO: turn gravity on later
- Model%doGrav = .false.
+ Model%doGrav = .true.
if (Model%doGrav) then
!Force spherical gravity (zap non-radial components)
-! Model%doSphGrav = .true.
-! Model%Phi => PhiGrav
+ Model%doSphGrav = .true.
+ Model%Phi => PhiGrav
endif
!Change console output pointer
! don't use for now
-! timeString => helioTime
+ timeString => helioTime
if (Model%isLoud) then
write(*,*) '---------------'
write(*,*) 'Heliosphere normalization'
write(*,*) 'T0 [hr] = ', gT0/3600.
write(*,*) 'x0 [Rsolar] = ', gx0
+ write(*,*) 'D0 [cm-3] = ' , gD0
write(*,*) 'v0 [km/s] = ' , gv0*1.e-5
write(*,*) 'P0 [erg/cm3] = ', gP0
write(*,*) 'B0 [nT] = ' , gB0*1.e5
@@ -90,19 +94,26 @@ module helioutils
! without setting the scaling below, it defaults to 1.
!Add normalization/labels to output slicing
- ! Model%gamOut%tScl = gT0 !/3600.
- ! Model%gamOut%dScl = gD0
- ! Model%gamOut%vScl = gv0 !*1.0e-5 !km/s
- ! Model%gamOut%pScl = gP0
- ! Model%gamOut%bScl = gB0 !*1.e5
+ Model%gamOut%tScl = gT0 !/3600.
+ Model%gamOut%dScl = gD0
+ Model%gamOut%vScl = gv0*1.0e-5 !km/s
+ Model%gamOut%pScl = gP0
+ Model%gamOut%bScl = gB0*1.e5 !nT
- ! Model%gamOut%tID = 'Helio'
- ! Model%gamOut%tID = 's' !'hr'
- ! Model%gamOut%dID = '#/cc'
- ! Model%gamOut%vID = 'km/s'
- ! Model%gamOut%pID = 'erg/cm3'
- ! Model%gamOut%bID = 'nT'
+ Model%gamOut%uID = 'Helio'
+ Model%gamOut%tID = 's'
+ Model%gamOut%dID = '#/cc'
+ Model%gamOut%vID = 'km/s'
+ Model%gamOut%pID = 'erg/cm3'
+ Model%gamOut%bID = 'nT'
+ ! finally rescale the relevant time constants
+ ! note, assume xml file specifies them in [hr]
+ Model%tFin = Model%tFin*3600./gT0
+
+ ! using IOSync from base/ioclock.F90 to sync the other time contants
+ clockScl = Model%IO
+ call IOSync(clockScl,Model%IO,3600./gT0)
end subroutine setHeliosphere
subroutine helioTime(T,tStr)
@@ -112,6 +123,19 @@ module helioutils
write(tStr,'(f9.3,a)' ) T*gT0/3600.0, ' [hr]'
end subroutine helioTime
+ ! slightly different version of PhiGrav in msphutils.F90
+ ! TODO: make this generic (use gravitational constant rather than little g for planets)
+ subroutine PhiGrav(x,y,z,pot)
+ real(rp), intent(in) :: x,y,z
+ real(rp), intent(out) :: pot
+
+ real(rp) :: rad
+ rad = sqrt(x**2.0 + y**2.0 + z**2.0)
+ ! (Msolar*1.D3) converts Msolar into g
+ ! G*M has the units of length * speed^2, thus the gx0*gv0**2 conversion
+ ! the result is the gravitational potential in code units
+ pot = - G_cgs*(Msolar*1.D3)/(gx0*gv0**2)/rad
+ end subroutine PhiGrav
subroutine helio_ibcJ(bc,Model,Grid,State)
! improved versions of Kareems zeroGrad_(i,o)bcJ