aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetLib
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetLib')
-rw-r--r--Carpet/CarpetLib/src/defs.cc2
-rw-r--r--Carpet/CarpetLib/src/ggf.cc69
-rw-r--r--Carpet/CarpetLib/src/ggf.hh7
-rw-r--r--Carpet/CarpetLib/src/gh.cc4
-rw-r--r--Carpet/CarpetLib/src/th.cc140
-rw-r--r--Carpet/CarpetLib/src/th.hh75
6 files changed, 208 insertions, 89 deletions
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 3bfdcd585..a2b366ea2 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -354,6 +354,7 @@ template size_t memoryof (vector<vector<ibbox> > const & v);
template size_t memoryof (vector<vector<dh::dboxes> > const & v);
template size_t memoryof (vector<vector<dh::fast_dboxes> > const & v);
template size_t memoryof (vector<vector<region_t> > const & v);
+template size_t memoryof (vector<vector<vector<CCTK_REAL> > > const & v);
template size_t memoryof (vector<vector<vector<dh::dboxes> > > const & v);
template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v);
template size_t memoryof (vector<vector<vector<region_t> > > const & v);
@@ -394,6 +395,7 @@ template ostream& output (ostream& os, const vector<CCTK_REAL>& v);
template ostream& output (ostream& os, const vector<ibbox>& v);
template ostream& output (ostream& os, const vector<rbbox>& v);
template ostream& output (ostream& os, const vector<ivect>& v);
+template ostream& output (ostream& os, const vector<rvect>& v);
template ostream& output (ostream& os, const vector<bbvect>& v);
template ostream& output (ostream& os, const vector<i2vect>& v);
template ostream& output (ostream& os, const vector<dh::dboxes> & v);
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 7d947fae5..8a7c643da 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -38,7 +38,7 @@ ggf::ggf (const int varindex_, const operator_type transport_operator_,
vectorlength(vectorlength_), vectorindex(vectorindex_),
vectorleader(vectorleader_)
{
- assert (&t.h == &d.h);
+ // assert (&t.h == &d.h);
assert (vectorlength >= 1);
assert (vectorindex >= 0 and vectorindex < vectorlength);
@@ -176,18 +176,6 @@ void ggf::recompose_fill (comm_state & state, int const rl,
assert (d.fast_boxes.AT(ml).AT(rl).do_init);
- vector <int> tls;
- if (do_prolongate and rl > 0 and
- transport_operator != op_none and transport_operator != op_sync and
- transport_operator != op_restrict)
- {
- int const numtl = timelevels (ml, rl);
- tls.resize (numtl);
- for (int tl = 0; tl < numtl; ++ tl) {
- tls.AT(tl) = tl;
- }
- }
-
// Initialise from the same level of the old hierarchy, where
// possible
if (rl < (int)oldstorage.AT(ml).size()) {
@@ -207,12 +195,18 @@ void ggf::recompose_fill (comm_state & state, int const rl,
if (transport_operator != op_none and transport_operator != op_sync and
transport_operator != op_restrict)
{
+ int const numtl = timelevels (ml, rl);
+ vector <int> tls (numtl);
+ for (int tl = 0; tl < numtl; ++ tl) {
+ tls.AT(tl) = tl;
+ }
+
for (int tl = 0; tl < timelevels (ml, rl); ++tl) {
transfer_from_all (state,
tl, rl, ml,
& dh::fast_dboxes::fast_old2new_ref_prol_sendrecv,
tls, rl - 1, ml,
- t.time (tl, rl, ml));
+ t.get_time (ml, rl, tl));
} // for tl
} // if transport_operator
} // if rl
@@ -277,6 +271,22 @@ void ggf::cycle_all (int const rl, int const ml) {
}
}
+// Uncycle the time levels by rotating the data sets
+void ggf::uncycle_all (int const rl, int const ml) {
+ assert (rl>=0 and rl<h.reflevels());
+ assert (ml>=0 and ml<h.mglevels());
+ int const ntl = timelevels(ml,rl);
+ assert (ntl > 0);
+ for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) {
+ fdata & fdatas = storage.AT(ml).AT(rl).AT(lc);
+ gdata * const tmpdata = fdatas.AT(0);
+ for (int tl=0; tl<ntl-1; ++tl) {
+ fdatas.AT(tl) = fdatas.AT(tl+1);
+ }
+ fdatas.AT(ntl-1) = tmpdata;
+ }
+}
+
// Flip the time levels by exchanging the data sets
void ggf::flip_all (int const rl, int const ml) {
assert (rl>=0 and rl<h.reflevels());
@@ -285,9 +295,9 @@ void ggf::flip_all (int const rl, int const ml) {
assert (ntl > 0);
for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) {
fdata & fdatas = storage.AT(ml).AT(rl).AT(lc);
- for (int tl=0; tl<ntl/2; ++tl) {
- const int tl1 = tl;
- const int tl2 = ntl-1 - tl;
+ for (int tl=1; tl<(ntl+1)/2; ++tl) {
+ const int tl1 = tl;
+ const int tl2 = ntl - tl;
assert (tl1 < tl2);
gdata * const tmpdata = fdatas.AT(tl1);
fdatas.AT(tl1) = fdatas.AT(tl2);
@@ -383,15 +393,15 @@ ref_bnd_prolongate_all (comm_state & state,
void
ggf::
mg_restrict_all (comm_state & state,
- int const tl, int const rl, int const ml,
+ int const tl, int const rl, int const ml,
CCTK_REAL const time)
{
static Timer timer ("mg_restrict_all");
timer.start ();
// Require same times
static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
- assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml-1))
- <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+ assert (abs(t.get_time(ml,rl,0) - t.get_time(ml-1,rl,0))
+ <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
vector<int> const tl2s(1,tl);
transfer_from_all (state,
tl ,rl,ml,
@@ -414,8 +424,8 @@ mg_prolongate_all (comm_state & state,
timer.start ();
// Require same times
static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
- assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml+1))
- <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+ assert (abs(t.get_time(ml,rl,0) - t.get_time(ml+1,rl,0))
+ <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
vector<int> const tl2s(1,tl);
transfer_from_all (state,
tl ,rl,ml,
@@ -431,22 +441,19 @@ mg_prolongate_all (comm_state & state,
void
ggf::
ref_restrict_all (comm_state & state,
- int const tl, int const rl, int const ml,
- CCTK_REAL const time)
+ int const tl, int const rl, int const ml)
{
// Require same times
static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
- assert (abs(t.get_time(rl,ml) - t.get_time(rl+1,ml))
- <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+ assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
+ 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
if (transport_operator == op_none or transport_operator == op_sync) return;
static Timer timer ("ref_restrict_all");
timer.start ();
- vector<int> const tl2s(1,tl);
transfer_from_all (state,
- tl ,rl ,ml,
+ tl,rl ,ml,
& dh::fast_dboxes::fast_ref_rest_sendrecv,
- tl2s,rl+1,ml,
- time);
+ tl,rl+1,ml);
timer.stop (0);
}
@@ -522,7 +529,7 @@ transfer_from_all (comm_state & state,
for (size_t j=0; j<i; ++j) {
assert (tl2s.AT(i) != tl2s.AT(j));
}
- times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
+ times.AT(i) = t.get_time(ml2,rl2,tl2s.AT(i));
}
// Interpolation orders
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 43c5c2e93..5538f3941 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -105,6 +105,9 @@ public:
// Cycle the time levels by rotating the data sets
void cycle_all (int rl, int ml);
+
+ // Uncycle the time levels by rotating the data sets
+ void uncycle_all (int rl, int ml);
// Flip the time levels by exchanging the data sets
void flip_all (int rl, int ml);
@@ -135,7 +138,7 @@ public:
// Restrict a refinement level
void ref_restrict_all (comm_state& state,
- int tl, int rl, int ml, CCTK_REAL time);
+ int tl, int rl, int ml);
// Prolongate a refinement level
void ref_prolongate_all (comm_state& state,
@@ -171,7 +174,7 @@ protected:
{
vector <int> tl2s(1);
tl2s.AT(0) = tl2;
- CCTK_REAL const time = t.time (tl2,rl2,ml2);
+ CCTK_REAL const time = t.get_time (ml2,rl2,tl2);
transfer_from_all (state,
tl1, rl1, ml1,
sendrecvs,
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 21cb399a6..05c1b55ea 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -549,7 +549,7 @@ gh::
output (ostream & os)
const
{
- os << "gh:"
+ os << "gh:{"
<< "reffacts=" << reffacts << ",refcentering=" << refcent << ","
<< "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
<< "baseextents=" << baseextents << ","
@@ -565,6 +565,6 @@ output (ostream & os)
os << *d;
}
}
- os << "}";
+ os << "}}";
return os;
}
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index dd441e455..bed473245 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -19,8 +19,10 @@ list<th*> th::allth;
// Constructors
-th::th (gh& h_, const vector<int> & reffacts_, const CCTK_REAL basedelta)
- : h(h_), reffacts(reffacts_), delta(basedelta)
+th::th (gh& h_, int const timelevels_, vector<int> const& reffacts_,
+ bool const time_interpolation_during_regridding_)
+ : h(h_), timelevels(timelevels_), reffacts(reffacts_),
+ time_interpolation_during_regridding (time_interpolation_during_regridding_)
{
assert (reffacts.size() >= 1);
assert (reffacts.front() == 1);
@@ -42,6 +44,9 @@ th::~th ()
// Modifiers
void th::regrid ()
{
+ CCTK_REAL const basetime = 0.0;
+ CCTK_REAL const basedelta = 1.0;
+
const int old_mglevels = times.size();
times.resize(h.mglevels());
deltas.resize(h.mglevels());
@@ -50,17 +55,64 @@ void th::regrid ()
times.AT(ml).resize(h.reflevels());
deltas.AT(ml).resize(h.reflevels());
for (int rl=0; rl<h.reflevels(); ++rl) {
+ if (ml==0) {
+ deltas.AT(ml).AT(rl) = basedelta / reffacts.AT(rl);
+ } else {
+ deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact;
+ }
if (old_mglevels==0) {
- times.AT(ml).AT(rl) = 0;
+ times.AT(ml).AT(rl).resize(timelevels);
+ for (int tl=0; tl<timelevels; ++tl) {
+ times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl);
+ }
} else if (rl < old_reflevels) {
// do nothing
+ } else if (rl == 0) {
+ assert (0);
+ times.AT(ml).AT(rl).resize(timelevels);
+ for (int tl=0; tl<timelevels; ++tl) {
+ times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl);
+ }
} else {
- times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1);
+ if (time_interpolation_during_regridding) {
+#warning "We probably don't want to do this, but it is nice for compatibility"
+ times.AT(ml).AT(rl).resize(timelevels);
+ for (int tl=0; tl<timelevels; ++tl) {
+ // linear interpolation between the two surrounding coarse
+ // grid times
+ assert (reffacts.AT(rl) % reffacts.AT(rl-1) == 0);
+ int const rf = reffacts.AT(rl) / reffacts.AT(rl-1);
+ int const ctl = tl / rf;
+ int const mtl = tl % rf;
+ if (mtl == 0) {
+ times.AT(ml).AT(rl).AT(tl) = times.AT(ml).AT(rl-1).AT(ctl);
+ } else {
+ assert (ctl+1<timelevels);
+ CCTK_REAL const alpha = (CCTK_REAL)mtl / rf;
+ assert (alpha>0 and alpha<1);
+ times.AT(ml).AT(rl).AT(tl) =
+ (1-alpha) * times.AT(ml).AT(rl-1).AT(ctl ) +
+ ( alpha) * times.AT(ml).AT(rl-1).AT(ctl+1);
+ }
+ }
+ } else {
+ times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1);
+ }
}
- if (ml==0) {
- deltas.AT(ml).AT(rl) = delta / reffacts.AT(rl);
- } else {
- deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact;
+ }
+ }
+ for (int ml=0; ml<h.mglevels(); ++ml) {
+ for (int rl=0; rl<h.reflevels(); ++rl) {
+ for (int tl=1; tl<timelevels; ++tl) {
+ assert (times.AT(ml).AT(rl).AT(tl) < times.AT(ml).AT(rl).AT(tl-1));
+ }
+ }
+ }
+ for (int ml=0; ml<h.mglevels(); ++ml) {
+ for (int rl=0; rl<h.reflevels(); ++rl) {
+ // assert (isfinite(deltas.AT(ml).AT(rl)));
+ for (int tl=0; tl<timelevels; ++tl) {
+ // assert (isfinite(times.AT(ml).AT(rl).AT(tl)));
}
}
}
@@ -72,6 +124,37 @@ void th::regrid_free ()
+void th::advance_time (int const ml, int const rl)
+{
+ for (int tl=timelevels-1; tl>0; --tl) {
+ set_time (ml,rl,tl, get_time(ml,rl,tl-1));
+ }
+ set_time (ml,rl,0, get_time(ml,rl,0) + get_delta(ml,rl));
+}
+
+void th::retreat_time (int const ml, int const rl)
+{
+ CCTK_REAL const t = get_time(ml,rl,0);
+ for (int tl=0; tl<timelevels-1; ++tl) {
+ set_time (ml,rl,tl, get_time(ml,rl,tl+1));
+ }
+ set_time (ml,rl,timelevels-1, t);
+}
+
+void th::flip_timelevels (int const ml, int const rl)
+{
+ for (int tl=1; tl<(timelevels+1)/2; ++tl) {
+ int const tl2 = timelevels - tl;
+ CCTK_REAL const t = get_time(ml,rl,tl );
+ CCTK_REAL const t2 = get_time(ml,rl,tl2);
+ set_time (ml,rl,tl ,t2);
+ set_time (ml,rl,tl2,t );
+ }
+ set_delta (ml,rl, -get_delta(ml,rl));
+}
+
+
+
// Memory usage
size_t
th::
@@ -80,7 +163,6 @@ memory ()
{
return
memoryof (reffacts) +
- memoryof (delta) +
memoryof (times) +
memoryof (deltas);
}
@@ -100,19 +182,33 @@ allmemory ()
+// Input
+istream& th::input (istream& is)
+{
+ skipws (is);
+ consume (is, "th:{");
+ consume (is, "timelevels=");
+ is >> timelevels;
+ consume (is, ",");
+ consume (is, "reffacts=");
+ is >> reffacts;
+ consume (is, ",");
+ consume (is, "times=");
+ is >> times;
+ consume (is, ",");
+ consume (is, "deltas=");
+ is >> deltas;
+ consume (is, "}");
+ return is;
+}
+
// Output
-void th::output (ostream& os) const
+ostream& th::output (ostream& os) const
{
- os << "th:"
- << "times={";
- const char * sep = "";
- for (int ml=0; ml<h.mglevels(); ++ml) {
- for (int rl=0; rl<h.reflevels(); ++rl) {
- os << sep
- << ml << ":" << rl << ":"
- << times.AT(ml).AT(rl) << "(" << deltas.AT(ml).AT(rl) << ")";
- sep = ",";
- }
- }
- os << "}";
+ os << "th:{"
+ << "timelevels=" << timelevels << ","
+ << "reffacts=" << reffacts << ","
+ << "times=" << times << ","
+ << "deltas=" << deltas << "}";
+ return os;
}
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index f1a7b72b7..c1453cf4d 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -2,10 +2,11 @@
#define TH_HH
#include <cassert>
+#include <cmath>
#include <iostream>
#include <vector>
-#include "cctk.h"
+#include <cctk.h>
#include "defs.hh"
#include "gh.hh"
@@ -17,6 +18,8 @@ using namespace std;
// Forward declaration
class th;
+// Input
+istream& operator>> (istream& is, th& t);
// Output
ostream& operator<< (ostream& os, const th& t);
@@ -34,18 +37,22 @@ public: // should be readonly
gh& h; // hierarchy
gh::th_handle gh_handle;
+ int timelevels; // const
+
private:
- const vector<int> reffacts;
-
- CCTK_REAL delta; // time step
- vector<vector<CCTK_REAL> > times; // current times [ml][rl]
- vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
+ vector<int> reffacts; // const
+
+ bool const time_interpolation_during_regridding;
+
+ vector<vector<vector<CCTK_REAL> > > times; // current times [ml][rl][tl]
+ vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
public:
// Constructors
- th (gh& h, const vector<int> & reffacts, const CCTK_REAL basedelta);
+ th (gh& h, int timelevels, vector<int> const& reffacts,
+ bool time_interpolation_during_regridding);
// Destructors
~th ();
@@ -55,50 +62,52 @@ public:
void regrid_free ();
// Time management
- CCTK_REAL get_time (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+ void set_time (int const ml, int const rl, int const tl, CCTK_REAL const& t)
{
- assert (rl>=0 and rl<h.reflevels());
assert (ml>=0 and ml<h.mglevels());
- return times.AT(ml).AT(rl);
- }
-
- void set_time (const int rl, const int ml, const CCTK_REAL t)
- {
assert (rl>=0 and rl<h.reflevels());
- assert (ml>=0 and ml<h.mglevels());
- times.AT(ml).AT(rl) = t;
- }
-
- void advance_time (const int rl, const int ml)
- {
- set_time(rl,ml, get_time(rl,ml) + get_delta(rl,ml));
+ assert (tl>=0 and tl<timelevels);
+ // assert (isfinite(t));
+ times.AT(ml).AT(rl).AT(tl) = t;
}
- CCTK_REAL get_delta (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+ CCTK_REAL get_time (int const ml, int const rl, int const tl)
+ const CCTK_ATTRIBUTE_PURE
{
- assert (rl>=0 and rl<h.reflevels());
assert (ml>=0 and ml<h.mglevels());
- return deltas.AT(ml).AT(rl);
+ assert (rl>=0 and rl<h.reflevels());
+ assert (tl>=0 and tl<timelevels);
+ CCTK_REAL const t = times.AT(ml).AT(rl).AT(tl);
+ // assert (isfinite(t));
+ return t;
}
- void set_delta (const int rl, const int ml, const CCTK_REAL dt)
+ void set_delta (int const ml, int const rl, CCTK_REAL const& dt)
{
- assert (rl>=0 and rl<h.reflevels());
assert (ml>=0 and ml<h.mglevels());
+ assert (rl>=0 and rl<h.reflevels());
+ // assert (isfinite(dt));
deltas.AT(ml).AT(rl) = dt;
}
- CCTK_REAL time (const int tl, const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+ CCTK_REAL get_delta (int const ml, int const rl) const CCTK_ATTRIBUTE_PURE
{
- assert (rl>=0 and rl<h.reflevels());
assert (ml>=0 and ml<h.mglevels());
- return get_time(rl, ml) - tl * get_delta(rl, ml);
+ assert (rl>=0 and rl<h.reflevels());
+ CCTK_REAL const dt = deltas.AT(ml).AT(rl);
+ // assert (isfinite(dt));
+ return dt;
}
+ void advance_time (int const ml, int const rl);
+ void retreat_time (int const ml, int const rl);
+ void flip_timelevels (int const ml, int const rl);
+
// Output
size_t memory () const CCTK_ATTRIBUTE_PURE;
static size_t allmemory () CCTK_ATTRIBUTE_PURE;
- void output (ostream& os) const;
+ istream& input (istream& is);
+ ostream& output (ostream& os) const;
};
@@ -108,9 +117,11 @@ inline size_t memoryof (th const & t)
{
return t.memory ();
}
+inline istream& operator>> (istream& is, th& t) {
+ return t.input(is);
+}
inline ostream& operator<< (ostream& os, const th& t) {
- t.output(os);
- return os;
+ return t.output(os);
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 01:19:50 +0000