NewRadX: Relaxed requirements of specific patch system setup

Author: lucass-carneiro

NewRadX/param.ccl

@@ -1 +1,5 @@
 # Parameter definitions for thorn NewRadX
+
+CCTK_BOOLEAN apply_inner_boundary "Apply radiative BC at inner boundary of a multipatch systems, if such boundary exists (e.g. Thornburg06 inner sphere)" STEERABLE=always
+{
+} "no"

NewRadX/src/newradx.cxx

@@ -1,4 +1,5 @@
 #include <cctk.h>
+#include <cctk_Parameters.h>
 
 #include <loop_device.hxx>
 #include <driver.hxx>
@@ -207,6 +208,7 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
   using namespace CapyrX::MultiPatch::GlobalDerivatives;
 
   DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
 
   const auto symmetries = CarpetX::ghext->patchdata.at(cctk_patch).symmetries;
   const vect<vect<bool, Loop::dim>, 2> is_sym_bnd{
@@ -221,35 +223,41 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
   grid.loop_outermost_int_device<0, 0, 0>(
       grid.nghostzones, is_sym_bnd,
       [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
-        // Make sure we are always on the outer boundary of a patch system
-        assert(p.patch != 0);
-        assert(p.NI[2] >= 0);
+        // Skip inner radial boundary unless requested
+        if (!apply_inner_boundary && p.NI[2] == -1) {
+          return;
+        }
+
+        // At the outer boundary absorb outgoing waves; at the inner
+        // boundary absorb incoming waves (flip sign of v0).
+        const CCTK_REAL sign = (p.NI[2] == -1) ? -1.0 : 1.0;
+        const auto sv0 = sign * v0;
 
         // Find local wave speeds at radiative boundary point p.I
         const auto x = vcoordx(p.I);
         const auto y = vcoordy(p.I);
         const auto z = vcoordz(p.I);
         const auto r = sqrt(pow2(x) + pow2(y) + pow2(z));
 
-        const auto vx = v0 * vcoordx(p.I) / r;
-        const auto vy = v0 * vcoordy(p.I) / r;
-        const auto vz = v0 * vcoordz(p.I) / r;
+        const auto vx = sv0 * x / r;
+        const auto vy = sv0 * y / r;
+        const auto vz = sv0 * z / r;
 
-        // Local derivatives
-        const LocalFirstDerivatives l_dvar{.da = r2o<0>(p, p.I, var),
-                                           .db = r2o<1>(p, p.I, var),
-                                           .dc = r2o<2>(p, p.I, var)};
+        // Local derivatives (stencil auto-selected by boundary location)
+        const LocalFirstDerivatives l_dvar{.da = calc_deriv<0>(p, var),
+                                           .db = calc_deriv<1>(p, var),
+                                           .dc = calc_deriv<2>(p, var)};
 
         // Global derivatives
         const Jacobians jac{VERTEX_JACOBIANS(p)};
         const auto g_dvar{project_first(l_dvar, jac)};
 
         // radiative rhs
         rhs(p.I) = -vx * g_dvar.dx - vy * g_dvar.dy - vz * g_dvar.dz -
-                   v0 * (var(p.I) - var0) / r;
+                   sv0 * (var(p.I) - var0) / r;
 
         if (radpower > 0.0) {
-          // TODO
+          // TODO: Coulomb correction (port from Cartesian overload)
         }
       });
 }

NewRadX/src/newradx.hxx

@@ -1,3 +1,6 @@
+#ifndef NEWRADX_HXX
+#define NEWRADX_HXX
+
 #include <cctk.h>
 #include <loop_device.hxx>
 
@@ -27,22 +30,8 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
       vJ_db_dz, vJ_dc_dx, vJ_dc_dy, vJ_dc_dz
 
 /**
- * @brief Applies radiative boundary condition to the RHS of a state variable.
- * Assumes that:
- *   1. Using multiple patches
- *   2. Patch 0 is cartesian
- *   3. Patches != 0 are spherical-like
- *   4. The local c coordinate is radial and points outward
- *
- * @param cctkGH Pointer to Cactus grid hierarchy struct.
- * @param var State variable which will have boundary conditions applied to it.
- * @param rhs RHS of the evolution equation for @param var
- * @param vcoordx x coordinates grid function.
- * @param vcoordy y coordinates grid function.
- * @param vcoordz z coordinates grid function.
- * @param var0 Value at infinity.
- * @param v0 Propagation speed.
- * @param radpower Radial fall-off exponent
+ * @brief Applies radiative boundary condition to the RHS of a state variable
+ * on multipatch grids.
  */
 void NewRadX_Apply(const cGH *restrict const cctkGH,
                    const Loop::GF3D2<const CCTK_REAL> &var,
@@ -63,3 +52,5 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
                    const CCTK_REAL radpower);
 
 } // namespace NewRadX
+
+#endif // NEWRADX_HXX