Skip to content

Commit cc0df57

Browse files
lucass-carneirolucass-carneiro
committed
NewRadX: Use multipatch projections
1 parent a1931d7 commit cc0df57

File tree

3 files changed

+56
-112
lines changed

3 files changed

+56
-112
lines changed

‎NewRadX/interface.ccl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,5 +4,6 @@ IMPLEMENTS: NewRadX
44

55
USES INCLUDE HEADER: driver.hxx
66
USES INCLUDE HEADER: loop_device.hxx
7+
USES INCLUDE HEADER: global_derivatives.hxx
78

89
INCLUDE HEADER: newradx.hxx in newradx.hxx

‎NewRadX/src/newradx.cxx‎

Lines changed: 42 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@
33
#include <loop_device.hxx>
44
#include <driver.hxx>
55

6+
#include <global_derivatives.hxx>
7+
68
#include "newradx.hxx"
79

810
#include <cmath>
@@ -77,31 +79,6 @@ static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL calc_deriv(
7779
}
7880
}
7981

80-
template <std::size_t dir>
81-
static inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_REAL
82-
calc_deriv_interior_upper(const Loop::PointDesc &p,
83-
const Loop::GF3D2<const CCTK_REAL> &gf) noexcept {
84-
if (p.NI[dir] == 0) {
85-
/* interior
86-
* in direction parallel to face/edge, apply symmetric stencil
87-
* currently second-order accurate
88-
* TODO: apply same finite-difference order as interior
89-
*/
90-
return c2o<dir>(p, p.I, gf);
91-
92-
} else if (p.NI[dir] == +1) {
93-
/* upper boundary
94-
* in direction normal to face/edge/corner, apply asymmetric stencil
95-
* currently second-order accurate
96-
* TODO: apply same finite-difference order as interior
97-
*/
98-
return r2o<dir>(p, p.I, gf);
99-
100-
} else {
101-
assert(0);
102-
}
103-
}
104-
10582
void NewRadX_Apply(const cGH *restrict const cctkGH,
10683
const Loop::GF3D2<const CCTK_REAL> &var,
10784
const Loop::GF3D2<CCTK_REAL> &rhs, const CCTK_REAL var0,
@@ -216,8 +193,19 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
216193
const Loop::GF3D2<const CCTK_REAL> &vcoordx,
217194
const Loop::GF3D2<const CCTK_REAL> &vcoordy,
218195
const Loop::GF3D2<const CCTK_REAL> &vcoordz,
196+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dx,
197+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dy,
198+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dz,
199+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dx,
200+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dy,
201+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dz,
202+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dx,
203+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dy,
204+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dz,
219205
const CCTK_REAL var0, const CCTK_REAL v0,
220206
const CCTK_REAL radpower) {
207+
using namespace CapyrX::MultiPatch::GlobalDerivatives;
208+
221209
DECLARE_CCTK_ARGUMENTS;
222210

223211
const auto symmetries = CarpetX::ghext->patchdata.at(cctk_patch).symmetries;
@@ -233,93 +221,35 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
233221
grid.loop_outermost_int_device<0, 0, 0>(
234222
grid.nghostzones, is_sym_bnd,
235223
[=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
236-
if (p.patch != 0) {
237-
/*
238-
* The main part of the boundary condition assumes that we have an
239-
* outgoing radial wave with some speed v0:
240-
*
241-
* var = var0 + u(r-v0*t)/r
242-
*
243-
* This implies the following differential equation:
244-
*
245-
* d_t var = - v^i d_i var - v0 (var - var0) / r
246-
*
247-
* where vi = v0 xi/r
248-
*/
249-
250-
// coordinate radius at p.I
251-
const auto x = vcoordx(p.I);
252-
const auto y = vcoordy(p.I);
253-
const auto z = vcoordz(p.I);
254-
const auto r = sqrt(pow2(x) + pow2(y) + pow2(z));
255-
256-
// Find local wave speeds at radiative boundary point p.I
257-
const auto vx = v0 * x / r;
258-
const auto vy = v0 * y / r;
259-
const auto vz = v0 * z / r;
260-
const auto vr = sqrt(pow2(vx) + pow2(vy) + pow2(vz));
261-
262-
// Derivatives
263-
const auto derivz = calc_deriv_interior_upper<2>(p, var);
264-
265-
// Radiative rhs
266-
rhs(p.I) = -vr * derivz - v0 * (var(p.I) - var0) / r;
267-
268-
if (radpower > 0.0) {
269-
// When solution is known to have a Coulomb-like component
270-
// asymptotically, this is estimated and extrapolated from
271-
// physical interior points to the radiative boundary. I.e.:
272-
//
273-
// var = var0 + u(r-v0*t)/r + h(t)/r^n
274-
//
275-
// This implies the following differential equation:
276-
//
277-
// d_t var = - v^i d_i var - v0 (var - var0) / r
278-
// + v0 (1 - n) h / r^n+1 + d_t h / r^n
279-
// ~ - v^i d_i var - v0 (var - var0) / r
280-
// + d_t h / r^n + O(1/r^n+1)
281-
//
282-
// where vi = v0 xi/r
283-
//
284-
// the Coulomb term, d_t h / r ^n , is estimated by extrapolating
285-
// from the nearest interior points
286-
//
287-
// Displacement to get from p.I to interior point placed
288-
// nghostpoints away
289-
const vect<int, dim> displacement{0, 0, grid.nghostzones[2]};
290-
const vect<int, dim> intp = p.I - displacement;
291-
292-
assert(intp[0] >= grid.nghostzones[0]);
293-
assert(intp[1] >= grid.nghostzones[1]);
294-
assert(intp[2] >= grid.nghostzones[2]);
295-
assert(intp[0] <= grid.lsh[0] - grid.nghostzones[0] - 1);
296-
assert(intp[1] <= grid.lsh[1] - grid.nghostzones[1] - 1);
297-
assert(intp[2] <= grid.lsh[2] - grid.nghostzones[2] - 1);
298-
299-
// coordinates at p.I-displacement
300-
const auto xint = vcoordx(intp);
301-
const auto yint = vcoordy(intp);
302-
const auto zint = vcoordz(intp);
303-
const auto rint = sqrt(pow2(xint) + pow2(yint) + pow2(zint));
304-
305-
// Find local wave speeds at physical point p.I-displacement
306-
const auto vxint = v0 * xint / rint;
307-
const auto vyint = v0 * yint / rint;
308-
const auto vzint = v0 * zint / rint;
309-
const auto vrint = sqrt(pow2(vxint) + pow2(vyint) + pow2(vzint));
310-
311-
// Derivative at physical point p.I-displacement
312-
const auto derivzint = c2o<2>(p, intp, var);
313-
314-
// Extrapolate Coulomb component, rescale to account for radial
315-
// fall-off
316-
const auto rad =
317-
-vrint * derivzint - v0 * (var(intp) - var0) / rint;
318-
const auto aux = (rhs(intp) - rad) * pow(rint / r, radpower);
319-
320-
// Radiative rhs with extrapolated Coulomb correction
321-
rhs(p.I) += aux;
322-
}
224+
// Make sure we are always on the outer boundary of a patch system
225+
assert(p.patch != 0);
226+
assert(p.NI[2] >= 0);
227+
228+
// Find local wave speeds at radiative boundary point p.I
229+
const auto x = vcoordx(p.I);
230+
const auto y = vcoordy(p.I);
231+
const auto z = vcoordz(p.I);
232+
const auto r = sqrt(pow2(x) + pow2(y) + pow2(z));
233+
234+
const auto vx = v0 * vcoordx(p.I) / r;
235+
const auto vy = v0 * vcoordy(p.I) / r;
236+
const auto vz = v0 * vcoordz(p.I) / r;
237+
238+
// Local derivatives
239+
const LocalFirstDerivatives l_dvar{.da = r2o<0>(p, p.I, var),
240+
.db = r2o<1>(p, p.I, var),
241+
.dc = r2o<2>(p, p.I, var)};
242+
243+
// Global derivatives
244+
const Jacobians jac{VERTEX_JACOBIANS(p)};
245+
const auto g_dvar{project_first(l_dvar, jac)};
246+
247+
// radiative rhs
248+
rhs(p.I) = -vx * g_dvar.dx - vy * g_dvar.dy - vz * g_dvar.dz -
249+
v0 * (var(p.I) - var0) / r;
250+
251+
if (radpower > 0.0) {
252+
// TODO
323253
}
324254
});
325255
}

‎NewRadX/src/newradx.hxx‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,10 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
2222
const Loop::GF3D2<CCTK_REAL> &rhs, const CCTK_REAL var0,
2323
const CCTK_REAL v0, const CCTK_REAL radpower);
2424

25+
#define NEWRADX_MULTIPATCH_QUANTITIES \
26+
vcoordx, vcoordy, vcoordz, vJ_da_dx, vJ_da_dy, vJ_da_dz, vJ_db_dx, vJ_db_dy, \
27+
vJ_db_dz, vJ_dc_dx, vJ_dc_dy, vJ_dc_dz
28+
2529
/**
2630
* @brief Applies radiative boundary condition to the RHS of a state variable.
2731
* Assumes that:
@@ -46,6 +50,15 @@ void NewRadX_Apply(const cGH *restrict const cctkGH,
4650
const Loop::GF3D2<const CCTK_REAL> &vcoordx,
4751
const Loop::GF3D2<const CCTK_REAL> &vcoordy,
4852
const Loop::GF3D2<const CCTK_REAL> &vcoordz,
53+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dx,
54+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dy,
55+
const Loop::GF3D2<const CCTK_REAL> &vJ_da_dz,
56+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dx,
57+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dy,
58+
const Loop::GF3D2<const CCTK_REAL> &vJ_db_dz,
59+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dx,
60+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dy,
61+
const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dz,
4962
const CCTK_REAL var0, const CCTK_REAL v0,
5063
const CCTK_REAL radpower);
5164

0 commit comments

Comments
 (0)