From c84b35535d5f241f3d4e415004c15c7b2c861823 Mon Sep 17 00:00:00 2001
From: Ben Ryan <bryan10@illinois.edu>
Date: Tue, 12 Dec 2023 08:27:00 -0700
Subject: [PATCH] Fix v^r, add v_r

---
 scripts/python/get_torus_radial_profiles.py | 11 +++++++----
 scripts/python/phoedf.py                    | 17 ++++++++++++++++-
 2 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/scripts/python/get_torus_radial_profiles.py b/scripts/python/get_torus_radial_profiles.py
index f0c15389..1f887348 100644
--- a/scripts/python/get_torus_radial_profiles.py
+++ b/scripts/python/get_torus_radial_profiles.py
@@ -46,6 +46,7 @@ def get_torus_radial_profiles(dfile):
     ucon = dfile.Getucon()
     ucov = dfile.Getucov()
     vpcon = dfile.GetVpcon()
+    vpcov = dfile.GetVpcov()
     bsq = dfile.GetPm() * 2
     sigma = bsq / rho
     alpha = dfile.alpha
@@ -54,7 +55,6 @@ def get_torus_radial_profiles(dfile):
     # Initialize profiles
     Volume = np.zeros(Nx1_profile) # For volume averages
     Mass = np.zeros(Nx1_profile) # For density-weighted volume averages
-    #Mdot = np.zeros(Nx1_profile) # Total accretion rate
     F_M_in = np.zeros(Nx1_profile) # Inflow mass flux
     F_M_out = np.zeros(Nx1_profile) # Outflow mass flux
     F_Eg = np.zeros(Nx1_profile) # Gas (MHD) energy flux
@@ -65,7 +65,8 @@ def get_torus_radial_profiles(dfile):
     JetP_m = np.zeros(Nx1_profile) # MHD jet power (sigma > 1)
     JetP_g = np.zeros(Nx1_profile) # hydro jet power (sigma > 1)
     beta = np.zeros(Nx1_profile) # Plasma beta
-    vr = np.zeros(Nx1_profile) # SADW radial velocity
+    vconr = np.zeros(Nx1_profile) # SADW contravariant radial three-velocity
+    vcovr = np.zeros(Nx1_profile) # SADW covariant radial three-velocity
     # TODO(BRR) Precalculate some azimuthal averages like <u u_g> for Bernoulli. Store
     # these 2D arrays as well?
 
@@ -110,7 +111,8 @@ def get_torus_radial_profiles(dfile):
                             JetP_m[ip] += dA[b,k,j,i]*Tmunu_concov[1,0]
                             JetP_g[ip] += dA[b,k,j,i]*rho[b,k,j,i]*ucon[b,1,k,j,i]
 
-                vr[ip] += np.sum(dV[b,:,:,i] * vpcon[b,0,:,:,i] / Gamma[b,:,:,i])
+                vconr[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * vpcon[b,0,:,:,i] / Gamma[b,:,:,i])
+                vcovr[ip] += np.sum(dV[b,:,:,i]*rho[b,:,:,i] * vpcov[b,0,:,:,i] / Gamma[b,:,:,i])
 
     beta = Pg_sadw / Pm_sadw
 
@@ -127,7 +129,8 @@ def get_torus_radial_profiles(dfile):
     profiles['beta'] = beta
     profiles['JetP_m'] = JetP_m
     profiles['JetP_g'] = JetP_g
-    profiles['vr'] = vr
+    profiles['vconr'] = vconr / Mass
+    profiles['vcovr'] = vcovr / Mass
 
     return profiles
 
diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py
index 464d1de4..ca16d71e 100644
--- a/scripts/python/phoedf.py
+++ b/scripts/python/phoedf.py
@@ -28,8 +28,11 @@
 
 # ---------------------------------------------------------------------------- #
 # Phoebus-specific derived class from Parthenon's phdf datafile reader
+#
+# geomfile : optional geometry file to read many files more efficiently
+# clip     : Whether to restrict output variables to the range [-1e100, 1e100]
 class phoedf(phdf.phdf):
-    def __init__(self, filename, geomfile=None):
+    def __init__(self, filename, geomfile=None, clip=True):
         super().__init__(filename)
 
         try:
@@ -173,6 +176,7 @@ def flatten_indices(mu, nu):
         self.Pr = None
         self.tau = None
         self.Tg = None
+        self.vpcov = None
         self.ucon = None
         self.ucov = None
         self.bcon = None
@@ -398,6 +402,17 @@ def GetGamma(self):
 
         return self.Gamma
 
+    def GetVpcov(self):
+        if self.vpcov is None:
+            self.vpcov = np.zeros(self.ThreeVectorField)
+
+            vpcon = self.GetVpCon()
+            for ii in range(3):
+                for jj in range(3):
+                    vpcov[:,ii,:,:,:] += self.gcov[:, ii + 1,jj + 1,:,:,:]*vpcon[:,jj,:,:,:]
+
+        return self.vpcov
+
     def GetVpcon(self):
         if self.vpcon is None:
             self.vpcon = np.clip(