5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC
SII5E2TZYKWQOLI6QKMQ2RTOT5OC32VNLIC6KWNWVNW6YAHCNP4QC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
QN2UTSQP4IMCMVXZNR24J22N24ASGXF6EWXQ2P3VWVQERUPFAFDQC
RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC
3KQ5ZQE45TNGATMW4R4NHXVBULGOXHF7JYC53BZL2LNDQ2V7C2CAC
UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC
VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC
2XYZZL42IEZHGDJA6NDKGSQKGJP24LOTLFJ6RNHOKWHHSUYIHGKQC
TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC
KG47IF4CPCUT3BHS34WDRHTH5HYMBTY4OSTB3X7APR2E5ZJ32IYQC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
UUGQGVC4WEKN64WAP7F5QPS2UHGQB5ZLMFRIYNWKMIEBDO3QRX4AC
IVHURSHY4636OGIF3PNDO5CWOVRLJ75M4LP65J6I2E6KAM4QKF4AC
33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC
5IAXY3XZJTRMMVT2OVIJ6OXQJI6OJPTPCHHA4IVLVMHANCCC5NKAC
TOBGHRPKEPSXDN56WGSZNWOMCBVJ4KUSLWYWI56MC2RR3MM3KLZAC
NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC
4UAVM7QNVW2KOPZCOWOCKMBIPDNDDIUOLNO5EFIWWHKKJ53ALSDAC
HNQLJR6ZWME6VBJ2Y7PSENLJPXC7INSS7NC2CKIWQAC776CQ74TQC
ostringstream buf;
buf << out_dir << "/" << groupname;
buf << ".it" << setw(6) << setfill('0') << cctk_iteration;
string filename = buf.str();
string filename = [&]() {
ostringstream buf;
buf << out_dir << "/" << groupname << ".it" << setw(6) << setfill('0')
<< cctk_iteration;
return buf.str();
}();
const bool is_root = CCTK_MyProc(nullptr) == 0;
if (is_root) {
string visitname = [&]() {
ostringstream buf;
buf << out_dir << "/" << groupname << ".visit";
return buf.str();
}();
ofstream visit(visitname, ios::app);
assert(visit.good());
// visit << filename << "/Header\n";
visit << groupname << ".it" << setw(6) << setfill('0') << cctk_iteration
<< "/Header\n";
visit.close();
}
#else
reduction<CCTK_REAL> red;
for (auto &restrict leveldata : ghext->leveldata) {
auto &restrict groupdata = leveldata.groupdata.at(gi);
MultiFab &mfab = *groupdata.mfab.at(tl);
reduction<CCTK_REAL> red1;
red1.min = mfab.min(vi);
red1.max = mfab.max(vi);
red1.sum = mfab.sum(vi);
// red1.sum2 = mfab.sum2(vi);
// red1.vol = mfab.vol(vi);
red1.maxabs = mfab.norminf(vi);
red1.sumabs = mfab.norm1(vi, mfab.fb_period);
red1.sum2abs = pow(mfab.norm2(vi), 2);
red += red1;
}
#endif
<< ": maxabs=" << red.maxabs << "\n";
// CCTK_REAL maxabs = 0.0;
// for (auto &restrict leveldata : ghext->leveldata) {
// auto &restrict groupdata = leveldata.groupdata.at(gi);
// MultiFab &mfab = *groupdata.mfab.at(tl);
// maxabs = fmax(maxabs, mfab.norminf(vi));
// }
// CCTK_VINFO(" %g", maxabs);
<< ": maxabs=" << red.maxabs << " vol=" << red.vol << "\n";
namespace {
template <typename T>
void mpi_reduce_typed(const void *restrict x0, void *restrict y0,
const int *restrict length) {
const T *restrict x = static_cast<const T *>(x0);
T *restrict y = static_cast<T *>(y0);
#pragma omp simd
for (int i = 0; i < *length; ++i)
y[i] += x[i];
}
void mpi_reduce(void *restrict x, void *restrict y, int *restrict length,
MPI_Datatype *restrict datatype) {
if (*datatype == MPI_FLOAT)
return mpi_reduce_typed<float>(x, y, length);
if (*datatype == MPI_DOUBLE)
return mpi_reduce_typed<double>(x, y, length);
if (*datatype == MPI_LONG_DOUBLE)
return mpi_reduce_typed<long double>(x, y, length);
abort();
}
} // namespace
MPI_Op reduction_mpi_op() {
static MPI_Op op = MPI_OP_NULL;
if (op == MPI_OP_NULL)
MPI_Op_create(mpi_reduce, 1, &op);
return op;
}
////////////////////////////////////////////////////////////////////////////////
#pragma omp parallel reduction(reduction : red)
{
for (auto &restrict leveldata : ghext->leveldata) {
const auto &restrict geom = ghext->amrcore->Geom(leveldata.level);
auto &restrict groupdata = leveldata.groupdata.at(gi);
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
for (auto &restrict leveldata : ghext->leveldata) {
const auto &restrict geom = ghext->amrcore->Geom(leveldata.level);
const auto &restrict groupdata = leveldata.groupdata.at(gi);
const MultiFab &mfab = *groupdata.mfab.at(tl);
unique_ptr<iMultiFab> finemask;
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
const CCTK_REAL *restrict ptr = grid.ptr(vars, vi);
const int fine_level = leveldata.level + 1;
if (fine_level < int(ghext->leveldata.size())) {
const auto &restrict fine_leveldata = ghext->leveldata.at(fine_level);
const auto &restrict fine_groupdata = fine_leveldata.groupdata.at(gi);
const MultiFab &fine_mfab = *fine_groupdata.mfab.at(tl);
const IntVect reffact{2, 2, 2};
finemask = make_unique<iMultiFab>(
makeFineMask(mfab, fine_mfab.boxArray(), reffact));
}
#warning "TODO: Skip refined regions"
red += reduce_array(geom, groupdata, ptr, grid);
}
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel reduction(reduction : red)
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
unique_ptr<Array4<const int> > finemask_array4;
if (finemask)
finemask_array4 = make_unique<Array4<const int> >(finemask->array(mfi));
const Array4<const CCTK_REAL> &vars = mfab.array(mfi);
const CCTK_REAL *restrict ptr = grid.ptr(vars, vi);
red += reduce_array(geom, groupdata, grid, ptr, finemask_array4.get());
#include "timer.hxx"
namespace AMReX {
Timer::Timer(const string &name)
: name(name), handle(CCTK_TimerCreate(name.c_str())) {}
Interval::Interval(const Timer &timer) : timer(timer) {
CCTK_TimerStartI(timer.handle);
}
Interval::~Interval() { CCTK_TimerStopI(timer.handle); }
} // namespace AMReX
#ifndef TIMER_HXX
#define TIMER_HXX
#include <cctk.h>
#undef copysign
#undef fpclassify
#undef isfinite
#undef isinf
#undef isnan
#undef isnormal
#undef signbit
#include <string>
namespace AMReX {
using namespace std;
class Interval;
class Timer {
friend class Interval;
string name;
int handle;
public:
Timer() = delete;
Timer(const string &name);
};
class Interval {
const Timer &timer;
public:
Interval() = delete;
Interval(const Timer &timer);
~Interval();
};
} // namespace AMReX
#endif // #ifndef TIMER_HXX
# *phi,t + d*A = 0 phi,t + A^i,i = 0
# *phi,t + d*A = 0 phi,t + A^i,i = 0
#
# E = -dphi - A,t E_i = - phi,i - A_i,t
# B = dA B_ij = A_j,i - A_i,j
# dE + B,t = 0 E_j,i - E_i,j + B_ij,t = 0
# dB = 0 B_jk,i + B_ki,j + B_ij,k = 0
# d*E = rho/epsilon E^i,i = rho/epsilon
# d*B - dt*E = mu J - B^ij,i - E^j,t = mu J^j
ActiveThorns = "
AMReX
Formaline
IOUtil
MaxwellToyAMReX
SystemTopology
# TerminationTrigger
TimerReport
"
$nlevels = 5
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 2.0
# AMReX::verbose = yes
AMReX::xmin = -1.0
AMReX::ymin = -1.0
AMReX::zmin = -1.0
AMReX::xmax = +1.0
AMReX::ymax = +1.0
AMReX::zmax = +1.0
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.01
AMReX::dtfac = 0.5
MaxwellToyAMReX::evolve_b = no
MaxwellToyAMReX::initial_condition = "Gaussian wave"
MaxwellToyAMReX::spatial_frequency_x = 0.5
MaxwellToyAMReX::spatial_frequency_y = 0.0
MaxwellToyAMReX::spatial_frequency_z = 0.0
MaxwellToyAMReX::analyse_every = $ncells * 2 ** ($nlevels - 1) / 32
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32
# AMReX::out_tsv = yes
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50
# TerminationTrigger::termination_from_file = yes
# TerminationTrigger::create_termination_file = yes
# TerminationTrigger::termination_file = "../TERMINATE"
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50
# TerminationTrigger::termination_from_file = yes
# TerminationTrigger::create_termination_file = yes
# TerminationTrigger::termination_file = "../TERMINATE"
CCTK_INT analyse_every "How often to calculate analysis quantities" STEERABLE=always
{
0 :: "never"
1:* :: "every that many iterations"
} 1
SCHEDULE MaxwellToyAMReX_Initialize AT initial
{
LANG: C
} "Set up initial conditions for the Maxwell equations"
if (evolve_b) {
SCHEDULE MaxwellToyAMReX_Initialize AT initial
{
LANG: C
SYNC: potential_phi
SYNC: potential_ax potential_ay potential_az
SYNC: field_ex field_ey field_ez
SYNC: field_bx field_by field_bz
SYNC: current_rho
SYNC: current_jx current_jy current_jz
} "Set up initial conditions for the Maxwell equations"
} else {
SCHEDULE MaxwellToyAMReX_Initialize AT initial
{
LANG: C
SYNC: potential_phi
SYNC: potential_ax potential_ay potential_az
SYNC: field_ex field_ey field_ez
SYNC: current_rho
SYNC: current_jx current_jy current_jz
} "Set up initial conditions for the Maxwell equations"
SCHEDULE MaxwellToyAMReX_Dependents1 AT initial AFTER MaxwellToyAMReX_Initialize
{
LANG: C
SYNC: field_bx field_by field_bz
} "Calculate the magnetic field"
}
if (evolve_b) {
SCHEDULE MaxwellToyAMReX_Evolve1 AT evol
{
LANG: C
SYNC: potential_ax potential_ay potential_az
SYNC: field_bx field_by field_bz
SYNC: current_jx current_jy current_jz
} "Evolve the Maxwell equations, step 1"
} else {
SCHEDULE MaxwellToyAMReX_Evolve1 AT evol
{
LANG: C
SYNC: potential_ax potential_ay potential_az
SYNC: current_jx current_jy current_jz
} "Evolve the Maxwell equations, step 1"
SCHEDULE MaxwellToyAMReX_Evolve1 AT evol
{
LANG: C
SYNC: potential_ax potential_ay potential_az
SYNC: field_bx field_by field_bz
SYNC: current_jx current_jy current_jz
} "Evolve the Maxwell equations, step 1"
SCHEDULE MaxwellToyAMReX_Dependents1 AT evol AFTER MaxwellToyAMReX_Evolve1
{
LANG: C
SYNC: field_bx field_by field_bz
} "Calculate the magnetic field"
}
# SCHEDULE MaxwellToyAMReX_Constraints AT analysis
# {
# LANG: C
# } "Calculate constraints of the Maxwell equations"
SCHEDULE MaxwellToyAMReX_Constraints AT analysis
{
LANG: C
SYNC: constraints_dive
SYNC: constraints_divb
# TRIGGERS: constraints_dive
# TRIGGERS: constraints_divb
} "Calculate constraints of the Maxwell equations"
# TRIGGERS: errors_potential_phi
# TRIGGERS: errors_potential_ax errors_potential_ay errors_potential_az
# TRIGGERS: errors_field_ex errors_field_ey errors_field_ez
# TRIGGERS: errors_field_bx errors_field_by errors_field_bz
# TRIGGERS: errors_current_rho
# TRIGGERS: errors_current_jx errors_current_jy errors_current_jz
if (CCTK_EQUALS(initial_condition, "plane wave")) {
Loop::loop_all<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
phi_(p.i, p.j, p.k) = linear(t + dt / 2, p.x, p.y, p.z).phi;
});
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ax_(p.i, p.j, p.k) = linear(t, p.x, p.y, p.z).ax;
});
Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ay_(p.i, p.j, p.k) = linear(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
az_(p.i, p.j, p.k) = linear(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ex_(p.i, p.j, p.k) =
-xderiv(linear<CCTK_REAL>, dx)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ax;
});
Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ey_(p.i, p.j, p.k) =
-yderiv(linear<CCTK_REAL>, dy)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ay;
});
Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
ez_(p.i, p.j, p.k) =
-zderiv(linear<CCTK_REAL>, dz)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).az;
});
if (evolve_b) {
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
yderiv(linear<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(linear<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
zderiv(linear<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(linear<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
xderiv(linear<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(linear<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax;
});
}
Loop::loop_all<0, 0, 0>(
cctkGH, [&](const Loop::PointDesc &p) { rho_(p.i, p.j, p.k) = 0.0; });
Loop::loop_all<1, 0, 0>(
cctkGH, [&](const Loop::PointDesc &p) { jx_(p.i, p.j, p.k) = 0.0; });
Loop::loop_all<0, 1, 0>(
cctkGH, [&](const Loop::PointDesc &p) { jy_(p.i, p.j, p.k) = 0.0; });
Loop::loop_all<0, 0, 1>(
cctkGH, [&](const Loop::PointDesc &p) { jz_(p.i, p.j, p.k) = 0.0; });
} else if (CCTK_EQUALS(initial_condition, "plane wave")) {
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
yderiv(plane_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(plane_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
zderiv(plane_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(plane_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
xderiv(plane_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(plane_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax;
});
if (evolve_b) {
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
yderiv(plane_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(plane_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
zderiv(plane_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(plane_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
xderiv(plane_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(plane_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax;
});
}
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
yderiv(gaussian_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(gaussian_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
zderiv(gaussian_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(gaussian_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
xderiv(gaussian_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(gaussian_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax;
});
if (evolve_b) {
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
yderiv(gaussian_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(gaussian_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
zderiv(gaussian_wave<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(gaussian_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
xderiv(gaussian_wave<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(gaussian_wave<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax;
});
}
}
extern "C" void MaxwellToyAMReX_Dependents1(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);
const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);
const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_(cctkGH, az);
const Loop::GF3D<CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);
const Loop::GF3D<CCTK_REAL, 1, 0, 1> by_(cctkGH, by);
const Loop::GF3D<CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);
Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) = (az_(p.i, p.j + 1, p.k) - az_(p.i, p.j, p.k)) / dy -
(ay_(p.i, p.j, p.k + 1) - ay_(p.i, p.j, p.k)) / dz;
});
Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) = (ax_(p.i, p.j, p.k + 1) - ax_(p.i, p.j, p.k)) / dz -
(az_(p.i + 1, p.j, p.k) - az_(p.i, p.j, p.k)) / dx;
});
Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) = (ay_(p.i + 1, p.j, p.k) - ay_(p.i, p.j, p.k)) / dx -
(ax_(p.i, p.j + 1, p.k) - ax_(p.i, p.j, p.k)) / dy;
});
Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
bx_p_(p.i, p.j, p.k) +
dt * ((ey_p_(p.i, p.j, p.k + 1) - ey_p_(p.i, p.j, p.k)) / dz -
(ez_p_(p.i, p.j + 1, p.k) - ez_p_(p.i, p.j, p.k)) / dy);
});
Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
by_p_(p.i, p.j, p.k) +
dt * ((ez_p_(p.i + 1, p.j, p.k) - ez_p_(p.i, p.j, p.k)) / dx -
(ex_p_(p.i, p.j, p.k + 1) - ex_p_(p.i, p.j, p.k)) / dz);
});
Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
bz_p_(p.i, p.j, p.k) +
dt * ((ex_p_(p.i, p.j + 1, p.k) - ex_p_(p.i, p.j, p.k)) / dy -
(ey_p_(p.i + 1, p.j, p.k) - ey_p_(p.i, p.j, p.k)) / dx);
});
if (evolve_b) {
Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) =
bx_p_(p.i, p.j, p.k) +
dt * ((ey_p_(p.i, p.j, p.k + 1) - ey_p_(p.i, p.j, p.k)) / dz -
(ez_p_(p.i, p.j + 1, p.k) - ez_p_(p.i, p.j, p.k)) / dy);
});
Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) =
by_p_(p.i, p.j, p.k) +
dt * ((ez_p_(p.i + 1, p.j, p.k) - ez_p_(p.i, p.j, p.k)) / dx -
(ex_p_(p.i, p.j, p.k + 1) - ex_p_(p.i, p.j, p.k)) / dz);
});
Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) =
bz_p_(p.i, p.j, p.k) +
dt * ((ex_p_(p.i, p.j + 1, p.k) - ex_p_(p.i, p.j, p.k)) / dy -
(ey_p_(p.i + 1, p.j, p.k) - ey_p_(p.i, p.j, p.k)) / dx);
});
}
});
#else
Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
jx_(p.i, p.j, p.k) = jx_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
jy_(p.i, p.j, p.k) = jy_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
jz_(p.i, p.j, p.k) = jz_p_(p.i, p.j, p.k);
if (evolve_b) {
Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bx_(p.i, p.j, p.k) = bx_p_(p.i, p.j, p.k);
});
Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
by_(p.i, p.j, p.k) = by_p_(p.i, p.j, p.k);
});
Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bz_(p.i, p.j, p.k) = bz_p_(p.i, p.j, p.k);
});
}
Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ax_(p.i, p.j, p.k) = ax_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ay_(p.i, p.j, p.k) = ay_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
az_(p.i, p.j, p.k) = az_p_(p.i, p.j, p.k);
});
#endif
(ay_(p.i, p.j, p.k) - az_(p.i, p.j, p.k - 1)) / dz);
(az_(p.i, p.j, p.k) - az_(p.i, p.j, p.k - 1)) / dz);
});
#else
Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
rho_(p.i, p.j, p.k) = rho_p_(p.i, p.j, p.k);
});
Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ex_(p.i, p.j, p.k) = ex_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ey_(p.i, p.j, p.k) = ey_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
ez_(p.i, p.j, p.k) = ez_p_(p.i, p.j, p.k);
});
Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
phi_(p.i, p.j, p.k) = phi_p_(p.i, p.j, p.k);
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);
const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);
const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> by_(cctkGH, by);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);
// Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
// regrid_error[p.idx] = fabs(p.x) < 0.5 && fabs(p.y) < 0.5 && fabs(p.z) <
// 0.5;
// });
// Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
// CCTK_REAL err = 0;
// for (int b = 0; b < 2; ++b)
// for (int a = 0; a < 2; ++a)
// err += fabs(ax_(p.i - 1, p.j + a, p.k + b) -
// 2 * ax_(p.i, p.j + a, p.k + b) +
// ax_(p.i + 1, p.j + a, p.k + b)) +
// fabs(ay_(p.i + a, p.j - 1, p.k + b) -
// 2 * ay_(p.i + a, p.j, p.k + b) +
// ay_(p.i + a, p.j + 1, p.k + b)) +
// fabs(az_(p.i + a, p.j + b, p.k - 1) -
// 2 * az_(p.i + a, p.j + b, p.k) +
// az_(p.i + a, p.j + b, p.k + 1));
// regrid_error[p.idx] = err;
// });
for (int b = 0; b < 2; ++b)
for (int a = 0; a < 2; ++a)
err += fabs(ax_(p.i - 1, p.j + a, p.k + b) -
2 * ax_(p.i, p.j + a, p.k + b) +
ax_(p.i + 1, p.j + a, p.k + b)) +
fabs(ay_(p.i + a, p.j - 1, p.k + b) -
2 * ay_(p.i + a, p.j, p.k + b) +
ay_(p.i + a, p.j + 1, p.k + b)) +
fabs(az_(p.i + a, p.j + b, p.k - 1) -
2 * az_(p.i + a, p.j + b, p.k) +
az_(p.i + a, p.j + b, p.k + 1));
const auto accum = [&](CCTK_REAL x) { err += fabs(x); };
const auto diffx = [&](const auto &var_, int i, int j, int k) {
for (int b = 0; b < 2; ++b)
for (int a = 0; a < 2; ++a)
accum(var_(i + 1, j + a, k + b) - var_(i, j + a, k + b));
};
const auto diffy = [&](const auto &var_, int i, int j, int k) {
for (int b = 0; b < 2; ++b)
for (int a = 0; a < 2; ++a)
accum(var_(i + a, j + 1, k + b) - var_(i + a, j, k + b));
};
const auto diffz = [&](const auto &var_, int i, int j, int k) {
for (int b = 0; b < 2; ++b)
for (int a = 0; a < 2; ++a)
accum(var_(i + a, j + b, k + 1) - var_(i + a, j + b, k));
};
const auto diff000 = [&](const auto &var_) {
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k);
};
const auto diff100 = [&](const auto &var_) {
diffx(var_, p.i - 1, p.j, p.k);
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k);
};
const auto diff010 = [&](const auto &var_) {
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j - 1, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k);
};
const auto diff001 = [&](const auto &var_) {
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k - 1);
diffz(var_, p.i, p.j, p.k);
};
const auto diff011 = [&](const auto &var_) {
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j - 1, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k - 1);
diffz(var_, p.i, p.j, p.k);
};
const auto diff101 = [&](const auto &var_) {
diffx(var_, p.i - 1, p.j, p.k);
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k - 1);
diffz(var_, p.i, p.j, p.k);
};
const auto diff110 = [&](const auto &var_) {
diffx(var_, p.i - 1, p.j, p.k);
diffx(var_, p.i, p.j, p.k);
diffy(var_, p.i, p.j - 1, p.k);
diffy(var_, p.i, p.j, p.k);
diffz(var_, p.i, p.j, p.k);
};
diff000(phi_);
diff100(ax_);
diff010(ay_);
diff001(az_);
diff100(ex_);
diff010(ey_);
diff001(ez_);
diff011(bx_);
diff101(by_);
diff110(bz_);
if (CCTK_EQUALS(initial_condition, "linear")) {
Loop::loop_all<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
phierr_(p.i, p.j, p.k) =
phi_(p.i, p.j, p.k) - linear(t + dt / 2, p.x, p.y, p.z).phi;
});
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
axerr_(p.i, p.j, p.k) = ax_(p.i, p.j, p.k) - linear(t, p.x, p.y, p.z).ax;
});
Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
ayerr_(p.i, p.j, p.k) = ay_(p.i, p.j, p.k) - linear(t, p.x, p.y, p.z).ay;
});
Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
azerr_(p.i, p.j, p.k) = az_(p.i, p.j, p.k) - linear(t, p.x, p.y, p.z).az;
});
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
exerr_(p.i, p.j, p.k) =
ex_(p.i, p.j, p.k) -
(-xderiv(linear<CCTK_REAL>, dx)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ax);
});
Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
eyerr_(p.i, p.j, p.k) =
ey_(p.i, p.j, p.k) -
(-yderiv(linear<CCTK_REAL>, dy)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ay);
});
Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
ezerr_(p.i, p.j, p.k) =
ez_(p.i, p.j, p.k) -
(-zderiv(linear<CCTK_REAL>, dz)(t + dt / 2, p.x, p.y, p.z).phi -
tderiv(linear<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).az);
});
if (CCTK_EQUALS(initial_condition, "plane wave")) {
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
bxerr_(p.i, p.j, p.k) =
bx_(p.i, p.j, p.k) -
(yderiv(linear<CCTK_REAL>, dy)(t, p.x, p.y, p.z).az -
zderiv(linear<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ay);
});
Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
byerr_(p.i, p.j, p.k) =
by_(p.i, p.j, p.k) -
(zderiv(linear<CCTK_REAL>, dz)(t, p.x, p.y, p.z).ax -
xderiv(linear<CCTK_REAL>, dx)(t, p.x, p.y, p.z).az);
});
Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
bzerr_(p.i, p.j, p.k) =
bz_(p.i, p.j, p.k) -
(xderiv(linear<CCTK_REAL>, dx)(t, p.x, p.y, p.z).ay -
yderiv(linear<CCTK_REAL>, dy)(t, p.x, p.y, p.z).ax);
});
Loop::loop_all<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
rhoerr_(p.i, p.j, p.k) = 0.0;
});
Loop::loop_all<1, 0, 0>(
cctkGH, [&](const Loop::PointDesc &p) { jxerr_(p.i, p.j, p.k) = 0.0; });
Loop::loop_all<0, 1, 0>(
cctkGH, [&](const Loop::PointDesc &p) { jyerr_(p.i, p.j, p.k) = 0.0; });
Loop::loop_all<0, 0, 1>(
cctkGH, [&](const Loop::PointDesc &p) { jzerr_(p.i, p.j, p.k) = 0.0; });
} else if (CCTK_EQUALS(initial_condition, "plane wave")) {
if (analyse_every <= 0 || cctk_iteration % analyse_every != 0)
return;
const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);
const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);
const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> by_(cctkGH, by);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);
const Loop::GF3D<const CCTK_REAL, 0, 0, 0> rho_(cctkGH, rho);
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> jx_(cctkGH, jx);
const Loop::GF3D<const CCTK_REAL, 0, 1, 0> jy_(cctkGH, jy);
const Loop::GF3D<const CCTK_REAL, 0, 0, 1> jz_(cctkGH, jz);
const Loop::GF3D<CCTK_REAL, 0, 0, 0> dive_(cctkGH, dive);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> divb_(cctkGH, divb);
Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
dive_(p.i, p.j, p.k) = (ex_(p.i, p.j, p.k) - ex_(p.i - 1, p.j, p.k)) / dx +
(ey_(p.i, p.j, p.k) - ey_(p.i, p.j - 1, p.k)) / dy +
(ez_(p.i, p.j, p.k) - ez_(p.i, p.j, p.k - 1)) / dz -
rho_(p.i, p.j, p.k);
});
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
divb_(p.i, p.j, p.k) = (bx_(p.i + 1, p.j, p.k) - bx_(p.i, p.j, p.k)) / dx +
(by_(p.i, p.j + 1, p.k) - by_(p.i, p.j, p.k)) / dy +
(bz_(p.i, p.j, p.k + 1) - bz_(p.i, p.j, p.k)) / dz;
});
}