ORAAQXS3UNKXYDJBCTHNEIDJIZDHWMM5EVCZKVMFRKQK2NIQNJGAC
}
namespace {
enum class mode_t { unknown, local, level, global, meta };
mode_t decode_mode(const cFunctionData *attribute) {
bool local_mode = attribute->local;
bool level_mode = attribute->level;
bool global_mode = attribute->global;
bool meta_mode = attribute->meta;
assert(int(local_mode) + int(level_mode) + int(global_mode) +
int(meta_mode) <=
1);
if (attribute->local)
return mode_t::local;
if (attribute->level)
return mode_t::level;
if (attribute->global)
return mode_t::global;
if (attribute->meta)
return mode_t::meta;
return mode_t::local; // default
// Decode mode
enum class mode_t { unknown, local, level, global, meta };
const auto decode_mode = [&]() {
bool local_mode = attribute->local;
bool level_mode = attribute->level;
bool global_mode = attribute->global;
bool meta_mode = attribute->meta;
assert(int(local_mode) + int(level_mode) + int(global_mode) +
int(meta_mode) <=
1);
if (attribute->local)
return mode_t::local;
if (attribute->level)
return mode_t::level;
if (attribute->global)
return mode_t::global;
if (attribute->meta)
return mode_t::meta;
return mode_t::local; // default
};
const mode_t mode = decode_mode();
CCTK_CallFunction(function, attribute, data);
mfis.resize(omp_get_max_threads(), nullptr);
#pragma omp parallel
for (MFIter mfi(ghext->mfab, MFItInfo().SetDynamic(true).EnableTiling(
{1024000, 16, 32}));
mfi.isValid(); ++mfi) {
mfis.at(omp_get_thread_num()) = &mfi;
CCTK_CallFunction(function, attribute, data);
}
mfis.clear();
# CCTK_REAL phi TYPE=gf TIMELEVELS=3
# {
# phi
# } "Scalar potential for wave equation"
# CCTK_REAL phi TYPE=gf
# {
# phi
# phi_p
# phi_p_p
# } "Scalar potential for wave equation"
const CCTK_REAL *restrict const x0 = ghext->geom.ProbLo();
const CCTK_REAL *restrict const x1 = ghext->geom.ProbHi();
for (int d = 0; d < 3; ++d) {
CCTK_REAL dx = (x1[d] - x0[d]) / ghext->ncells;
cctkGH->cctk_origin_space[d] = x0[d] + 0.5 * dx;
cctkGH->cctk_delta_space[d] = dx;
}
cctkGH->cctk_time = 0.0;
cctkGH->cctk_delta_time = 0.5 / ghext->ncells;
}
// Initialize phi
#pragma omp parallel
for (MFIter mfi(ghext->mfab,
MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32}));
mfi.isValid(); ++mfi) {
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.growntilebox();
// Initialize phi
const MFIter &mfi = *mfis.at(omp_get_thread_num());
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.growntilebox();
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
phi[idx] = standing(t0, x, y, z);
phi_p[idx] = standing(t0 - dt, x, y, z);
}
}
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
phi[idx] = standing(t0, x, y, z);
phi_p[idx] = standing(t0 - dt, x, y, z);
}
}
extern "C" void WaveToyAMReX_Cycle(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
// Cycle time levels
MultiFab::Copy(ghext->mfab, ghext->mfab, 1, 2, 1, ghext->nghostzones);
MultiFab::Copy(ghext->mfab, ghext->mfab, 0, 1, 1, ghext->nghostzones);
// Step time
cctkGH->cctk_time += cctkGH->cctk_delta_time;
#pragma omp parallel
for (MFIter mfi(ghext->mfab,
MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32}));
mfi.isValid(); ++mfi) {
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.tilebox();
const MFIter &mfi = *mfis.at(omp_get_thread_num());
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.tilebox();
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
const CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
const CCTK_REAL *restrict const phi_p_p = vars.ptr(0, 0, 0, 2);
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
const CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
const CCTK_REAL *restrict const phi_p_p = vars.ptr(0, 0, 0, 2);
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL ddx_phi =
(phi_p[idx - di] - 2 * phi_p[idx] + phi_p[idx + di]) / sqr(dx[0]);
CCTK_REAL ddy_phi =
(phi_p[idx - dj] - 2 * phi_p[idx] + phi_p[idx + dj]) / sqr(dx[1]);
CCTK_REAL ddz_phi =
(phi_p[idx - dk] - 2 * phi_p[idx] + phi_p[idx + dk]) / sqr(dx[2]);
phi[idx] = -phi_p_p[idx] + 2 * phi_p[idx] +
sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
}
}
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL ddx_phi =
(phi_p[idx - di] - 2 * phi_p[idx] + phi_p[idx + di]) / sqr(dx[0]);
CCTK_REAL ddy_phi =
(phi_p[idx - dj] - 2 * phi_p[idx] + phi_p[idx + dj]) / sqr(dx[1]);
CCTK_REAL ddz_phi =
(phi_p[idx - dk] - 2 * phi_p[idx] + phi_p[idx + dk]) / sqr(dx[2]);
phi[idx] = -phi_p_p[idx] + 2 * phi_p[idx] +
sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
}
}
extern "C" void WaveToyAMReX_Sync(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
// Initialize phi
#pragma omp parallel
for (MFIter mfi(ghext->mfab,
MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32}));
mfi.isValid(); ++mfi) {
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.growntilebox();
const MFIter &mfi = *mfis.at(omp_get_thread_num());
const Box &fbx = mfi.fabbox();
const Box &bx = mfi.growntilebox();
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const err = vars.ptr(0, 0, 0, 3);
const CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
CCTK_REAL *restrict const err = vars.ptr(0, 0, 0, 3);
const CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
err[idx] = phi[idx] - standing(t0, x, y, z);
}
}
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
err[idx] = phi[idx] - standing(t0, x, y, z);
}