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 modeenum 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 parallelfor (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 parallelfor (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 phiconst 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 levelsMultiFab::Copy(ghext->mfab, ghext->mfab, 1, 2, 1, ghext->nghostzones);MultiFab::Copy(ghext->mfab, ghext->mfab, 0, 1, 1, ghext->nghostzones);// Step timecctkGH->cctk_time += cctkGH->cctk_delta_time;
#pragma omp parallelfor (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 parallelfor (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);}