JYHNXXG3OJTYBJ2U3WMKUS6N5WHXO773WXVRNSJLBC5FGEYLJP2AC KCIWCVZOHG44WBOLKI2XK33WPHPRI5FWCETF4AOGTPZISKCW3CLQC NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC 5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC NE2O3IMQBR3OAPLJC22SPJ7AX2ABULA7XNSIZYFHCQLWC45T6GJAC UD2O4RHXGNMJII7BYM75WDELMJFIAEJMMZM72VJMY34FGAB2DT4AC PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC Solve the Maxwell equations.2. Derivatives, and constraint preservationWe store variables and calculate derivatives in a way that is verysimilar to discrete differential forms. Scalars live on vertices,vectors on edges, fluxes on faces, and densities in cells. Thishappens in 4D, leading to a staggered evolution scheme. We have thefollowing variables:- electric potential phi: vertex, staggered in time (equivalent to atimelike edge)- magnetic vector potential A^i: (spacelike) edges- electric field E^i: (spacelike) edges, staggered in time
Solve the Maxwell equations
- magnetic field B^i: faces- electric charge density rho: vertex (equivalent to dual cell),staggered in time- electric current density j^i: edges (equivalent to dual face)- electric constraint div E: vertices (equivalent to dual cell),staggered in time- magnetic constraint div B: cellWhen calculating derivatives, we divide by the grid spacing (unlikefrom differential forms). This is for convenience only.In terms of differential forms, our variables correspond to thefollowing. The notation [phi, A] means that phi and A^i are the timeand space components of a 4-vector.- A: [phi, A]: 1-form- F: [E^i, B^i]: 2-form- j: *[rho, j^i]: 3-form living on dual grid- *[div E]: 4-form living on dual grid- *[div B]: 4-form4d Equations:- F = dA- dF = 0- d*F = 03+1 Equations:- E = dphi- B = dA- *d*E = *rho- *d*B = *j(TODO: Correct and complete equations)3. How to calculate consistent derivatives with mesh refinementProlongation opterators need to be conservative in cell centreddirections, and interpolating in vertex centred directions. We wantlinear accuracy for vertex centred and constant accuracy for cellcentred directions.
CCTK_REAL constraints_dive TYPE=gf TAGS='index={0 0 0} checkpoint="no"'{dive} "Constraints"CCTK_REAL constraints_divb TYPE=gf TAGS='index={1 1 1} checkpoint="no"'{divb} "Constraints"
WRITES: potential_phi(everywhere)WRITES: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)WRITES: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)WRITES: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)WRITES: current_rho(everywhere)WRITES: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
WRITES: potential_phi(interior)WRITES: potential_ax(interior) potential_ay(interior) potential_az(interior)WRITES: field_ex(interior) field_ey(interior) field_ez(interior)WRITES: field_bx(interior) field_by(interior) field_bz(interior)WRITES: current_rho(interior)WRITES: current_jx(interior) current_jy(interior) current_jz(interior)
WRITES: potential_phi(everywhere)WRITES: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)WRITES: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)WRITES: current_rho(everywhere)WRITES: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
WRITES: potential_phi(interior)WRITES: potential_ax(interior) potential_ay(interior) potential_az(interior)WRITES: field_ex(interior) field_ey(interior) field_ez(interior)WRITES: current_rho(interior)WRITES: current_jx(interior) current_jy(interior) current_jz(interior)
SCHEDULE MaxwellToyAMReX_EstimateError AT poststep{LANG: CREADS: potential_phi(everywhere)READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)WRITES: AMReX::regrid_error(interior)} "Estimate local error for regridding during evolution"
READS: potential_phi(everywhere)READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)WRITES: AMReX::regrid_error(interior)} "Estimate local error for regridding during evolution"
READS: potential_phi(everywhere)READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)READS: current_rho(everywhere)READS: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)} "Check for NaNs"
# SCHEDULE MaxwellToyAMReX_NaNCheck AT analysis AS zzz_MaxwellToyAMReX_NaNCheck# {# LANG: C## READS: potential_phi(everywhere)# READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)# READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)# READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)# READS: current_rho(everywhere)# READS: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)# READS: dive(everywhere)# READS: divb(everywhere)# } "Check for NaNs"
extern "C" void MaxwellToyAMReX_EstimateError(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_MaxwellToyAMReX_EstimateError;DECLARE_CCTK_PARAMETERS;const Loop::GF3D<const CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);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<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);#if 0Loop::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;});#endif#if 0Loop::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;});#endif#if 0Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {auto closeto = [&](CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL r) {return fabs(p.x - x) <= r && fabs(p.y - y) <= r && fabs(p.z - z) <= r;};constexpr CCTK_REAL x = 0.25;constexpr CCTK_REAL y = 0.25;constexpr CCTK_REAL z = 0.25;constexpr CCTK_REAL r = 0.2;regrid_error[p.idx] = // closeto(-x, -y, -z, r) ||closeto(+x, -y, -z, r) || closeto(-x, +y, -z, r) ||closeto(+x, +y, -z, r) || closeto(-x, -y, +z, r) ||closeto(+x, -y, +z, r) || closeto(-x, +y, +z, r) ||closeto(+x, +y, +z, r);});#endif
#if 1Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {CCTK_REAL err = 0;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_);regrid_error[p.idx] = err;});#endif}extern "C" void MaxwellToyAMReX_Boundaries(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_MaxwellToyAMReX_Boundaries;DECLARE_CCTK_PARAMETERS;// Do nothing; boundary conditions consist of synchronization only}
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 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> 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<const CCTK_REAL, 0, 0, 0> dive_(cctkGH, dive);const Loop::GF3D<const CCTK_REAL, 1, 1, 1> divb_(cctkGH, divb);const auto nancheck = [&](const CCTK_REAL *restrict var,const Loop::PointDesc &p) {if (CCTK_isnan(var[p.idx]))CCTK_VERROR("Found nan at [%d,%d,%d]", p.i, p.j, p.k);};Loop::loop_all<0, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(phi, p); });Loop::loop_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ax, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ay, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(az, p); });Loop::loop_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ex, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ey, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ez, p); });Loop::loop_all<0, 1, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(bx, p); });Loop::loop_all<1, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(by, p); });Loop::loop_all<1, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(bz, p); });
Loop::loop_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jx, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jy, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jz, p); });
const Loop::GF3D<CCTK_REAL, 1, 1, 1> divb_(cctkGH, divb);
// Loop::loop_all<0, 0, 0>(cctkGH,// [&](const Loop::PointDesc &p) { nancheck(dive, p);// });// Loop::loop_all<1, 1, 1>(cctkGH,// [&](const Loop::PointDesc &p) { nancheck(divb, p);// });}extern "C" void MaxwellToyAMReX_EstimateError(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;const Loop::GF3D<const CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);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<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);#if 0Loop::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;});#endif#if 0Loop::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;
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);
auto closeto = [&](CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL r) {return fabs(p.x - x) <= r && fabs(p.y - y) <= r && fabs(p.z - z) <= r;};constexpr CCTK_REAL x = 0.25;constexpr CCTK_REAL y = 0.25;constexpr CCTK_REAL z = 0.25;constexpr CCTK_REAL r = 0.2;regrid_error[p.idx] = // closeto(-x, -y, -z, r) ||closeto(+x, -y, -z, r) || closeto(-x, +y, -z, r) ||closeto(+x, +y, -z, r) || closeto(-x, -y, +z, r) ||closeto(+x, -y, +z, r) || closeto(-x, +y, +z, r) ||closeto(+x, +y, +z, r);});#endif#if 1Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {CCTK_REAL err = 0;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_);regrid_error[p.idx] = err;
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;
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, 0, 0> dive_(cctkGH, dive);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> divb_(cctkGH, divb);const auto nancheck = [&](const CCTK_REAL *restrict var,const Loop::PointDesc &p) {if (CCTK_isnan(var[p.idx]))CCTK_VERROR("Found nan at [%d,%d,%d]", p.i, p.j, p.k);};
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_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ax, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ay, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(az, p); });
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;});
Loop::loop_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ex, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ey, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(ez, p); });Loop::loop_all<0, 1, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(bx, p); });Loop::loop_all<1, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(by, p); });Loop::loop_all<1, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(bz, p); });Loop::loop_all<0, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(rho, p); });Loop::loop_all<1, 0, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jx, p); });Loop::loop_all<0, 1, 0>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jy, p); });Loop::loop_all<0, 0, 1>(cctkGH,[&](const Loop::PointDesc &p) { nancheck(jz, p); });// Loop::loop_all<0, 0, 0>(cctkGH,// [&](const Loop::PointDesc &p) { nancheck(dive, p);// });// Loop::loop_all<1, 1, 1>(cctkGH,// [&](const Loop::PointDesc &p) { nancheck(divb, p);// });