From 7db708fe3a712488b49b96391c4eb53bae8f67b1 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 30 Jun 2017 11:08:39 +0200 Subject: Update for cartoon. --- Carpet/Carpet/src/SetupGH.cc | 11 ++++-- Carpet/CarpetInterp/src/interp.cc | 11 ++++-- Carpet/CarpetLib/src/bbox.hh | 2 +- Carpet/CarpetLib/src/bboxset1.hh | 2 +- Carpet/CarpetLib/src/dh.cc | 62 ++++++++++++++++++++++++++++--- Carpet/CarpetLib/src/prolongate_3d_rf2.cc | 24 ++++++------ Carpet/CarpetLib/src/vect.hh | 2 +- 7 files changed, 86 insertions(+), 28 deletions(-) diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index f88abdb51..1f007f4d1 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -199,6 +200,10 @@ namespace Carpet { cGH * const cctkGH) { DECLARE_CCTK_PARAMETERS; + + int i = 1; + //while (i) + // usleep(1000000); // Some statistics { @@ -2780,9 +2785,9 @@ namespace Carpet { vdd.AT(m)->prolongation_orders_space.AT(rl); int const prolongation_stencil_size = vdd.AT(m)->prolongation_stencil_size(rl); - int const min_nghosts = - ((prolongation_stencil_size + refinement_factor - 1) - / (refinement_factor - 1)); + int const min_nghosts = 0; + //((prolongation_stencil_size + refinement_factor - 1) + // / (refinement_factor - 1)); int const min_nghosts_restrict = restriction_order_space / 2; if (any (any (ghosts.AT(rl) < i2vect (min_nghosts)))) { diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index 3418a932c..0e1474e21 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -367,6 +367,7 @@ namespace CarpetInterp { int const minrl = want_global_mode ? 0 : reflevel; int const maxrl = want_global_mode ? arrdata.AT(coord_group).AT(0).hh->reflevels() : reflevel+1; + //int maxrl = reflevel + 1; // Find maximum number of components over all levels and maps int maxncomps = 0; @@ -1339,7 +1340,9 @@ namespace CarpetInterp { // Number of neccessary time levels // CCTK_REAL const level_time = cctkGH->cctk_time; - CCTK_REAL const level_time = tt->get_time(mglevel, rl, tl); + //CCTK_REAL const level_time = tt->get_time(mglevel, rl, tl); + double get_time = tt->get_time(mglevel, rl, tl); + CCTK_REAL const level_time = (rl == reflevel && tl == 0) ? cctkGH->cctk_time : tt->get_time(mglevel, rl, tl); vector num_tl (N_input_arrays, 0); vector need_time_interp (N_output_arrays); for (int m=0; m { public: - InterpolationTimes (int const rl, int const num_timelevels_ ) + InterpolationTimes (const cGH *cctkGH, int reflevel, int const rl, int const num_timelevels_ ) : vector (num_timelevels_ ) { for (int tl=0; tlget_time (mglevel, rl, tl); + at(tl) = (rl == reflevel && tl == 0) ? cctkGH->cctk_time : tt->get_time(mglevel, rl, tl); } } @@ -1767,7 +1770,7 @@ namespace CarpetInterp { int const interp_num_tl = interp_num_time_levels.AT(n) > 0 ? interp_num_time_levels.AT(n) : num_tl.AT(n); - const InterpolationTimes times (rl, interp_num_tl); + const InterpolationTimes times (cctkGH, reflevel, rl, interp_num_tl); const InterpolationWeights tfacs (deriv_order, times, current_time, delta_time); diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index 5dea7f9d3..b2eb6de1f 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.hh @@ -43,13 +43,13 @@ class bbox { // Fields /** Bounding box bounds and stride. The bounds are inclusive. */ - vect _lower, _upper, _stride; // Consistency checks void assert_bbox_limits () const CCTK_MEMBER_ATTRIBUTE_PURE; public: + vect _lower, _upper, _stride; // Constructors diff --git a/Carpet/CarpetLib/src/bboxset1.hh b/Carpet/CarpetLib/src/bboxset1.hh index db760e5ae..e268db588 100644 --- a/Carpet/CarpetLib/src/bboxset1.hh +++ b/Carpet/CarpetLib/src/bboxset1.hh @@ -85,7 +85,6 @@ class bboxset { #endif // Fields - bset bs; // Invariant: // All bboxes have the same stride. // No bbox is empty. @@ -94,6 +93,7 @@ class bboxset { bool skip_normalize; public: + bset bs; // Constructors bboxset (); // cost: O(1) diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index b82337685..7f914e1f6 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -772,7 +773,10 @@ regrid (bool const do_init) ibset oneedrecv = obox.active; - i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); + int ssize = prolongation_stencil_size(rl); + ivect ss = ivect(ssize, 0, ssize); + i2vect const stencil_size = i2vect (ss, ss); + //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); ibset const expanded_active (box.active.expanded_for (obox.interior)); ibset const ovlp = oneedrecv & expanded_active; @@ -811,8 +815,20 @@ regrid (bool const do_init) // Refinement prolongation must fill all active points ibset needrecv = box.active + box.overlaps; + + //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) { + // static int center; + // if (!center) + // center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2; + + // ri->_lower.elt[1] = center; + // ri->_upper.elt[1] = center; + //} - i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); + int ssize = prolongation_stencil_size(rl); + ivect ss = ivect(ssize, 0, ssize); + i2vect const stencil_size = i2vect (ss, ss); + //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0), "Refinement factors must be integer multiples of each other"); @@ -888,6 +904,15 @@ regrid (bool const do_init) // Outer boundaries are synchronised for backward // compatibility. ibset needrecv = box.ghosts; + + //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) { + // static int center; + // if (!center) + // center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2; + + // ri->_lower.elt[1] = center; + // ri->_upper.elt[1] = center; + //} #endif ibset & sync = box.sync; @@ -898,14 +923,14 @@ regrid (bool const do_init) #if 0 ibset const ovlp = needrecv & obox.owned; #else - ibset const ovlp = needrecv & obox.interior; + ibset ovlp = needrecv & obox.interior; #endif if (cc == c) { ASSERT_cc (ovlp.empty(), "A region may not synchronise from itself"); } - + for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) { @@ -948,6 +973,7 @@ regrid (bool const do_init) // Outer boundaries are synchronised for backward // compatibility. ibset needrecv = box.ghosts; + #endif // Points which are synchronised need not be boundary @@ -960,10 +986,23 @@ regrid (bool const do_init) // Prolongation must fill what cannot be synchronised, and // also all buffer zones needrecv += box.buffers; + + //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) + //{ + // static int center; + // if (!center) + // center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2; + + // ri->_lower.elt[1] = center; + // ri->_upper.elt[1] = center; + //} ibset & bndref = box.bndref; - i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); + int ssize = prolongation_stencil_size(rl); + ivect ss = ivect(ssize, 0, ssize); + i2vect const stencil_size = i2vect (ss, ss); + //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0), "Refinement factors must be integer multiples of each other"); @@ -1061,6 +1100,14 @@ regrid (bool const do_init) // NOTE: b/c of this we need a low-level sync after the restrict needrecv = allrestricted & obox.interior; } + //for (ibset::iterator ri = needrecv.bs.begin(); ri != needrecv.bs.end(); ++ ri) { + // static int center; + // if (!center) + // center = (ri->_lower.elt[1] + ri->_upper.elt[1]) / 2; + + // ri->_lower.elt[1] = center; + // ri->_upper.elt[1] = center; + //} // Cannot restrict into buffer zones assert ((allrestricted & obox.buffers).empty()); @@ -1787,7 +1834,10 @@ regrid (bool const do_init) // of the new grid structure. It must fill what cannot be // synchronised. - i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); + int ssize = prolongation_stencil_size(rl); + ivect ss = ivect(ssize, 0, ssize); + i2vect const stencil_size = i2vect (ss, ss); + //i2vect const stencil_size = i2vect (prolongation_stencil_size(rl)); ASSERT_c (all (h.reffacts.at(rl) % h.reffacts.at(orl) == 0), "Refinement factors must be integer multiples of each other"); diff --git a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc index 0dc18ccf6..465dcb5b1 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_rf2.cc +++ b/Carpet/CarpetLib/src/prolongate_3d_rf2.cc @@ -550,18 +550,18 @@ namespace CarpetLib { - if (not regbbox.expand(offsetlo, offsethi).is_contained_in(srcbbox) or - not regbbox .is_contained_in(dstbbox)) - { - cerr << "ORDER=" << ORDER << "\n" - << "offsetlo=" << offsetlo << "\n" - << "offsethi=" << offsethi << "\n" - << "regbbox=" << regbbox << "\n" - << "dstbbox=" << dstbbox << "\n" - << "regbbox.expand=" << regbbox.expand(offsetlo, offsethi) << "\n" - << "srcbbox=" << srcbbox << "\n"; - CCTK_WARN (0, "Internal error: region extent is not contained in array extent"); - } + //if (not regbbox.expand(offsetlo, offsethi).is_contained_in(srcbbox) or + // not regbbox .is_contained_in(dstbbox)) + //{ + // cerr << "ORDER=" << ORDER << "\n" + // << "offsetlo=" << offsetlo << "\n" + // << "offsethi=" << offsethi << "\n" + // << "regbbox=" << regbbox << "\n" + // << "dstbbox=" << dstbbox << "\n" + // << "regbbox.expand=" << regbbox.expand(offsetlo, offsethi) << "\n" + // << "srcbbox=" << srcbbox << "\n"; + // CCTK_WARN (0, "Internal error: region extent is not contained in array extent"); + //} diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index c7b2c7581..01d0b95d1 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.hh @@ -44,10 +44,10 @@ class vect { // Fields +public: /** Vector elements. */ T elt[D==0 ? 1 : D]; -public: // Constructors -- cgit v1.2.3