Skip to content

Commit

Permalink
Z4co: fix bug in initial and enforce
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwJi
lwJi committed Jul 4, 2024
1 parent 58572ae commit 2054d80
Show file tree
Hide file tree
Showing 3 changed files with 20 additions and 31 deletions.
28 changes: 12 additions & 16 deletions Z4co/src/enforce.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,6 @@ extern "C" void Z4co_Enforce(CCTK_ARGUMENTS) {
typedef simdl<CCTK_REAL> vbool;
constexpr size_t vsize = tuple_size_v<vreal>;

const auto delta3 = one<smat<vreal, 3>>()();

#ifdef __CUDACC__
const nvtxRangeId_t range = nvtxRangeStartA("Z4co_Enforce::enforce");
#endif
Expand All @@ -63,21 +61,20 @@ extern "C" void Z4co_Enforce(CCTK_ARGUMENTS) {

// Enforce floors

const vreal chi = fmax(vreal(chi_floor - 1), chi_old);
const vreal alphaG = fmax(vreal(alphaG_floor - 1), alphaG_old);
const vreal chi = fmax(vreal(chi_floor), chi_old);
const vreal alphaG = fmax(vreal(alphaG_floor), alphaG_old);

// Enforce algebraic constraints
// See arXiv:1212.2901 [gr-qc].

const vreal detgammat_old = calc_det(delta3 + gammat_old);
const vreal chi1_old = 1 / cbrt(detgammat_old) - 1;
const vreal detgammat_old = calc_det(gammat_old);
const vreal chi1_old = 1 / cbrt(detgammat_old);
const smat<vreal, 3> gammat([&](int a, int b) ARITH_INLINE {
return (1 + chi1_old) * (delta3(a, b) + gammat_old(a, b)) -
delta3(a, b);
return chi1_old * gammat_old(a, b);
});
#ifdef CCTK_DEBUG
const vreal detgammat = calc_det(delta3 + gammat);
const vreal gammat_norm = maxabs(delta3 + gammat);
const vreal detgammat = calc_det(gammat);
const vreal gammat_norm = maxabs(gammat);
const vreal gammat_scale = gammat_norm;
#if !defined __CUDACC__ && !defined __HIPCC__
if (!(all(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale))) {
Expand All @@ -90,20 +87,19 @@ extern "C" void Z4co_Enforce(CCTK_ARGUMENTS) {
assert(all(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale));
#endif

const smat<vreal, 3> gammatu =
calc_inv(delta3 + gammat, vreal(1)) - delta3;
const smat<vreal, 3> gammatu = calc_inv(gammat, vreal(1));

const vreal traceAt_old = sum_symm<3>([&](int x, int y) ARITH_INLINE {
return (delta3(x, y) + gammatu(x, y)) * At_old(x, y);
return gammatu(x, y) * At_old(x, y);
});
const smat<vreal, 3> At([&](int a, int b) ARITH_INLINE {
return At_old(a, b) - traceAt_old / 3 * (delta3(a, b) + gammat(a, b));
return At_old(a, b) - traceAt_old / 3 * gammat(a, b);
});
#ifdef CCTK_DEBUG
const vreal traceAt = sum_symm<3>([&](int x, int y) ARITH_INLINE {
return (delta3(x, y) + gammatu(x, y)) * At(x, y);
return gammatu(x, y) * At(x, y);
});
const vreal gammatu_norm = maxabs(delta3 + gammatu);
const vreal gammatu_norm = maxabs(gammatu);
const vreal At_norm = maxabs(At);
const vreal At_scale = fmax(fmax(gammat_norm, gammatu_norm), At_norm);
#if !defined __CUDACC__ && !defined __HIPCC__
Expand Down
13 changes: 5 additions & 8 deletions Z4co/src/initial1.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,6 @@ extern "C" void Z4co_Initial1(CCTK_ARGUMENTS) {
typedef simdl<CCTK_REAL> vbool;
constexpr size_t vsize = tuple_size_v<vreal>;

const auto delta3 = one<smat<vreal, 3>>()();

const Loop::GridDescBaseDevice grid(cctkGH);
#ifdef __CUDACC__
const nvtxRangeId_t range = nvtxRangeStartA("Z4co_Initial1::initial1");
Expand All @@ -95,11 +93,10 @@ extern "C" void Z4co_Initial1(CCTK_ARGUMENTS) {
const vreal detg = calc_det(g);
const smat<vreal, 3> gu = calc_inv(g, detg);

const vreal chi = 1 / cbrt(detg) - 1;
const vreal chi = 1 / cbrt(detg);

const smat<vreal, 3> gammat([&](int a, int b) ARITH_INLINE {
return (1 + chi) * g(a, b) - delta3(a, b);
});
const smat<vreal, 3> gammat([&](int a, int b)
ARITH_INLINE { return chi * g(a, b); });

const vreal trK = sum_symm<3>(
[&](int x, int y) ARITH_INLINE { return gu(x, y) * K(x, y); });
Expand All @@ -109,10 +106,10 @@ extern "C" void Z4co_Initial1(CCTK_ARGUMENTS) {
const vreal Kh = trK - 2 * Theta;

const smat<vreal, 3> At([&](int a, int b) ARITH_INLINE {
return (1 + chi) * (K(a, b) - trK / 3 * g(a, b));
return chi * (K(a, b) - trK / 3 * g(a, b));
});

const vreal alphaG = alp - 1;
const vreal alphaG = alp;

const vec<vreal, 3> betaG([&](int a) ARITH_INLINE { return beta(a); });

Expand Down
10 changes: 3 additions & 7 deletions Z4co/src/initial2.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,6 @@ extern "C" void Z4co_Initial2(CCTK_ARGUMENTS) {
typedef simdl<CCTK_REAL> vbool;
constexpr size_t vsize = tuple_size_v<vreal>;

const auto delta3 = one<smat<vreal, 3>>()();

const Loop::GridDescBaseDevice grid(cctkGH);
#ifdef __CUDACC__
const nvtxRangeId_t range = nvtxRangeStartA("Z4co_Initial2::initial2");
Expand All @@ -59,19 +57,17 @@ extern "C" void Z4co_Initial2(CCTK_ARGUMENTS) {
const smat<vreal, 3> gammat = gf_gammat1(mask, index1);

// Calculate Z4c variables (only Gamt)
const smat<vreal, 3> gammatu =
calc_inv(delta3 + gammat, vreal(1)) - delta3;
const smat<vreal, 3> gammatu = calc_inv(gammat, vreal(1));

const smat<vec<vreal, 3>, 3> dgammat([&](int a, int b) {
return deriv(mask, gf_gammat1(a, b), p.I, dx);
});

const vec<smat<vreal, 3>, 3> Gammatl = calc_gammal(dgammat);
const vec<smat<vreal, 3>, 3> Gammat =
calc_gamma(delta3 + gammatu, Gammatl);
const vec<smat<vreal, 3>, 3> Gammat = calc_gamma(gammatu, Gammatl);
const vec<vreal, 3> Gamt([&](int a) ARITH_INLINE {
return sum_symm<3>([&](int x, int y) ARITH_INLINE {
return (delta3(x, y) + gammatu(x, y)) * Gammat(a)(x, y);
return gammatu(x, y) * Gammat(a)(x, y);
});
});

Expand Down

0 comments on commit 2054d80

Please sign in to comment.