33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC
7Y37CMWGUQY3IOMFZW46UINY2W5DIEE3TXGN7RDKYZO6FHCLHDBQC
MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC
RTNZAS3UPI6GG3KY4Z5WVXJ4R2YF5427BB6WAV3GHRS5W7XPOSUQC
2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC
QN2UTSQP4IMCMVXZNR24J22N24ASGXF6EWXQ2P3VWVQERUPFAFDQC
KG47IF4CPCUT3BHS34WDRHTH5HYMBTY4OSTB3X7APR2E5ZJ32IYQC
UUGQGVC4WEKN64WAP7F5QPS2UHGQB5ZLMFRIYNWKMIEBDO3QRX4AC
2DD222JSYRPHTXKSRXLSOMSCQPZUORNFLLO2P3GMIDELAAMD5MEQC
UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC
TOBGHRPKEPSXDN56WGSZNWOMCBVJ4KUSLWYWI56MC2RR3MM3KLZAC
PQB3EKQ6MBCXPTW4HB7SGMSTOTYMB3EFZX2573OFCQI6PSOEKCSQC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
2XYZZL42IEZHGDJA6NDKGSQKGJP24LOTLFJ6RNHOKWHHSUYIHGKQC
ORAAQXS3UNKXYDJBCTHNEIDJIZDHWMM5EVCZKVMFRKQK2NIQNJGAC
AEVGZIZEUIC52MCK3J4V547YEV2R4YQL3JUJW7FSP4R34PSZ43DAC
ActiveThorns = "
AMReX
IOUtil
WaveToyAMReX
"
$nlevels = 4
$ncells = 16
Cactus::cctk_show_schedule = no
Cactus::cctk_itlast = $ncells * 2 ** ($nlevels - 1) * 2
AMReX::verbose = yes
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells
AMReX::max_num_levels = $nlevels
AMReX::regrid_every = 16
AMReX::regrid_error_threshold = 0.02 # 0.05 # 0.1
AMReX::dtfac = 0.5
WaveToyAMReX::initial_condition = "periodic Gaussian"
WaveToyAMReX::spatial_frequency_x = 0.5
WaveToyAMReX::spatial_frequency_y = 0.0
WaveToyAMReX::spatial_frequency_z = 0.0
IO::out_dir = "planewave"
IO::out_every = 1 #TODO $ncells * 2 ** ($nlevels - 1) / 2
// // refine everywhere
// tags.setVal(boxArray(level), TagBox::SET);
// // refine centre
// const Box &dom = Geom(level).Domain();
// Box nbx;
// for (int d = 0; d < dim; ++d) {
// int md = (dom.bigEnd(d) + dom.smallEnd(d) + 1) / 2;
// int rd = (dom.bigEnd(d) - dom.smallEnd(d) + 1) / 2;
// // mark one fewer cells; AMReX seems to add one cell
// nbx.setSmall(d, md - rd / (1 << (level + 1)) + 1);
// nbx.setBig(d, md + rd / (1 << (level + 1)) - 2);
// }
// cout << "EE nbx: " << nbx << "\n";
// const BoxArray &ba = boxArray(level);
// tags.setVal(intersect(ba, nbx), TagBox::SET);
const int gi = CCTK_GroupIndex("AMReX::regrid_tag");
const int gi = CCTK_GroupIndex("AMReX::regrid_error");
void SetupLevel(int level, const BoxArray &ba, const DistributionMapping &dm) {
DECLARE_CCTK_PARAMETERS;
if (verbose)
CCTK_VINFO("SetupLevel level %d", level);
assert(level == int(ghext->leveldata.size()));
ghext->leveldata.resize(level + 1);
GHExt::LevelData &leveldata = ghext->leveldata.at(level);
leveldata.level = level;
void CactusAmrCore::MakeNewLevelFromScratch(int lev, Real time,
const int numgroups = CCTK_NumGroups();
leveldata.groupdata.resize(numgroups);
for (int gi = 0; gi < numgroups; ++gi) {
cGroup group;
int ierr = CCTK_GroupData(gi, &group);
assert(!ierr);
assert(group.grouptype == CCTK_GF);
assert(group.vartype == CCTK_VARIABLE_REAL);
assert(group.disttype == CCTK_DISTRIB_DEFAULT);
assert(group.dim == dim);
GHExt::LevelData::GroupData &groupdata = leveldata.groupdata.at(gi);
groupdata.firstvarindex = CCTK_FirstVarIndexI(gi);
groupdata.numvars = group.numvars;
// Allocate grid hierarchies
groupdata.mfab.resize(group.numtimelevels);
for (int tl = 0; tl < int(groupdata.mfab.size()); ++tl) {
groupdata.mfab.at(tl) =
make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);
}
}
}
void CactusAmrCore::MakeNewLevelFromScratch(int level, Real time,
CCTK_VINFO("MakeNewLevelFromCoarse level %d", lev);
CCTK_VINFO("MakeNewLevelFromCoarse level %d", level);
assert(level > 0);
SetupLevel(level, ba, dm);
// Prolongate
auto &leveldata = ghext->leveldata.at(level);
auto &coarseleveldata = ghext->leveldata.at(level - 1);
const int num_groups = CCTK_NumGroups();
for (int gi = 0; gi < num_groups; ++gi) {
auto &restrict groupdata = leveldata.groupdata.at(gi);
auto &restrict coarsegroupdata = coarseleveldata.groupdata.at(gi);
assert(coarsegroupdata.numvars == groupdata.numvars);
PhysBCFunctNoOp cphysbc;
PhysBCFunctNoOp fphysbc;
const IntVect reffact{2, 2, 2};
// periodic boundaries
const BCRec bcrec(BCType::int_dir, BCType::int_dir, BCType::int_dir,
BCType::int_dir, BCType::int_dir, BCType::int_dir);
const Vector<BCRec> bcs(groupdata.numvars, bcrec);
// If there is more than one time level, then we don't
// prolongate the oldest.
int ntls = groupdata.mfab.size();
int prolongate_tl = ntls > 1 ? ntls - 1 : ntls;
for (int tl = 0; tl < prolongate_tl; ++tl)
InterpFromCoarseLevel(*groupdata.mfab.at(tl), 0.0,
*coarsegroupdata.mfab.at(tl), 0, 0,
groupdata.numvars, ghext->amrcore->Geom(level - 1),
ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0,
reffact, &cell_bilinear_interp, bcs, 0);
}
CCTK_VINFO("RemakeLevel level %d", lev);
CCTK_VINFO("RemakeLevel level %d", level);
// Copy or prolongate
auto &leveldata = ghext->leveldata.at(level);
const int num_groups = CCTK_NumGroups();
for (int gi = 0; gi < num_groups; ++gi) {
auto &restrict groupdata = leveldata.groupdata.at(gi);
// If there is more than one time level, then we don't
// prolongate the oldest.
int ntls = groupdata.mfab.size();
int prolongate_tl = ntls > 1 ? ntls - 1 : ntls;
if (leveldata.level == 0) {
// Coarsest level: Copy from same level
// This should not happen
assert(0);
PhysBCFunctNoOp fphysbc;
for (int tl = 0; tl < ntls; ++tl) {
auto mfab =
make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);
if (tl < prolongate_tl)
FillPatchSingleLevel(*mfab, 0.0, {&*groupdata.mfab.at(tl)}, {0.0}, 0,
0, groupdata.numvars,
ghext->amrcore->Geom(level), fphysbc, 0);
groupdata.mfab.at(tl) = move(mfab);
}
} else {
// Refined level: copy from same level and/or prolongate from
// next coarser level
auto &coarseleveldata = ghext->leveldata.at(level - 1);
auto &restrict coarsegroupdata = coarseleveldata.groupdata.at(gi);
assert(coarsegroupdata.numvars == groupdata.numvars);
PhysBCFunctNoOp cphysbc;
PhysBCFunctNoOp fphysbc;
const IntVect reffact{2, 2, 2};
// periodic boundaries
const BCRec bcrec(BCType::int_dir, BCType::int_dir, BCType::int_dir,
BCType::int_dir, BCType::int_dir, BCType::int_dir);
const Vector<BCRec> bcs(groupdata.numvars, bcrec);
for (int tl = 0; tl < ntls; ++tl) {
auto mfab =
make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);
if (tl < prolongate_tl)
FillPatchTwoLevels(*mfab, 0.0, {&*coarsegroupdata.mfab.at(tl)}, {0.0},
{&*groupdata.mfab.at(tl)}, {0.0}, 0, 0,
groupdata.numvars, ghext->amrcore->Geom(level - 1),
ghext->amrcore->Geom(level), cphysbc, 0, fphysbc,
0, reffact, &cell_bilinear_interp, bcs, 0);
groupdata.mfab.at(tl) = move(mfab);
}
}
}
const int numgroups = CCTK_NumGroups();
leveldata.groupdata.resize(numgroups);
for (int gi = 0; gi < numgroups; ++gi) {
cGroup group;
int ierr = CCTK_GroupData(gi, &group);
assert(!ierr);
assert(group.grouptype == CCTK_GF);
assert(group.vartype == CCTK_VARIABLE_REAL);
assert(group.disttype == CCTK_DISTRIB_DEFAULT);
assert(group.dim == dim);
GHExt::LevelData::GroupData &groupdata = leveldata.groupdata.at(gi);
groupdata.firstvarindex = CCTK_FirstVarIndexI(gi);
groupdata.numvars = group.numvars;
// Allocate grid hierarchies
groupdata.mfab.resize(group.numtimelevels);
for (auto &groupdata : leveldata.groupdata) {
groupdata.mfab.at(tl) = make_unique<MultiFab>(
ghext->amrcore->boxArray(leveldata.level),
ghext->amrcore->DistributionMap(leveldata.level), groupdata.numvars,
ghost_size);
#warning "TODO"
assert(0);
groupdata.mfab.at(tl) =
make_unique<MultiFab>(ghext->amrcore->boxArray(level),
ghext->amrcore->DistributionMap(level),
groupdata.numvars, ghost_size);
current_level = ghext->amrcore->finestLevel();
CCTK_VINFO("Initializing level %d...", current_level);
// TODO: Remove "current_level" mechanism
// current_level = ghext->amrcore->finestLevel();
const int level = ghext->amrcore->finestLevel();
CCTK_VINFO("Initializing level %d...", level);
if (regrid_every > 0 && cctkGH->cctk_iteration % regrid_every == 0) {
CCTK_VINFO("Regridding...");
CCTK_REAL time = 0.0; // dummy time
const int old_numlevels = ghext->amrcore->finestLevel() + 1;
ghext->amrcore->regrid(0, time);
const int new_numlevels = ghext->amrcore->finestLevel() + 1;
CCTK_VINFO(" old levels %d, new levels %d", old_numlevels,
new_numlevels);
CCTK_Traverse(cctkGH, "CCTK_BASEGRID");
}
if (current_level < 0) {
// Loop over all levels
// TODO: parallelize this loop
for (auto &restrict leveldata : ghext->leveldata)
callfunc(leveldata);
} else {
// Loop over a single level
auto &restrict leveldata = ghext->leveldata.at(current_level);
// Loop over all levels
// TODO: parallelize this loop
for (auto &restrict leveldata : ghext->leveldata)
KEYWORD initial_condition "Initial condition"
{
"standing wave" :: ""
"periodic Gaussian" :: ""
} "standing wave"
CCTK_REAL spatial_frequency_x "spatial frequency"
{
*:* :: ""
} 0.5
CCTK_REAL spatial_frequency_y "spatial frequency"
{
*:* :: ""
} 0.5
CCTK_REAL spatial_frequency_z "spatial frequency"
{
*:* :: ""
} 0.5
CCTK_REAL width "width"
{
(0.0:* :: ""
} 0.1
CCTK_REAL phi_abs "Typical value of phi to determine absolute error" STEERABLE=always
{
(0.0:* :: ""
} 1.0
// Standing wave
template <typename T> T gaussian(T t, T x, T y, T z) {
DECLARE_CCTK_PARAMETERS;
T kx = M_PI * spatial_frequency_x;
T ky = M_PI * spatial_frequency_y;
T kz = M_PI * spatial_frequency_z;
T omega = sqrt(sqr(kx) + sqr(ky) + sqr(kz));
return exp(-0.5 * sqr(sin(kx * x + ky * y + kz * z - omega * t) / width));
}
for (int i = imin[0]; i < imax[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phi[idx] = standing(t, x, y, z);
phi_p[idx] = standing(t - dt, x, y, z);
// psi[idx] = timederiv(standing, dt)(t, x, y, z);
for (int i = imin[0]; i < imax[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phi[idx] = standing(t, x, y, z);
// phi_p[idx] = standing(t - dt, x, y, z);
psi[idx] = timederiv(standing, dt)(t, x, y, z);
}
}
}
} else if (CCTK_EQUALS(initial_condition, "periodic Gaussian")) {
for (int k = imin[2]; k < imax[2]; ++k) {
for (int j = imin[1]; j < imax[1]; ++j) {
#pragma omp simd
for (int i = imin[0]; i < imax[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phi[idx] = gaussian(t, x, y, z);
// phi_p[idx] = gaussian(t - dt, x, y, z);
psi[idx] = timederiv(gaussian, dt)(t, x, y, z);
}
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
CCTK_REAL r = sqrt(sqr(x) + sqr(y) + sqr(z));
// CCTK_REAL r = fabs(x);
tag[idx] = r <= 0.5 / cctk_levfac[0] - 0.5 * dx[0];
CCTK_REAL base_phi = fabs(phi[idx]) + fabs(phi_abs);
CCTK_REAL errx_phi =
fabs(phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / base_phi;
CCTK_REAL erry_phi =
fabs(phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / base_phi;
CCTK_REAL errz_phi =
fabs(phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / base_phi;
regrid_error[idx] = errx_phi + erry_phi + errz_phi;
// CCTK_REAL ddx_psi =
// (psi_p[idx - di] - 2 * psi_p[idx] + psi_p[idx + di]) /
// sqr(dx[0]);
// CCTK_REAL ddy_psi =
// (psi_p[idx - dj] - 2 * psi_p[idx] + psi_p[idx + dj]) /
// sqr(dx[1]);
// CCTK_REAL ddz_psi =
// (psi_p[idx - dk] - 2 * psi_p[idx] + psi_p[idx + dk]) /
// sqr(dx[2]);
// phi[idx] = phi_p[idx] + dt * psi_p[idx];
// psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi) +
// sqr(dt) * (ddx_psi + ddy_psi + ddz_psi);
// phi[idx] = phi_p[idx] + dt * psi_p[idx];
// psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
phi[idx] = 2 * phi_p[idx] - phi_p_p[idx] +
sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
// phi[idx] = 2 * phi_p[idx] - phi_p_p[idx] +
// sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
phi[idx] = phi_p[idx] + dt * psi[idx];
CCTK_REAL dt_phi = (phi[idx] - phi_p[idx]) / dt;
CCTK_REAL dx_phi = (phi[idx + di] - phi[idx - di]) / (2 * dx[0]);
CCTK_REAL dy_phi = (phi[idx + dj] - phi[idx - dj]) / (2 * dx[1]);
CCTK_REAL dz_phi = (phi[idx + dk] - phi[idx - dk]) / (2 * dx[2]);
eps[idx] = (sqr(dt_phi) + sqr(dx_phi) + sqr(dy_phi) + sqr(dz_phi)) / 2;
CCTK_REAL ddx_phi =
(phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / sqr(dx[0]);
CCTK_REAL ddy_phi =
(phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / sqr(dx[1]);
CCTK_REAL ddz_phi =
(phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / sqr(dx[2]);
CCTK_REAL psi_n = psi[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
CCTK_REAL dt_phi_p = psi_n;
CCTK_REAL dt_phi_m = psi_p[idx];
CCTK_REAL dx_phi_p = (phi[idx + di] - phi[idx]) / dx[0];
CCTK_REAL dx_phi_m = (phi[idx] - phi[idx - di]) / dx[0];
CCTK_REAL dy_phi_p = (phi[idx + dj] - phi[idx]) / dx[1];
CCTK_REAL dy_phi_m = (phi[idx] - phi[idx - dj]) / dx[1];
CCTK_REAL dz_phi_p = (phi[idx + dk] - phi[idx]) / dx[2];
CCTK_REAL dz_phi_m = (phi[idx] - phi[idx - dk]) / dx[2];
eps[idx] =
(sqr(dt_phi_p) + sqr(dt_phi_m) + sqr(dx_phi_p) + sqr(dx_phi_m) +
sqr(dy_phi_p) + sqr(dy_phi_m) + sqr(dz_phi_p) + sqr(dz_phi_m)) /
4;
for (int i = 0; i < cctk_lsh[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phierr[idx] = phi[idx] - standing(t, x, y, z);
// psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);
for (int i = 0; i < cctk_lsh[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phierr[idx] = phi[idx] - standing(t, x, y, z);
psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);
}
} else if (CCTK_EQUALS(initial_condition, "periodic Gaussian")) {
for (int k = 0; k < cctk_lsh[2]; ++k) {
for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
for (int i = 0; i < cctk_lsh[0]; ++i) {
const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
phierr[idx] = phi[idx] - gaussian(t, x, y, z);
psierr[idx] = psi[idx] - timederiv(gaussian, dt)(t, x, y, z);
}
}
}