diff options
author | schnetter <> | 2001-12-09 15:43:00 +0000 |
---|---|---|
committer | schnetter <> | 2001-12-09 15:43:00 +0000 |
commit | 3319d0c8afdfe7e599b005be6f76656019f9b833 (patch) | |
tree | 602721e707302a87bbde0a05138409129a55e97b | |
parent | 2082bb5099476279b53877b27d9df295ab62ca2f (diff) |
Changed handling of interpolation orders; they are now stored in the
Changed handling of interpolation orders; they are now stored in the
grid functions and don't have to be passed around.
Speeded up some prolongation operators. They are still slow.
darcs-hash:20011209154308-07bb3-caa74a89a87c290852f2b59160ed26d9361f3cf1.gz
23 files changed, 252 insertions, 223 deletions
diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77 index e00dd8d54..8a07051dc 100644 --- a/Carpet/CarpetLib/src/copy_3d_real8.F77 +++ b/Carpet/CarpetLib/src/copy_3d_real8.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.4 2001/12/09 16:43:08 schnetter Exp $ #include "cctk.h" @@ -39,7 +39,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: strides disagree") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -47,7 +47,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index d21eb0bdd..efca6c699 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.13 2001/07/04 12:29:51 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.14 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -202,7 +202,8 @@ void data<T,D> ::interpolate_from_innerloop (const vector<const generic_data<D>*> gsrcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space) + const int order_space, + const int order_time) { assert (has_storage()); assert (all(box.lower()>=extent().lower())); @@ -223,14 +224,15 @@ void data<T,D> MPI_Comm_rank (dist::comm, &rank); assert (rank == proc()); - assert (order_space == 1); + assert (order_space <= 1); - vector<T> src_fac(srcs.size()); + assert (srcs.size() > order_time); + vector<T> src_fac(order_time+1); for (int t=0; t<(int)src_fac.size(); ++t) { src_fac[t] = 1; for (int tt=0; tt<(int)src_fac.size(); ++tt) { if (tt!=t) { - src_fac[t] *= (T)(t - tls[tt]) / (T)(tls[t] - tls[tt]); + src_fac[t] *= (T)(tl - tls[tt]) / (T)(tls[t] - tls[tt]); } } } @@ -413,7 +415,8 @@ void data<CCTK_REAL8,3> ::interpolate_from_innerloop (const vector<const generic_data<3>*> gsrcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space) + const int order_space, + const int order_time) { assert (has_storage()); assert (all(box.lower()>=extent().lower())); @@ -471,11 +474,12 @@ void data<CCTK_REAL8,3> } else if (all(sext.stride() > dext.stride())) { - switch (srcs.size()) { + switch (order_time) { - case 1: - // One timelevel + case 0: + assert (srcs.size()>=1); switch (order_space) { + case 0: case 1: CCTK_FNAME(prolongate_3d_real8) ((const CCTK_REAL8*)srcs[0]->storage(), @@ -484,6 +488,7 @@ void data<CCTK_REAL8,3> dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); break; + case 2: case 3: CCTK_FNAME(prolongate_3d_real8_o3) ((const CCTK_REAL8*)srcs[0]->storage(), @@ -497,9 +502,10 @@ void data<CCTK_REAL8,3> } break; - case 2: - // Two timelevels + case 1: + assert (srcs.size()>=2); switch (order_space) { + case 0: case 1: CCTK_FNAME(prolongate_3d_real8_2tl) ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], @@ -509,6 +515,7 @@ void data<CCTK_REAL8,3> dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); break; + case 2: case 3: CCTK_FNAME(prolongate_3d_real8_2tl_o3) ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], @@ -523,9 +530,10 @@ void data<CCTK_REAL8,3> } break; - case 3: - // Three timelevels + case 2: + assert (srcs.size()>=3); switch (order_space) { + case 0: case 1: CCTK_FNAME(prolongate_3d_real8_3tl) ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], @@ -536,6 +544,7 @@ void data<CCTK_REAL8,3> dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); break; + case 2: case 3: CCTK_FNAME(prolongate_3d_real8_3tl_o3) ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index 8d49998a7..00e239e22 100644 --- a/Carpet/CarpetLib/src/data.hh +++ b/Carpet/CarpetLib/src/data.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.7 2001/03/30 00:50:21 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.8 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -96,7 +96,8 @@ public: void interpolate_from_innerloop (const vector<const generic_data<D>*> gsrcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space); + const int order_space, + const int order_time); void write_ascii_output_element (ostream& os, const ivect& index) const; // void write_ieee (const string name, const int time, diff --git a/Carpet/CarpetLib/src/dgdh.cc b/Carpet/CarpetLib/src/dgdh.cc index 188ca8885..4287621c9 100644 --- a/Carpet/CarpetLib/src/dgdh.cc +++ b/Carpet/CarpetLib/src/dgdh.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.cc,v 1.1 2001/06/12 14:56:57 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.cc,v 1.2 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -36,10 +36,10 @@ using namespace std; // Constructors -dimgeneric_dh::dimgeneric_dh (int prolongation_order) - : prolongation_order(prolongation_order) +dimgeneric_dh::dimgeneric_dh (int prolongation_order_space) + : prolongation_order_space(prolongation_order_space) { - assert (prolongation_order>=0); + assert (prolongation_order_space>=0); } // Destructors @@ -49,8 +49,8 @@ dimgeneric_dh::~dimgeneric_dh () // Helpers int dimgeneric_dh::prolongation_stencil_size () const { - assert (prolongation_order>=0); - return prolongation_order/2; + assert (prolongation_order_space>=0); + return prolongation_order_space/2; } diff --git a/Carpet/CarpetLib/src/dgdh.hh b/Carpet/CarpetLib/src/dgdh.hh index ff7ed8b1e..222e2cfab 100644 --- a/Carpet/CarpetLib/src/dgdh.hh +++ b/Carpet/CarpetLib/src/dgdh.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.hh,v 1.1 2001/06/12 14:56:58 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.hh,v 1.2 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -48,12 +48,12 @@ class dimgeneric_dh { public: // should be readonly // Fields - int prolongation_order; // order of spatial prolongation operator + int prolongation_order_space; // order of spatial prolongation operator public: // Constructors - dimgeneric_dh (int prolongation_order); + dimgeneric_dh (int prolongation_order_space); // Destructors virtual ~dimgeneric_dh (); diff --git a/Carpet/CarpetLib/src/dggf.cc b/Carpet/CarpetLib/src/dggf.cc index 43c72a1df..998c5b90e 100644 --- a/Carpet/CarpetLib/src/dggf.cc +++ b/Carpet/CarpetLib/src/dggf.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.cc,v 1.1 2001/06/12 14:56:58 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.cc,v 1.2 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -39,8 +39,10 @@ using namespace std; // Constructors dimgeneric_gf::dimgeneric_gf (const string name, th& t, - const int tmin, const int tmax) - : name(name), t(t), tmin(tmin), tmax(tmax) + const int tmin, const int tmax, + const int prolongation_order_time) + : name(name), t(t), + tmin(tmin), tmax(tmax), prolongation_order_time(prolongation_order_time) { assert (tmin<=tmax+1); } diff --git a/Carpet/CarpetLib/src/dggf.hh b/Carpet/CarpetLib/src/dggf.hh index 76727afee..28004463e 100644 --- a/Carpet/CarpetLib/src/dggf.hh +++ b/Carpet/CarpetLib/src/dggf.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.hh,v 1.2 2001/07/02 13:22:13 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.hh,v 1.3 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -53,12 +53,14 @@ public: // should be readonly th &t; // time hierarchy int tmin, tmax; // timelevels + int prolongation_order_time; // order of temporal prolongation operator public: // Constructors dimgeneric_gf (const string name, th& t, - const int tmin, const int tmax); + const int tmin, const int tmax, + const int prolongation_order_time); // Destructors virtual ~dimgeneric_gf (); @@ -99,24 +101,19 @@ public: virtual void sync (int tl, int rl, int c, int ml) = 0; // Prolongate the boundaries of a component - virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml, - int order_space=1, int order_time=1) = 0; + virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml) = 0; // Restrict a multigrid level - virtual void mg_restrict (int tl, int rl, int c, int ml, int order_space=1) - = 0; + virtual void mg_restrict (int tl, int rl, int c, int ml) = 0; // Prolongate a multigrid level - virtual void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1) - = 0; + virtual void mg_prolongate (int tl, int rl, int c, int ml) = 0; // Restrict a refinement level - virtual void ref_restrict (int tl, int rl, int c, int ml, int order_space=1) - = 0; + virtual void ref_restrict (int tl, int rl, int c, int ml) = 0; // Prolongate a refinement level - virtual void ref_prolongate (int tl, int rl, int c, int ml, - int order_space=1) = 0; + virtual void ref_prolongate (int tl, int rl, int c, int ml) = 0; diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index f1512be42..2f3863dac 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.16 2001/07/04 12:29:51 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.17 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -39,8 +39,8 @@ using namespace std; // Constructors template<int D> dh<D>::dh (gh<D>& h, const ivect& lghosts, const ivect& ughosts, - int prolongation_order) - : dimgeneric_dh(prolongation_order), + int prolongation_order_space) + : dimgeneric_dh(prolongation_order_space), h(h), lghosts(lghosts), ughosts(ughosts) { diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index 7ba680691..73c2a54d6 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.8 2001/06/12 14:56:59 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.9 2001/12/09 16:43:09 schnetter Exp $ ***************************************************************************/ @@ -109,7 +109,7 @@ public: // Constructors dh (gh<D>& h, const ivect& lghosts, const ivect& ughosts, - int prolongation_order); + int prolongation_order_space); // Destructors virtual ~dh (); diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index 347bfb52c..41e265238 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.14 2001/07/04 12:29:52 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.15 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -90,7 +90,8 @@ void generic_data<D> ::interpolate_from (const vector<const generic_data*> srcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space) + const int order_space, + const int order_time) { assert (has_storage()); assert (all(box.lower()>=extent().lower())); @@ -110,7 +111,7 @@ void generic_data<D> int rank; MPI_Comm_rank (dist::comm, &rank); if (rank == proc()) { - interpolate_from_innerloop (srcs, tls, box, tl, order_space); + interpolate_from_innerloop (srcs, tls, box, tl, order_space, order_time); } } else { @@ -118,7 +119,7 @@ void generic_data<D> generic_data* const tmp = make_typed(); tmp->allocate (box, srcs[0]->proc()); - tmp->interpolate_from (srcs, tls, box, tl, order_space); + tmp->interpolate_from (srcs, tls, box, tl, order_space, order_time); tmp->change_processor (proc()); copy_from (tmp, box); delete tmp; diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh index 8e561c073..e827737f8 100644 --- a/Carpet/CarpetLib/src/gdata.hh +++ b/Carpet/CarpetLib/src/gdata.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.9 2001/06/12 14:56:59 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.10 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -102,7 +102,8 @@ public: void interpolate_from (const vector<const generic_data*> srcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space); + const int order_space, + const int order_time); protected: virtual void copy_from_innerloop (const generic_data* src, const ibbox& box) = 0; @@ -110,7 +111,8 @@ protected: interpolate_from_innerloop (const vector<const generic_data*> srcs, const vector<int> tls, const ibbox& box, const int tl, - const int order_space) = 0; + const int order_space, + const int order_time) = 0; public: diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc index 30eb4a5b4..4a026a574 100644 --- a/Carpet/CarpetLib/src/gf.cc +++ b/Carpet/CarpetLib/src/gf.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.7 2001/07/04 12:29:52 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.8 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -34,8 +34,8 @@ using namespace std; // Constructors template<class T,int D> gf<T,D>::gf (const string name, th& t, dh<D>& d, - const int tmin, const int tmax) - : generic_gf<D>(name, t, d, tmin, tmax) + const int tmin, const int tmax, const int prolongation_order_time) + : generic_gf<D>(name, t, d, tmin, tmax, prolongation_order_time) { recompose(); } diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh index 4a51bc98d..df500e074 100644 --- a/Carpet/CarpetLib/src/gf.hh +++ b/Carpet/CarpetLib/src/gf.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.4 2001/06/12 14:56:59 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.5 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -61,7 +61,8 @@ class gf: public generic_gf<D> { public: // Constructors - gf (const string name, th& t, dh<D>& d, const int tmin, const int tmax); + gf (const string name, th& t, dh<D>& d, + const int tmin, const int tmax, const int prolongation_order_time); // Destructors virtual ~gf (); diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index a81bb9bb2..0408417d4 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.12 2001/07/04 12:29:52 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.13 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -40,8 +40,9 @@ using namespace std; // Constructors template<int D> generic_gf<D>::generic_gf (const string name, th& t, dh<D>& d, - const int tmin, const int tmax) - : dimgeneric_gf(name, t, tmin, tmax), + const int tmin, const int tmax, + const int prolongation_order_time) + : dimgeneric_gf(name, t, tmin, tmax, prolongation_order_time), h(d.h), d(d), storage(tmax-tmin+1) { @@ -70,7 +71,9 @@ bool generic_gf<D>::operator== (const generic_gf<D>& f) const { template<int D> void generic_gf<D>::recompose () { // Retain storage that might be needed - fdata oldstorage(tmax-tmin+1); + fdata oldstorage; + + oldstorage.resize(tmax-tmin+1); for (int tl=tmin; tl<=tmax; ++tl) { oldstorage[tl-tmin].resize (min(h.reflevels(), (int)storage[tl-tmin].size())); @@ -88,9 +91,7 @@ void generic_gf<D>::recompose () { } // for rl } // for tl - // Delete storage - storage.clear(); - + // Resize structure storage.resize(tmax-tmin+1); for (int tl=tmin; tl<=tmax; ++tl) { storage[tl-tmin].resize(h.reflevels()); @@ -98,6 +99,13 @@ void generic_gf<D>::recompose () { storage[tl-tmin][rl].resize(h.components(rl)); for (int c=0; c<h.components(rl); ++c) { storage[tl-tmin][rl][c].resize(h.mglevels(rl,c)); + } // for c + } // for rl + } // for tl + + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int tl=tmin; tl<=tmax; ++tl) { + for (int c=0; c<h.components(rl); ++c) { for (int ml=0; ml<h.mglevels(rl,c); ++ml) { storage[tl-tmin][rl][c][ml] = typed_data(); @@ -258,8 +266,7 @@ template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, - const ibbox dh<D>::dboxes::* send_list, - int order_space) + const ibbox dh<D>::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); @@ -270,11 +277,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, } assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; - assert (ml2<h.mglevels(rl2,c2)); + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); vector<const generic_data<D>*> gsrcs(tl2s.size()); vector<int> tls(tl2s.size()); for (int i=0; i<(int)gsrcs.size(); ++i) { + assert (rl2<(int)storage[tl2s[i]-tmin].size()); + assert (c2<(int)storage[tl2s[i]-tmin][rl2].size()); + assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size()); gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; tls[i] = t.time(tl2s[i],rl2,ml2); } @@ -286,7 +296,8 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, // interpolate the content assert (recv==send); storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, tls, recv, tl, order_space); + (gsrcs, tls, recv, tl, + d.prolongation_order_space, prolongation_order_time); } // Interpolate regions @@ -294,8 +305,7 @@ template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, - const iblist dh<D>::dboxes::* send_list, - int order_space) + const iblist dh<D>::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); @@ -306,11 +316,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, } assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; - assert (ml2<h.mglevels(rl2,c2)); + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); vector<const generic_data<D>*> gsrcs(tl2s.size()); vector<int> tls(tl2s.size()); for (int i=0; i<(int)gsrcs.size(); ++i) { + assert (rl2<(int)storage[tl2s[i]-tmin].size()); + assert (c2<(int)storage[tl2s[i]-tmin][rl2].size()); + assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size()); gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; tls[i] = t.time(tl2s[i],rl2,ml2); } @@ -325,7 +338,8 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // interpolate the content storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, tls, *r, tl, order_space); + (gsrcs, tls, *r, tl, + d.prolongation_order_space, prolongation_order_time); } } @@ -334,8 +348,7 @@ template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, const vector<int> tl2s, int rl2, int ml2, - const iblistvect dh<D>::dboxes::* send_listvect, - int order_space) + const iblistvect dh<D>::dboxes::* send_listvect) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); @@ -347,11 +360,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, assert (rl2>=0 && rl2<h.reflevels()); // walk all components for (int c2=0; c2<h.components(rl2); ++c2) { - assert (ml2<h.mglevels(rl2,c2)); + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); vector<const generic_data<D>*> gsrcs(tl2s.size()); vector<int> tls(tl2s.size()); for (int i=0; i<(int)gsrcs.size(); ++i) { + assert (rl2<(int)storage[tl2s[i]-tmin].size()); + assert (c2<(int)storage[tl2s[i]-tmin][rl2].size()); + assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size()); gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; tls[i] = t.time(tl2s[i],rl2,ml2); } @@ -366,7 +382,8 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // interpolate the content storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, tls, *r, tl, order_space); + (gsrcs, tls, *r, tl, + d.prolongation_order_space, prolongation_order_time); } } } @@ -391,75 +408,83 @@ void generic_gf<D>::sync (int tl, int rl, int c, int ml) { // Prolongate the boundaries of a component template<int D> -void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml, - int order_space, int order_time) { +void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml) { // Interpolate assert (rl>=1); - const int tmod - = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1, ml) - + t.get_delta(rl-1, ml)) % t.get_delta(rl-1, ml); +// const int tmod +// = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml) +// + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml); vector<int> tl2s; - if (tmod == 0) { - // No interpolation in time - tl2s.resize(1); - tl2s[0] = tl; - } else { +// if (tmod == 0) { +// // No interpolation in time +// tl2s.resize(1); +// tl2s[0] = tl; +// } else { // Interpolation in time - assert (tmax-tmin+1 >= order_time+1); - tl2s.resize(order_time+1); - for (int i=0; i<=order_time; ++i) tl2s[i] = tmax - i; - } + assert (tmax-tmin+1 >= prolongation_order_time+1); + tl2s.resize(prolongation_order_time+1); + for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i; +// } intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, - tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine, - order_space); + tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine); } // Restrict a multigrid level template<int D> -void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml, - int order_space) { +void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml) +{ // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml-1)); const vector<int> tl2s(1,tl); intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, - tl2s,rl, ml-1, &dh<D>::dboxes::send_mg_fine, - order_space); + tl2s,rl, ml-1, &dh<D>::dboxes::send_mg_fine); } // Prolongate a multigrid level template<int D> -void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml, - int order_space) { +void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml) +{ // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml+1)); const vector<int> tl2s(1,tl); intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, - tl2s,rl, ml+1, &dh<D>::dboxes::send_mg_fine, - order_space); + tl2s,rl, ml+1, &dh<D>::dboxes::send_mg_fine); } // Restrict a refinement level template<int D> -void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml, - int order_space) { +void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml) +{ // Require same times assert (t.get_time(rl,ml) == t.get_time(rl+1,ml)); const vector<int> tl2s(1,tl); intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine, - tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse, - order_space); + tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse); } // Prolongate a refinement level template<int D> -void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml, - int order_space) { +void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml) +{ // Require same times assert (t.get_time(rl,ml) == t.get_time(rl-1,ml)); - const vector<int> tl2s(1,tl); + assert (rl>=1); +// const int tmod +// = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml) +// + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml); + vector<int> tl2s; +// if (tmod == 0) { +// // No interpolation in time +// tl2s.resize(1); +// tl2s[0] = tl; +// } else { + // Interpolation in time + assert (tmax-tmin+1 >= prolongation_order_time+1); + tl2s.resize(prolongation_order_time+1); + for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i; +// } intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse, - tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine, - order_space); + tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine); } diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 15d817b1b..f269d1772 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.7 2001/07/02 13:22:14 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.8 2001/12/09 16:43:10 schnetter Exp $ ***************************************************************************/ @@ -80,7 +80,8 @@ public: // Constructors generic_gf (const string name, th& t, dh<D>& d, - const int tmin, const int tmax); + const int tmin, const int tmax, + const int prolongation_order_time); // Destructors virtual ~generic_gf (); @@ -132,22 +133,19 @@ protected: void intercat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, - const ibbox dh<D>::dboxes::* send_list, - int order_space); + const ibbox dh<D>::dboxes::* send_list); // Interpolate regions void intercat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, - const iblist dh<D>::dboxes::* send_list, - int order_space); + const iblist dh<D>::dboxes::* send_list); // Interpolate regions void intercat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, const vector<int> tl2s, int rl2, int ml2, - const iblistvect dh<D>::dboxes::* send_listvect, - int order_space); + const iblistvect dh<D>::dboxes::* send_listvect); @@ -166,20 +164,19 @@ public: void sync (int tl, int rl, int c, int ml); // Prolongate the boundaries of a component - void ref_bnd_prolongate (int tl, int rl, int c, int ml, - int order_space=1, int order_time=1); + void ref_bnd_prolongate (int tl, int rl, int c, int ml); // Restrict a multigrid level - void mg_restrict (int tl, int rl, int c, int ml, int order_space=1); + void mg_restrict (int tl, int rl, int c, int ml); // Prolongate a multigrid level - void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1); + void mg_prolongate (int tl, int rl, int c, int ml); // Restrict a refinement level - void ref_restrict (int tl, int rl, int c, int ml, int order_space=1); + void ref_restrict (int tl, int rl, int c, int ml); // Prolongate a refinement level - void ref_prolongate (int tl, int rl, int c, int ml, int order_space=1); + void ref_prolongate (int tl, int rl, int c, int ml); diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 index 81939404a..cff34ddfb 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.4 2001/03/26 02:27:46 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.5 2001/12/09 16:43:10 schnetter Exp $ #include "cctk.h" @@ -29,12 +29,13 @@ c bbox(:,3) is stride integer srcioff, srcjoff, srckoff integer dstioff, dstjoff, dstkoff + CCTK_REAL8 dstdiv integer i, j, k integer i0, j0, k0 - CCTK_REAL8 fi, fj, fk - CCTK_REAL8 ifac(2), jfac(2), kfac(2) + integer fi, fj, fk + integer ifac(2), jfac(2), kfac(2) integer ii, jj, kk - CCTK_REAL8 fac + integer fac CCTK_REAL8 res integer d @@ -53,7 +54,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -61,7 +62,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do @@ -95,26 +96,25 @@ c bbox(:,3) is stride c Loop over fine region + dstdiv = one / (dstifac * dstjfac * dstkfac) + do k = 0, regkext-1 + k0 = (srckoff + k) / dstkfac + fk = mod(srckoff + k, dstkfac) + kfac(1) = (fk-dstkfac) * (-1) + kfac(2) = (fk ) * 1 + do j = 0, regjext-1 + j0 = (srcjoff + j) / dstjfac + fj = mod(srcjoff + j, dstjfac) + jfac(1) = (fj-dstjfac) * (-1) + jfac(2) = (fj ) * 1 + do i = 0, regiext-1 - i0 = (srcioff + i) / dstifac - j0 = (srcjoff + j) / dstjfac - k0 = (srckoff + k) / dstkfac - - fi = mod(srcioff + i, dstifac) * one / dstifac - fj = mod(srcjoff + j, dstjfac) * one / dstjfac - fk = mod(srckoff + k, dstkfac) * one / dstkfac - - ifac(1) = (fi-1) / (-1) - ifac(2) = (fi ) / 1 - - jfac(1) = (fj-1) / (-1) - jfac(2) = (fj ) / 1 - - kfac(1) = (fk-1) / (-1) - kfac(2) = (fk ) / 1 + fi = mod(srcioff + i, dstifac) + ifac(1) = (fi-dstifac) * (-1) + ifac(2) = (fi ) * 1 res = 0 @@ -132,7 +132,7 @@ c Loop over fine region end do end do - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 index 0cb837b3d..7b6e43da4 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -35,9 +35,11 @@ c bbox(:,3) is stride CCTK_REAL8 s1fac, s2fac + CCTK_REAL8 dstdiv integer i, j, k integer i0, j0, k0 integer fi, fj, fk + integer ifac(2), jfac(2), kfac(2) integer ii, jj, kk integer fac CCTK_REAL8 res @@ -58,7 +60,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -66,7 +68,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do @@ -110,41 +112,33 @@ c Linear (first order) interpolation c Loop over fine region + dstdiv = one / (dstifac * dstjfac * dstkfac) + do k = 0, regkext-1 + k0 = (srckoff + k) / dstkfac + fk = mod(srckoff + k, dstkfac) + kfac(1) = (fk-dstkfac) * (-1) + kfac(2) = (fk ) * 1 + do j = 0, regjext-1 + j0 = (srcjoff + j) / dstjfac + fj = mod(srcjoff + j, dstjfac) + jfac(1) = (fj-dstjfac) * (-1) + jfac(2) = (fj ) * 1 + do i = 0, regiext-1 - - i0 = (srcioff + i) / dstifac + 1 - j0 = (srcjoff + j) / dstjfac + 1 - k0 = (srckoff + k) / dstkfac + 1 - + i0 = (srcioff + i) / dstifac fi = mod(srcioff + i, dstifac) - fj = mod(srcjoff + j, dstjfac) - fk = mod(srckoff + k, dstkfac) + ifac(1) = (fi-dstifac) * (-1) + ifac(2) = (fi ) * 1 res = 0 - do kk=0,1 - do jj=0,1 - do ii=0,1 - - fac = 1 + do kk=1,2 + do jj=1,2 + do ii=1,2 - if (ii.eq.0) then - fac = fac * (dstifac - fi) - else - fac = fac * fi - end if - if (jj.eq.0) then - fac = fac * (dstjfac - fj) - else - fac = fac * fj - end if - if (kk.eq.0) then - fac = fac * (dstkfac - fk) - else - fac = fac * fk - end if + fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then res = res @@ -156,8 +150,7 @@ c Loop over fine region end do end do - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) - $ = res / (dstifac * dstjfac * dstkfac) + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 index 4f17e67c0..081150e87 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -58,7 +58,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -66,7 +66,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 index 17834ddd1..0ec7eb6ce 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.1 2001/03/24 22:38:48 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.2 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -60,7 +60,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -68,7 +68,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 index 7c16e8d35..eb4a006ec 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -60,7 +60,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -68,7 +68,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 index 5ee417885..906b67dd6 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -29,12 +29,13 @@ c bbox(:,3) is stride integer srcioff, srcjoff, srckoff integer dstioff, dstjoff, dstkoff + CCTK_REAL8 dstdiv integer i, j, k integer i0, j0, k0 - CCTK_REAL8 fi, fj, fk - CCTK_REAL8 ifac(4), jfac(4), kfac(4) + integer fi, fj, fk + integer ifac(4), jfac(4), kfac(4) integer ii, jj, kk - CCTK_REAL8 fac + integer fac CCTK_REAL8 res integer d @@ -53,7 +54,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -61,7 +62,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do @@ -95,32 +96,31 @@ c bbox(:,3) is stride c Loop over fine region + dstdiv = one / (6*dstifac * 6*dstjfac * 6*dstkfac) + do k = 0, regkext-1 + k0 = (srckoff + k) / dstkfac + fk = mod(srckoff + k, dstkfac) + kfac(1) = (fk ) * (fk-dstkfac) * (fk-2*dstkfac) * (-1) + kfac(2) = (fk+dstkfac) * (fk-dstkfac) * (fk-2*dstkfac) * 3 + kfac(3) = (fk+dstkfac) * (fk ) * (fk-2*dstkfac) * (-3) + kfac(4) = (fk+dstkfac) * (fk ) * (fk- dstkfac) * 1 + do j = 0, regjext-1 + j0 = (srcjoff + j) / dstjfac + fj = mod(srcjoff + j, dstjfac) + jfac(1) = (fj ) * (fj-dstjfac) * (fj-2*dstjfac) * (-1) + jfac(2) = (fj+dstjfac) * (fj-dstjfac) * (fj-2*dstjfac) * 3 + jfac(3) = (fj+dstjfac) * (fj ) * (fj-2*dstjfac) * (-3) + jfac(4) = (fj+dstjfac) * (fj ) * (fj- dstjfac) * 1 + do i = 0, regiext-1 - i0 = (srcioff + i) / dstifac - j0 = (srcjoff + j) / dstjfac - k0 = (srckoff + k) / dstkfac - - fi = mod(srcioff + i, dstifac) * one / dstifac - fj = mod(srcjoff + j, dstjfac) * one / dstjfac - fk = mod(srckoff + k, dstkfac) * one / dstkfac - - ifac(1) = (fi ) * (fi-1) * (fi-2) / (-6) - ifac(2) = (fi+1) * (fi-1) * (fi-2) / 2 - ifac(3) = (fi+1) * (fi ) * (fi-2) / (-2) - ifac(4) = (fi+1) * (fi ) * (fi-1) / 6 - - jfac(1) = (fj ) * (fj-1) * (fj-2) / (-6) - jfac(2) = (fj+1) * (fj-1) * (fj-2) / 2 - jfac(3) = (fj+1) * (fj ) * (fj-2) / (-2) - jfac(4) = (fj+1) * (fj ) * (fj-1) / 6 - - kfac(1) = (fk ) * (fk-1) * (fk-2) / (-6) - kfac(2) = (fk+1) * (fk-1) * (fk-2) / 2 - kfac(3) = (fk+1) * (fk ) * (fk-2) / (-2) - kfac(4) = (fk+1) * (fk ) * (fk-1) / 6 + fi = mod(srcioff + i, dstifac) + ifac(1) = (fi ) * (fi-dstifac) * (fi-2*dstifac) * (-1) + ifac(2) = (fi+dstifac) * (fi-dstifac) * (fi-2*dstifac) * 3 + ifac(3) = (fi+dstifac) * (fi ) * (fi-2*dstifac) * (-3) + ifac(4) = (fi+dstifac) * (fi ) * (fi- dstifac) * 1 res = 0 @@ -138,7 +138,7 @@ c Loop over fine region end do end do - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do end do diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77 index d580f423c..a8ba6c64c 100644 --- a/Carpet/CarpetLib/src/restrict_3d_real8.F77 +++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $ #include "cctk.h" @@ -44,7 +44,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: source strides are not integer multiples of the destination strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 + $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") end if @@ -52,7 +52,7 @@ c bbox(:,3) is stride $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") end if end do diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc index 18e26de1e..1e01d4787 100644 --- a/Carpet/CarpetLib/src/th.cc +++ b/Carpet/CarpetLib/src/th.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.5 2001/06/12 14:57:00 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.6 2001/12/09 16:43:11 schnetter Exp $ ***************************************************************************/ @@ -35,7 +35,8 @@ using namespace std; // Constructors -th::th (dimgeneric_gh* h, const int basedelta) : h(h), delta(basedelta) { +th::th (dimgeneric_gh* h, const int basedelta) + : h(h), delta(basedelta) { h->add(this); } |