In sph geom, (r|z)_coord contains radial position

They were previously set to r_coord - rcmb.

Author: amorison

stagpy/_step.py

@@ -83,12 +83,12 @@ def __init__(self, header, par):
                 self._coords[0] = np.array(np.pi / 2)
             elif self.twod_xz:
                 self._coords[1] = np.array(0)
-            r_coord = self.z_coord + self.rcmb
+            self._coords[2] += self.rcmb
             if par['magma_oceans_in']['magma_oceans_mode']:
-                r_coord += header['mo_lambda']
-                r_coord *= header['mo_thick_sol']
+                self._coord[2] += header['mo_lambda']
+                self._coord[2] *= header['mo_thick_sol']
             t_mesh, p_mesh, r_mesh = np.meshgrid(
-                self.x_coord, self.y_coord, r_coord, indexing='ij')
+                self.t_coord, self.p_coord, self.r_coord, indexing='ij')
             # compute cartesian coordinates
             # z along rotation axis at theta=0
             # x at th=90, phi=0
@@ -166,8 +166,10 @@ def threed(self):
     def at_z(self, zval):
         """Return iz closest to given zval position.
 
-        This is different than radial position in spherical geometry.
+        In spherical geometry, the bottom boundary is considered to be at z=0.
         """
+        if self.curvilinear:
+            zval += self.rcmb
         return np.argmin(np.abs(self.z_coord - zval))
 
     def __getattr__(self, attr):

stagpy/plates.py

@@ -18,7 +18,7 @@ def detect_plates_vzcheck(step, seuil_memz):
     n_z = step.geom.nztot
     nphi = step.geom.nptot  # -1? should be OK, ghost not included
     rcmb = max(0, step.geom.rcmb)
-    radius = step.geom.r_coord + rcmb
+    radius = step.geom.r_coord
     radiusgrid = step.geom.rgeom[:, 0] + rcmb
     dphi = 1 / nphi
 
@@ -93,7 +93,8 @@ def detect_plates(step, vrms_surface, fids, time):
         dsa = step.sdat.par['boundaries']['air_thickness']
         # we are a bit below the surface; should check if you are in the
         # thermal boundary layer
-        indsurf = np.argmin(abs((1 - dsa) - step.geom.r_coord)) - 4
+        indsurf = np.argmin(
+            np.abs(1 - dsa - step.geom.r_coord + step.geom.rcmb)) - 4
     else:
         indsurf = -1
 
@@ -211,9 +212,11 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
         # to be just below
         # the surface (that is considered plane here); should check if you are
         # in the thermal boundary layer
-        indsurf = np.argmin(abs((1 - dsa) - step.geom.r_coord)) - 4
+        indsurf = np.argmin(
+            np.abs(1 - dsa - step.geom.r_coord + step.geom.rcmb)) - 4
         # depth to detect the continents
-        indcont = np.argmin(abs((1 - dsa) - np.array(step.geom.r_coord))) - 10
+        indcont = np.argmin(
+            np.abs(1 - dsa - step.geom.r_coord + step.geom.rcmb)) - 10
     else:
         indsurf = -1
         indcont = -1  # depth to detect continents
@@ -514,7 +517,8 @@ def lithospheric_stress(step, trench, ridge, time):
         # to be just below
         dsa = step.sdat.par['boundaries']['air_thickness']
         # depth to detect the continents
-        indcont = np.argmin(abs((1 - dsa) - step.geom.r_coord)) - 10
+        indcont = np.argmin(
+            np.abs(1 - dsa - step.geom.r_coord + step.geom.rcmb)) - 10
     else:
         # depth to detect continents
         indcont = -1

stagpy/processing.py

@@ -349,7 +349,7 @@ def stream_function(step):
     psi = np.zeros_like(v_x)
     if step.geom.spherical:  # YZ annulus
         # positions
-        r_nc = step.geom.r_coord + step.geom.rcmb  # numerical centers
+        r_nc = step.geom.r_coord  # numerical centers
         r_pc = step.geom.r_mesh[0, 0, :]  # physical centers
         r_nw = r_edges(step)[0][:2]  # numerical walls of first cell
         # vz at center of bottom cells

stagpy/stagyyparsers.py

@@ -651,7 +651,6 @@ def _read_coord_h5(files, shapes, header, twod):
     if header['rcmb'] == 0:
         header['rcmb'] = -1
     else:
-        # could make the difference between r_coord and z_coord
         header['e3_coord'] = header['e3_coord'] - header['rcmb']
     if twod is None or 'X' in twod:
         header['e1_coord'] = header['e1_coord'][:-1]