|
4 | 4 | #include <loop_device.hxx> |
5 | 5 | #include <driver.hxx> |
6 | 6 |
|
7 | | -#include <global_derivatives.hxx> |
8 | | - |
9 | 7 | #include "newradx.hxx" |
10 | 8 |
|
11 | 9 | #include <cmath> |
@@ -205,8 +203,6 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, |
205 | 203 | const Loop::GF3D2<const CCTK_REAL> &vJ_dc_dz, |
206 | 204 | const CCTK_REAL var0, const CCTK_REAL v0, |
207 | 205 | const CCTK_REAL radpower) { |
208 | | - using namespace CapyrX::MultiPatch::GlobalDerivatives; |
209 | | - |
210 | 206 | DECLARE_CCTK_ARGUMENTS; |
211 | 207 | DECLARE_CCTK_PARAMETERS; |
212 | 208 |
|
@@ -243,17 +239,29 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, |
243 | 239 | const auto vy = sv0 * y / r; |
244 | 240 | const auto vz = sv0 * z / r; |
245 | 241 |
|
246 | | - // Local derivatives (stencil auto-selected by boundary location) |
247 | | - const LocalFirstDerivatives l_dvar{.da = calc_deriv<0>(p, var), |
248 | | - .db = calc_deriv<1>(p, var), |
249 | | - .dc = calc_deriv<2>(p, var)}; |
| 242 | + // Local derivatives |
| 243 | + const auto dgf_da = calc_deriv<0>(p, var); |
| 244 | + const auto dgf_db = calc_deriv<1>(p, var); |
| 245 | + const auto dgf_dc = calc_deriv<2>(p, var); |
| 246 | + |
| 247 | + // Jacobians |
| 248 | + const auto da_dx = vJ_da_dx(p.I); |
| 249 | + const auto da_dy = vJ_da_dy(p.I); |
| 250 | + const auto da_dz = vJ_da_dz(p.I); |
| 251 | + const auto db_dx = vJ_db_dx(p.I); |
| 252 | + const auto db_dy = vJ_db_dy(p.I); |
| 253 | + const auto db_dz = vJ_db_dz(p.I); |
| 254 | + const auto dc_dx = vJ_dc_dx(p.I); |
| 255 | + const auto dc_dy = vJ_dc_dy(p.I); |
| 256 | + const auto dc_dz = vJ_dc_dz(p.I); |
250 | 257 |
|
251 | 258 | // Global derivatives |
252 | | - const Jacobians jac{VERTEX_JACOBIANS(p)}; |
253 | | - const auto g_dvar{project_first(l_dvar, jac)}; |
| 259 | + const auto dgf_dx = dgf_db * db_dx + dgf_dc * dc_dx + da_dx * dgf_da; |
| 260 | + const auto dgf_dy = dgf_dc * dc_dy + db_dy * dgf_db + da_dy * dgf_da; |
| 261 | + const auto dgf_dz = dc_dz * dgf_dc + db_dz * dgf_db + da_dz * dgf_da; |
254 | 262 |
|
255 | 263 | // radiative rhs |
256 | | - rhs(p.I) = -vx * g_dvar.dx - vy * g_dvar.dy - vz * g_dvar.dz - |
| 264 | + rhs(p.I) = -vx * dgf_dx - vy * dgf_dy - vz * dgf_dz - |
257 | 265 | sv0 * (var(p.I) - var0) / r; |
258 | 266 |
|
259 | 267 | // Coulomb correction: estimate and extrapolate the 1/r^n |
@@ -287,23 +295,29 @@ void NewRadX_Apply(const cGH *restrict const cctkGH, |
287 | 295 | const auto vyint = sv0 * yint / rint; |
288 | 296 | const auto vzint = sv0 * zint / rint; |
289 | 297 |
|
290 | | - // Local derivatives at interior point (centered stencils) |
291 | | - const LocalFirstDerivatives l_dvar_int{.da = c2o<0>(p, intp, var), |
292 | | - .db = c2o<1>(p, intp, var), |
293 | | - .dc = c2o<2>(p, intp, var)}; |
| 298 | + // Local derivatives at interior point |
| 299 | + const auto dgf_da = c2o<0>(p, intp, var); |
| 300 | + const auto dgf_db = c2o<1>(p, intp, var); |
| 301 | + const auto dgf_dc = c2o<2>(p, intp, var); |
294 | 302 |
|
295 | 303 | // Jacobians at interior point |
296 | | - const Jacobians jac_int{ |
297 | | - vJ_da_dx(intp), vJ_da_dy(intp), vJ_da_dz(intp), |
298 | | - vJ_db_dx(intp), vJ_db_dy(intp), vJ_db_dz(intp), |
299 | | - vJ_dc_dx(intp), vJ_dc_dy(intp), vJ_dc_dz(intp)}; |
| 304 | + const auto da_dx = vJ_da_dx(intp); |
| 305 | + const auto da_dy = vJ_da_dy(intp); |
| 306 | + const auto da_dz = vJ_da_dz(intp); |
| 307 | + const auto db_dx = vJ_db_dx(intp); |
| 308 | + const auto db_dy = vJ_db_dy(intp); |
| 309 | + const auto db_dz = vJ_db_dz(intp); |
| 310 | + const auto dc_dx = vJ_dc_dx(intp); |
| 311 | + const auto dc_dy = vJ_dc_dy(intp); |
| 312 | + const auto dc_dz = vJ_dc_dz(intp); |
300 | 313 |
|
301 | 314 | // Global derivatives at interior point |
302 | | - const auto g_dvar_int{project_first(l_dvar_int, jac_int)}; |
| 315 | + const auto dgf_dx = dgf_db * db_dx + dgf_dc * dc_dx + da_dx * dgf_da; |
| 316 | + const auto dgf_dy = dgf_dc * dc_dy + db_dy * dgf_db + da_dy * dgf_da; |
| 317 | + const auto dgf_dz = dc_dz * dgf_dc + db_dz * dgf_db + da_dz * dgf_da; |
303 | 318 |
|
304 | 319 | // Radiative part at interior point |
305 | | - const auto rad = -vxint * g_dvar_int.dx - vyint * g_dvar_int.dy - |
306 | | - vzint * g_dvar_int.dz - |
| 320 | + const auto rad = -vxint * dgf_dx - vyint * dgf_dy - vzint * dgf_dz - |
307 | 321 | sv0 * (var(intp) - var0) / rint; |
308 | 322 |
|
309 | 323 | // Extrapolate Coulomb component, rescale to account for radial |
|
0 commit comments