From a1931d7011704cf62c794a357556b6d817c662cd Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Tue, 30 Sep 2025 11:59:31 -0500 Subject: [PATCH 1/8] NewRadX: Added the ability to use CapyrX Multipatch --- NewRadX/src/newradx.cxx | 494 ++++++++++++++++++++++++---------------- NewRadX/src/newradx.hxx | 57 ++++- 2 files changed, 344 insertions(+), 207 deletions(-) diff --git a/NewRadX/src/newradx.cxx b/NewRadX/src/newradx.cxx index 3134f8cc..22a790d9 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -1,225 +1,327 @@ #include -#include -#include -#include "loop.hxx" -#include "loop_device.hxx" -#include "driver.hxx" + +#include +#include + #include "newradx.hxx" +#include + namespace NewRadX { using namespace Loop; using namespace std; -namespace { -template constexpr T pow2(const T x) { return x * x; } -} // namespace +template static inline constexpr T pow2(const T x) { + return x * x; +} -// Adapted from NewRad thorn by E. Schnetter, used with Carpet. -// Original code adapted from BSSN_MoL's files NewRad.F and newrad.h. -// This code was probably originally written by Miguel Alcubierre. -void NewRadX_Apply(const cGH *restrict const cctkGH, - const Loop::GF3D2 var, - Loop::GF3D2 rhs, - const CCTK_REAL var0, //!< value at infinity - const CCTK_REAL v0, //!< propagation speed - const CCTK_REAL radpower //!< exponent in radial fall-off -) { - DECLARE_CCTK_ARGUMENTS; +template +static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL +c2o(const Loop::PointDesc &p, const vect &pI, + const Loop::GF3D2 &gf) noexcept { + const auto num{gf(pI + p.DI[dir]) - gf(pI - p.DI[dir])}; + const auto den{1.0 / (2.0 * p.DX[dir])}; + return num * den; +} - constexpr vect DI{1, 0, 0}; - constexpr vect DJ{0, 1, 0}; - constexpr vect DK{0, 0, 1}; - - const CCTK_REAL dx = CCTK_DELTA_SPACE(0); - const CCTK_REAL dy = CCTK_DELTA_SPACE(1); - const CCTK_REAL dz = CCTK_DELTA_SPACE(2); - - const auto derivx = [=] CCTK_DEVICE CCTK_HOST( - const GF3D2 &u_, - const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto I = p.I; - if (p.NI[0] == 0) - // interior - // in direction parallel to face/edge, apply symmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return (u_(I + DI) - u_(I - DI)) / (2 * dx); - if (p.NI[0] == +1) - // upper boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return +(3 * u_(I) - 4 * u_(I - DI) + u_(I - 2 * DI)) / (2 * dx); - if (p.NI[0] == -1) - // lower boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return -(3 * u_(I) - 4 * u_(I + DI) + u_(I + 2 * DI)) / (2 * dx); - assert(0); - }; - - const auto derivy = [=] CCTK_DEVICE CCTK_HOST( - const GF3D2 &u_, - const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto I = p.I; - if (p.NI[1] == 0) - // interior - // in direction parallel to face/edge, apply symmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return (u_(I + DJ) - u_(I - DJ)) / (2 * dy); - if (p.NI[1] == +1) - // upper boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return +(3 * u_(I) - 4 * u_(I - DJ) + u_(I - 2 * DJ)) / (2 * dy); - if (p.NI[1] == -1) - // lower boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return -(3 * u_(I) - 4 * u_(I + DJ) + u_(I + 2 * DJ)) / (2 * dy); +template +static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL +l2o(const Loop::PointDesc &p, const vect &pI, + const Loop::GF3D2 &gf) noexcept { + const auto num{-3.0 * gf(pI) + 4.0 * gf(pI + p.DI[dir]) - + gf(pI + 2 * p.DI[dir])}; + const auto den{1.0 / (2.0 * p.DX[dir])}; + return num * den; +} + +template +static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL +r2o(const Loop::PointDesc &p, const vect &pI, + const Loop::GF3D2 &gf) noexcept { + const auto num{3.0 * gf(pI) - 4.0 * gf(pI - p.DI[dir]) + + gf(pI - 2 * p.DI[dir])}; + const auto den{1.0 / (2.0 * p.DX[dir])}; + return num * den; +} + +template +static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL calc_deriv( + const Loop::PointDesc &p, const Loop::GF3D2 &gf) noexcept { + if (p.NI[dir] == 0) { + /* interior + * in direction parallel to face/edge, apply symmetric stencil + * currently second-order accurate + * TODO: apply same finite-difference order as interior + */ + return c2o(p, p.I, gf); + + } else if (p.NI[dir] == +1) { + /* upper boundary + * in direction normal to face/edge/corner, apply asymmetric stencil + * currently second-order accurate + * TODO: apply same finite-difference order as interior + */ + return r2o(p, p.I, gf); + + } else if (p.NI[dir] == -1) { + /* lower boundary + * in direction normal to face/edge/corner, apply asymmetric stencil + * currently second-order accurate + * TODO: apply same finite-difference order as interior + */ + return l2o(p, p.I, gf); + + } else { assert(0); - }; - - const auto derivz = [=] CCTK_DEVICE CCTK_HOST( - const GF3D2 &u_, - const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - const auto I = p.I; - if (p.NI[2] == 0) - // interior - // in direction parallel to face/edge, apply symmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return (u_(I + DK) - u_(I - DK)) / (2 * dz); - if (p.NI[2] == +1) - // upper boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return +(3 * u_(I) - 4 * u_(I - DK) + u_(I - 2 * DK)) / (2 * dz); - if (p.NI[2] == -1) - // lower boundary - // in direction normal to face/edge/corner, apply asymmetric stencil - // currently second-order accurate - // TODO: apply same finite-difference order as interior - return -(3 * u_(I) - 4 * u_(I + DK) + u_(I + 2 * DK)) / (2 * dz); + } +} + +template +static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL +calc_deriv_interior_upper(const Loop::PointDesc &p, + const Loop::GF3D2 &gf) noexcept { + if (p.NI[dir] == 0) { + /* interior + * in direction parallel to face/edge, apply symmetric stencil + * currently second-order accurate + * TODO: apply same finite-difference order as interior + */ + return c2o(p, p.I, gf); + + } else if (p.NI[dir] == +1) { + /* upper boundary + * in direction normal to face/edge/corner, apply asymmetric stencil + * currently second-order accurate + * TODO: apply same finite-difference order as interior + */ + return r2o(p, p.I, gf); + + } else { assert(0); - }; + } +} + +void NewRadX_Apply(const cGH *restrict const cctkGH, + const Loop::GF3D2 &var, + const Loop::GF3D2 &rhs, const CCTK_REAL var0, + const CCTK_REAL v0, const CCTK_REAL radpower) { + DECLARE_CCTK_ARGUMENTS; + + const auto symmetries = CarpetX::ghext->patchdata.at(cctk_patch).symmetries; + const vect, 2> is_sym_bnd{ + {symmetries[0][0] != CarpetX::symmetry_t::none, + symmetries[0][1] != CarpetX::symmetry_t::none, + symmetries[0][2] != CarpetX::symmetry_t::none}, + {symmetries[1][0] != CarpetX::symmetry_t::none, + symmetries[1][1] != CarpetX::symmetry_t::none, + symmetries[1][2] != CarpetX::symmetry_t::none}}; + + const Loop::GridDescBaseDevice grid(cctkGH); + grid.loop_outermost_int_device<0, 0, 0>( + grid.nghostzones, is_sym_bnd, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // The main part of the boundary condition assumes that we have an + // outgoing radial wave with some speed v0: + // + // var = var0 + u(r-v0*t)/r + // + // This implies the following differential equation: + // + // d_t var = - v^i d_i var - v0 (var - var0) / r + // + // where vi = v0 xi/r + + // coordinate radius at p.I + const auto r = sqrt(pow2(p.x) + pow2(p.y) + pow2(p.z)); + + // Find local wave speeds at radiative boundary point p.I + const auto vx = v0 * p.x / r; + const auto vy = v0 * p.y / r; + const auto vz = v0 * p.z / r; + + // Derivatives + const auto varx = calc_deriv<0>(p, var); + const auto vary = calc_deriv<1>(p, var); + const auto varz = calc_deriv<2>(p, var); + + // radiative rhs + rhs(p.I) = + -vx * varx - vy * vary - vz * varz - v0 * (var(p.I) - var0) / r; + + if (radpower > 0) { + // When solution is known to have a Coulomb-like component + // asymptotically, this is estimated and extrapolated from + // physical interior points to the radiative boundary. I.e.: + // + // var = var0 + u(r-v0*t)/r + h(t)/r^n + // + // This implies the following differential equation: + // + // d_t var = - v^i d_i var - v0 (var - var0) / r + // + v0 (1 - n) h / r^n+1 + d_t h / r^n + // ~ - v^i d_i var - v0 (var - var0) / r + // + d_t h / r^n + O(1/r^n+1) + // + // where vi = v0 xi/r + // + // the Coulomb term, d_t h / r ^n , is estimated by extrapolating + // from the nearest interior points + // + // Displacement to get from p.I to interior point placed + // nghostpoints away + const vect displacement{grid.nghostzones[0] * p.NI[0], + grid.nghostzones[1] * p.NI[1], + grid.nghostzones[2] * p.NI[2]}; + const vect intp = p.I - displacement; + + assert(intp[0] >= grid.nghostzones[0]); + assert(intp[1] >= grid.nghostzones[1]); + assert(intp[2] >= grid.nghostzones[2]); + assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1); + assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1); + assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1); + + // coordinates at p.I-displacement + const auto xint = p.x - displacement[0] * p.DX[0]; + const auto yint = p.y - displacement[1] * p.DX[1]; + const auto zint = p.z - displacement[2] * p.DX[2]; + const auto rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint)); + + // Find local wave speeds at physical point p.I-displacement + const auto vxint = v0 * xint / rint; + const auto vyint = v0 * yint / rint; + const auto vzint = v0 * zint / rint; + + // Derivatives at physical point p.I-displacement + const auto varxint = c2o<0>(p, intp, var); + const auto varyint = c2o<1>(p, intp, var); + const auto varzint = c2o<2>(p, intp, var); + + // Eextrapolate Coulomb component, rescale to account for radial + // fall-off + const auto rad = -vxint * varxint - vyint * varyint - + vzint * varzint - v0 * (var(intp) - var0) / rint; + const auto aux = (rhs(intp) - rad) * pow(rint / r, radpower); + + // Radiative rhs with extrapolated Coulomb correction + rhs(p.I) += aux; + } + }); +} + +void NewRadX_Apply(const cGH *restrict const cctkGH, + const Loop::GF3D2 &var, + const Loop::GF3D2 &rhs, + const Loop::GF3D2 &vcoordx, + const Loop::GF3D2 &vcoordy, + const Loop::GF3D2 &vcoordz, + const CCTK_REAL var0, const CCTK_REAL v0, + const CCTK_REAL radpower) { + DECLARE_CCTK_ARGUMENTS; const auto symmetries = CarpetX::ghext->patchdata.at(cctk_patch).symmetries; - const vect, 2> is_sym_bnd { - { - symmetries[0][0] != CarpetX::symmetry_t::none, - symmetries[0][1] != CarpetX::symmetry_t::none, - symmetries[0][2] != CarpetX::symmetry_t::none - }, - { - symmetries[1][0] != CarpetX::symmetry_t::none, - symmetries[1][1] != CarpetX::symmetry_t::none, - symmetries[1][2] != CarpetX::symmetry_t::none - } - }; + const vect, 2> is_sym_bnd{ + {symmetries[0][0] != CarpetX::symmetry_t::none, + symmetries[0][1] != CarpetX::symmetry_t::none, + symmetries[0][2] != CarpetX::symmetry_t::none}, + {symmetries[1][0] != CarpetX::symmetry_t::none, + symmetries[1][1] != CarpetX::symmetry_t::none, + symmetries[1][2] != CarpetX::symmetry_t::none}}; const Loop::GridDescBaseDevice grid(cctkGH); grid.loop_outermost_int_device<0, 0, 0>( grid.nghostzones, is_sym_bnd, - [=] CCTK_DEVICE(const Loop::PointDesc &p) - CCTK_ATTRIBUTE_ALWAYS_INLINE { - // The main part of the boundary condition assumes that we have an - // outgoing radial wave with some speed v0: + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + if (p.patch != 0) { + /* + * The main part of the boundary condition assumes that we have an + * outgoing radial wave with some speed v0: + * + * var = var0 + u(r-v0*t)/r + * + * This implies the following differential equation: + * + * d_t var = - v^i d_i var - v0 (var - var0) / r + * + * where vi = v0 xi/r + */ + + // coordinate radius at 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)); + + // Find local wave speeds at radiative boundary point p.I + const auto vx = v0 * x / r; + const auto vy = v0 * y / r; + const auto vz = v0 * z / r; + const auto vr = sqrt(pow2(vx) + pow2(vy) + pow2(vz)); + + // Derivatives + const auto derivz = calc_deriv_interior_upper<2>(p, var); + + // Radiative rhs + rhs(p.I) = -vr * derivz - v0 * (var(p.I) - var0) / r; + + if (radpower > 0.0) { + // When solution is known to have a Coulomb-like component + // asymptotically, this is estimated and extrapolated from + // physical interior points to the radiative boundary. I.e.: // - // var = var0 + u(r-v0*t)/r + // var = var0 + u(r-v0*t)/r + h(t)/r^n // // This implies the following differential equation: // // d_t var = - v^i d_i var - v0 (var - var0) / r + // + v0 (1 - n) h / r^n+1 + d_t h / r^n + // ~ - v^i d_i var - v0 (var - var0) / r + // + d_t h / r^n + O(1/r^n+1) // // where vi = v0 xi/r + // + // the Coulomb term, d_t h / r ^n , is estimated by extrapolating + // from the nearest interior points + // + // Displacement to get from p.I to interior point placed + // nghostpoints away + const vect displacement{0, 0, grid.nghostzones[2]}; + const vect intp = p.I - displacement; + + assert(intp[0] >= grid.nghostzones[0]); + assert(intp[1] >= grid.nghostzones[1]); + assert(intp[2] >= grid.nghostzones[2]); + assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1); + assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1); + assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1); + + // coordinates at p.I-displacement + const auto xint = vcoordx(intp); + const auto yint = vcoordy(intp); + const auto zint = vcoordz(intp); + const auto rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint)); + + // Find local wave speeds at physical point p.I-displacement + const auto vxint = v0 * xint / rint; + const auto vyint = v0 * yint / rint; + const auto vzint = v0 * zint / rint; + const auto vrint = sqrt(pow2(vxint) + pow2(vyint) + pow2(vzint)); + + // Derivative at physical point p.I-displacement + const auto derivzint = c2o<2>(p, intp, var); + + // Extrapolate Coulomb component, rescale to account for radial + // fall-off + const auto rad = + -vrint * derivzint - v0 * (var(intp) - var0) / rint; + const auto aux = (rhs(intp) - rad) * pow(rint / r, radpower); - // coordinate radius at p.I - const CCTK_REAL r = sqrt(pow2(p.x) + pow2(p.y) + pow2(p.z)); - - // Find local wave speeds at radiative boundary point p.I - const CCTK_REAL vx = v0 * p.x / r; - const CCTK_REAL vy = v0 * p.y / r; - const CCTK_REAL vz = v0 * p.z / r; - // CCTK_REAL const vr = sqrt(pow2(vx) + pow2(vy) + pow2(vz)); - - // Derivatives - const CCTK_REAL varx = derivx(var, p); - const CCTK_REAL vary = derivy(var, p); - const CCTK_REAL varz = derivz(var, p); - - // radiative rhs - rhs(p.I) = - -vx * varx - vy * vary - vz * varz - v0 * (var(p.I) - var0) / r; - - if (radpower > 0) { - // When solution is known to have a Coulomb-like component - // asymptotically, this is estimated and extrapolated from - // physical interior points to the radiative boundary. I.e.: - // - // var = var0 + u(r-v0*t)/r + h(t)/r^n - // - // This implies the following differential equation: - // - // d_t var = - v^i d_i var - v0 (var - var0) / r - // + v0 (1 - n) h / r^n+1 + d_t h / r^n - // ~ - v^i d_i var - v0 (var - var0) / r - // + d_t h / r^n + O(1/r^n+1) - // - // where vi = v0 xi/r - // - // the Coulomb term, d_t h / r ^n , is estimated by extrapolating - // from the nearest interior points - // - // Displacement to get from p.I to interior point placed - // nghostpoints away - vect displacement{grid.nghostzones[0] * p.NI[0], - grid.nghostzones[1] * p.NI[1], - grid.nghostzones[2] * p.NI[2]}; - vect intp = p.I - displacement; - - assert(intp[0] >= grid.nghostzones[0]); - assert(intp[1] >= grid.nghostzones[1]); - assert(intp[2] >= grid.nghostzones[2]); - assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1); - assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1); - assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1); - - // coordinates at p.I-displacement - const CCTK_REAL xint = p.x - displacement[0] * p.DX[0]; - const CCTK_REAL yint = p.y - displacement[1] * p.DX[1]; - const CCTK_REAL zint = p.z - displacement[2] * p.DX[2]; - const CCTK_REAL rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint)); - - // Find local wave speeds at physical point p.I-displacement - const CCTK_REAL vxint = v0 * xint / rint; - const CCTK_REAL vyint = v0 * yint / rint; - const CCTK_REAL vzint = v0 * zint / rint; - - // Derivatives at physical point p.I-displacement - const CCTK_REAL varxint = - (var(intp + p.DI[0]) - var(intp - p.DI[0])) / (2 * p.DX[0]); - const CCTK_REAL varyint = - (var(intp + p.DI[1]) - var(intp - p.DI[1])) / (2 * p.DX[1]); - const CCTK_REAL varzint = - (var(intp + p.DI[2]) - var(intp - p.DI[2])) / (2 * p.DX[2]); - - // extrapolate Coulomb component, rescale to account for radial - // fall-off - CCTK_REAL rad = -vxint * varxint - vyint * varyint - - vzint * varzint - v0 * (var(intp) - var0) / rint; - CCTK_REAL aux = (rhs(intp) - rad) * pow(rint / r, radpower); - - // radiative rhs with extrapolated Coulomb correction - rhs(p.I) += aux; - } - }); + // Radiative rhs with extrapolated Coulomb correction + rhs(p.I) += aux; + } + } + }); } } // namespace NewRadX diff --git a/NewRadX/src/newradx.hxx b/NewRadX/src/newradx.hxx index 657e310e..551606e2 100644 --- a/NewRadX/src/newradx.hxx +++ b/NewRadX/src/newradx.hxx @@ -1,17 +1,52 @@ #include -#include -#include -#include "loop.hxx" -#include "loop_device.hxx" +#include namespace NewRadX { -// Adapted from BSSN_MoL's files NewRad.F and newrad.h. This code was probably -// originally written by Miguel Alcubierre. +/** + * @brief Applies radiative boundary condition to the RHS of a state variable + * + * Adapted from NewRad thorn by E. Schnetter, used with Carpet. Original code + * adapted from BSSN_MoL's files NewRad.F and newrad.h. This code was probably + * originally written by Miguel Alcubierre. + * + * @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 var0 Value at infinity. + * @param v0 Propagation speed. + * @param radpower Radial fall-off exponent + */ void NewRadX_Apply(const cGH *restrict const cctkGH, - const Loop::GF3D2 var, Loop::GF3D2 rhs, - const CCTK_REAL var0, //!< value at infinity - const CCTK_REAL v0, //!< propagation speed - const CCTK_REAL radpower //!< radial fall-off exponent -); + const Loop::GF3D2 &var, + const Loop::GF3D2 &rhs, const CCTK_REAL var0, + const CCTK_REAL v0, const CCTK_REAL radpower); + +/** + * @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 + */ +void NewRadX_Apply(const cGH *restrict const cctkGH, + const Loop::GF3D2 &var, + const Loop::GF3D2 &rhs, + const Loop::GF3D2 &vcoordx, + const Loop::GF3D2 &vcoordy, + const Loop::GF3D2 &vcoordz, + const CCTK_REAL var0, const CCTK_REAL v0, + const CCTK_REAL radpower); + } // namespace NewRadX From cc0df575467707dc6d3fa75ba18f4e2ed718104b Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Wed, 1 Oct 2025 00:36:48 -0500 Subject: [PATCH 2/8] NewRadX: Use multipatch projections --- NewRadX/interface.ccl | 1 + NewRadX/src/newradx.cxx | 154 +++++++++++----------------------------- NewRadX/src/newradx.hxx | 13 ++++ 3 files changed, 56 insertions(+), 112 deletions(-) diff --git a/NewRadX/interface.ccl b/NewRadX/interface.ccl index 7a62e44f..7d4329f4 100644 --- a/NewRadX/interface.ccl +++ b/NewRadX/interface.ccl @@ -4,5 +4,6 @@ IMPLEMENTS: NewRadX USES INCLUDE HEADER: driver.hxx USES INCLUDE HEADER: loop_device.hxx +USES INCLUDE HEADER: global_derivatives.hxx INCLUDE HEADER: newradx.hxx in newradx.hxx diff --git a/NewRadX/src/newradx.cxx b/NewRadX/src/newradx.cxx index 22a790d9..d04dada2 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -3,6 +3,8 @@ #include #include +#include + #include "newradx.hxx" #include @@ -77,31 +79,6 @@ static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL calc_deriv( } } -template -static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL -calc_deriv_interior_upper(const Loop::PointDesc &p, - const Loop::GF3D2 &gf) noexcept { - if (p.NI[dir] == 0) { - /* interior - * in direction parallel to face/edge, apply symmetric stencil - * currently second-order accurate - * TODO: apply same finite-difference order as interior - */ - return c2o(p, p.I, gf); - - } else if (p.NI[dir] == +1) { - /* upper boundary - * in direction normal to face/edge/corner, apply asymmetric stencil - * currently second-order accurate - * TODO: apply same finite-difference order as interior - */ - return r2o(p, p.I, gf); - - } else { - assert(0); - } -} - void NewRadX_Apply(const cGH *restrict const cctkGH, const Loop::GF3D2 &var, const Loop::GF3D2 &rhs, const CCTK_REAL var0, @@ -216,8 +193,19 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const Loop::GF3D2 &vcoordx, const Loop::GF3D2 &vcoordy, const Loop::GF3D2 &vcoordz, + const Loop::GF3D2 &vJ_da_dx, + const Loop::GF3D2 &vJ_da_dy, + const Loop::GF3D2 &vJ_da_dz, + const Loop::GF3D2 &vJ_db_dx, + const Loop::GF3D2 &vJ_db_dy, + const Loop::GF3D2 &vJ_db_dz, + const Loop::GF3D2 &vJ_dc_dx, + const Loop::GF3D2 &vJ_dc_dy, + const Loop::GF3D2 &vJ_dc_dz, const CCTK_REAL var0, const CCTK_REAL v0, const CCTK_REAL radpower) { + using namespace CapyrX::MultiPatch::GlobalDerivatives; + DECLARE_CCTK_ARGUMENTS; const auto symmetries = CarpetX::ghext->patchdata.at(cctk_patch).symmetries; @@ -233,93 +221,35 @@ 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 { - if (p.patch != 0) { - /* - * The main part of the boundary condition assumes that we have an - * outgoing radial wave with some speed v0: - * - * var = var0 + u(r-v0*t)/r - * - * This implies the following differential equation: - * - * d_t var = - v^i d_i var - v0 (var - var0) / r - * - * where vi = v0 xi/r - */ - - // coordinate radius at 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)); - - // Find local wave speeds at radiative boundary point p.I - const auto vx = v0 * x / r; - const auto vy = v0 * y / r; - const auto vz = v0 * z / r; - const auto vr = sqrt(pow2(vx) + pow2(vy) + pow2(vz)); - - // Derivatives - const auto derivz = calc_deriv_interior_upper<2>(p, var); - - // Radiative rhs - rhs(p.I) = -vr * derivz - v0 * (var(p.I) - var0) / r; - - if (radpower > 0.0) { - // When solution is known to have a Coulomb-like component - // asymptotically, this is estimated and extrapolated from - // physical interior points to the radiative boundary. I.e.: - // - // var = var0 + u(r-v0*t)/r + h(t)/r^n - // - // This implies the following differential equation: - // - // d_t var = - v^i d_i var - v0 (var - var0) / r - // + v0 (1 - n) h / r^n+1 + d_t h / r^n - // ~ - v^i d_i var - v0 (var - var0) / r - // + d_t h / r^n + O(1/r^n+1) - // - // where vi = v0 xi/r - // - // the Coulomb term, d_t h / r ^n , is estimated by extrapolating - // from the nearest interior points - // - // Displacement to get from p.I to interior point placed - // nghostpoints away - const vect displacement{0, 0, grid.nghostzones[2]}; - const vect intp = p.I - displacement; - - assert(intp[0] >= grid.nghostzones[0]); - assert(intp[1] >= grid.nghostzones[1]); - assert(intp[2] >= grid.nghostzones[2]); - assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1); - assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1); - assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1); - - // coordinates at p.I-displacement - const auto xint = vcoordx(intp); - const auto yint = vcoordy(intp); - const auto zint = vcoordz(intp); - const auto rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint)); - - // Find local wave speeds at physical point p.I-displacement - const auto vxint = v0 * xint / rint; - const auto vyint = v0 * yint / rint; - const auto vzint = v0 * zint / rint; - const auto vrint = sqrt(pow2(vxint) + pow2(vyint) + pow2(vzint)); - - // Derivative at physical point p.I-displacement - const auto derivzint = c2o<2>(p, intp, var); - - // Extrapolate Coulomb component, rescale to account for radial - // fall-off - const auto rad = - -vrint * derivzint - v0 * (var(intp) - var0) / rint; - const auto aux = (rhs(intp) - rad) * pow(rint / r, radpower); - - // Radiative rhs with extrapolated Coulomb correction - rhs(p.I) += aux; - } + // Make sure we are always on the outer boundary of a patch system + assert(p.patch != 0); + assert(p.NI[2] >= 0); + + // 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; + + // 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)}; + + // 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; + + if (radpower > 0.0) { + // TODO } }); } diff --git a/NewRadX/src/newradx.hxx b/NewRadX/src/newradx.hxx index 551606e2..5f694cf9 100644 --- a/NewRadX/src/newradx.hxx +++ b/NewRadX/src/newradx.hxx @@ -22,6 +22,10 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const Loop::GF3D2 &rhs, const CCTK_REAL var0, const CCTK_REAL v0, const CCTK_REAL radpower); +#define NEWRADX_MULTIPATCH_QUANTITIES \ + vcoordx, vcoordy, vcoordz, vJ_da_dx, vJ_da_dy, vJ_da_dz, vJ_db_dx, vJ_db_dy, \ + 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: @@ -46,6 +50,15 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const Loop::GF3D2 &vcoordx, const Loop::GF3D2 &vcoordy, const Loop::GF3D2 &vcoordz, + const Loop::GF3D2 &vJ_da_dx, + const Loop::GF3D2 &vJ_da_dy, + const Loop::GF3D2 &vJ_da_dz, + const Loop::GF3D2 &vJ_db_dx, + const Loop::GF3D2 &vJ_db_dy, + const Loop::GF3D2 &vJ_db_dz, + const Loop::GF3D2 &vJ_dc_dx, + const Loop::GF3D2 &vJ_dc_dy, + const Loop::GF3D2 &vJ_dc_dz, const CCTK_REAL var0, const CCTK_REAL v0, const CCTK_REAL radpower); From dc2a421e07f331ee5c725755625f364df1c1f293 Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Wed, 11 Feb 2026 23:39:07 -0600 Subject: [PATCH 3/8] NewRadX: Relaxed requirements of specific patch system setup --- NewRadX/param.ccl | 4 ++++ NewRadX/src/newradx.cxx | 32 ++++++++++++++++++++------------ NewRadX/src/newradx.hxx | 23 +++++++---------------- 3 files changed, 31 insertions(+), 28 deletions(-) diff --git a/NewRadX/param.ccl b/NewRadX/param.ccl index 7ba1a5f0..7f70ee40 100644 --- a/NewRadX/param.ccl +++ b/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" diff --git a/NewRadX/src/newradx.cxx b/NewRadX/src/newradx.cxx index d04dada2..db8b0961 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -1,4 +1,5 @@ #include +#include #include #include @@ -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, 2> is_sym_bnd{ @@ -221,9 +223,15 @@ 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); @@ -231,14 +239,14 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, 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)}; @@ -246,10 +254,10 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, // 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) } }); } diff --git a/NewRadX/src/newradx.hxx b/NewRadX/src/newradx.hxx index 5f694cf9..a4986c75 100644 --- a/NewRadX/src/newradx.hxx +++ b/NewRadX/src/newradx.hxx @@ -1,3 +1,6 @@ +#ifndef NEWRADX_HXX +#define NEWRADX_HXX + #include #include @@ -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 &var, @@ -63,3 +52,5 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const CCTK_REAL radpower); } // namespace NewRadX + +#endif // NEWRADX_HXX \ No newline at end of file From 353ca7df25f713ee055931fdba18a211b006fb7a Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Thu, 12 Feb 2026 00:43:06 -0600 Subject: [PATCH 4/8] NewRadX: Added radpower correction for multipatch --- NewRadX/src/newradx.cxx | 58 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 56 insertions(+), 2 deletions(-) diff --git a/NewRadX/src/newradx.cxx b/NewRadX/src/newradx.cxx index db8b0961..442aceae 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -256,8 +256,62 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, rhs(p.I) = -vx * g_dvar.dx - vy * g_dvar.dy - vz * g_dvar.dz - sv0 * (var(p.I) - var0) / r; - if (radpower > 0.0) { - // TODO: Coulomb correction (port from Cartesian overload) + // Coulomb correction: estimate and extrapolate the 1/r^n + // component from interior points to the radiative boundary. + // Only applied at the outer boundary — at the inner boundary + // the asymptotic 1/r^n assumption breaks down (r is small, + // rint > r, so the rescaling amplifies rather than attenuates). + if (radpower > 0.0 && p.NI[2] != -1) { + // Displacement to get from p.I to interior point placed + // nghostpoints away + const vect displacement{grid.nghostzones[0] * p.NI[0], + grid.nghostzones[1] * p.NI[1], + grid.nghostzones[2] * p.NI[2]}; + const vect intp = p.I - displacement; + + assert(intp[0] >= grid.nghostzones[0]); + assert(intp[1] >= grid.nghostzones[1]); + assert(intp[2] >= grid.nghostzones[2]); + assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1); + assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1); + assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1); + + // Global coordinates at interior point + const auto xint = vcoordx(intp); + const auto yint = vcoordy(intp); + const auto zint = vcoordz(intp); + const auto rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint)); + + // Find local wave speeds at interior point + const auto vxint = sv0 * xint / rint; + const auto vyint = sv0 * yint / rint; + const auto vzint = sv0 * zint / rint; + + // Local derivatives at interior point (centered stencils) + const LocalFirstDerivatives l_dvar_int{.da = c2o<0>(p, intp, var), + .db = c2o<1>(p, intp, var), + .dc = c2o<2>(p, intp, var)}; + + // Jacobians at interior point + const Jacobians jac_int{ + vJ_da_dx(intp), vJ_da_dy(intp), vJ_da_dz(intp), + vJ_db_dx(intp), vJ_db_dy(intp), vJ_db_dz(intp), + vJ_dc_dx(intp), vJ_dc_dy(intp), vJ_dc_dz(intp)}; + + // Global derivatives at interior point + const auto g_dvar_int{project_first(l_dvar_int, jac_int)}; + + // Radiative part at interior point + const auto rad = -vxint * g_dvar_int.dx - vyint * g_dvar_int.dy - + vzint * g_dvar_int.dz - + sv0 * (var(intp) - var0) / rint; + + // Extrapolate Coulomb component, rescale to account for radial + // fall-off + const auto aux = (rhs(intp) - rad) * pow(rint / r, radpower); + + // Radiative rhs with extrapolated Coulomb correction + rhs(p.I) += aux; } }); } From 8097aed0c5b18785dea28fd23900242b790c0d23 Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Thu, 12 Feb 2026 14:03:39 -0600 Subject: [PATCH 5/8] Add CapyrX thorns to spacetimex.th Added CapyrX thorns configuration to spacetimex.th --- scripts/spacetimex.th | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/scripts/spacetimex.th b/scripts/spacetimex.th index 11cfe2a4..1d2e8ebe 100644 --- a/scripts/spacetimex.th +++ b/scripts/spacetimex.th @@ -161,6 +161,18 @@ CarpetX/Loop CarpetX/ODESolvers CarpetX/TmunuBaseX +# CapyrX thorns +!TARGET = $ARR +!TYPE = git +!URL = https://github.com/lucass-carneiro/CapyrX +!REPO_PATH= $2 +!CHECKOUT = +CapyrX_GlobalDerivatives +CapyrX_MultiPatch +CapyrX_TestMultiPatch +CapyrX_WaveToy +CapyrX_WaveToy_SecondOrder + # SpacetimeX thorns !TARGET = $ARR !TYPE = ignore From c5c8c0ded0636a99fd2b4ac5b0542a8450d6343e Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Thu, 12 Feb 2026 15:12:16 -0600 Subject: [PATCH 6/8] Update repository paths for CapyrX thorns --- scripts/spacetimex.th | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/spacetimex.th b/scripts/spacetimex.th index 1d2e8ebe..c5dfd41b 100644 --- a/scripts/spacetimex.th +++ b/scripts/spacetimex.th @@ -167,11 +167,11 @@ CarpetX/TmunuBaseX !URL = https://github.com/lucass-carneiro/CapyrX !REPO_PATH= $2 !CHECKOUT = -CapyrX_GlobalDerivatives -CapyrX_MultiPatch -CapyrX_TestMultiPatch -CapyrX_WaveToy -CapyrX_WaveToy_SecondOrder +CapyrX/CapyrX_GlobalDerivatives +CapyrX/CapyrX_MultiPatch +CapyrX/CapyrX_TestMultiPatch +CapyrX/CapyrX_WaveToy +CapyrX/CapyrX_WaveToy_SecondOrder # SpacetimeX thorns !TARGET = $ARR From cac008af18a697ab223b7704a5543d774b703e8d Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Thu, 19 Feb 2026 09:05:34 -0600 Subject: [PATCH 7/8] NewRadX: Removed explicit references to CapyrX --- NewRadX/interface.ccl | 1 - NewRadX/src/newradx.cxx | 58 +++++++++++++++++++++++++---------------- 2 files changed, 36 insertions(+), 23 deletions(-) diff --git a/NewRadX/interface.ccl b/NewRadX/interface.ccl index 7d4329f4..7a62e44f 100644 --- a/NewRadX/interface.ccl +++ b/NewRadX/interface.ccl @@ -4,6 +4,5 @@ IMPLEMENTS: NewRadX USES INCLUDE HEADER: driver.hxx USES INCLUDE HEADER: loop_device.hxx -USES INCLUDE HEADER: global_derivatives.hxx INCLUDE HEADER: newradx.hxx in newradx.hxx diff --git a/NewRadX/src/newradx.cxx b/NewRadX/src/newradx.cxx index 442aceae..601ce614 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -4,8 +4,6 @@ #include #include -#include - #include "newradx.hxx" #include @@ -205,8 +203,6 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const Loop::GF3D2 &vJ_dc_dz, const CCTK_REAL var0, const CCTK_REAL v0, const CCTK_REAL radpower) { - using namespace CapyrX::MultiPatch::GlobalDerivatives; - DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; @@ -243,17 +239,29 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const auto vy = sv0 * y / r; const auto vz = sv0 * z / r; - // 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)}; + // Local derivatives + const auto dgf_da = calc_deriv<0>(p, var); + const auto dgf_db = calc_deriv<1>(p, var); + const auto dgf_dc = calc_deriv<2>(p, var); + + // Jacobians + const auto da_dx = vJ_da_dx(p.I); + const auto da_dy = vJ_da_dy(p.I); + const auto da_dz = vJ_da_dz(p.I); + const auto db_dx = vJ_db_dx(p.I); + const auto db_dy = vJ_db_dy(p.I); + const auto db_dz = vJ_db_dz(p.I); + const auto dc_dx = vJ_dc_dx(p.I); + const auto dc_dy = vJ_dc_dy(p.I); + const auto dc_dz = vJ_dc_dz(p.I); // Global derivatives - const Jacobians jac{VERTEX_JACOBIANS(p)}; - const auto g_dvar{project_first(l_dvar, jac)}; + const auto dgf_dx = dgf_db * db_dx + dgf_dc * dc_dx + da_dx * dgf_da; + const auto dgf_dy = dgf_dc * dc_dy + db_dy * dgf_db + da_dy * dgf_da; + const auto dgf_dz = dc_dz * dgf_dc + db_dz * dgf_db + da_dz * dgf_da; // radiative rhs - rhs(p.I) = -vx * g_dvar.dx - vy * g_dvar.dy - vz * g_dvar.dz - + rhs(p.I) = -vx * dgf_dx - vy * dgf_dy - vz * dgf_dz - sv0 * (var(p.I) - var0) / r; // Coulomb correction: estimate and extrapolate the 1/r^n @@ -287,23 +295,29 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, const auto vyint = sv0 * yint / rint; const auto vzint = sv0 * zint / rint; - // Local derivatives at interior point (centered stencils) - const LocalFirstDerivatives l_dvar_int{.da = c2o<0>(p, intp, var), - .db = c2o<1>(p, intp, var), - .dc = c2o<2>(p, intp, var)}; + // Local derivatives at interior point + const auto dgf_da = c2o<0>(p, intp, var); + const auto dgf_db = c2o<1>(p, intp, var); + const auto dgf_dc = c2o<2>(p, intp, var); // Jacobians at interior point - const Jacobians jac_int{ - vJ_da_dx(intp), vJ_da_dy(intp), vJ_da_dz(intp), - vJ_db_dx(intp), vJ_db_dy(intp), vJ_db_dz(intp), - vJ_dc_dx(intp), vJ_dc_dy(intp), vJ_dc_dz(intp)}; + const auto da_dx = vJ_da_dx(intp); + const auto da_dy = vJ_da_dy(intp); + const auto da_dz = vJ_da_dz(intp); + const auto db_dx = vJ_db_dx(intp); + const auto db_dy = vJ_db_dy(intp); + const auto db_dz = vJ_db_dz(intp); + const auto dc_dx = vJ_dc_dx(intp); + const auto dc_dy = vJ_dc_dy(intp); + const auto dc_dz = vJ_dc_dz(intp); // Global derivatives at interior point - const auto g_dvar_int{project_first(l_dvar_int, jac_int)}; + const auto dgf_dx = dgf_db * db_dx + dgf_dc * dc_dx + da_dx * dgf_da; + const auto dgf_dy = dgf_dc * dc_dy + db_dy * dgf_db + da_dy * dgf_da; + const auto dgf_dz = dc_dz * dgf_dc + db_dz * dgf_db + da_dz * dgf_da; // Radiative part at interior point - const auto rad = -vxint * g_dvar_int.dx - vyint * g_dvar_int.dy - - vzint * g_dvar_int.dz - + const auto rad = -vxint * dgf_dx - vyint * dgf_dy - vzint * dgf_dz - sv0 * (var(intp) - var0) / rint; // Extrapolate Coulomb component, rescale to account for radial From 0bed7590f2f63c886be5f5e7017d9238e392936a Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Thu, 19 Feb 2026 10:03:44 -0600 Subject: [PATCH 8/8] Updated CI thornlist --- scripts/spacetimex.th | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/scripts/spacetimex.th b/scripts/spacetimex.th index 3e368b7d..742abba7 100644 --- a/scripts/spacetimex.th +++ b/scripts/spacetimex.th @@ -161,18 +161,6 @@ CarpetX/Loop CarpetX/ODESolvers CarpetX/TmunuBaseX -# CapyrX thorns -!TARGET = $ARR -!TYPE = git -!URL = https://github.com/lucass-carneiro/CapyrX -!REPO_PATH= $2 -!CHECKOUT = -CapyrX/CapyrX_GlobalDerivatives -CapyrX/CapyrX_MultiPatch -CapyrX/CapyrX_TestMultiPatch -CapyrX/CapyrX_WaveToy -CapyrX/CapyrX_WaveToy_SecondOrder - # SpacetimeX thorns !TARGET = $ARR !TYPE = ignore @@ -194,3 +182,4 @@ SpacetimeX/TwoPuncturesX SpacetimeX/Weyl SpacetimeX/WeylScal4 SpacetimeX/Z4c +