NUOLOGCKMF5UOBGBYEOX4O7NQ5AEVVLCH6KRBQRJQXIRDNJ2C2ZQC RZPXXQUUAG7JIADVCLYQOY3I4CQDR23PF57L6XDGHNCT4EBPSO6QC FWX77KPLT4HVGASOPYOJAMG7DT5HLCNMZ42UMG2J7LYX5SJJSSDQC UNLFJ6OBGZOKEVGFS3DEF75WM2BY4QKYVNIT443DINI2JUU4WYBQC 36DLQFITHIMQZEQHPJIEHLIWRVTSD7UZSWCQAIT5ZY3JU2W66CFQC O4RDTKKLC3OK5MSXNMTMOXOONLV3XSMH4OPZMANRBVYRFVJJ5KGQC ES2T6C6XURRCKVEFML2YPDBARPNZZNR3E77CT2LWU5LWPSSTCIQQC LJ65UCZG376GW736Z4KJUPCD6WNZB6DH2W6BCFNBGQZWZOOOOUVAC JUWPJ45RRUIBK7GC3TQGJEDV4P5NUK3PDMQTLAWSN6PMSM5V7V5AC M7FZ5YH4U6TXWMCSM6A5GDXXUH636DROLQJPZV5SG52IR7227VJQC D3M3FHH64RTQZV75QAOSF7BCT2VOFK6GN4GBHBTYZMRR2TKUO5BAC JCY6DNGV3U62ISIUZDX3QDRIESX2RZ4A4MEH3QHETQSAKXI345TAC OXDROTMRIJ42WSAJSZAZSDM2LDRDIBQH62VUGNAIID5BMWN3VETQC RTOD77HMOIHRJW4EYT2X4BXSLRLI5TTECWHZPDLQO766WOR7VS6QC VWXEYEDCJ3R7JN4GEZ23FAAZVOG5GOSF2Z6IFVKMLCFEBMBEJS4QC 67UXOK4QSZFAPVZXU2MJMD7I2232JCK4QDMIOMVMEQNQUGN7632AC 2CPXGSD5DDJJQ3LXW5LEW5UJYAQYX3ISVCRIMGDQJ5IQYYPCJNEAC 3BFBGMNXEPJUUTSJDJ222PQATNDNX4GGRJA2PZBPMWFMG2ZRGH2AC FHE6M7GL6TJSS4VHGV6XY5PQSZBCU6SIXAGZG67Z76NRGKZZAGTQC GP6UYE5IVVG25EHEMT4EPKFUUWI2QUZSJCRSHOIVTM4AX4IVXBLAC X7A7HYGU2NPWCZCG5SYFD4NINYDP6UYWFFE3WP6YNLSDIABG2KJAC HVV73EHEMEFVCPS75XZAFUGWOY7LBYOJOIDNHJUK4MB5QTS3CTJQC 722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC 2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC BVO33OXTG3QDBJ5YLBBVMUDHZGVFE777UOGZKVAHQK6MSDVR5RRAC BJDGFYBMECTJG7BHLNHLSCUCBVYHAY6OGY37FIJP6JDGNDXQNQVAC 5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC QNV4LD7UGYSSNYDXGJC6SRP6Y3USUVDQLEERBMBWRC7LILWREB7AC UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC 6PYEITKIMEZZLF2UZMYD5QJWQYCXKGWTRBE7OENW5EEAXAQSKHPQC N42E473LAYAM5NC4UXPM6WD37FLCUM26JCCKCARQEMXM7V6GWYNQC DJPITUE36JN2YUFH5ZO6M33SEROO5IDTTXIIL4CI7IHLMK6V5STQC VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC FS7Q6TUHBK5WSRDC3TM6KV2BPGWATRBLDHFGEJ72BR3FRDEOC3WAC M5R6KQLXLGYSVKHVAX5AJKD6NYE6IM5Z6WVTR3BTKPJDNNKF3ARAC 5UMVJBOXMEKW2H5JG7T4APWP4QSLSPQSXW6RJBZS6FQQ3BG2XXQAC FUMQRGIGTXORIGJACBDFANU6ZOHBT45EVI4IZTPR7IE2NW2XE77AC TUDUMVD5MTJJJGVYVOCALOOKOUYNB7LOFJRCYCZDVGN536JAO2OQC OHAARXLBTFOIV5FOTY7IDMOIAE275CKXCE6XQMWV3EC63YCSMIWAC U4IDI3M3RBNIQIY43PFKMPBHVCKG37FI5ARJ4GOKY7KKCD2DILWAC L2PGFWL2UYCKFV7XQTNLUESMHCP2TBNE4IWFBDU4E4UNBNEMSKXAC QU4IBVXLLNKNFHT25ZP23INGTVGISDTTLPUBD5S7PCDBC3ZXF3QAC KOOPA77MHDMNHAKH7YXUD2POVQCT3R2HHDHEFBW3G7MVM7QJQVYAC KK3NMMOX7Z77JZ4D6UDD44VVQOGYRUN4SSFCSHJDELRFZMQE7HXAC T3TZRPPAIA24I3YL3JFB4XEAYCWU3HJAJUCF7NNIFMP4I5X4SM5QC DTKQAPJB6JCJZYLQ4MOFA5QHDCCIQ3KG5ED6PVF5PGEH33YSJFBAC NH72X2CBW6LZIU354XZFSV7TSVZFDVESLCUSAIUZGDMSQMU4BQYAC DGEGTDYVGRHEGNFDHCCVTD6AVLRJ2XHDXQPWR5MMR62V5NU2XGIQC RKGV7HOLE7Q4PVL3B5NWS7X4BBIT7WCKTBANPIGR6KMFKGZENTVAC EK7F2IQUVSSHTDV7MPOZXX34UHJA2ZNDHIZX4G4HCSUJYXKFWGRAC DQZRKU4B6C3WWWXFMBBSVKYD5AXDGUMHHPACVWKERGKB5WLHZGYAC C7Y67X3DJBDM4QW3NAVSVJF6GSSLD466S4FFJSIKUEJANLPC6FIAC I4P6OKQG6OD3JX2KV23QLKMG7DSNYMTYUZCIOKNPTHOFHOFLMSNAC KG47IF4CPCUT3BHS34WDRHTH5HYMBTY4OSTB3X7APR2E5ZJ32IYQC BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC JD6PQOJ6YYNQYEEWEXO2NM7NVYNBUI6V7ZU6Q3FNHGAT2VYOF5WAC BSMJ4V7GV3EOGY4KCSTOJQUOFE2OOCIKQETE4WC2WRNLWBQIBW3QC XU5HOJREK4XY4NBCJINLZPKQNKSYOLUDTWR47REFSNQKDOSNDXLQC VMCDMDXKME66ESRMB3PYSUZZH2XG2GIQMEOEKRH33WGCEBPXTWUQC GKKJ75HX2ERLVBZVE2CUB6T3J2SUT7R3UKEKTEYNOG43ZKX6X5MQC P2M2LDDITHKC7HC6RMWQDTS5CT6CCH4LSIN37HQUA3NTC6SFNZOQC V37ACC7MAD57VQOKAT5FBO2ITNX73AWD2DUEIHMF7UPKV63YRTQAC LAPWK2M55JGEUK5ZGFN5DKUOAV3HJYUBWHO7ZWVU7FBYMCQXHQZAC 7XNNFKOUEK2MQEI67B3U3ZLGLTTVZPAEJWFHAE6JGMMRT2B7WKBAC ULG227G2VMOOSDSKHY36QQQGEABGXZIX2B3JA4YIBS6HCVDVJ7IAC AUQBSIC45Q6FWYPSNFO3SM35WZ7DJZHXBANWY7RM72UHTV4IHGEQC VJOYOPMITY6PJCBSZ7HCYOSRDMAEB2IJGYYXVDOOS6WLUE5XZQSAC TKGR5VICXNXS4LWLDKKNRQ6BVOEZTJOL66D6YFSPFTXLLSFDLXZAC GLSGWVXRSIWTKLLTKJTSTDL6IISOKX77Z23MQRJUZ5TCSMA7MF6AC HNQLJR6ZWME6VBJ2Y7PSENLJPXC7INSS7NC2CKIWQAC776CQ74TQC 2YLGRIZKWTRG6Q3BU3W6ZJICHNEQBWJDEKBWOMM7YLGZREVV5GUQC CK7YISOYBIBXPO7OHYSA7GKFX3QA3IB7V2K2VIB5FYEBRRU2THZQC TWBYQDSGSRLSHPYATCLQBLSWBVPMKRBJYCPZ4ADC2LDV5X64PITAC W5IGLN4UQIR5MF5MIQT25JQZ5F2DGRRHJFLN3WDL3T6PLSWRS5DAC K5LKSBGB5T27IFKIIMA6COTUL3MXPZBMAWVZQFKJSD52X3AIUWAQC CarpetX/examples
CarpetX::ghost_size = 2
CarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yes
CarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yes
ActiveThorns = "ADMBaseBaikalXBrillLindquistCarpetXErrorEstimatorIOUtilODESolvers"$nlevels = 1$ncells = 64Cactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 1.0CarpetX::verbose = yesCarpetX::xmin = -1.6 #TODO -16.0CarpetX::ymin = -1.6 #TODO -16.0CarpetX::zmin = -1.6 #TODO -16.0CarpetX::xmax = +1.6 #TODO +16.0CarpetX::ymax = +1.6 #TODO +16.0CarpetX::zmax = +1.6 #TODO +16.0CarpetX::ncells_x = $ncellsCarpetX::ncells_y = $ncellsCarpetX::ncells_z = $ncellsCarpetX::ghost_size = 2CarpetX::max_num_levels = $nlevelsCarpetX::regrid_every = 16CarpetX::regrid_error_threshold = 1.0 #TODO 1.0 / 16.0ErrorEstimator::region_shape = "cube"ErrorEstimator::scale_by_resolution = no #TODO yesCarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 3CarpetX::dtfac = 0.25ADMBase::initial_data = "Brill-Lindquist"ADMBase::initial_lapse = "Brill-Lindquist"IO::out_dir = $parfileIO::out_every = 1 #TODO $ncells * 2 ** ($nlevels - 1) / 32CarpetX::out_tsv = yes
constexpr array<int, 3> CarpetX_widths{16, 5, 5};
template <typename T> constexpr T lapse(const T x, const T y, const T z) {constexpr T m = 1.0;const T r = smooth(sqrt(pow2(expand(x)) + pow2(expand(y)) + pow2(expand(z))));const T lapse1 = (1 - m / (2 * r)) / (1 + m / (2 * r));return (1 + lapse1) / 2; // average}
// width: 4 + 5 + 5, height: 5, XZ" CCC AAA RRRR ", //"C A A R R", //"C AAAAA RRRR ", //"C A A R R ", //" CCC A A R RR", //
// width: 5 + 5 + 5, height: 5, XZ" CCCC AAA RRRR ", //"C A A R R", //"C AAAAA RRRR ", //"C A A R R ", //" CCCC A A R RR", //
USES KEYWORD out_modeUSES CCTK_INT out_proc_every
////////////////////////////////////////////////////////////////////////////////#ifdef HAVE_CAPABILITY_Silostruct mesh_props_t {IntVect ngrow;auto to_tuple() const { return make_tuple(ngrow); }friend bool operator==(const mesh_props_t &p, const mesh_props_t &q) {return p.to_tuple() == q.to_tuple();}friend bool operator<(const mesh_props_t &p, const mesh_props_t &q) {return p.to_tuple() < q.to_tuple();}};void OutputSilo(const cGH *restrict const cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;int ierr;// Set up timersstatic Timer timer("OutputSilo");Interval interval(timer);const int myproc = CCTK_MyProc(cctkGH);const int nprocs = CCTK_nProcs(cctkGH);// Notes for laterout_mode;out_proc_every;// Configure Silo libraryDBShowErrors(DB_ALL_AND_DRVR, nullptr);// DBSetAllowEmptyObjects(1);DBSetCompression("METHOD=GZIP");DBSetEnableChecksums(1);
// Create output fileconst string parfilename = [&]() {string buf(1024, '\0');const int len =CCTK_ParameterFilename(buf.length(), const_cast<char *>(buf.data()));buf.resize(len);return buf;}();const string name = [&]() {string name = parfilename;const size_t last_slash = name.rfind('/');if (last_slash != string::npos && last_slash < name.length())name = name.substr(last_slash + 1);const size_t last_dot = name.rfind('.');if (last_dot != string::npos && last_dot > 0)name = name.substr(0, last_dot);return name;}();const string filename = [&]() {ostringstream buf;buf << out_dir << "/" << name << ".it" << setw(6) << setfill('0')<< cctk_iteration << ".p" << setw(6) << setfill('0') << myproc<< ".silo";return buf.str();}();const DB::ptr<DBfile> file = DB::make(DBCreate(filename.c_str(), DB_CLOBBER, DB_LOCAL, name.c_str(), DB_HDF5));assert(file);// TODO: directories instead of carefully chosen names// Find output groupsconst vector<bool> group_enabled = [&]() {vector<bool> enabled(CCTK_NumGroups(), false);const auto callback{[](const int index, const char *const optstring, void *const arg) {vector<bool> &enabled = *static_cast<vector<bool> *>(arg);enabled.at(CCTK_GroupIndexFromVarI(index)) = true;}};CCTK_TraverseString(out_silo_vars, callback, &enabled, CCTK_GROUP_OR_VAR);return enabled;}();// Iterate over all groups that should be outputmap<mesh_props_t, string> all_multimeshnames;for (int gi = 0; gi < CCTK_NumGroups(); ++gi) {if (!group_enabled.at(gi))continue;if (CCTK_GroupTypeI(gi) != CCTK_GF)continue;// Write groupconst int v0 = CCTK_FirstVarIndexI(gi);const int nv = CCTK_NumVarsInGroupI(gi);for (int vi = 0; vi < nv; ++vi) {const auto &leveldata0 = *ghext->leveldata.begin();// const auto &groupdata0 = *leveldata0.groupdata.at(gi);const int tl = 0;const MultiFab &mfab0 = *leveldata0.groupdata.at(gi)->mfab[tl];const IndexType &indextype = mfab0.ixType();const IntVect &ngrow = mfab0.nGrowVect();const mesh_props_t mesh_props{ngrow};const bool have_multimesh = all_multimeshnames.count(mesh_props);constexpr int ndims = dim;vector<string> meshnames;vector<array<array<double, ndims>, 2> > meshextents;vector<int> meshzonecounts;vector<string> varnames;for (const auto &leveldata : ghext->leveldata) {const auto &groupdata = *leveldata.groupdata.at(gi);const int tl = 0;const MultiFab &mfab = *groupdata.mfab[tl];assert(mfab.ixType() == indextype);assert(mfab.nGrowVect() == ngrow);for (MFIter mfi(mfab); mfi.isValid(); ++mfi) {// const Box &box = mfi.validbox(); // interiorconst Box &fabbox = mfi.fabbox(); // exteriorconst FArrayBox &fab = mfab[mfi];assert(fab.box() == fabbox);const string meshname = [&]() {ostringstream buf;buf << "box.rl" << setw(2) << setfill('0') << leveldata.level<< ".c" << setw(8) << setfill('0') << mfi.index();return DB::legalize_name(buf.str());}();if (!have_multimesh) {meshnames.push_back(meshname);array<int, ndims> dims_vc;for (int d = 0; d < ndims; ++d)dims_vc[d] = fabbox.length(d) + int(indextype.cellCentered(d));const Geometry &geom = ghext->amrcore->Geom(leveldata.level);const double *const x0 = geom.ProbLo();const double *const dx = geom.CellSize();array<vector<double>, ndims> coords;for (int d = 0; d < ndims; ++d) {coords[d].resize(dims_vc[d]);for (int i = 0; i < dims_vc[d]; ++i)coords[d][i] = x0[d] + (fabbox.smallEnd(d) + i) * dx[d];}array<void *, ndims> coord_ptrs;for (int d = 0; d < ndims; ++d)coord_ptrs[d] = coords[d].data();array<array<double, ndims>, 2> extent;for (int d = 0; d < ndims; ++d) {extent[0][d] = *coords[d].begin();extent[1][d] = *coords[d].rbegin();}meshextents.push_back(extent);int zonecount = 1;for (int d = 0; d < ndims; ++d)zonecount *= dims_vc[d];meshzonecounts.push_back(zonecount);const DB::ptr<DBoptlist> optlist = DB::make(DBMakeOptlist(10));assert(optlist);int cartesian = DB_CARTESIAN;ierr = DBAddOption(optlist.get(), DBOPT_COORDSYS, &cartesian);assert(!ierr);int cycle = cctk_iteration;ierr = DBAddOption(optlist.get(), DBOPT_CYCLE, &cycle);assert(!ierr);array<int, ndims> min_index, max_index;for (int d = 0; d < ndims; ++d) {min_index[d] = ngrow[d];max_index[d] = ngrow[d];}ierr =DBAddOption(optlist.get(), DBOPT_LO_OFFSET, min_index.data());assert(!ierr);ierr =DBAddOption(optlist.get(), DBOPT_HI_OFFSET, max_index.data());assert(!ierr);int column_major = 1;ierr = DBAddOption(optlist.get(), DBOPT_MAJORORDER, &column_major);assert(!ierr);// float time = cctk_time;// ierr = DBAddOption(optlist.get(), DBOPT_TIME, &time);// assert(!ierr);double dtime = cctk_time;ierr = DBAddOption(optlist.get(), DBOPT_DTIME, &dtime);assert(!ierr);int hide_from_gui = 1;ierr =DBAddOption(optlist.get(), DBOPT_HIDE_FROM_GUI, &hide_from_gui);assert(!ierr);ierr = DBPutQuadmesh(file.get(), meshname.c_str(), nullptr,coord_ptrs.data(), dims_vc.data(), ndims,DB_DOUBLE, DB_COLLINEAR, optlist.get());assert(!ierr);}const string varname = [&]() {ostringstream buf;buf << CCTK_FullVarName(v0 + vi) << ".rl" << setw(2) << setfill('0')<< leveldata.level << ".c" << setw(8) << setfill('0')<< mfi.index();return DB::legalize_name(buf.str());}();varnames.push_back(varname);array<int, ndims> dims;for (int d = 0; d < ndims; ++d)dims[d] = fabbox.length(d);const int centering = [&]() {if (indextype.nodeCentered())return DB_NODECENT;if (indextype.cellCentered())return DB_ZONECENT;assert(0);}();const DB::ptr<DBoptlist> optlist = DB::make(DBMakeOptlist(10));assert(optlist);int cartesian = DB_CARTESIAN;ierr = DBAddOption(optlist.get(), DBOPT_COORDSYS, &cartesian);assert(!ierr);int cycle = cctk_iteration;ierr = DBAddOption(optlist.get(), DBOPT_CYCLE, &cycle);assert(!ierr);int column_major = 1;ierr = DBAddOption(optlist.get(), DBOPT_MAJORORDER, &column_major);assert(!ierr);// float time = cctk_time;// ierr = DBAddOption(optlist.get(), DBOPT_TIME, &time);// assert(!ierr);double dtime = cctk_time;ierr = DBAddOption(optlist.get(), DBOPT_DTIME, &dtime);assert(!ierr);// int hide_from_gui = 1;// ierr =// DBAddOption(optlist.get(), DBOPT_HIDE_FROM_GUI,// &hide_from_gui);// assert(!ierr);ierr = DBPutQuadvar1(file.get(), varname.c_str(), meshname.c_str(),fab.dataPtr(vi), dims.data(), ndims, nullptr, 0,DB_DOUBLE, centering, optlist.get());assert(!ierr);} // box} // levelconst string multimeshname = [&]() {ostringstream buf;buf << "multimesh";return DB::legalize_name(buf.str());}();if (!have_multimesh) {all_multimeshnames[mesh_props] = multimeshname;vector<const char *> meshname_ptrs;meshname_ptrs.reserve(meshnames.size());for (const auto &meshname : meshnames)meshname_ptrs.push_back(meshname.c_str());const DB::ptr<DBoptlist> optlist = DB::make(DBMakeOptlist(10));assert(optlist);int cycle = cctk_iteration;ierr = DBAddOption(optlist.get(), DBOPT_CYCLE, &cycle);assert(!ierr);// float time = cctk_time;// ierr = DBAddOption(optlist.get(), DBOPT_TIME, &time);// assert(!ierr);double dtime = cctk_time;ierr = DBAddOption(optlist.get(), DBOPT_DTIME, &dtime);assert(!ierr);int quadmesh = DB_QUADMESH;ierr = DBAddOption(optlist.get(), DBOPT_MB_BLOCK_TYPE, &quadmesh);assert(!ierr);int extents_size = 2 * ndims;assert(sizeof(*meshextents.data()) == extents_size * sizeof(double));ierr = DBAddOption(optlist.get(), DBOPT_EXTENTS_SIZE, &extents_size);assert(!ierr);assert(meshextents.size() == meshname_ptrs.size());ierr = DBAddOption(optlist.get(), DBOPT_EXTENTS, meshextents.data());assert(!ierr);assert(meshzonecounts.size() == meshname_ptrs.size());ierr =DBAddOption(optlist.get(), DBOPT_ZONECOUNTS, meshzonecounts.data());assert(!ierr);ierr = DBPutMultimesh(file.get(), multimeshname.c_str(),meshname_ptrs.size(), meshname_ptrs.data(),nullptr, optlist.get());assert(!ierr);}const string multivarname = [&]() {ostringstream buf;buf << CCTK_FullVarName(v0 + vi);return DB::legalize_name(buf.str());}();vector<const char *> varname_ptrs;varname_ptrs.reserve(varnames.size());for (const auto &varname : varnames)varname_ptrs.push_back(varname.c_str());const DB::ptr<DBoptlist> optlist = DB::make(DBMakeOptlist(10));assert(optlist);int cycle = cctk_iteration;ierr = DBAddOption(optlist.get(), DBOPT_CYCLE, &cycle);assert(!ierr);// float time = cctk_time;// ierr = DBAddOption(optlist.get(), DBOPT_TIME, &time);// assert(!ierr);double dtime = cctk_time;ierr = DBAddOption(optlist.get(), DBOPT_DTIME, &dtime);assert(!ierr);ierr = DBAddOption(optlist.get(), DBOPT_MMESH_NAME,const_cast<char *>(multimeshname.c_str()));assert(!ierr);int vartype_scalar = DB_VARTYPE_SCALAR;ierr = DBAddOption(optlist.get(), DBOPT_TENSOR_RANK, &vartype_scalar);assert(!ierr);int quadvar = DB_QUADVAR;ierr = DBAddOption(optlist.get(), DBOPT_MB_BLOCK_TYPE, &quadvar);assert(!ierr);ierr =DBPutMultivar(file.get(), multivarname.c_str(), varname_ptrs.size(),varname_ptrs.data(), nullptr, optlist.get());assert(!ierr);} // vi} // gi}#endif////////////////////////////////////////////////////////////////////////////////
file << "# 1:iteration" << sep << "2:time" << sep << "3:varname" << sep<< "4:min" << sep << "5:max" << sep << "6:sum" << sep << "7:avg" << sep<< "8:stddev" << sep << "9:volume" << sep << "10:L1norm" << sep<< "11:L2norm" << sep << "12:maxabs"<< "\n";
int col = 0;file << "# " << ++col << ":iteration";file << sep << ++col << ":time";file << sep << ++col << ":varname";file << sep << ++col << ":min";file << sep << ++col << ":max";file << sep << ++col << ":sum";file << sep << ++col << ":avg";file << sep << ++col << ":stddev";file << sep << ++col << ":volume";file << sep << ++col << ":L1norm";file << sep << ++col << ":L2norm";file << sep << ++col << ":maxabs";for (int d = 0; d < dim; ++d)file << sep << ++col << ":minloc[" << d << "]";for (int d = 0; d < dim; ++d)file << sep << ++col << ":maxloc[" << d << "]";file << "\n";
#if 1reduction<CCTK_REAL> red = reduce(gi, vi, tl);#elsereduction<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
const reduction<CCTK_REAL, dim> red = reduce(gi, vi, tl);
constexpr int dim = 3;////////////////////////////////////////////////////////////////////////////////// TODO: This is broken (at least on MacOS, with gcc 9.2, with -O0).// namespace {// template <typename T, size_t N, typename... Ts>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array<T, N>{xs...};// }//// template <typename T, size_t N, typename... Ts>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array_from_initializer_list<T, N>(beg + 1, end, xs...,// beg < end ? *beg : T{});// }// } // namespace//// template <typename T, size_t N>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr array<T, N>// array_from_initializer_list(initializer_list<T> l) {// return array_from_initializer_list<T, N>(l.begin(), l.end());// }// namespace {// template <typename T, size_t N, typename... Ts>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array<T, N>{xs...};// }// template <typename T, size_t N, typename... Ts>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array_from_initializer_list<T, N>(beg, end - 1,// beg < end ? *(end - 1) : T{},// xs...);// }// } // namespace//// template <typename T, size_t N>// /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr array<T, N>// array_from_initializer_list(initializer_list<T> l) {// return array_from_initializer_list<T, N>(l.begin(), l.end());// }template <typename T, size_t N>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr enable_if_t<(N == 3), array<T, N> >array_from_initializer_list(initializer_list<T> l) {assert(l.size() >= N);const T *restrict const p = l.begin();return {p[0], p[1], p[2]};}template <typename T, size_t N>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr enable_if_t<(N == 6), array<T, N> >array_from_initializer_list(initializer_list<T> l) {assert(l.size() >= N);const T *restrict const p = l.begin();return {p[0], p[1], p[2], p[3], p[4], p[5]};}static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[0] == 0,"");static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[1] == 1,"");static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[2] == 2,"");template <typename T, int D> struct vect {array<T, D> elts;// initializes all elts to zero/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect() : elts() {}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect(const array<T, D> &arr): elts(arr) {}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect(initializer_list<T> lst): elts(array_from_initializer_list<T, D>(lst)) {#ifdef CCTK_DEBUGassert(lst.size() == D);#endif}static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect pure(T a) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = a;return r;}static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect unit(int dir) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = d == dir;return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr const T &operator[](int d) const {#ifdef CCTK_DEBUGassert(d >= 0 && d < D);#endifreturn elts[d];}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T &operator[](int d) {#ifdef CCTK_DEBUGassert(d >= 0 && d < D);#endifreturn elts[d];}
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator+(const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = +x.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator-(const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = -x.elts[d];return r;}
using CarpetX::dim;using CarpetX::vect;
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator+(const vect &x, const vect &y) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] + y.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator-(const vect &x, const vect &y) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] - y.elts[d];return r;}// friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect operator+(const vect// &x, const array<T, D> &y) {// return x + vect(y);// }// friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect operator-(const vect// &x, const array<T, D> &y) {// return x - vect(y);// }// friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect operator+(const vect// &x, T a) {// vect r;// for (int d = 0; d < D; ++d)// r.elts[d] = x.elts[d] + a;// return r;// }// friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect operator-(const vect// &x, T a) {// vect r;// for (int d = 0; d < D; ++d)// r.elts[d] = x.elts[d] - a;// return r;// }friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator*(T a, const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = a * x.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vectoperator*(const vect &x, T a) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] * a;return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D> operator!(const vect &x) {vect<bool, D> r;for (int d = 0; d < dim; ++d)r.elts[d] = !x.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator&&(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] && y.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator||(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] || y.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator==(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] == y.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator!=(const vect &x, const vect &y) {return !(x == y);}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator<(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] < y.elts[d];return r;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator>(const vect &x, const vect &y) {return y < x;}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator<=(const vect &x, const vect &y) {return !(x > y);}friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<bool, D>operator>=(const vect &x, const vect &y) {return !(x < y);}template <typename U>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<U, D>ifelse(const Loop::vect<U, D> &x, const Loop::vect<U, D> &y) const {vect<U, D> r;for (int d = 0; d < D; ++d)r.elts[d] = elts[d] ? x.elts[d] : y.elts[d];return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool all() const {bool r = true;for (int d = 0; d < D; ++d)r &= elts[d];return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool any() const {bool r = false;for (int d = 0; d < D; ++d)r |= elts[d];return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T maxabs() const {T r = 0;for (int d = 0; d < D; ++d)r = fmax(r, fabs(elts[d]));return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T prod() const {T r = 1;for (int d = 0; d < D; ++d)r *= elts[d];return r;}/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T sum() const {T r = 0;for (int d = 0; d < D; ++d)r += elts[d];return r;}friend ostream &operator<<(ostream &os, const vect &x) {os << "[";for (int d = 0; d < D; ++d) {if (d > 0)os << ",";os << x.elts[d];}os << "]";return os;}};} // namespace Loopnamespace std {template <typename T, int D> struct equal_to<Loop::vect<T, D> > {/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr booloperator()(const Loop::vect<T, D> &lhs, const Loop::vect<T, D> &rhs) const {return equal_to<array<T, D> >()(lhs.elts, rhs.elts);}};template <typename T, int D> struct less<Loop::vect<T, D> > {/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr booloperator()(const Loop::vect<T, D> &lhs, const Loop::vect<T, D> &rhs) const {return less<array<T, D> >(lhs.elts, rhs.elts);}};template <typename T, int D>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr Loop::vect<T, D>abs(const Loop::vect<T, D> &x) {Loop::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = abs(x.elts[d]);return r;}template <typename T, int D>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr Loop::vect<T, D>max(const Loop::vect<T, D> &x, const Loop::vect<T, D> &y) {Loop::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = max(x.elts[d], y.elts[d]);return r;}template <typename T, int D>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr Loop::vect<T, D>min(const Loop::vect<T, D> &x, const Loop::vect<T, D> &y) {Loop::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = min(x.elts[d], y.elts[d]);return r;}} // namespace stdnamespace Loop {namespace static_test {constexpr vect<int, 3> vect1{0, 1, 2};static_assert(vect1[0] == 0, "");static_assert(vect1[1] == 1, "");static_assert(vect1[2] == 2, "");constexpr vect<vect<int, 3>, 3> vect2{{100, 101, 102},{110, 111, 112},{120, 121, 122},};static_assert(vect2[0][0] == 100, "");static_assert(vect2[0][1] == 101, "");static_assert(vect2[0][2] == 102, "");static_assert(vect2[1][0] == 110, "");static_assert(vect2[1][1] == 111, "");static_assert(vect2[1][2] == 112, "");static_assert(vect2[2][0] == 120, "");static_assert(vect2[2][1] == 121, "");static_assert(vect2[2][2] == 122, "");} // namespace static_test////////////////////////////////////////////////////////////////////////////////
void loop_box(const F &f, const array<int, dim> &restrict imin,const array<int, dim> &restrict imax,const array<int, dim> &restrict inormal) const {
CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_box(const F &f, const array<int, dim> &restrict imin,const array<int, dim> &restrict imax,const array<int, dim> &restrict inormal) const {
void loop_all(const array<int, dim> &group_nghostzones, const F &f) const {constexpr array<int, dim> offset{!CI, !CJ, !CK};
inline CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_all(const array<int, dim> &group_nghostzones, const F &f) const {const array<int, dim> offset{!CI, !CJ, !CK};
void loop_int(const array<int, dim> &group_nghostzones, const F &f) const {constexpr array<int, dim> offset{!CI, !CJ, !CK};
inline CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_int(const array<int, dim> &group_nghostzones, const F &f) const {const array<int, dim> offset{!CI, !CJ, !CK};
void loop_bnd(const array<int, dim> &group_nghostzones, const F &f) const {constexpr array<int, dim> offset{!CI, !CJ, !CK};
inline CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_bnd(const array<int, dim> &group_nghostzones, const F &f) const {const array<int, dim> offset{!CI, !CJ, !CK};
void loop_ghosts_inclusive(const array<int, dim> &group_nghostzones,const F &f) const {constexpr array<int, dim> offset{!CI, !CJ, !CK};
inline CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_ghosts_inclusive(const array<int, dim> &group_nghostzones,const F &f) const {const array<int, dim> offset{!CI, !CJ, !CK};
void loop_ghosts(const array<int, dim> &group_nghostzones, const F &f) const {constexpr array<int, dim> offset{!CI, !CJ, !CK};
inline CCTK_ATTRIBUTE_ALWAYS_INLINE voidloop_ghosts(const array<int, dim> &group_nghostzones, const F &f) const {const array<int, dim> offset{!CI, !CJ, !CK};
int dj, dk;int ni, nj, nk;static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr array<int, dim>indextype() {
const int dj, dk, np;const int ni, nj, nk;vector<remove_cv_t<T> > data;T *restrict const ptr;static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE array<int, dim> indextype() {
: ptr(ptr), dj(di * (cctkGH->cctk_ash[0] + 1 - CI)),
: dj(di * (cctkGH->cctk_ash[0] + 1 - CI)),dk(dj * (cctkGH->cctk_ash[1] + 1 - CJ)),np(dk * (cctkGH->cctk_ash[2] + 1 - CK)),ni(cctkGH->cctk_lsh[0] + 1 - CI), nj(cctkGH->cctk_lsh[1] + 1 - CJ),nk(cctkGH->cctk_lsh[2] + 1 - CK), ptr(ptr) {}inline GF3D(const cGH *restrict cctkGH, allocate): dj(di * (cctkGH->cctk_ash[0] + 1 - CI)),
constexpr int i0min = (ORDER + 1) / 2 - 1;constexpr int i0max = (ORDER + 1) / 2;constexpr int imin = 0;constexpr int imax = (ORDER + 1) / 2 - 1;constexpr int i1min = ORDER - imax;constexpr int i1max = ORDER - imin;static_assert(abs(imin - i0min) <= required_ghosts, "");static_assert(abs(imin - i0max) <= required_ghosts, "");static_assert(abs(imax - i0min) <= required_ghosts, "");static_assert(abs(imax - i0max) <= required_ghosts, "");static_assert(abs(i1min - i0min) <= required_ghosts, "");static_assert(abs(i1min - i0max) <= required_ghosts, "");static_assert(abs(i1max - i0min) <= required_ghosts, "");static_assert(abs(i1max - i0max) <= required_ghosts, "");
for (int k = targetbox.loVect()[2]; k <= targetbox.hiVect()[2]; ++k) {for (int j = targetbox.loVect()[1]; j <= targetbox.hiVect()[1]; ++j) {
constexpr int required_ghosts =interp1d<CENTERING, CONSERVATIVE, ORDER>::required_ghosts;{const IntVect fineind(targetbox.loVect());IntVect crseind = fineind;crseind.getVect()[D] = coarsen(fineind.getVect()[D], 2) - required_ghosts;for (int d = 0; d < 3; ++d)assert(crseind.getVect()[d] >= crsebox.loVect()[d]);for (int d = 0; d < 3; ++d)assert(targetbox.loVect()[d] >= finebox.loVect()[d]);}{const IntVect fineind(targetbox.hiVect());IntVect crseind = fineind;crseind.getVect()[D] = coarsen(fineind.getVect()[D], 2) + required_ghosts;for (int d = 0; d < 3; ++d)assert(crseind.getVect()[d] <= crsebox.hiVect()[d]);for (int d = 0; d < 3; ++d)assert(targetbox.hiVect()[d] <= finebox.hiVect()[d]);}const array<int, 3> imin{targetbox.loVect()[0],targetbox.loVect()[1],targetbox.loVect()[2],};const array<int, 3> imax{targetbox.hiVect()[0] + 1,targetbox.hiVect()[1] + 1,targetbox.hiVect()[2] + 1,};const ptrdiff_t fined0 = finebox.index(IntVect(0, 0, 0));constexpr ptrdiff_t finedi = 1;assert(finebox.index(IntVect(1, 0, 0)) - fined0 == finedi);const ptrdiff_t finedj = finebox.index(IntVect(0, 1, 0)) - fined0;const ptrdiff_t finedk = finebox.index(IntVect(0, 0, 1)) - fined0;const ptrdiff_t crsed0 = crsebox.index(IntVect(0, 0, 0));constexpr ptrdiff_t crsedi = 1;assert(crsebox.index(IntVect(1, 0, 0)) - crsed0 == crsedi);const ptrdiff_t crsedj = crsebox.index(IntVect(0, 1, 0)) - crsed0;const ptrdiff_t crsedk = crsebox.index(IntVect(0, 0, 1)) - crsed0;for (int k = imin[2]; k < imax[2]; ++k) {for (int j = imin[1]; j < imax[1]; ++j) {
#endif#endif// Note: fineind = 2 * coarseind + offconst int ci = D == 0 ? i >> 1 : i;const int cj = D == 1 ? j >> 1 : j;const int ck = D == 2 ? k >> 1 : k;const int off = (D == 0 ? i : D == 1 ? j : k) & 0x1;if (D == 0) {// allow vectorizationconst T *restrict const ptr =&crseptr[crsed0 + ck * crsedk + cj * crsedj + ci * crsedi];const T res0 = interp1d<CENTERING, CONSERVATIVE, ORDER>()(ptr, di, 0);const T res1 = interp1d<CENTERING, CONSERVATIVE, ORDER>()(ptr, di, 1);const T res = off == 0 ? res0 : res1;fineptr[fined0 + k * finedk + j * finedj + i * finedi] = res;} else {fineptr[fined0 + k * finedk + j * finedj + i * finedi] =interp1d<CENTERING, CONSERVATIVE, ORDER>()(&crseptr[crsed0 + ck * crsedk + cj * crsedj + ci * crsedi],di, off);}#ifdef CCTK_DEBUGassert(CCTK_isfinite(fineptr[fined0 + i * finedi + j * finedj + k * finedk]));
static vector<Timer> timers;static bool have_timers = false;const int thread_num = omp_get_thread_num();bool my_have_timers;#pragma omp atomic readmy_have_timers = have_timers;if (!my_have_timers) {#pragma omp critical{#pragma omp atomic readmy_have_timers = have_timers;if (!my_have_timers) {const int num_threads = omp_get_num_threads();timers.reserve(num_threads);for (int i = 0; i < num_threads; ++i) {ostringstream buf;buf << "prolongate_3d_rf2<CENT=" << CENTI << CENTJ << CENTK<< ",CONS=" << CONSI << CONSJ << CONSK << ",ORDER=" << ORDERI<< ORDERJ << ORDERK << ">[thread=" << i << "]";timers.emplace_back(buf.str());}#pragma omp atomic writehave_timers = true;}}}const Timer &timer = timers.at(thread_num);Interval interval(timer);
reduction<T> reduce_array(const Array4<const T> &restrict vars, int n,const array<int, dim> &imin,const array<int, dim> &imax,const Array4<const int> *restrict finemask, T dV) {reduction<T> red;
reduction<T, dim> reduce_array(const Array4<const T> &restrict vars,const int n, const array<int, dim> &imin,const array<int, dim> &imax,const Array4<const int> *restrict const finemask,const vect<T, dim> &x0, const vect<T, dim> &dx) {CCTK_REAL dV = 1.0;for (int d = 0; d < dim; ++d)dV *= dx[d];reduction<T, dim> red;
bool is_masked = finemask && (*finemask)(i, j, k);if (!is_masked)red += reduction<T>(dV, vars(i, j, k, n));
const bool is_masked = finemask && (*finemask)(i, j, k);if (!is_masked) {const vect<T, dim> x{x0[0] + i * dx[0], x0[1] + j * dx[1],x0[2] + k * dx[2]};red += reduction<T, dim>(x, dV, vars(i, j, k, n));}
const CCTK_REAL *restrict dx = geom.CellSize();CCTK_REAL dV = 1.0;for (int d = 0; d < dim; ++d)dV *= dx[d];
const CCTK_REAL *restrict const x01 = geom.ProbLo();const CCTK_REAL *restrict const dx1 = geom.CellSize();const vect<CCTK_REAL, dim> x0{x01[0], x01[1], x01[2]};const vect<CCTK_REAL, dim> dx{dx1[0], dx1[1], dx1[2]};
friend ostream &operator<<(ostream &os, const reduction &red) {return os << "reduction{min:" << red.min << ",max:" << red.max<< ",sum:" << red.sum << ",sum2:" << red.sum2<< ",vol:" << red.vol << ",maxabs:" << red.maxabs<< ",sumabs:" << red.sumabs << ",sum2abs:" << red.sum2abs << "}";}
template <typename T1, int D1>friend ostream &operator<<(ostream &os, const reduction<T1, D1> &red);
sum2abs(x.sum2abs + y.sum2abs) {}
sum2abs(x.sum2abs + y.sum2abs),minloc(x.min <= y.min ? x.minloc : y.minloc),maxloc(x.max >= y.max ? x.maxloc : y.maxloc) {}template <typename T, int D>ostream &operator<<(ostream &os, const reduction<T, D> &red) {return os << "reduction{min:" << red.min << ",max:" << red.max<< ",sum:" << red.sum << ",sum2:" << red.sum2 << ",vol:" << red.vol<< ",maxabs:" << red.maxabs << ",sumabs:" << red.sumabs<< ",sum2abs:" << red.sum2abs << ",minloc:" << red.minloc<< ",maxloc:" << red.maxloc << "}";}
FillPatchSingleLevel(*groupdata.mfab.at(tl), 0.0, {&*groupdata.mfab.at(tl)}, {0.0}, 0,0, groupdata.numvars, ghext->amrcore->Geom(level), physbc, 0);
{static Timer timer("Sync::FillPatchSingleLevel");Interval interval(timer);FillPatchSingleLevel(*groupdata.mfab.at(tl), 0.0,{&*groupdata.mfab.at(tl)}, {0.0}, 0, 0,groupdata.numvars, ghext->amrcore->Geom(level),physbc, 0);}
FillPatchTwoLevels(*groupdata.mfab.at(tl), 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),physbc, 0, physbc, 0, reffact, interpolator, bcs, 0);
{static Timer timer("Sync::FillPatchTwoLevels");Interval interval(timer);FillPatchTwoLevels(*groupdata.mfab.at(tl), 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), physbc, 0, physbc, 0, reffact,interpolator, bcs, 0);}
#ifndef VECT_HXX#define VECT_HXX#include <cctk.h>#undef copysign#undef fpclassify#undef isfinite#undef isinf#undef isnan#undef isnormal#undef signbit#include <array>#include <initializer_list>#include <type_traits>#include <utility>#include <vector>namespace CarpetX {using namespace std;constexpr int dim = 3;////////////////////////////////////////////////////////////////////////////////// TODO: This is broken (at least on MacOS, with gcc 9.2, with -O0).// namespace {// template <typename T, size_t N, typename... Ts>// CCTK_ATTRIBUTE_ALWAYS_INLINE// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array<T, N>{xs...};// }//// template <typename T, size_t N, typename... Ts>// CCTK_ATTRIBUTE_ALWAYS_INLINE// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array_from_initializer_list<T, N>(beg + 1, end, xs...,// beg < end ? *beg : T{});// }// } // namespace//// template <typename T, size_t N>// constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE array<T, N>// array_from_initializer_list(initializer_list<T> l) {// return array_from_initializer_list<T, N>(l.begin(), l.end());// }// namespace {// template <typename T, size_t N, typename... Ts>// CCTK_ATTRIBUTE_ALWAYS_INLINE// constexpr enable_if_t<(sizeof...(Ts) == N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array<T, N>{xs...};// }// template <typename T, size_t N, typename... Ts>// CCTK_ATTRIBUTE_ALWAYS_INLINE// constexpr enable_if_t<(sizeof...(Ts) < N), array<T, N> >// array_from_initializer_list(const T *const beg, const T *const end,// const Ts &... xs) {// return array_from_initializer_list<T, N>(beg, end - 1,// beg < end ? *(end - 1) : T{},// xs...);// }// } // namespace//// template <typename T, size_t N>// constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE array<T, N>// array_from_initializer_list(initializer_list<T> l) {// return array_from_initializer_list<T, N>(l.begin(), l.end());// }template <typename T, size_t N>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 3), array<T, N> >array_from_initializer_list(initializer_list<T> l) {assert(l.size() >= N);const T *restrict const p = l.begin();return {p[0], p[1], p[2]};}template <typename T, size_t N>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 6), array<T, N> >array_from_initializer_list(initializer_list<T> l) {assert(l.size() >= N);const T *restrict const p = l.begin();return {p[0], p[1], p[2], p[3], p[4], p[5]};}template <typename T, size_t N>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 3), array<T, N> >array_from_vector(vector<T> &&v) {assert(v.size() >= N);typename vector<T>::iterator const p = v.begin();return {move(p[0]), move(p[1]), move(p[2])};}template <typename T, size_t N>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(N == 6), array<T, N> >array_from_vector(vector<T> &&v) {assert(v.size() >= N);typename vector<T>::iterator const p = v.begin();return {move(p[0]), move(p[1]), move(p[2]),move(p[3]), move(p[4]), move(p[5])};}static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[0] == 0,"");static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[1] == 1,"");static_assert(static_cast<const array<int, 3> &>(array_from_initializer_list<int, 3>({0, 1, 2}))[2] == 2,"");////////////////////////////////////////////////////////////////////////////////template <typename T, int D> struct vect {array<T, D> elts;// initializes all elts to zeroconstexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect() : elts() {}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const array<T, D> &arr): elts(arr) {}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(initializer_list<T> lst): elts(array_from_initializer_list<T, D>(lst)) {#ifdef CCTK_DEBUGassert(lst.size() == D);#endif}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(const vector<T> &vec): elts(array_from_vector<T, D>(vec)) {#ifdef CCTK_DEBUGassert(vec.size() == D);#endif}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect(vector<T> &&vec): elts(array_from_vector<T, D>(move(vec))) {}static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect pure(T a) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = a;return r;}static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect unit(int dir) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = d == dir;return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE const T &operator[](int d) const {#ifdef CCTK_DEBUGassert(d >= 0 && d < D);#endifreturn elts[d];}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T &operator[](int d) {#ifdef CCTK_DEBUGassert(d >= 0 && d < D);#endifreturn elts[d];}template <typename F,typename R = remove_cv_t<remove_reference_t<result_of_t<F(T)> > > >constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map(const F &f) const {vect<R, D> r;for (int d = 0; d < D; ++d)r.elts[d] = f(elts[d]);return r;}template <typename F, typename U,typename R = remove_cv_t<remove_reference_t<result_of_t<F(T, U)> > > >constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R map2(const F &f,const vect<U, D> &y) const {vect<R, D> r;for (int d = 0; d < D; ++d)r.elts[d] = f(elts[d], y.elts[d]);return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = +x.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = -x.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect &x,const vect &y) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] + y.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect &x,const vect &y) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] - y.elts[d];return r;}// friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect// &x, const array<T, D> &y) {// return x + vect(y);// }// friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect// &x, const array<T, D> &y) {// return x - vect(y);// }// friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator+(const vect// &x, T a) {// vect r;// for (int d = 0; d < D; ++d)// r.elts[d] = x.elts[d] + a;// return r;// }// friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator-(const vect// &x, T a) {// vect r;// for (int d = 0; d < D; ++d)// r.elts[d] = x.elts[d] - a;// return r;// }friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator*(T a,const vect &x) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = a * x.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect operator*(const vect &x,T a) {vect r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] * a;return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D> operator!(const vect &x) {vect<bool, D> r;for (int d = 0; d < dim; ++d)r.elts[d] = !x.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator&&(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] && y.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator||(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] || y.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator==(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] == y.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator!=(const vect &x, const vect &y) {return !(x == y);}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator<(const vect &x, const vect &y) {vect<bool, D> r;for (int d = 0; d < D; ++d)r.elts[d] = x.elts[d] < y.elts[d];return r;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator>(const vect &x, const vect &y) {return y < x;}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator<=(const vect &x, const vect &y) {return !(x > y);}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<bool, D>operator>=(const vect &x, const vect &y) {return !(x < y);}template <typename U>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vect<U, D>ifelse(const CarpetX::vect<U, D> &x, const CarpetX::vect<U, D> &y) const {vect<U, D> r;for (int d = 0; d < D; ++d)r.elts[d] = elts[d] ? x.elts[d] : y.elts[d];return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool all() const {bool r = true;for (int d = 0; d < D; ++d)r &= elts[d];return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE bool any() const {bool r = false;for (int d = 0; d < D; ++d)r |= elts[d];return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T maxabs() const {T r = 0;for (int d = 0; d < D; ++d)r = fmax(r, fabs(elts[d]));return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T prod() const {T r = 1;for (int d = 0; d < D; ++d)r *= elts[d];return r;}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE T sum() const {T r = 0;for (int d = 0; d < D; ++d)r += elts[d];return r;}friend ostream &operator<<(ostream &os, const vect &x) {os << "[";for (int d = 0; d < D; ++d) {if (d > 0)os << ",";os << x.elts[d];}os << "]";return os;}};} // namespace CarpetXnamespace std {template <typename T, int D> struct equal_to<CarpetX::vect<T, D> > {constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE booloperator()(const CarpetX::vect<T, D> &lhs,const CarpetX::vect<T, D> &rhs) const {return equal_to<array<T, D> >()(lhs.elts, rhs.elts);}};template <typename T, int D> struct less<CarpetX::vect<T, D> > {constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE booloperator()(const CarpetX::vect<T, D> &lhs,const CarpetX::vect<T, D> &rhs) const {return less<array<T, D> >(lhs.elts, rhs.elts);}};template <typename T, int D>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>abs(const CarpetX::vect<T, D> &x) {CarpetX::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = abs(x.elts[d]);return r;}template <typename T, int D>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>max(const CarpetX::vect<T, D> &x, const CarpetX::vect<T, D> &y) {CarpetX::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = max(x.elts[d], y.elts[d]);return r;}template <typename T, int D>constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CarpetX::vect<T, D>min(const CarpetX::vect<T, D> &x, const CarpetX::vect<T, D> &y) {CarpetX::vect<T, D> r;for (int d = 0; d < D; ++d)r.elts[d] = min(x.elts[d], y.elts[d]);return r;}} // namespace stdnamespace CarpetX {namespace static_test {constexpr vect<int, 3> vect1{0, 1, 2};static_assert(vect1[0] == 0, "");static_assert(vect1[1] == 1, "");static_assert(vect1[2] == 2, "");constexpr vect<vect<int, 3>, 3> vect2{{100, 101, 102},{110, 111, 112},{120, 121, 122},};static_assert(vect2[0][0] == 100, "");static_assert(vect2[0][1] == 101, "");static_assert(vect2[0][2] == 102, "");static_assert(vect2[1][0] == 110, "");static_assert(vect2[1][1] == 111, "");static_assert(vect2[1][2] == 112, "");static_assert(vect2[2][0] == 120, "");static_assert(vect2[2][1] == 121, "");static_assert(vect2[2][2] == 122, "");} // namespace static_test} // namespace CarpetX#endif // #ifndef VECT_HXX
Cactus Code Thorn SiloAuthor(s) : Erik SchnetterMaintainer(s): Cactus teamLicence : ?--------------------------------------------------------------------------1. PurposeConfigure the Silo library; see<https://wci.llnl.gov/simulation/computer-codes/silo>From the web site:A Mesh and Field I/O Library and Scientific DatabaseSilo is a library for reading and writing a wide variety of scientificdata to binary, disk files. The files Silo produces and the datawithin them can be easily shared and exchanged between whollyindependently developed applications running on disparate computingplatforms. Consequently, Silo facilitates the development of generalpurpose tools for processing scientific data. One of the more populartools that process Silo data files is the VisIt visualization tool.
# Configuration definitions for thorn SiloREQUIRES HDF5PROVIDES Silo{SCRIPT configure.shLANG bashOPTIONS SILO_DIR}REQUIRES Silo
#! /bin/bash################################################################################# Prepare################################################################################# Set up shellif [ "$(echo ${VERBOSE} | tr '[:upper:]' '[:lower:]')" = 'yes' ]; thenset -x # Output commandsfiset -e # Abort on errors################################################################################# Configure Cactus################################################################################if [ -z "${SILO_DIR}" ]; thenecho "BEGIN ERROR"echo "Configuration variable SILO_DIR is not set"echo "END ERROR"exit 1fi# Set options: ${SILO_INC_DIRS="${SILO_DIR}/include"}: ${SILO_LIB_DIRS="${SILO_DIR}/lib"}: ${SILO_LIBS="siloh5"}SILO_INC_DIRS="$(${CCTK_HOME}/lib/sbin/strip-incdirs.sh ${SILO_INC_DIRS})"SILO_LIB_DIRS="$(${CCTK_HOME}/lib/sbin/strip-libdirs.sh ${SILO_LIB_DIRS})"# Pass options to Cactusecho "BEGIN MAKE_DEFINITION"echo "SILO_DIR = ${SILO_DIR}"echo "SILO_INC_DIRS = ${SILO_INC_DIRS}"echo "SILO_LIB_DIRS = ${SILO_LIB_DIRS}"echo "SILO_LIBS = ${SILO_LIBS}"echo "END MAKE_DEFINITION"echo 'INCLUDE_DIRECTORY $(SILO_INC_DIRS)'echo 'LIBRARY_DIRECTORY $(SILO_LIB_DIRS)'echo 'LIBRARY $(SILO_LIBS)'
# Interface definition for thorn SiloIMPLEMENTS: SiloINCLUDES HEADER: silo.hxx IN silo.hxx
# Parameter definitions for thorn Silo
# Schedule definitions for thorn Silo
# Main make.code.defn file for thorn Silo# Source files in this directorySRCS = silo.cxx# Subdirectories containing source filesSUBDIRS =
#include "silo.hxx"#include <cctype>#include <sstream>namespace DB {using namespace std;string legalize_name(const string &name) {ostringstream buf;for (const char c : name) {if (isalnum(c) || c == '_')buf << c;elsebuf << '_';}return buf.str();}} // namespace DB
#ifndef SILO_HXX#define SILO_HXX#include <silo.h>#include <cstdlib>#include <memory>#include <string>#include <type_traits>namespace DB {namespace detail {inline void free(DBfile *const obj) { DBClose(obj); }inline void free(DBmultimesh *const obj) { DBFreeMultimesh(obj); }inline void free(DBmultivar *const obj) { DBFreeMultivar(obj); }inline void free(DBoptlist *const obj) { DBFreeOptlist(obj); }inline void free(DBquadmesh *const obj) { DBFreeQuadmesh(obj); }inline void free(DBquadvar *const obj) { DBFreeQuadvar(obj); }inline void free(DBtoc *const obj) {}// inline void free(char *const obj) { std::free(obj); }} // namespace detailtemplate <typename T> using ptr = std::shared_ptr<T>;template <typename T> ptr<T> make(T *const obj) {return std::shared_ptr<T>(obj,static_cast<void (*)(typename std::remove_cv<T>::type *)>(detail::free));}std::string legalize_name(const std::string &name);} // namespace DB#endif // #ifndef SILO_HXX
ActiveThorns = "ADMBaseBrillLindquistCarpetXErrorEstimatorFormalineIOUtilODESolversTimerReportZ4c"# domain_size = 16.0# domain_spacing = 1.0# fine_box_size = 1.0# fine_box_spacing = 0.0625# dtfac = 0.25Cactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 0.0 #TODO 100.0CarpetX::verbose = noCarpetX::poison_undefined_values = noCarpetX::xmin = -16.0CarpetX::ymin = -16.0CarpetX::zmin = -16.0CarpetX::xmax = 16.0CarpetX::ymax = 16.0CarpetX::zmax = 16.0CarpetX::ncells_x = 32CarpetX::ncells_y = 32CarpetX::ncells_z = 32CarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yesCarpetX::ghost_size = 3CarpetX::max_num_levels = 2CarpetX::regrid_every = 0 #TODO 16CarpetX::regrid_error_threshold = 0.5ErrorEstimator::region_shape = "cube"ErrorEstimator::scale_by_resolution = yesCarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 5ODESolvers::method = "RK4"CarpetX::dtfac = 0.25ADMBase::initial_data = "Brill-Lindquist"ADMBase::initial_lapse = "Brill-Lindquist"Z4c::epsdiss = 0.32IO::out_dir = $parfileIO::out_every = 8CarpetX::out_plotfile_groups = "#TODO ADMBase::metric#TODO ADMBase::curv#TODO ADMBase::lapse#TODO ADMBase::shift#TODO Z4c::chi#TODO Z4c::gamma_tilde#TODO Z4c::K_hat#TODO Z4c::A_tilde#TODO Z4c::Gam_tilde#TODO Z4c::Theta#TODO Z4c::alphaG#TODO Z4c::betaG#TODO Z4c::ZtC#TODO Z4c::HC#TODO Z4c::MtC#TODO Z4c::allC"CarpetX::out_tsv = noTimerReport::out_every = 128TimerReport::out_filename = "TimerReport"TimerReport::output_schedule_timers = noTimerReport::n_top_timers = 100
#!/usr/bin/env pythonfrom math import *import refrom string import Templateimport sys################################################################################# domain_size = 16.0# domain_spacing = 1/2.0# fine_box_size = 1.0# fine_box_spacing = 1/16.0#GOOD fine_box_size = 1.0#GOOD fine_box_spacing = 1/16.0#GOOD fine_box_size = 1/2.0#GOOD fine_box_spacing = 1/32.0#GOOD fine_box_size = 1/4.0#GOOD fine_box_spacing = 1/64.0#GOOD fine_box_size = 1/4.0#GOOD fine_box_spacing = 1/128.0# domain_size = 128.0# domain_spacing = 1.0# fine_box_size = 2.0# fine_box_spacing = 1/64.0domain_size = 16.0 #TODO 128.0domain_spacing = 1.0fine_box_size = 1.0fine_box_spacing = 1/16.0dtfac = 1/4.0################################################################################xmin = -domain_sizexmax = +domain_sizenlevels = 2 #TODO round(log(1.0 * domain_spacing / fine_box_spacing, 2)) + 1ncells = round(2.0 * domain_size / domain_spacing)large_box_size = fine_box_size * 2 ** (nlevels - 1)error_threshold = 1.0 / large_box_sizeout_every = round(1.0 / dtfac * 2 ** (nlevels - 1))parfile = """ActiveThorns = "ADMBaseBrillLindquistCarpetXErrorEstimatorFormalineIOUtilODESolversTimerReportZ4c"# domain_size = $domain_size# domain_spacing = $domain_spacing# fine_box_size = $fine_box_size# fine_box_spacing = $fine_box_spacing# dtfac = $dtfacCactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 0.0 #TODO 100.0CarpetX::verbose = noCarpetX::poison_undefined_values = noCarpetX::xmin = $xminCarpetX::ymin = $xminCarpetX::zmin = $xminCarpetX::xmax = $xmaxCarpetX::ymax = $xmaxCarpetX::zmax = $xmaxCarpetX::ncells_x = $ncellsCarpetX::ncells_y = $ncellsCarpetX::ncells_z = $ncellsCarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yesCarpetX::ghost_size = 3CarpetX::max_num_levels = $nlevelsCarpetX::regrid_every = 0 #TODO 16CarpetX::regrid_error_threshold = $error_thresholdErrorEstimator::region_shape = "cube"ErrorEstimator::scale_by_resolution = yesCarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 5ODESolvers::method = "RK4"CarpetX::dtfac = $dtfacADMBase::initial_data = "Brill-Lindquist"ADMBase::initial_lapse = "Brill-Lindquist"Z4c::epsdiss = 0.32IO::out_dir = $$parfileIO::out_every = $out_everyCarpetX::out_plotfile_groups = "#TODO ADMBase::metric#TODO ADMBase::curv#TODO ADMBase::lapse#TODO ADMBase::shift#TODO Z4c::chi#TODO Z4c::gamma_tilde#TODO Z4c::K_hat#TODO Z4c::A_tilde#TODO Z4c::Gam_tilde#TODO Z4c::Theta#TODO Z4c::alphaG#TODO Z4c::betaG#TODO Z4c::ZtC#TODO Z4c::HC#TODO Z4c::MtC#TODO Z4c::allC"CarpetX::out_silo_vars = "ADMBase::lapse"CarpetX::out_tsv = noTimerReport::out_every = 128TimerReport::out_filename = "TimerReport"TimerReport::output_schedule_timers = noTimerReport::n_top_timers = 100"""open(re.sub(r'(.*)\.rpar$', r'\1.par', sys.argv[0]), 'w').write(re.sub(r'\n *',r'\n', Template(parfile).substitute(locals())))
# CarpetX::periodic_x = no# CarpetX::periodic_y = no# CarpetX::periodic_z = no# CarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yes
loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Load and calculateconst z4c_vars_noderivs<CCTK_REAL> vars(kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,gf_Atyy_, gf_Atyz_, gf_Atzz_,//// gf_Gamtx_, gf_Gamty_, gf_Gamtz_,gf_chi_, gf_chi_, gf_chi_,//gf_Theta_, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_, //p.I);
loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {// Load and calculateconst z4c_vars_noderivs<CCTK_REAL> vars(kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,gf_Atyy_, gf_Atyz_, gf_Atzz_,//// gf_Gamtx_, gf_Gamty_, gf_Gamtz_,gf_chi_, gf_chi_, gf_chi_,//gf_Theta_, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_, //p.I);
// Storevars.g.store(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_, gf_gzz_, p);vars.k.store(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_, gf_kzz_, p);gf_alp_(p.I) = vars.alp;vars.beta.store(gf_betax_, gf_betay_, gf_betaz_, p);});
// Storevars.g.store(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_, gf_gzz_, p.I);vars.k.store(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_, gf_kzz_, p.I);gf_alp_(p.I) = vars.alp;vars.beta.store(gf_betax_, gf_betay_, gf_betaz_, p.I);});
#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>namespace Z4c {using namespace Loop;extern "C" void Z4c_ConstraintBoundaries(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_ConstraintBoundaries;const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCx_(cctkGH, ZtCx);const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCy_(cctkGH, ZtCy);const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCz_(cctkGH, ZtCz);const GF3D<CCTK_REAL, 0, 0, 0> gf_HC_(cctkGH, HC);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCx_(cctkGH, MtCx);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCy_(cctkGH, MtCy);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCz_(cctkGH, MtCz);const GF3D<CCTK_REAL, 0, 0, 0> gf_allC_(cctkGH, allC);loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_ZtCx_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_ZtCy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_ZtCz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_HC_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_MtCx_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_MtCy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_MtCz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_allC_(p.I) = 0; });}} // namespace Z4c
const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);
const mat3<GF3D<const CCTK_REAL, 0, 0, 0>, DN, DN> gf_gammat_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatyy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatyz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatzz));
const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);
const mat3<GF3D<const CCTK_REAL, 0, 0, 0>, DN, DN> gf_At_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atyy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atyz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atzz));
const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);
const vec3<GF3D<const CCTK_REAL, 0, 0, 0>, UP> gf_Gamt_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamtx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamty),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamtz));
//const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dchi_(cctkGH, allocate());const mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN> gf_ddchi_(cctkGH, allocate());calc_derivs2(cctkGH, gf_chi_, gf_dchi_, gf_ddchi_);const mat3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, DN, DN> gf_dgammat_(cctkGH, allocate());const mat3<mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN>, DN, DN> gf_ddgammat_(cctkGH, allocate());calc_derivs2(cctkGH, gf_gammat_, gf_dgammat_, gf_ddgammat_);const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dKh_(cctkGH, allocate());calc_derivs(cctkGH, gf_Kh_, gf_dKh_);const mat3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, DN, DN> gf_dAt_(cctkGH,allocate());calc_derivs(cctkGH, gf_At_, gf_dAt_);
const CCTK_REAL alphaG{1};const vec3<CCTK_REAL, DN> dalphaG{0, 0, 0};const mat3<CCTK_REAL, DN, DN> ddalphaG{0, 0, 0, 0, 0, 0};const vec3<CCTK_REAL, UP> betaG{0, 0, 0};const vec3<vec3<CCTK_REAL, DN>, UP> dbetaG{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};const vec3<mat3<CCTK_REAL, DN, DN>, UP> ddbetaG{{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}};const CCTK_REAL rho{0};const vec3<CCTK_REAL, DN> Si{0, 0, 0};const mat3<CCTK_REAL, DN, DN> Sij{0, 0, 0, 0, 0, 0};
kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,gf_Atyy_, gf_Atyz_, gf_Atzz_, gf_Gamtx_, gf_Gamty_, gf_Gamtz_,gf_Theta_,//// gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_chi_, gf_chi_, gf_chi_, gf_chi_,//p.I, dx);
kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_(p.I), gf_dchi_(p.I), gf_ddchi_(p.I), //gf_gammat_(p.I), gf_dgammat_(p.I), gf_ddgammat_(p.I), //gf_Kh_(p.I), gf_dKh_(p.I), //gf_At_(p.I), gf_dAt_(p.I), //gf_Gamt_(p.I), gf_dGamt_(p.I), //gf_Theta_(p.I), gf_dTheta_(p.I), //alphaG, dalphaG, ddalphaG, //betaG, dbetaG, ddbetaG, //rho, //Si, //Sij);
#ifndef DERIVS_HXX#define DERIVS_HXX#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>#include <algorithm>#include <cassert>#include <cmath>#include <initializer_list>#include <ostream>#include <sstream>#include <type_traits>namespace Z4c {using namespace Loop;using namespace std;constexpr int deriv_order = 4;////////////////////////////////////////////////////////////////////////////////template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T deriv1d(const T *restrict const var,const ptrdiff_t di, const T dx) {switch (deriv_order) {case 2:return -1 / T(2) * (var[-di] - var[+di]) / dx;case 4:return (1 / T(12) * (var[-2 * di] - var[+2 * di]) //- 2 / T(3) * (var[-di] - var[+di])) /dx;default:assert(0);}}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE Tderiv1d_upwind(const T *restrict const var, const ptrdiff_t di, const bool sign,const T dx) {// arXiv:1111.2177 [gr-qc], (71)switch (deriv_order) {case 2:if (sign)// + [ 0 -1 +1 0 0]// + 1/2 [+1 -2 +1 0 0]return (1 / T(2) * var[-2 * di] //- 2 * var[-di] //+ 3 / T(2) * var[0]) /dx;else// + [ 0 0 -1 +1 0 ]// - 1/2 [ 0 0 +1 -2 +1 ]// [ 0 0 -3/2 +2 -1/2]return (-3 / T(2) * var[0] //+ 2 * var[+di] //- 1 / T(2) * var[+2 * di]) /dx;case 4:// A fourth order stencil for a first derivative, shifted by one grid pointif (sign)return (-1 / T(12) * var[-3 * di] //+ 1 / T(2) * var[-2 * di] //- 3 / T(2) * var[-di] //+ 5 / T(6) * var[0] //+ 1 / T(4) * var[+di]) /dx;elsereturn (-1 / T(4) * var[-di] //- 5 / T(6) * var[0] //+ 3 / T(2) * var[+di] //- 1 / T(2) * var[+2 * di] //+ 1 / T(12) * var[+3 * di]) /dx;default:assert(0);}}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T deriv2_1d(const T *restrict const var,const ptrdiff_t di,const T dx) {switch (deriv_order) {case 2:return ((var[-di] + var[+di]) //- 2 * var[0]) /pow2(dx);case 4:return (-1 / T(12) * (var[-2 * di] + var[+2 * di]) //+ 4 / T(3) * (var[-di] + var[+di]) //- 5 / T(2) * var[0]) /pow2(dx);default:assert(0);}}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T deriv2_2d(const T *restrict const var,const ptrdiff_t di,const ptrdiff_t dj, const T dx,const T dy) {array<T, deriv_order + 1> arrx;T *const varx = &arrx[arrx.size() / 2];for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j)varx[j] = deriv1d(&var[j * dj], di, dx);return deriv1d(varx, 1, dy);}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T deriv1d_diss(const T *restrict const var,const ptrdiff_t di,const T dx) {switch (deriv_order) {case 2:return ((var[-2 * di] + var[+2 * di]) //- 4 * (var[-di] + var[+di]) //+ 6 * var[0]) /dx;case 4:return ((var[-3 * di] + var[+3 * di]) //- 6 * (var[-2 * di] + var[+2 * di]) //+ 15 * (var[-di] + var[+di]) //- 20 * var[0]) /dx;default:assert(0);}}////////////////////////////////////////////////////////////////////////////////template <int dir, typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T deriv(const GF3D<const T, 0, 0, 0> &gf_,const vect<int, dim> &I,const vec3<T, UP> &dx) {const auto &DI = vect<int, dim>::unit;const ptrdiff_t di = gf_.offset(DI(dir));return deriv1d(&gf_(I), di, dx(dir));}template <int dir, typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE Tderiv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const bool sign, const vec3<T, UP> &dx) {const auto &DI = vect<int, dim>::unit;const ptrdiff_t di = gf_.offset(DI(dir));return deriv1d_upwind(&gf_(I), di, sign, dx(dir));}template <int dir1, int dir2, typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(dir1 == dir2), T>deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {const auto &DI = vect<int, dim>::unit;const ptrdiff_t di = gf_.offset(DI(dir1));return deriv2_1d(&gf_(I), di, dx(dir1));}template <int dir1, int dir2, typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE enable_if_t<(dir1 != dir2), T>deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {const auto &DI = vect<int, dim>::unit;const ptrdiff_t di = gf_.offset(DI(dir1));const ptrdiff_t dj = gf_.offset(DI(dir2));return deriv2_2d(&gf_(I), di, dj, dx(dir1), dx(dir2));}template <int dir, typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE Tderiv_diss(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {const auto &DI = vect<int, dim>::unit;const ptrdiff_t di = gf_.offset(DI(dir));return deriv1d_diss(&gf_(I), di, dx(dir));}////////////////////////////////////////////////////////////////////////////////template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE vec3<T, DN>deriv(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {return {deriv<0>(gf_, I, dx),deriv<1>(gf_, I, dx),deriv<2>(gf_, I, dx),};}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE vec3<T, DN>deriv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dir, const vec3<T, UP> &dx) {return {deriv_upwind<0>(gf_, I, signbit(dir(0)), dx),deriv_upwind<1>(gf_, I, signbit(dir(1)), dx),deriv_upwind<2>(gf_, I, signbit(dir(2)), dx),};}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE mat3<T, DN, DN>deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {return {deriv2<0, 0>(gf_, I, dx), deriv2<0, 1>(gf_, I, dx),deriv2<0, 2>(gf_, I, dx), deriv2<1, 1>(gf_, I, dx),deriv2<1, 2>(gf_, I, dx), deriv2<2, 2>(gf_, I, dx),};}template <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T diss(const GF3D<const T, 0, 0, 0> &gf_,const vect<int, dim> &I,const vec3<T, UP> &dx) {// arXiv:gr-qc/0610128, (63), with r=2constexpr int diss_order = deriv_order + 2;constexpr int sign = diss_order % 4 == 0 ? -1 : +1;return sign / T(pown(2, deriv_order + 2)) *(deriv_diss<0>(gf_, I, dx) //+ deriv_diss<1>(gf_, I, dx) //+ deriv_diss<2>(gf_, I, dx));}////////////////////////////////////////////////////////////////////////////////template <typename T>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs(const cGH *restrict const cctkGH, const GF3D<const T, 0, 0, 0> &gf_,const vec3<GF3D<T, 0, 0, 0>, DN> &dgf_) {DECLARE_CCTK_ARGUMENTS;const vec3<CCTK_REAL, UP> dx([&](int a) { return CCTK_DELTA_SPACE(a); });loop_int<0, 0, 0>(cctkGH,[&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {const auto dval = deriv(gf_, p.I, dx);for (int a = 0; a < 3; ++a)dgf_(a)(p.I) = dval(a);});}template <typename T>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs2(const cGH *restrict const cctkGH,const GF3D<const T, 0, 0, 0> &gf_,const vec3<GF3D<T, 0, 0, 0>, DN> &dgf_,const mat3<GF3D<T, 0, 0, 0>, DN, DN> &ddgf_) {DECLARE_CCTK_ARGUMENTS;const vec3<CCTK_REAL, UP> dx([&](int a) { return CCTK_DELTA_SPACE(a); });loop_int<0, 0, 0>(cctkGH,[&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {const auto dval = deriv(gf_, p.I, dx);for (int a = 0; a < 3; ++a)dgf_(a)(p.I) = dval(a);const auto ddval = deriv2(gf_, p.I, dx);for (int a = 0; a < 3; ++a)for (int b = a; b < 3; ++b)ddgf_(a, b)(p.I) = ddval(a, b);});}template <typename T, dnup_t dnup>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs(const cGH *restrict const cctkGH,const vec3<GF3D<const T, 0, 0, 0>, dnup> &gf_,const vec3<vec3<GF3D<T, 0, 0, 0>, DN>, dnup> &dgf_) {for (int a = 0; a < 3; ++a)calc_derivs(cctkGH, gf_(a), dgf_(a));}template <typename T, dnup_t dnup>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs2(const cGH *restrict const cctkGH,const vec3<GF3D<const T, 0, 0, 0>, dnup> &gf_,const vec3<vec3<GF3D<T, 0, 0, 0>, DN>, dnup> &dgf_,const vec3<mat3<GF3D<T, 0, 0, 0>, DN, DN>, dnup> &ddgf_) {for (int a = 0; a < 3; ++a)calc_derivs2(cctkGH, gf_(a), dgf_(a), ddgf_(a));}template <typename T, dnup_t dnup1, dnup_t dnup2>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs(const cGH *restrict const cctkGH,const mat3<GF3D<const T, 0, 0, 0>, dnup1, dnup2> &gf_,const mat3<vec3<GF3D<T, 0, 0, 0>, DN>, dnup1, dnup2> &dgf_) {for (int a = 0; a < 3; ++a)for (int b = a; b < 3; ++b)calc_derivs(cctkGH, gf_(a, b), dgf_(a, b));}template <typename T, dnup_t dnup1, dnup_t dnup2>CCTK_ATTRIBUTE_NOINLINE voidcalc_derivs2(const cGH *restrict const cctkGH,const mat3<GF3D<const T, 0, 0, 0>, dnup1, dnup2> &gf_,const mat3<vec3<GF3D<T, 0, 0, 0>, DN>, dnup1, dnup2> &dgf_,const mat3<mat3<GF3D<T, 0, 0, 0>, DN, DN>, dnup1, dnup2> &ddgf_) {for (int a = 0; a < 3; ++a)for (int b = a; b < 3; ++b)calc_derivs2(cctkGH, gf_(a, b), dgf_(a, b), ddgf_(a, b));}template <typename T>CCTK_ATTRIBUTE_NOINLINE voidapply_upwind(const cGH *restrict const cctkGH,const GF3D<const T, 0, 0, 0> &gf_,const vec3<GF3D<const T, 0, 0, 0>, UP> &gf_betaG_,const GF3D<T, 0, 0, 0> &gf_rhs_) {DECLARE_CCTK_ARGUMENTS;const vec3<CCTK_REAL, UP> dx([&](int a) { return CCTK_DELTA_SPACE(a); });loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {const vec3<CCTK_REAL, UP> betaG = gf_betaG_(p.I);const vec3<CCTK_REAL, DN> dgf_upwind(deriv_upwind(gf_, p.I, betaG, dx));gf_rhs_(p.I) += sum1([&](int x) { return betaG(x) * dgf_upwind(x); });});}template <typename T>CCTK_ATTRIBUTE_NOINLINE void apply_diss(const cGH *restrict const cctkGH,const GF3D<const T, 0, 0, 0> &gf_,const GF3D<T, 0, 0, 0> &gf_rhs_) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;if (epsdiss == 0)return;const vec3<CCTK_REAL, UP> dx([&](int a) { return CCTK_DELTA_SPACE(a); });loop_int<0, 0, 0>(cctkGH,[&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {gf_rhs_(p.I) += epsdiss * diss(gf_, p.I, dx);});}} // namespace Z4c#endif // #ifndef DERIVS_HXX
loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst mat3<CCTK_REAL, DN, DN> gammat_old(gf_gammatxx_, gf_gammatxy_,gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p.I);const mat3<CCTK_REAL, DN, DN> At_old(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_,gf_Atyz_, gf_Atzz_, p.I);
loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {// Loadconst mat3<CCTK_REAL, DN, DN> gammat_old(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p.I);const mat3<CCTK_REAL, DN, DN> At_old(gf_Atxx_, gf_Atxy_, gf_Atxz_,gf_Atyy_, gf_Atyz_, gf_Atzz_, p.I);
const CCTK_REAL detgammat_old = gammat_old.det();const CCTK_REAL chi_old = 1 / cbrt(detgammat_old);const mat3<CCTK_REAL, DN, DN> gammat([&](int a, int b) { return chi_old * gammat_old(a, b); });
const CCTK_REAL detgammat_old = gammat_old.det();const CCTK_REAL chi_old = 1 / cbrt(detgammat_old);const mat3<CCTK_REAL, DN, DN> gammat([&](int a, int b) { return chi_old * gammat_old(a, b); });
const CCTK_REAL detgammat = gammat.det();const CCTK_REAL gammat_norm = gammat.maxabs();const CCTK_REAL gammat_scale = gammat_norm;if (!(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale)) {ostringstream buf;buf << "det gammat is not one: gammat=" << gammat<< " det(gammat)=" << detgammat;CCTK_VERROR("%s", buf.str().c_str());}assert(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale);
const CCTK_REAL detgammat = gammat.det();const CCTK_REAL gammat_norm = gammat.maxabs();const CCTK_REAL gammat_scale = gammat_norm;if (!(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale)) {ostringstream buf;buf << "det gammat is not one: gammat=" << gammat<< " det(gammat)=" << detgammat;CCTK_VERROR("%s", buf.str().c_str());}assert(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale);
const CCTK_REAL traceAt_old =sum2([&](int x, int y) { return gammatu(x, y) * At_old(x, y); });const mat3<CCTK_REAL, DN, DN> At([&](int a, int b) {return At_old(a, b) - traceAt_old / 3 * gammat(a, b);});
const CCTK_REAL traceAt_old =sum2([&](int x, int y) { return gammatu(x, y) * At_old(x, y); });const mat3<CCTK_REAL, DN, DN> At([&](int a, int b) {return At_old(a, b) - traceAt_old / 3 * gammat(a, b);});
const CCTK_REAL traceAt =sum2([&](int x, int y) { return gammatu(x, y) * At(x, y); });const CCTK_REAL gammatu_norm = gammatu.maxabs();const CCTK_REAL At_norm = At.maxabs();const CCTK_REAL At_scale = fmax(fmax(gammat_norm, gammatu_norm), At_norm);if (!(fabs(traceAt) <= 1.0e-12 * At_scale)) {ostringstream buf;buf << "tr At: At=" << At << " tr(At)=" << traceAt;CCTK_VERROR("%s", buf.str().c_str());}assert(fabs(traceAt) <= 1.0e-12 * At_scale);
const CCTK_REAL traceAt =sum2([&](int x, int y) { return gammatu(x, y) * At(x, y); });const CCTK_REAL gammatu_norm = gammatu.maxabs();const CCTK_REAL At_norm = At.maxabs();const CCTK_REAL At_scale =fmax(fmax(gammat_norm, gammatu_norm), At_norm);if (!(fabs(traceAt) <= 1.0e-12 * At_scale)) {ostringstream buf;buf << "tr At: At=" << At << " tr(At)=" << traceAt;CCTK_VERROR("%s", buf.str().c_str());}assert(fabs(traceAt) <= 1.0e-12 * At_scale);
// Storegammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p);At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_, p);});
// Storegammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p.I);At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_,p.I);});
return mat3<vec3<T, DN>, UP, UP>([&](int a, int b) {return vec3<T, DN>([&](int c) {return sum2([&](int x, int y) { return -gu(a, x) * gu(b, y) * dg(x, y)(c); });});});
return mat3<vec3<T, DN>, UP, UP>([&](int a, int b) CCTK_ATTRIBUTE_ALWAYS_INLINE {return vec3<T, DN>([&](int c) CCTK_ATTRIBUTE_ALWAYS_INLINE {return sum2([&](int x, int y) { return -gu(a, x) * gu(b, y) * dg(x, y)(c); });});});
return mat3<vec3<T, DN>, UP, UP>([&](int a, int b) {return vec3<T, DN>([&](int c) {return sum2([&](int x, int y) {return dgu(a, x)(c) * gu(b, y) * A(x, y) //+ gu(a, x) * dgu(b, y)(c) * A(x, y) //+ gu(a, x) * gu(b, y) * dA(x, y)(c);
return mat3<vec3<T, DN>, UP, UP>([&](int a, int b) CCTK_ATTRIBUTE_ALWAYS_INLINE {return vec3<T, DN>([&](int c) CCTK_ATTRIBUTE_ALWAYS_INLINE {return sum2([&](int x, int y) CCTK_ATTRIBUTE_ALWAYS_INLINE {return dgu(a, x)(c) * gu(b, y) * A(x, y) //+ gu(a, x) * dgu(b, y)(c) * A(x, y) //+ gu(a, x) * gu(b, y) * dA(x, y)(c);});});
return vec3<mat3<T, DN, DN>, UP>([&](int a) {return mat3<T, DN, DN>([&](int b, int c) {return sum1([&](int x) { return gu(a, x) * Gammal(x)(b, c); });
return vec3<mat3<T, DN, DN>, UP>([&](int a) CCTK_ATTRIBUTE_ALWAYS_INLINE {return mat3<T, DN, DN>([&](int b, int c) CCTK_ATTRIBUTE_ALWAYS_INLINE {return sum1([&](int x) CCTK_ATTRIBUTE_ALWAYS_INLINE {return gu(a, x) * Gammal(x)(b, c);});
template <typename T>void apply_upwind(const cGH *restrict const cctkGH,const GF3D<const T, 0, 0, 0> gf_,const GF3D<const T, 0, 0, 0> gf_betaGx_,const GF3D<const T, 0, 0, 0> gf_betaGy_,const GF3D<const T, 0, 0, 0> gf_betaGz_,const GF3D<T, 0, 0, 0> gf_rhs_) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;const vec3<CCTK_REAL, UP> dx{CCTK_DELTA_SPACE(0),CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),};loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {const vec3<CCTK_REAL, UP> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p.I);const vec3<CCTK_REAL, DN> dgf_upwind(deriv_upwind(gf_, p.I, betaG, dx));gf_rhs_(p.I) += sum1([&](int x) { return betaG(x) * dgf_upwind(x); });});}template <typename T>void apply_diss(const cGH *restrict const cctkGH,const GF3D<const T, 0, 0, 0> gf_,const GF3D<T, 0, 0, 0> gf_rhs_) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;const vec3<CCTK_REAL, UP> dx{CCTK_DELTA_SPACE(0),CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),};
const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);
const mat3<GF3D<const CCTK_REAL, 0, 0, 0>, DN, DN> gf_gammat_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatxz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatyy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatyz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, gammatzz));
const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);
const mat3<GF3D<const CCTK_REAL, 0, 0, 0>, DN, DN> gf_At_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atxz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atyy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atyz),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Atzz));
const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);
const vec3<GF3D<const CCTK_REAL, 0, 0, 0>, UP> gf_Gamt_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamtx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamty),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, Gamtz));
const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGx_(cctkGH, betaGx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGy_(cctkGH, betaGy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGz_(cctkGH, betaGz);
const vec3<GF3D<const CCTK_REAL, 0, 0, 0>, UP> gf_betaG_(GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, betaGx),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, betaGy),GF3D<const CCTK_REAL, 0, 0, 0>(cctkGH, betaGz));//const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dchi_(cctkGH, allocate());const mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN> gf_ddchi_(cctkGH, allocate());calc_derivs2(cctkGH, gf_chi_, gf_dchi_, gf_ddchi_);const mat3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, DN, DN> gf_dgammat_(cctkGH, allocate());const mat3<mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN>, DN, DN> gf_ddgammat_(cctkGH, allocate());calc_derivs2(cctkGH, gf_gammat_, gf_dgammat_, gf_ddgammat_);const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dKh_(cctkGH, allocate());calc_derivs(cctkGH, gf_Kh_, gf_dKh_);const mat3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, DN, DN> gf_dAt_(cctkGH,allocate());calc_derivs(cctkGH, gf_At_, gf_dAt_);const vec3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, UP> gf_dGamt_(cctkGH,allocate());calc_derivs(cctkGH, gf_Gamt_, gf_dGamt_);const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dTheta_(cctkGH, allocate());calc_derivs(cctkGH, gf_Theta_, gf_dTheta_);const vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN> gf_dalphaG_(cctkGH, allocate());const mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN> gf_ddalphaG_(cctkGH, allocate());calc_derivs2(cctkGH, gf_alphaG_, gf_dalphaG_, gf_ddalphaG_);const vec3<vec3<GF3D<CCTK_REAL, 0, 0, 0>, DN>, UP> gf_dbetaG_(cctkGH,allocate());const vec3<mat3<GF3D<CCTK_REAL, 0, 0, 0>, DN, DN>, UP> gf_ddbetaG_(cctkGH, allocate());calc_derivs2(cctkGH, gf_betaG_, gf_dbetaG_, gf_ddbetaG_);//
loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Load and calculateconst z4c_vars<CCTK_REAL> vars(kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,gf_Atyy_, gf_Atyz_, gf_Atzz_, gf_Gamtx_, gf_Gamty_, gf_Gamtz_,gf_Theta_, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_, //p.I, dx);
//
// Storegf_chi_rhs_(p.I) = vars.chi_rhs;vars.gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_, gf_gammatxz_rhs_,gf_gammatyy_rhs_, gf_gammatyz_rhs_, gf_gammatzz_rhs_,p);gf_Kh_rhs_(p.I) = vars.Kh_rhs;vars.At_rhs.store(gf_Atxx_rhs_, gf_Atxy_rhs_, gf_Atxz_rhs_, gf_Atyy_rhs_,gf_Atyz_rhs_, gf_Atzz_rhs_, p);vars.Gamt_rhs.store(gf_Gamtx_rhs_, gf_Gamty_rhs_, gf_Gamtz_rhs_, p);gf_Theta_rhs_(p.I) = vars.Theta_rhs;gf_alphaG_rhs_(p.I) = vars.alphaG_rhs;vars.betaG_rhs.store(gf_betaGx_rhs_, gf_betaGy_rhs_, gf_betaGz_rhs_, p);});
loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {// Load and calculateconst CCTK_REAL rho{0};const vec3<CCTK_REAL, DN> Si{0, 0, 0};const mat3<CCTK_REAL, DN, DN> Sij{0, 0, 0, 0, 0, 0};const z4c_vars<CCTK_REAL> vars(kappa1, kappa2, f_mu_L, f_mu_S, eta, //gf_chi_(p.I), gf_dchi_(p.I), gf_ddchi_(p.I), //gf_gammat_(p.I), gf_dgammat_(p.I), gf_ddgammat_(p.I), //gf_Kh_(p.I), gf_dKh_(p.I), //gf_At_(p.I), gf_dAt_(p.I), //gf_Gamt_(p.I), gf_dGamt_(p.I), //gf_Theta_(p.I), gf_dTheta_(p.I), //gf_alphaG_(p.I), gf_dalphaG_(p.I), gf_ddalphaG_(p.I), //gf_betaG_(p.I), gf_dbetaG_(p.I), gf_ddbetaG_(p.I), //rho, //Si, //Sij);// Storegf_chi_rhs_(p.I) = vars.chi_rhs;vars.gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_,gf_gammatxz_rhs_, gf_gammatyy_rhs_,gf_gammatyz_rhs_, gf_gammatzz_rhs_, p.I);gf_Kh_rhs_(p.I) = vars.Kh_rhs;vars.At_rhs.store(gf_Atxx_rhs_, gf_Atxy_rhs_, gf_Atxz_rhs_,gf_Atyy_rhs_, gf_Atyz_rhs_, gf_Atzz_rhs_, p.I);vars.Gamt_rhs.store(gf_Gamtx_rhs_, gf_Gamty_rhs_, gf_Gamtz_rhs_, p.I);gf_Theta_rhs_(p.I) = vars.Theta_rhs;gf_alphaG_rhs_(p.I) = vars.alphaG_rhs;vars.betaG_rhs.store(gf_betaGx_rhs_, gf_betaGy_rhs_, gf_betaGz_rhs_,p.I);});
apply_upwind(cctkGH, gf_gammatxx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatxx_rhs_);apply_upwind(cctkGH, gf_gammatxy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatxy_rhs_);apply_upwind(cctkGH, gf_gammatxz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatxz_rhs_);apply_upwind(cctkGH, gf_gammatyy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatyy_rhs_);apply_upwind(cctkGH, gf_gammatyz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatyz_rhs_);apply_upwind(cctkGH, gf_gammatzz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_gammatzz_rhs_);
apply_upwind(cctkGH, gf_gammat_(0, 0), gf_betaG_, gf_gammatxx_rhs_);apply_upwind(cctkGH, gf_gammat_(0, 1), gf_betaG_, gf_gammatxy_rhs_);apply_upwind(cctkGH, gf_gammat_(0, 2), gf_betaG_, gf_gammatxz_rhs_);apply_upwind(cctkGH, gf_gammat_(1, 1), gf_betaG_, gf_gammatyy_rhs_);apply_upwind(cctkGH, gf_gammat_(1, 2), gf_betaG_, gf_gammatyz_rhs_);apply_upwind(cctkGH, gf_gammat_(2, 2), gf_betaG_, gf_gammatzz_rhs_);
apply_upwind(cctkGH, gf_Atxx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atxx_rhs_);apply_upwind(cctkGH, gf_Atxy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atxy_rhs_);apply_upwind(cctkGH, gf_Atxz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atxz_rhs_);apply_upwind(cctkGH, gf_Atyy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atyy_rhs_);apply_upwind(cctkGH, gf_Atyz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atyz_rhs_);apply_upwind(cctkGH, gf_Atzz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Atzz_rhs_);
apply_upwind(cctkGH, gf_At_(0, 0), gf_betaG_, gf_Atxx_rhs_);apply_upwind(cctkGH, gf_At_(0, 1), gf_betaG_, gf_Atxy_rhs_);apply_upwind(cctkGH, gf_At_(0, 2), gf_betaG_, gf_Atxz_rhs_);apply_upwind(cctkGH, gf_At_(1, 1), gf_betaG_, gf_Atyy_rhs_);apply_upwind(cctkGH, gf_At_(1, 2), gf_betaG_, gf_Atyz_rhs_);apply_upwind(cctkGH, gf_At_(2, 2), gf_betaG_, gf_Atzz_rhs_);
apply_upwind(cctkGH, gf_Gamtx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Gamtx_rhs_);apply_upwind(cctkGH, gf_Gamty_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Gamty_rhs_);apply_upwind(cctkGH, gf_Gamtz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_Gamtz_rhs_);
apply_upwind(cctkGH, gf_Gamt_(0), gf_betaG_, gf_Gamtx_rhs_);apply_upwind(cctkGH, gf_Gamt_(1), gf_betaG_, gf_Gamty_rhs_);apply_upwind(cctkGH, gf_Gamt_(2), gf_betaG_, gf_Gamtz_rhs_);
apply_upwind(cctkGH, gf_betaGx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_betaGx_rhs_);apply_upwind(cctkGH, gf_betaGy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_betaGy_rhs_);apply_upwind(cctkGH, gf_betaGz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,gf_betaGz_rhs_);
apply_upwind(cctkGH, gf_betaG_(0), gf_betaG_, gf_betaGx_rhs_);apply_upwind(cctkGH, gf_betaG_(1), gf_betaG_, gf_betaGy_rhs_);apply_upwind(cctkGH, gf_betaG_(2), gf_betaG_, gf_betaGz_rhs_);
apply_diss(cctkGH, gf_gammatxx_, gf_gammatxx_rhs_);apply_diss(cctkGH, gf_gammatxy_, gf_gammatxy_rhs_);apply_diss(cctkGH, gf_gammatxz_, gf_gammatxz_rhs_);apply_diss(cctkGH, gf_gammatyy_, gf_gammatyy_rhs_);apply_diss(cctkGH, gf_gammatyz_, gf_gammatyz_rhs_);apply_diss(cctkGH, gf_gammatzz_, gf_gammatzz_rhs_);
apply_diss(cctkGH, gf_gammat_(0, 0), gf_gammatxx_rhs_);apply_diss(cctkGH, gf_gammat_(0, 1), gf_gammatxy_rhs_);apply_diss(cctkGH, gf_gammat_(0, 2), gf_gammatxz_rhs_);apply_diss(cctkGH, gf_gammat_(1, 1), gf_gammatyy_rhs_);apply_diss(cctkGH, gf_gammat_(1, 2), gf_gammatyz_rhs_);apply_diss(cctkGH, gf_gammat_(2, 2), gf_gammatzz_rhs_);
apply_diss(cctkGH, gf_Atxx_, gf_Atxx_rhs_);apply_diss(cctkGH, gf_Atxy_, gf_Atxy_rhs_);apply_diss(cctkGH, gf_Atxz_, gf_Atxz_rhs_);apply_diss(cctkGH, gf_Atyy_, gf_Atyy_rhs_);apply_diss(cctkGH, gf_Atyz_, gf_Atyz_rhs_);apply_diss(cctkGH, gf_Atzz_, gf_Atzz_rhs_);
apply_diss(cctkGH, gf_At_(0, 0), gf_Atxx_rhs_);apply_diss(cctkGH, gf_At_(0, 1), gf_Atxy_rhs_);apply_diss(cctkGH, gf_At_(0, 2), gf_Atxz_rhs_);apply_diss(cctkGH, gf_At_(1, 1), gf_Atyy_rhs_);apply_diss(cctkGH, gf_At_(1, 2), gf_Atyz_rhs_);apply_diss(cctkGH, gf_At_(2, 2), gf_Atzz_rhs_);
apply_diss(cctkGH, gf_Gamtx_, gf_Gamtx_rhs_);apply_diss(cctkGH, gf_Gamty_, gf_Gamty_rhs_);apply_diss(cctkGH, gf_Gamtz_, gf_Gamtz_rhs_);
apply_diss(cctkGH, gf_Gamt_(0), gf_Gamtx_rhs_);apply_diss(cctkGH, gf_Gamt_(1), gf_Gamty_rhs_);apply_diss(cctkGH, gf_Gamt_(2), gf_Gamtz_rhs_);
apply_diss(cctkGH, gf_betaGx_, gf_betaGx_rhs_);apply_diss(cctkGH, gf_betaGy_, gf_betaGy_rhs_);apply_diss(cctkGH, gf_betaGz_, gf_betaGz_rhs_);}
apply_diss(cctkGH, gf_betaG_(0), gf_betaGx_rhs_);apply_diss(cctkGH, gf_betaG_(1), gf_betaGy_rhs_);apply_diss(cctkGH, gf_betaG_(2), gf_betaGz_rhs_);} // namespace Z4c
template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T pow4(const T x) {return pow2(pow2(x));
template <typename T> CCTK_ATTRIBUTE_ALWAYS_INLINE constexpr T pow4(const T x) {const T x2 = x * x;return x2 * x2;}template <typename T> CCTK_ATTRIBUTE_ALWAYS_INLINE constexpr T pow5(const T x) {const T x2 = x * x;const T x4 = x2 * x2;return x4 * x;}template <typename T> CCTK_ATTRIBUTE_ALWAYS_INLINE constexpr T pow6(const T x) {const T x2 = x * x;const T x4 = x2 * x2;return x4 * x2;}namespace detail {template <typename T> constexpr T pown(const T x, int n) {T r{1};T y{x};while (n) {if (n & 1)r *= y;y *= y;n >>= 1;}return r;}} // namespace detailtemplate <typename T> constexpr T pown(const T x, const int n) {return n >= 0 ? detail::pown(x, n) : 1 / detail::pown(x, -n);}constexpr int factorial(int n) {int r{1};while (n > 1) {r *= n;--n;}return r;
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3(initializer_list<T> v)
private:static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vector<T>make_vector(T vx, T vy, T vz) {vector<T> vec;vec.reserve(3);vec.push_back(move(vx));vec.push_back(move(vy));vec.push_back(move(vz));return vec;}public:explicit constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vec3(T vx, T vy, T vz): elts(make_vector(move(vx), move(vy), move(vz))) {}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vec3(initializer_list<T> v)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3(const GF3D<const T, 0, 0, 0> &gf_vx_,const GF3D<const T, 0, 0, 0> &gf_vy_,const GF3D<const T, 0, 0, 0> &gf_vz_,const vect<int, 3> &I)
CCTK_ATTRIBUTE_ALWAYS_INLINEvec3(const GF3D<add_const_t<T>, 0, 0, 0> &gf_vx_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_vy_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_vz_, const vect<int, 3> &I)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/vec3(const GF3D<T, 0, 0, 0> &gf_vx_, const GF3D<T, 0, 0, 0> &gf_vy_,const GF3D<T, 0, 0, 0> &gf_vz_, const vect<int, 3> &I)
CCTK_ATTRIBUTE_ALWAYS_INLINEvec3(const GF3D<remove_const_t<T>, 0, 0, 0> &gf_vx_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_vy_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_vz_, const vect<int, 3> &I)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ void store(const GF3D<T, 0, 0, 0> &gf_vx_,const GF3D<T, 0, 0, 0> &gf_vy_,const GF3D<T, 0, 0, 0> &gf_vz_,const PointDesc &p) const {
vec3(const cGH *const cctkGH, allocate): vec3(T(cctkGH, allocate()), T(cctkGH, allocate()),T(cctkGH, allocate())) {}CCTK_ATTRIBUTE_ALWAYS_INLINE void store(const GF3D<T, 0, 0, 0> &gf_vx_,const GF3D<T, 0, 0, 0> &gf_vy_,const GF3D<T, 0, 0, 0> &gf_vz_,const vect<int, 3> &I) const {
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T, dnup>
template <typename U = T>CCTK_ATTRIBUTE_ALWAYS_INLINEvec3<remove_cv_t<remove_reference_t<result_of_t<U(vect<int, 3>)> > >,dnup>operator()(const vect<int, 3> &I) const {return {elts[0](I), elts[1](I), elts[2](I)};}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vec3<T, dnup>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3(initializer_list<T> A)
private:static constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE vector<T>make_vector(T Axx, T Axy, T Axz, T Ayy, T Ayz, T Azz) {vector<T> vec;vec.reserve(6);vec.push_back(move(Axx));vec.push_back(move(Axy));vec.push_back(move(Axz));vec.push_back(move(Ayy));vec.push_back(move(Ayz));vec.push_back(move(Azz));return vec;}public:explicit constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE mat3(T Axx, T Axy, T Axz,T Ayy, T Ayz, T Azz): elts(make_vector(move(Axx), move(Axy), move(Axz), move(Ayy), move(Ayz),move(Azz))) {}constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE mat3(initializer_list<T> A)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ mat3(const GF3D<const T, 0, 0, 0> &gf_Axx_,const GF3D<const T, 0, 0, 0> &gf_Axy_,const GF3D<const T, 0, 0, 0> &gf_Axz_,const GF3D<const T, 0, 0, 0> &gf_Ayy_,const GF3D<const T, 0, 0, 0> &gf_Ayz_,const GF3D<const T, 0, 0, 0> &gf_Azz_,const vect<int, 3> &I)
CCTK_ATTRIBUTE_ALWAYS_INLINEmat3(const GF3D<add_const_t<T>, 0, 0, 0> &gf_Axx_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_Axy_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_Axz_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_Ayy_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_Ayz_,const GF3D<add_const_t<T>, 0, 0, 0> &gf_Azz_, const vect<int, 3> &I)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/mat3(const GF3D<T, 0, 0, 0> &gf_Axx_, const GF3D<T, 0, 0, 0> &gf_Axy_,const GF3D<T, 0, 0, 0> &gf_Axz_, const GF3D<T, 0, 0, 0> &gf_Ayy_,const GF3D<T, 0, 0, 0> &gf_Ayz_, const GF3D<T, 0, 0, 0> &gf_Azz_,const vect<int, 3> &I)
CCTK_ATTRIBUTE_ALWAYS_INLINEmat3(const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Axx_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Axy_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Axz_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Ayy_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Ayz_,const GF3D<remove_const_t<T>, 0, 0, 0> &gf_Azz_, const vect<int, 3> &I)
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ void
mat3(const cGH *const cctkGH, allocate): mat3(T(cctkGH, allocate()), T(cctkGH, allocate()),T(cctkGH, allocate()), T(cctkGH, allocate()),T(cctkGH, allocate()), T(cctkGH, allocate())) {}CCTK_ATTRIBUTE_ALWAYS_INLINE void
gf_Axx_(p.I) = A(0, 0);gf_Axy_(p.I) = A(0, 1);gf_Axz_(p.I) = A(0, 2);gf_Ayy_(p.I) = A(1, 1);gf_Ayz_(p.I) = A(1, 2);gf_Azz_(p.I) = A(2, 2);
gf_Axx_(I) = A(0, 0);gf_Axy_(I) = A(0, 1);gf_Axz_(I) = A(0, 2);gf_Ayy_(I) = A(1, 1);gf_Ayz_(I) = A(1, 2);gf_Azz_(I) = A(2, 2);
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T, dnup1, dnup2>
template <typename U = T>CCTK_ATTRIBUTE_ALWAYS_INLINEmat3<remove_cv_t<remove_reference_t<result_of_t<U(vect<int, 3>)> > >,dnup1, dnup2>operator()(const vect<int, 3> &I) const {return {elts[0](I), elts[1](I), elts[2](I),elts[3](I), elts[4](I), elts[5](I)};}friend constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE mat3<T, dnup1, dnup2>
return mat3<T, dnup1, dnup2>([&](int a, int b) {return sum1([&](int x) { return A(a, x) * B(x, b); });
return mat3<T, dnup1, dnup2>([&](int a, int b) CCTK_ATTRIBUTE_ALWAYS_INLINE {return sum1([&](int x)CCTK_ATTRIBUTE_ALWAYS_INLINE { return A(a, x) * B(x, b); });
template <typename F>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int)>sum1(const F &f) {result_of_t<F(int)> s{0};
template <typename F,typename R = remove_cv_t<remove_reference_t<result_of_t<F(int)> > > >constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R sum1(const F &f) {R s{0};
template <typename F>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int, int)>sum2(const F &f) {result_of_t<F(int, int)> s{0};
template <typename F, typename R = remove_cv_t<remove_reference_t<result_of_t<F(int, int)> > > >constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R sum2(const F &f) {R s{0};
template <typename F>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int, int, int)>sum3(const F &f) {result_of_t<F(int, int, int)> s{0};
template <typename F, typename R = remove_cv_t<remove_reference_t<result_of_t<F(int, int, int)> > > >constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE R sum3(const F &f) {R s{0};
}////////////////////////////////////////////////////////////////////////////////template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T deriv(const GF3D<const T, 0, 0, 0> &gf_,const vect<int, dim> &I,const vec3<T, UP> &dx, const int dir) {constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};return (gf_(I + DI[dir]) - gf_(I - DI[dir])) / (2 * dx(dir));}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ Tderiv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const bool sign, const vec3<T, UP> &dx, const int dir) {constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};// arXiv:1111.2177 [gr-qc], (71)if (sign)// + [ 0 -1 +1 0 0]// + 1/2 [+1 -2 +1 0 0]// [+1/2 -2 +3/2 0 0]return (+1 / T(2) * gf_(I - 2 * DI[dir]) //- 2 * gf_(I - DI[dir]) //+ 3 / T(2) * gf_(I)) /dx(dir);else// + [ 0 0 -1 +1 0 ]// - 1/2 [ 0 0 +1 -2 +1 ]// [ 0 0 -3/2 +2 -1/2]return (-3 / T(2) * gf_(I) //+ 2 * gf_(I + DI[dir]) //- 1 / T(2) * gf_(I + 2 * DI[dir])) /dx(dir);}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ Tderiv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx, const int dir1, const int dir2) {constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};if (dir1 == dir2)return ((gf_(I - DI[dir1]) + gf_(I + DI[dir1])) - 2 * gf_(I)) /pow2(dx(dir1));else// return (deriv(gf_, I + DI[dir2], dx, dir1) -// deriv(gf_, I - DI[dir2], dx, dir1)) /// (2 * dx(dir2));return ((gf_(I + DI[dir1] + DI[dir2]) //- gf_(I - DI[dir1] + DI[dir2])) //- (gf_(I + DI[dir1] - DI[dir2]) //- gf_(I - DI[dir1] - DI[dir2]))) /// (4 * dx(dir1) * dx(dir2));
template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ Tderiv4(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx, const int dir) {constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};return ((gf_(I - 2 * DI[dir]) + gf_(I + 2 * DI[dir])) //- 4 * (gf_(I - DI[dir]) + gf_(I + DI[dir])) //+ 6 * gf_(I)) /pow4(dx(dir));}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3<T, DN>deriv(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {return {deriv(gf_, I, dx, 0),deriv(gf_, I, dx, 1),deriv(gf_, I, dx, 2),};}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3<T, DN>deriv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dir, const vec3<T, UP> &dx) {return {deriv_upwind(gf_, I, signbit(dir(0)), dx, 0),deriv_upwind(gf_, I, signbit(dir(1)), dx, 1),deriv_upwind(gf_, I, signbit(dir(2)), dx, 2),};}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ mat3<T, DN, DN>deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T, UP> &dx) {return {deriv2(gf_, I, dx, 0, 0), deriv2(gf_, I, dx, 0, 1),deriv2(gf_, I, dx, 0, 2), deriv2(gf_, I, dx, 1, 1),deriv2(gf_, I, dx, 1, 2), deriv2(gf_, I, dx, 2, 2),};}template <typename T>/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T diss(const GF3D<const T, 0, 0, 0> &gf_,const vect<int, dim> &I,const vec3<T, UP> &dx) {// arXiv:gr-qc/0610128, (63), with r=2return -1 / T(16) *(pow3(dx(0)) * deriv4(gf_, I, dx, 0) //+ pow3(dx(1)) * deriv4(gf_, I, dx, 1) //+ pow3(dx(2)) * deriv4(gf_, I, dx, 2));}
}// Test derivativesstatic_assert(deriv_order % 2 == 0, "");constexpr int required_ghosts = deriv_order / 2 + 1;const double eps = 1.0e-12;// derivfor (int order = 0; order <= deriv_order; ++order) {// CCTK_VINFO("Testing deriv order=%d", order);array<double, 2 * required_ghosts + 7> arr;for (size_t i = 0; i < arr.size(); ++i)arr[i] = NAN;double *const var = &arr[arr.size() / 2];for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)var[i] = pown(i, order);const double expected = order == 1 ? 1 : 0;const double found = deriv1d(var, 1, 1.0);assert(fabs(found - expected) <= eps);
// deriv (upwind)for (int order = 0; order <= deriv_order; ++order) {for (int sign = 0; sign <= 1; ++sign) {// CCTK_VINFO("Testing deriv (upwind) order=%d sign=%d", order, sign);array<double, 2 * required_ghosts + 7> arr;for (size_t i = 0; i < arr.size(); ++i)arr[i] = NAN;double *const var = &arr[arr.size() / 2];if (sign)for (int i = -deriv_order / 2 - 1; i <= deriv_order / 2 - 1; ++i)var[i] = pown(i, order);elsefor (int i = -deriv_order / 2 + 1; i <= deriv_order / 2 + 1; ++i)var[i] = pown(i, order);const double expected = order == 1 ? 1 : 0;const double found = deriv1d_upwind(var, 1, sign, 1.0);assert(fabs(found - expected) <= eps);}}// deriv2for (int order = 0; order <= deriv_order; ++order) {// CCTK_VINFO("Testing deriv2 order=%d", order);array<double, 2 * required_ghosts + 7> arr;for (size_t i = 0; i < arr.size(); ++i)arr[i] = NAN;double *const var = &arr[arr.size() / 2];for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)var[i] = pown(i, order);const double expected = order == 2 ? 2 : 0;const double found = deriv2_1d(var, 1, 1.0);assert(fabs(found - expected) <= eps);}// deriv2 (mixed)for (int orderj = 0; orderj <= deriv_order; ++orderj) {for (int orderi = 0; orderi <= deriv_order; ++orderi) {// CCTK_VINFO("Testing deriv2 (mixed) order=%d,%d", orderi, orderj);array<array<double, 2 * required_ghosts + 7>, 2 * required_ghosts + 7>arr;for (size_t j = 0; j < arr.size(); ++j)for (size_t i = 0; i < arr[j].size(); ++i)arr[j][i] = NAN;const int di = 1;const int dj = arr[0].size();double *const var = &arr[arr.size() / 2][arr[0].size() / 2];for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j)for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)var[j * dj + i * di] = pown(i, orderi) * pown(j, orderj);const double expected = orderi == 1 && orderj == 1 ? 1 : 0;const double found = deriv2_2d(var, di, dj, 1.0, 1.0);assert(fabs(found - expected) <= eps);}}// deriv (dissipation)for (int order = 0; order <= deriv_order + 2; ++order) {// CCTK_VINFO("Testing deriv (dissipation) order=%d", order);array<double, 2 * required_ghosts + 7> arr;for (size_t i = 0; i < arr.size(); ++i)arr[i] = NAN;double *const var = &arr[arr.size() / 2];for (int i = -deriv_order / 2 - 1; i <= deriv_order / 2 + 1; ++i)var[i] = pown(i, order);const double expected = order == deriv_order + 2 ? factorial(order) : 0;const double found = deriv1d_diss(var, 1, 1.0);assert(fabs(found - expected) <= eps);}
cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DSIMD=AVX2 -DSIMD_OPTIONALS=FMA -DCMAKE_INSTALL_PREFIX=/cactus/nsimd .. && \
cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DSIMD=AVX2 -DSIMD_OPTIONALS=FMA -DCMAKE_INSTALL_PREFIX=/usr/local .. && \
true ) && \
true) && \rm -rf src# Install SiloRUN mkdir src && \(cd src && \wget https://wci.llnl.gov/content/assets/docs/simulation/computer-codes/silo/silo-4.10.2/silo-4.10.2-bsd.tar.gz && \tar xzf silo-4.10.2-bsd.tar.gz && \cd silo-4.10.2-bsd && \mkdir build && \cd build && \../configure --disable-fortran --enable-optimization --with-hdf5=/usr/lib/x86_64-linux-gnu/hdf5/serial/include,/usr/lib/x86_64-linux-gnu/hdf5/serial/lib --prefix=/usr/local && \make -j$(nproc) && \make -j$(nproc) install && \true) && \
( cd src && \wget https://github.com/AMReX-Codes/amrex/archive/20.01.tar.gz && \tar xzf 20.01.tar.gz && \cd amrex-20.01 && \
(cd src && \wget https://github.com/AMReX-Codes/amrex/archive/20.02.tar.gz && \tar xzf 20.02.tar.gz && \cd amrex-20.02 && \
cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_PARTICLES=ON -DENABLE_ASSERTIONS=ON -DENABLE_BACKTRACE=ON -DENABLE_OMP=ON -DCMAKE_INSTALL_PREFIX=/cactus/amrex .. && \
cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_PARTICLES=ON -DENABLE_ASSERTIONS=ON -DENABLE_BACKTRACE=ON -DENABLE_OMP=ON -DCMAKE_INSTALL_PREFIX=/usr-local .. && \
# How to build this Docker image:# docker build . --tag eschnett/cactus-amrex# docker push eschnett/cactus-amrex# Use UbuntuFROM ubuntu:19.10RUN mkdir /cactusWORKDIR /cactus# Install system packagesRUN apt-get update && \apt-get --yes --no-install-recommends install \build-essential \ca-certificates \clang-7 \flang-7 \git \libgsl-dev \libhwloc-dev \libmpich-dev \perl \pkg-config \python \python3 \rsync \subversion \wget \zlib1g-dev \&& \rm -rf /var/lib/apt/lists/*# Install cmake# (AMReX 20.01 requires at least cmake 3.14)RUN mkdir dist && \( cd dist && \wget https://github.com/Kitware/CMake/releases/download/v3.16.2/cmake-3.16.2-Linux-x86_64.tar.gz && \tar xzf cmake-3.16.2-Linux-x86_64.tar.gz -C /usr/local --strip-components=1 && \true ) && \rm -rf dist# Install NSIMD# Note: This assumes that the system has x86_64 CPUs with AVX2 and FMARUN mkdir src && \( cd src && \git clone -b eschnett/storeu_masked https://github.com/eschnett/nsimd.git && \cd nsimd && \python3 egg/hatch.py --all --force --disable-clang-format && \mkdir build && \cd build && \cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_COMPILER=clang-7 -DCMAKE_CXX_COMPILER=clang++-7 -DSIMD=AVX2 -DSIMD_OPTIONALS=FMA -DCMAKE_INSTALL_PREFIX=/cactus/nsimd .. && \make -j$(nproc) && \make -j$(nproc) install && \true ) && \rm -rf src# Install AMReXRUN mkdir src && \( cd src && \wget https://github.com/AMReX-Codes/amrex/archive/20.01.tar.gz && \tar xzf 20.01.tar.gz && \cd amrex-20.01 && \mkdir build && \cd build && \cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_C_COMPILER=clang-7 -DCMAKE_CXX_COMPILER=clang++-7 -DCMAKE_Fortran_COMPILER=flang-7 -DENABLE_PARTICLES=ON -DENABLE_ASSERTIONS=ON -DENABLE_BACKTRACE=ON -DENABLE_OMP=OFF -DCMAKE_INSTALL_PREFIX=/cactus/amrex .. && \make -j$(nproc) && \make -j$(nproc) install && \true ) && \rm -rf src# As documentationCOPY Dockerfile /Dockerfile
# How to build this Docker image:# docker build . --tag eschnett/cactus-amrex# docker push eschnett/cactus-amrex# Use NixOSFROM nixos/nix# RUN nix-channel --add https://nixos.org/channels/nixpkgs-unstable nixpkgs# RUN nix-channel --update# RUN nix-build -A pythonFull '<nixpkgs>'RUN nix-env -i \cmake \gfortran-wrapper \git \gnumake \gsl \hwloc \mpich \perl \pkg-config \python \python3 \rsync \subversion \wget \zlibRUN mkdir /cactusWORKDIR /cactus# # Install system packages# RUN apt-get update && \# apt-get --yes --no-install-recommends install \# build-essential \# ca-certificates \# g++ \# gfortran \# git \# libgsl-dev \# libhwloc-dev \# libmpich-dev \# perl \# pkg-config \# python \# python3 \# rsync \# subversion \# wget \# zlib1g-dev \# && \# rm -rf /var/lib/apt/lists/*# Install AMReXRUN mkdir srcWORKDIR /cactus/srcRUN wget https://github.com/AMReX-Codes/amrex/archive/20.01.tar.gzRUN tar xzf 20.01.tar.gzWORKDIR /cactus/src/amrex-20.01RUN mkdir buildWORKDIR /cactus/src/amrex-20.01/buildRUN cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_PARTICLES=ON -DENABLE_ASSERTIONS=ON -DENABLE_BACKTRACE=ON -DENABLE_OMP=ON -DCMAKE_INSTALL_PREFIX=/cactus/amrex ..RUN make -j$(nproc)RUN make -j$(nproc) installWORKDIR /cactus# RUN rm -rf src# As documentationCOPY Dockerfile /Dockerfile
# Configuration for an Ubuntu installation, assuming the following# list of packages is installed:## build-essential perl python gfortran g++ libmpich-dev## In addition, installing the following list of packaed will prevent# Cactus from compiling its own versions of these libraries:## libfftw3-dev libgsl-dev libatlas-base-dev libjpeg-dev libssl-dev# libhdf5-dev hdf5-tools libnuma-dev libltdl-dev libhwloc-dev zlib1g-dev## Tools like GetComponents and Simfactory like to have the following list# installed too## python subversion git
# Configuration for an Ubuntu installation