change distance function to use coord directly also fix diffusion_tes…

…t readme and analysis script

Exec/unit_tests/diffusion_test/analysis/profiles.py

@@ -12,6 +12,8 @@
 import matplotlib.pyplot as plt
 
 import yt
+from yt.frontends.boxlib.api import CastroDataset
+
 
 ## Define RGBA to HEX
 def rgba_to_hex(rgba):
@@ -22,15 +24,22 @@ def rgba_to_hex(rgba):
 
 def get_T_profile(plotfile):
 
-    ds = yt.load(plotfile)
+    ds = CastroDataset(plotfile)
 
     time = float(ds.current_time)
     ad = ds.all_data()
 
     # Sort the ray values by 'x' so there are no discontinuities
     # in the line plot
-    srt = np.argsort(ad['x'])
-    x_coord = np.array(ad['x'][srt])
+
+    coords = {"cartesian":"x",
+              "cylindrical":"z",
+              "spherical":"r"}
+
+    coord = coords[ds.geometry]
+
+    srt = np.argsort(ad[coord])
+    x_coord = np.array(ad[coord][srt])
     temp = np.array(ad['Temp'][srt])
 
     return time, x_coord, temp

Source/driver/Castro.cpp

@@ -3634,6 +3634,7 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time
 #endif
 
     auto geomdata = geom.data();
+    const auto coord = geomdata.Coord();
 
     // Allow the user to limit tagging outside of some distance from the problem center.
 #ifdef AMREX_USE_OMP
@@ -3661,7 +3662,7 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time
             loc[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
 #endif
 
-            Real r = distance(geomdata, loc);
+            Real r = distance(coord, loc);
 
             Real max_dist_lo = 0.0;
             Real max_dist_hi = 0.0;

Source/driver/Castro_util.H

@@ -158,14 +158,12 @@ void position(int i, int j, int k,
 
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-Real distance(GeometryData const& geomdata, GpuArray<Real, 3>& loc)
+Real distance(const int coord, GpuArray<Real, 3>& loc)
 {
     // Returns the distance from the center provided with loc, the position.
     // loc is the form of [x, y, z,] in Cartesian, [r, z, phi] in cylindrical
     // and [r, theta, phi] in spherical
 
-    const auto coord = geomdata.Coord();
-
     if (coord == CoordSys::cartesian) {
         return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1] + loc[2]*loc[2]);
     }

Source/driver/Derive.cpp

@@ -579,7 +579,7 @@ extern "C"
 
       auto dx = geom.CellSizeArray();
       auto problo = geom.ProbLoArray();
-      auto geomdata = geom.data();
+      auto coord = geom.Coord();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -608,7 +608,7 @@ extern "C"
           der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + dat(i,j,k,2)*loc[1]) / (dat(i,j,k,0)*r);
 #endif
         } else {
-          Real r = distance(geomdata, loc);
+          Real r = distance(coord, loc);
 
           der(i,j,k,0) = (dat(i,j,k,1)*loc[0] +
                           dat(i,j,k,2)*loc[1] +
@@ -631,7 +631,7 @@ extern "C"
 
       auto dx = geom.CellSizeArray();
       auto problo = geom.ProbLoArray();
-      auto geomdata = geom.data();
+      const auto coord = geom.Coord();
 
       amrex::ParallelFor(bx,
       [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -661,7 +661,7 @@ extern "C"
           der(i,j,k,0) = (-dat(i,j,k,1)*loc[1] + dat(i,j,k,2)*loc[0]) / (dat(i,j,k,0)*r);
 #endif
         } else {
-          Real r = distance(geomdata, loc);
+          Real r = distance(coord, loc);
 
           // we really mean just the velocity component that is
           // perpendicular to radial, and in general 3-d (e.g. a
@@ -951,7 +951,7 @@ extern "C"
 
       auto dx = geom.CellSizeArray();
 
-      const int coord_type = geom.Coord();
+      const auto coord_type = geom.Coord();
 
 #if AMREX_SPACEDIM == 2
       auto problo = geom.ProbLoArray();

Source/gravity/Gravity.cpp

@@ -1334,7 +1334,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab&
     const Real dr = dx[0] / static_cast<Real>(gravity::drdxfac);
 
     const auto problo = geom.ProbLoArray();
-    const auto geomdata = geom.data();
+    const auto coord_type = geom.Coord();
 
 #ifdef _OPENMP
 #pragma omp parallel
@@ -1368,7 +1368,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab&
             loc[2] = 0.0_rt;
 #endif
 
-            Real r = distance(geomdata, loc);
+            Real r = distance(coord_type, loc);
 
             int index = static_cast<int>(r / dr);
 
@@ -1457,8 +1457,7 @@ Gravity::compute_radial_mass(const Box& bx,
     Real dr = dx[0] / static_cast<Real>(gravity::drdxfac);
     Real drinv = 1.0_rt / dr;
 
-    const int coord_type = geom.Coord();
-    const auto geomdata = geom.data();
+    const auto coord_type = geom.Coord();
 
     AMREX_ALWAYS_ASSERT(coord_type >= 0 && coord_type <= 2);
 
@@ -1509,7 +1508,7 @@ Gravity::compute_radial_mass(const Box& bx,
         loc[2]= problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
         Real lo_k = problo[2] + static_cast<Real>(k) * dx[2] - problem::center[2];
 
-        Real r = distance(geomdata, loc);
+        Real r = distance(coord_type, loc);
         int index = static_cast<int>(r * drinv);
 
         // We may be coming in here with a masked out zone (in a zone on a coarse

Source/reactions/Castro_react.cpp

@@ -209,7 +209,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra
         const auto dx = geom.CellSizeArray();
 #ifdef MODEL_PARSER
         const auto problo = geom.ProbLoArray();
-        const auto geomdata = geom.data();
+        const auto coord = geom.Coord();
 #endif
 
 #if defined(AMREX_USE_GPU)
@@ -286,7 +286,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra
                 if (domain_is_plane_parallel) {
                     dist = rr[AMREX_SPACEDIM-1];
                 } else {
-                    dist = distance(geomdata, rr);
+                    dist = distance(coord, rr);
                 }
 
                 burn_state.T_fixed = interpolate(dist, model::itemp);
@@ -566,7 +566,7 @@ Castro::react_state(Real time, Real dt)
 
         const auto dx = geom.CellSizeArray();
         const auto problo = geom.ProbLoArray();
-        const auto geomdata = geom.data();
+        const auto coord = geom.Coord();
 
 #if defined(AMREX_USE_GPU)
         ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
@@ -635,7 +635,7 @@ Castro::react_state(Real time, Real dt)
                 if (domain_is_plane_parallel) {
                     dist = rr[AMREX_SPACEDIM-1];
                 } else {
-                    dist = distance(geomdata, rr);
+                    dist = distance(coord, rr);
                 }
 
                 burn_state.T_fixed = interpolate(dist, model::itemp);

Source/sources/Castro_sponge.cpp

@@ -81,7 +81,7 @@ Castro::apply_sponge(const Box& bx,
 
   auto dx = geom.CellSizeArray();
   auto problo = geom.ProbLoArray();
-  auto geomdata = geom.data();
+  const auto coord = geom.Coord();
 
   amrex::ParallelFor(bx,
   [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -125,7 +125,7 @@ Castro::apply_sponge(const Box& bx,
     Real sponge_factor = 0.0_rt;
 
     if (sponge_lower_radius >= 0.0_rt && sponge_upper_radius > sponge_lower_radius) {
-      Real rad = distance(geomdata, r);
+      Real rad = distance(coord, r);
 
       if (rad < sponge_lower_radius) {
         sponge_factor = sponge_lower_factor;