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

Merged in heliograv (pull request #60)

Browse Source 
Enabling physical units in helio output

Approved-by: Slava Merkin
Approved-by: Jeff
This commit is contained in:
Elena Provornikova
2022-05-20 02:55:02 +00:00
parent d8cc06189f 82cee4d35b
commit ea32907cce
14 changed files with 306 additions and 172 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
kaiju.sublime-project
View File

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

39
kaipy/gamhelio/ConfigScripts/startup.config
View File

@@ -1,4 +1,11 @@
#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
;nCS is the number density in the current sheet for pressure balance calculation
;TCS is the temperature in the current sheet for pressure balance calculation
[Gamera]
gameraGridFile = heliogrid.h5
@@ -9,28 +16,32 @@ 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
nCS = 1100. ; in [cm-3]
TCS = 1.e6 ; in [K]
[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]

60
kaipy/gamhelio/helioViz.py
View File

@@ -7,10 +7,9 @@ import numpy as np
import kaipy.kaiViz as kv
import kaipy.gamhelio.heliosphere as hsph
from kaipy.kdefs import *
import os
#Tsolar = 25.38
Tsolar = 1.e6
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

100
kaipy/gamhelio/heliosphere.py
View File

@@ -9,12 +9,12 @@ 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
MK = 1.e6 #MegaKelvin
#Adapted to helio grid
class GamsphPipe(GameraPipe):
@@ -28,9 +28,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 +46,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 +111,27 @@ class GamsphPipe(GameraPipe):
Qj[:,:] = 0.5*( Q[:,ja,:] + Q[:,jb,:] )
return Qj
#Radial profile thru cell centers
def RadialProfileGrid(self):
self.GetGrid(doVerbose=True)
#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 +235,53 @@ 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)
r = self.RadialProfileGrid()
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)
r = self.RadialProfileGrid()
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

8
kaipy/gamhelio/wsa2gamera/params.py
View File

@@ -12,16 +12,20 @@ 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.TCS = config.getfloat('Constants','TCS')
self.nCS = config.getfloat('Constants','nCS')
self.B0 = config.getfloat('Normalization','B0')
self.n0 = config.getfloat('Normalization','n0')
self.tMin = config.getfloat('Grid','tMin')
self.tMax = config.getfloat('Grid','tMax')
self.Rin = config.getfloat('Grid','Rin')

8
kaipy/kdefs.py
View File

@@ -43,5 +43,11 @@ NeptuneM0g = 0.142 # Gauss
RNeptuneXE = 3.860 # Rx = X*Re
#------
#Helio
Rsolar = 6.956E5 # [km] Solar radius
#------
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

1
scripts/preproc/wsa2TDgamera.py
View File

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

96
scripts/preproc/wsa2gamera.py
View File

@@ -9,80 +9,66 @@ import matplotlib.pyplot as plt
import kaipy.gamhelio.wsa2gamera.params as params
import kaipy.gamhelio.lib.wsa as wsa
from kaipy.kdefs import *
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 ---------#
# Read params from config file
prm = params.params(args.ConfigFileName)
Ng=prm.NO2
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
T0 = 3.44e6 #2.88e6
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
#grid parameters
# Grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
Rin = prm.Rin
Rout = prm.Rout
Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
Ni = prm.Ni
Nj = prm.Nj
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 +76,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 +96,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 +105,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*kbltz*TCS - (br*B0)**2/8./np.pi)/(rho*n0)/kbltz
# 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 +134,17 @@ 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]
# 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 +153,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()

20
scripts/quicklook/heliopic.py
View File

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

5
src/base/kdefs.F90
View File

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

10
src/chimp/chmpunits.F90
View File

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

8
src/chimp/gridloc.F90
View File

@@ -157,6 +157,9 @@ module gridloc
write(*,*) 'Cartesian not implemented'
stop
case(SPHGRID)
!Take Rin/Rout
DomR(1) = norm2(ebGr%xyz(ebGr%is,ebGr%js,ebGr%ks,XDIR:ZDIR))
DomR(2) = norm2(ebGr%xyz(ebGr%ie-1,ebGr%js,ebGr%ks,XDIR:ZDIR))
write(*,*) 'Initializing SPH locator'
locate=>Loc_SPH
locAux%isInit = .true.
@@ -381,7 +384,7 @@ module gridloc
real(rp) :: dphi,dtheta,dr, thetaMin, rMin
integer :: i0,j0,k0
!E: Add localization routine here
!Add localization routine here
ijk = 0
isInO = .false.
@@ -406,9 +409,6 @@ module gridloc
thetaMin = acos(ebGr%xyz(ebGr%is,ebGr%js,ebGr%ks,ZDIR)/norm2(ebGr%xyz(ebGr%is,ebGr%js,ebGr%ks,:)))
rMin = norm2(ebGr%xyz(ebGr%is,ebGr%js,ebGr%ks,:))
!write(*,*) 'dr, dtheta, dphi', dr, dtheta, dphi
!write(*,*) 'thetaMin, rMin', thetaMin, rMin
! pick the closest one
i0 = min(floor((helioC(IDIR)-rMin)/dr) + 1,ebGr%Nip) !Evenly spaced i

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

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

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

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