Manual merge from heliograv branch.

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\""
+		}
+	],
 }

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
+Rin  = 21.5
-Rout = 215.
+Rout = 220.
-Ni = 128
+Ni   = 128
-Nj = 64
+Nj   = 64
-Nk = 128
+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
+gauss_smooth_width         = 0 ; 8
-plots = yes
+normalized                 = no
-normalized = no
 
 [Constants]
-gamma = 1.5  ;1.05
+gamma  = 1.5 
-NO2   = 4
+Nghost = 4
+Tsolar = 25.38 
 
 [Normalization]
-B0 = 1.e-3  ; 100 nT
+B0 = 1.e-3  ; in [Gs] equals to 100 [nT]
-n0 = 200.   ; 200/cc
+n0 = 200.   ; in [cm-3]
+T0 = 1.e6   ; in [K]
 
 
 

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 = 25.38
-Tsolar = 1.e6
 
 VMax = 800.
 VMin = 300.
@@ -35,8 +34,8 @@ TMin = 0.2
 TMax = 2.
 TCM = "copper"
 
-T0Min = 0.05
+T0Min = 0.01
-T0Max = 0.15
+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 

kaipy/gamhelio/heliosphere.py

@@ -9,12 +9,13 @@ import timeit
 
 #Object to pull from MPI/Serial heliosphere runs (H5 data), extends base
 
-ffam =  "monospace"
+ffam   =  "monospace"
-dLabC = "black" #Default label color
+dLabC  = "black" #Default label color
 dLabFS = "medium" #Default label size
-dBoxC = "lightgrey" #Default box color
+dBoxC  = "lightgrey" #Default box color
-TINY = 1.0e-8
+TINY   = 1.0e-8
-rmStr = "mixtest"
+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.bScl   = 1.0          #->nT
-			self.vScl   = 1.0  #-> km/s
+			self.vScl   = 1.0          #-> km/s
-			self.tScl   = 1.0 #->Seconds
+			self.tScl   = 1.0          #-> Seconds
-			# [EP] added
+			self.dScl   = 1.0          #-> cm-3
-			self.dScl = 1.0
+			self.TScl   = 1.0/kbltz/MK #-> MKelvin
-			self.TScl = 1.0 
 
 		#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

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')

kaipy/kdefs.py

@@ -2,3 +2,4 @@
 import numpy as np
 
 EarthM0g = 0.31 #Gauss
+kbltz    = 1.38e-16 #Boltzmann constant [erg/K]

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]
+    # <HACK>
+    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]
+    # </HACK>
     scTrackName = os.path.join(fdir,scId+".sc.h5")
     with h5py.File(scTrackName,'w') as hf:
-        hf.create_dataset("T" ,data=elapsedSecs)
+        hf.create_dataset("T" ,data=elapsed)
-        hf.create_dataset("X" ,data=x)
+        # Reverse x for gamhelio frame.
-        hf.create_dataset("Y" ,data=y)
+        hf.create_dataset("X" ,data=-c.cartesian.x)
-        hf.create_dataset("Z" ,data=z)
+        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)

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]
+    loc = np.argwhere(gamT > 0.0)[0][0]
-    # t0 = gamUT[loc]  # First non-0 time
+    t0 = gamUT[loc]  # First non-0 time
-    t0 = gamUT[0]
+    # 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)
     # <HACK>
@@ -177,14 +177,14 @@ if __name__ == "__main__":
     # </HACK>
 
     # Save the (float) MJD of the first step.
-    # mjdFileStart = gamMJD[loc]
+    mjdFileStart = gamMJD[loc]
-    mjdFileStart = gamMJD[0]
+    # 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[loc]
-    secFileStart = gamT[0]
+    # secFileStart = gamT[0]
     if debug:
         print("secFileStart = %s" % secFileStart)
 

scripts/datamodel/sunpy.ipynb

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a442a301447dc28f4f41d6f88accb9bd3a5ac6f3340dcf4298832117ec05ec8b
+size 7775

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))

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
+# Generate spherical helio grid
-mp = 1.67e-24
+print("Generating gamera-helio grid Ni = %d, Nj  = %d, Nk = %d " % (Ni, Nj, Nk))
-
-#----------GENERATE HELIO GRID------
-
-print("Generating gamera-helio grid ...")
 
 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!")
+if os.path.exists(prm.gameraGridFile):
-
+    print("Grid file heliogrid.h5 is ready!")
-#----------GENERATE HELIO GRID------
 
 
-############### 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 #####################
+# Get GAMERA grid for further interpolation
-
-# GAMERA GRID
 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
+# radius of the inner boundary
-R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2)  # 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)  
+# Not interpolating temperature, but calculating from the total pressure balance
-#temp = f(Pc[:,0,0],Tc[0,:,0])
+# AFTER interpolating br and rho to the gamera grid
-temp =  1.*T0/rho + (1.**2-(br)**2)*V0**2 / 2e8/1.38 * 1.67/rho   # *****  
+# n_CS*k*T_CS = n*k*T + Br^2/8pi  
-temp_T = temp.T
+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
+#check
-pressure_therm = (n0*rho* temp)*1.38e-16 * 0.1
+#print ("Max and min of temperature in MK")
-pressure_B = (br)**2 *V0**2 / 2.*mp*n0 *0.1
+#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()

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")

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

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
+            in2cms = 1.0e+5  ! km/s -> cm/s
-            in2G   = 1.0e-3 !in [G]
+            in2G   = 1.0e-5 ! nT -> Gs
-            in2s   = L0/in2cms ! time in s 
+            in2s   = 1.0 ! already in s
-            M0g = 0.0 
+            M0g    = 0.0 
-            inPScl = 1.0e-6*1.0e+8/4/pi  !Pressure  unit B[G]^2/4pi *1.e8 in [nPa]
+            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

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

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

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 = .true.
-        Model%doGrav = .false.
         if (Model%doGrav) then
             !Force spherical gravity (zap non-radial components)
-!            Model%doSphGrav = .true.
+            Model%doSphGrav = .true.
-!            Model%Phi => PhiGrav
+            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%tScl = gT0   !/3600.
-        ! Model%gamOut%dScl = gD0
+        Model%gamOut%dScl = gD0
-        ! Model%gamOut%vScl = gv0   !*1.0e-5 !km/s
+        Model%gamOut%vScl = gv0*1.0e-5 !km/s
-        ! Model%gamOut%pScl = gP0
+        Model%gamOut%pScl = gP0
-        ! Model%gamOut%bScl = gB0   !*1.e5
+        Model%gamOut%bScl = gB0*1.e5 !nT
 
-        ! Model%gamOut%tID = 'Helio'
+        Model%gamOut%uID = 'Helio'
-        ! Model%gamOut%tID = 's'  !'hr'
+        Model%gamOut%tID = 's'  
-        ! Model%gamOut%dID = '#/cc'
+        Model%gamOut%dID = '#/cc'
-        ! Model%gamOut%vID = 'km/s'
+        Model%gamOut%vID = 'km/s'
-        ! Model%gamOut%pID = 'erg/cm3'
+        Model%gamOut%pID = 'erg/cm3'
-        ! Model%gamOut%bID = 'nT'
+        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