1-10 au helio all fixed

2026-01-09 15:17:56 -05:00 · 2021-11-05 07:58:54 -06:00
parent 717f5af79f
commit 994e394018
4 changed files with 62 additions and 63 deletions
--- a/kaipy/gamhelio/helioViz.py
+++ b/kaipy/gamhelio/helioViz.py
@@ -15,22 +15,34 @@ VMax = 800.
 VMin = 300.
 MagVCM = "inferno"

-DMax = 150.
-DMin = 2000.
+#inner helio
+#DMax = 150.
+#DMin = 2000.
+#1-10 au helio
+DMax = 25.
+DMin = 5.
 DCM = "copper_r"

-D0Max = 15.
-D0Min = 1.
+#limits for iSlice
+#21.5 R_S
 #D0Max = 1000.
 #D0Min = 300.
+#10 au
+#D0Max = 0.2
+#D0Min = 0.002
+#1 au
+D0Max = 1.
+D0Min = 25.
 D0CM = "copper_r"

-TMax = 0.12
 TMin = 0.01
+TMax = 0.12
 TCM = "copper"

 #BMax = 150.
 #BMin = -150.
+#BMax = 0.05
+#BMin = -0.05
 BMax = 5.
 BMin = -5.
 BCM = "coolwarm"
@@ -148,7 +160,15 @@ def PlotMerBrNorm(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
 	xr, zr, xl, zl, r = gsph.MeridGridHalfs()
 	Ax.pcolormesh(xr,zr,Br_r,cmap=BCM,norm=vB,shading='auto')
 	Ax.pcolormesh(xl,zl,Br_l,cmap=BCM,norm=vB,shading='auto')
-
+	#plot heliospheric current sheet
+	#cell-cent coords first
+	xr_c = 0.25*( xr[:-1,:-1]+xr[:-1,1:]+xr[1:,:-1]+xr[1:,1:] )
+	zr_c = 0.25*( zr[:-1,:-1]+zr[:-1,1:]+zr[1:,:-1]+zr[1:,1:] )
+	xl_c = 0.25*( xl[:-1,:-1]+xl[:-1,1:]+xl[1:,:-1]+xl[1:,1:] )
+	zl_c = 0.25*( zl[:-1,:-1]+zl[:-1,1:]+zl[1:,:-1]+zl[1:,1:] )
+	#plot Br=0
+	Ax.contour(xr_c,zr_c,Br_r,[0.],colors='black')
+	Ax.contour(xl_c,zl_c,Br_l,[0.],colors='black')
 	kv.SetAx(xyBds,Ax)

 	if (doDeco):
@@ -296,8 +316,6 @@ def PlotiSlD(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):

 #Plot Br and current sheet (Br=0) at 1 AU
 def PlotiSlBr(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
-	BMin = -5.
-	BMax = 5.
 	vB = kv.genNorm(BMin, BMax, doLog=False, midP=None)
 	if (AxCB is not None):
 		AxCB.clear()
@@ -381,9 +399,6 @@ def PlotiSlBrRotatingFrame(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True

 #Plot Temperature at 1 AU
 def PlotiSlTemp(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True):
-	#colorbar limits
-	TMin = 0.02
-	TMax = 0.12
 	vT = kv.genNorm(TMin, TMax, doLog=False, midP=None)

 	if (AxCB is not None):
--- a/kaipy/gamhelio/heliosphere.py
+++ b/kaipy/gamhelio/heliosphere.py
@@ -28,13 +28,14 @@ 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/2./1.e6 #in MK
-
+		#self.TScl = 1.e-6/4/np.pi/200./1.38e-16/1.e6 #in MK
+    
+                # [OHelio]
 		self.bScl = 5.    #->nT
 		self.vScl = 34.5  #-> km/s
 		self.tScl = 1.4e8/34.5
 		self.dScl = 10. #cm-3
-		self.TScl = 0.07 #in MK
+		self.TScl = 0.144 #in MK

 		#2D equatorial grid
 		self.xxi = [] ; self.yyi = [] #corners
@@ -45,7 +46,6 @@ class GamsphPipe(GameraPipe):

 		#inner boundary distance
 		self.R0 = self.xxc[0,0]
-		print (self.R0)
 		
 	def OpenPipe(self,doVerbose=True):
 		GameraPipe.OpenPipe(self,doVerbose)
@@ -205,9 +205,9 @@ class GamsphPipe(GameraPipe):
 		Q = self.GetVar(vID,sID,vScl,doVerb)

 		#cell centered values from the last cell
-		#Qi = Q[-1,:,:]
+		Qi = Q[-1,:,:]
                #cell centered values from the first cell
-		Qi = Q[0,:,:]
+		#Qi = Q[0,:,:]
 		#jd_c = self.MJDs[sID]
 		#print ('jd_c = ', jd_c)
 		return Qi
--- a/scripts/helio/ih2oh.py
+++ b/scripts/helio/ih2oh.py
@@ -25,20 +25,20 @@ gamma = prm.gamma
 mp = 1.67e-24
 kb = 1.38e-16

+#normalization in IH
 B0 = prm.B0
 n0 = prm.n0
 V0 = B0/np.sqrt(4*np.pi*mp*n0)
-T0 = B0*B0/4/np.pi/n0/kb/2. #in K
+T0 = B0*B0/4/np.pi/n0/kb #in K p = nkT

-print ("inner helio units")
+print ("inner helio normalization")
 print (B0, n0, V0, T0)

 #normalization in OH
-B0OH = 5.e-5 #Gs
-n0OH = 10 #cm-3
-V0OH = B0OH/np.sqrt(4*np.pi*mp*n0OH) #Alfven speed at 1 AU
-#V0OH = 400 #km/s
-T0OH = B0OH*B0OH/4/np.pi/n0OH/kb/2.#in K
+B0OH = 5.e-5 # [Gs] 5 nT = 5.e-5 Gs
+n0OH = 10 # [cm-3]
+V0OH = B0OH/np.sqrt(4*np.pi*mp*n0OH) #Alfven speed at 1 AU 34.5 [km/s]
+T0OH = B0OH*B0OH/4/np.pi/n0OH/kb #in K p = nkT

 print ("outer helio units")
 print (B0OH, n0OH, V0OH, T0OH)
@@ -46,6 +46,7 @@ print (B0OH, n0OH, V0OH, T0OH)

 ############### READ GAMERA solution at 1 AU #####################
 f = h5py.File(prm.wsaFile,'r')
+#the latest Step saved in inner helio solution wsa.h5
 step = 'Step#4'

 f[step].attrs.keys()
@@ -89,31 +90,14 @@ T_wsa = T[:,:,Nr-1]
 print ("1AU arrays")
 print (bi_wsa.shape, v_wsa.shape, n_wsa.shape, T_wsa.shape)

-#renormalize 
+#renormalize inner helio solution
 bi_wsa = bi_wsa * B0/B0OH
 n_wsa = n_wsa * n0/n0OH
-v_wsa = v_wsa *V0/V0OH
-T_wsa = T_wsa *T0/T0OH
-
-
-############### WSA STUFF #####################
-#jd_c,phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa,T_wsa = wsa.read(prm.wsaFile,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
-
-#bi_wsa /= B0
-#n_wsa  /= (mp*n0)
-#v_wsa /= V0
-#convert julian date from wsa fits into modified julian date
-#mjd_c = jd_c - 2400000.5
+v_wsa = v_wsa * V0/V0OH
+T_wsa = T_wsa * T0     
 # keep temperature in K
-############### WSA STUFF #####################
-
-############### GAMERA STUFF #####################

+#######Interpolate to GAMERA grid###########
 # GAMERA GRID
 with h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile),'r') as f:
    x=f['X'][:]
@@ -174,15 +158,8 @@ vr = fv(Pc[:,0,0],Tc[0,:,0])
 f      = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa,kx=1,ky=1)  
 rho = f(Pc[:,0,0],Tc[0,:,0])

-f      = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa,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
-
-#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
-
+fT      = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa,kx=1,ky=1)  
+temp = fT(Pc[:,0,0],Tc[0,:,0])

 fbi = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],br,kx=1,ky=1)
 fv  = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],vr,kx=1,ky=1)
@@ -192,10 +169,13 @@ 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)]
-rho*=(R0/Rc[0,0,:Ng])**2
-br*=(R0/Rc[0,0,:Ng])**2
-br_kface*=(R0/Rc[0,0,:Ng])**2
+rho *= (R0/Rc[0,0,:Ng])**2
+br *= (R0/Rc[0,0,:Ng])**2
+br_kface *= (R0/Rc[0,0,:Ng])**2

+#FIX: 
+#For now I use wsa.h5 which do not have mjd inside
+#so hardcoded using WSA fits value + 200*4637/60/60/24 [number of days]
 mjd_c = 58005.83415

 print ("writing out innerbc.h5...")
--- a/src/gamera/ICs/wsa.F90
+++ b/src/gamera/ICs/wsa.F90
@@ -37,7 +37,11 @@ module usergamic
    ! 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.
-    real(rp) :: Rbc = 21.5   
+
+    !for inner heliosphere
+    !real(rp) :: Rbc = 21.5
+    ![OHelio] for 1-10 au simualtion   
+     real(rp) :: Rbc = 1.

    real(rp) :: Tsolar ! Solar rotation period, defined in apps/helioutils.F90
    
@@ -95,11 +99,11 @@ module usergamic
        !Rho0  = 1.0    ! 200/cc
        !P0    = 1.0e-4*Rho0*Cs0**2.0/Model%gamma

-        !for 1-10 au helio
-        Cs0 = 0.78 !27 km/s
-        B0 = 1. ! 5nT
-        Rho0 = 1. ! 10/cc
-        Vslow = 8.5 !300 km/s 
+        ![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