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 3134f8cc..601ce614 100644 --- a/NewRadX/src/newradx.cxx +++ b/NewRadX/src/newradx.cxx @@ -1,225 +1,333 @@ #include -#include -#include -#include "loop.hxx" -#include "loop_device.hxx" -#include "driver.hxx" +#include + +#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; +} + +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; +} + +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); + } +} -// 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 -) { + const Loop::GF3D2 &var, + const Loop::GF3D2 &rhs, const CCTK_REAL var0, + const CCTK_REAL v0, const CCTK_REAL radpower) { DECLARE_CCTK_ARGUMENTS; - 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); - 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); - assert(0); - }; + 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 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) { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; 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: - // - // 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 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; - } - }); + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // 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 = sv0 * x / r; + const auto vy = sv0 * y / r; + const auto vz = sv0 * z / r; + + // 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 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 * dgf_dx - vy * dgf_dy - vz * dgf_dz - + sv0 * (var(p.I) - var0) / r; + + // 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 + 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 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 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 * dgf_dx - vyint * dgf_dy - vzint * dgf_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; + } + }); } } // namespace NewRadX diff --git a/NewRadX/src/newradx.hxx b/NewRadX/src/newradx.hxx index 657e310e..a4986c75 100644 --- a/NewRadX/src/newradx.hxx +++ b/NewRadX/src/newradx.hxx @@ -1,17 +1,56 @@ +#ifndef NEWRADX_HXX +#define NEWRADX_HXX + #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, + 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 + * on multipatch grids. + */ 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 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); + } // namespace NewRadX + +#endif // NEWRADX_HXX \ No newline at end of file diff --git a/scripts/spacetimex.th b/scripts/spacetimex.th index 3625e66c..742abba7 100644 --- a/scripts/spacetimex.th +++ b/scripts/spacetimex.th @@ -182,3 +182,4 @@ SpacetimeX/TwoPuncturesX SpacetimeX/Weyl SpacetimeX/WeylScal4 SpacetimeX/Z4c +