diff options
author | eschnett <> | 2001-03-24 21:38:00 +0000 |
---|---|---|
committer | eschnett <> | 2001-03-24 21:38:00 +0000 |
commit | 63aad58172f57c2d2d91af630ae3b13a7cda70eb (patch) | |
tree | 054fe9b47ceca6067cb3dbe982d3c87c3c5b5255 | |
parent | 4f27fd634e6772a8075f9737c0d5e2f9545109fe (diff) |
Added support for higher-order interpolation in space and time.
darcs-hash:20010324213842-f6438-3ccfdd7797b28055f08d28b77e33205c69c60e27.gz
22 files changed, 639 insertions, 696 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl index 1690c5430..59ee423d4 100644 --- a/Carpet/Carpet/param.ccl +++ b/Carpet/Carpet/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn Carpet -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.5 2001/03/17 16:04:28 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.6 2001/03/24 22:38:42 eschnett Exp $ shares: Cactus @@ -71,6 +71,19 @@ CCTK_INT multigrid_factor "Multigrid factor" 1:* :: "must be positive" } 2 +CCTK_INT prolongation_order_space "Order of prolongation operator in space" +{ + 1 :: "first order (linear)" + 3 :: "third order (cubic)" +} 1 + +CCTK_INT prolongation_order_time "Order of prolongation operator in time" +{ + 0 :: "zeroth order (constant)" + 1 :: "first order (linear)" + 2 :: "second order (quadratic)" +} 1 + BOOLEAN verbose "Display info on the screen" diff --git a/Carpet/Carpet/src/carpet.cc b/Carpet/Carpet/src/carpet.cc index 010e783d2..f39e26126 100644 --- a/Carpet/Carpet/src/carpet.cc +++ b/Carpet/Carpet/src/carpet.cc @@ -35,7 +35,7 @@ #include "carpet.hh" -static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.17 2001/03/22 18:42:05 eschnett Exp $"; +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.18 2001/03/24 22:38:48 eschnett Exp $"; @@ -131,6 +131,8 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; + assert (cgh->cctk_dim == dim); + // Not sure what to do with that assert (convLevel==0); @@ -531,6 +533,8 @@ namespace Carpet { int SyncGroup (cGH* cgh, const char* groupname) { + DECLARE_CCTK_PARAMETERS; + assert (component == -1); Checkpoint ("%*sSyncGroup %s", 2*reflevel, "", groupname); @@ -560,7 +564,8 @@ namespace Carpet { if (num_tl>1 && reflevel>0) { for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { arrdata[group].data[var]->ref_bnd_prolongate - (tl, reflevel, c, mglevel); + (tl, reflevel, c, mglevel, + prolongation_order_space, prolongation_order_time); } } for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { @@ -575,7 +580,8 @@ namespace Carpet { if (num_tl>1 && reflevel>0) { for (int c=0; c<hh->components(reflevel); ++c) { gfdata[group].data[var]->ref_bnd_prolongate - (tl, reflevel, c, mglevel); + (tl, reflevel, c, mglevel, + prolongation_order_space, prolongation_order_time); } } for (int c=0; c<hh->components(reflevel); ++c) { @@ -908,7 +914,6 @@ namespace Carpet { if (CCTK_QueryGroupStorageI(cgh, group)) { for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) { const int n = CCTK_FirstVarIndexI(group) + var; - const int num_tl = CCTK_NumTimeLevelsFromVarI(n); switch (CCTK_GroupTypeFromVarI(n)) { case CCTK_SCALAR: { assert (group<(int)scdata.size()); diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77 index ae5d1c34f..e00dd8d54 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.2 2001/03/22 18:42:05 eschnett Exp $ +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 $ #include "cctk.h" @@ -43,6 +43,12 @@ c bbox(:,3) is stride $ .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 + if (regbbox(d,1).lt.srcbbox(d,1) + $ .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") + end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 23cd28b57..b75c770ff 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.8 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.9 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -198,109 +198,64 @@ void data<T,D> template<class T, int D> void data<T,D> -::interpolate_from_innerloop (const generic_data<D>* gsrc, const ibbox& box) +::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 data* src = (const data*)gsrc; - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); assert (all(box.stride()==extent().stride())); assert (all((box.lower()-extent().lower())%box.stride() == 0)); - - assert (proc() == src->proc()); + vector<const data*> srcs(gsrcs.size()); + for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t]; + assert (srcs.size() == tls.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); + assert (proc() == srcs[t]->proc()); + } int rank; MPI_Comm_rank (dist::comm, &rank); assert (rank == proc()); - for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { - const ivect& pos = *posi; - - // get box around current position - const ibbox frombox - = ibbox(pos,pos,extent().stride()).expanded_for(src->extent()); - - // interpolate from box to position - T sum = 0; - for (ibbox::iterator fromposi=frombox.begin(); - fromposi!=frombox.end(); ++fromposi) - { - const ivect& frompos = *fromposi; - - // interpolation weight - const ivect str = src->extent().stride(); - const T f = prod(vect<T,D>(str - abs(pos - frompos)) - / vect<T,D>(str)); - sum += f * (*src)[frompos]; - } - (*this)[pos] = sum; - - } // for pos - -} - - - -template<class T, int D> -void data<T,D> -::interpolate_from_innerloop (const generic_data<D>* gsrc1, const int t1, - const generic_data<D>* gsrc2, const int t2, - const ibbox& box, const int t) -{ - const data* src1 = (const data*)gsrc1; - const data* src2 = (const data*)gsrc2; - assert (has_storage() && src1->has_storage() && src2->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src1->extent().lower() - && box.lower()>=src2->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src1->extent().upper() - && box.upper()<=src2->extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src1->extent().lower())%box.stride() == 0 - && (box.lower()-src2->extent().lower())%box.stride() == 0)); - assert (src1->proc() == src2->proc()); + assert (order_space == 1); - assert (proc() == src1->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const T one = 1; - const T src1_fac = (T)((t - t2) * one / (t1 - t2)); - const T src2_fac = (T)((t - t1) * one / (t2 - t1)); + vector<T> src_fac(srcs.size()); + 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]); + } + } + } for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { const ivect& pos = *posi; // get box around current position const ibbox frombox - = ibbox(pos,pos,extent().stride()).expanded_for(src1->extent()); + = ibbox(pos,pos,extent().stride()).expanded_for(srcs[0]->extent()); // interpolate from box to position - T src1_sum = 0; - T src2_sum = 0; + T sum = 0; for (ibbox::iterator fromposi=frombox.begin(); fromposi!=frombox.end(); ++fromposi) { const ivect& frompos = *fromposi; // interpolation weight - const ivect str = src1->extent().stride(); - const T f = prod(vect<T,D>(str - abs(pos - frompos)) - / vect<T,D>(str)); - src1_sum += f * (*src1)[frompos]; - src2_sum += f * (*src2)[frompos]; + const ivect str = srcs[0]->extent().stride(); + const T f = prod(vect<T,D>(str - abs(pos - frompos)) / vect<T,D>(str)); + for (int t=0; t<(int)src_fac.size(); ++t) { + sum += f * src_fac[t] * (*srcs[t])[frompos]; + } } - (*this)[pos] = src1_fac * src1_sum + src2_fac * src2_sum; + (*this)[pos] = sum; } // for pos @@ -381,6 +336,18 @@ void data<CCTK_REAL8,3> extern "C" { + + void CCTK_FCALL CCTK_FNAME(restrict_3d_real8) + (const CCTK_REAL8* src, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); + + + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8) (const CCTK_REAL8* src, const int& srciext, const int& srcjext, const int& srckext, @@ -389,7 +356,7 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(restrict_3d_real8) + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3) (const CCTK_REAL8* src, const int& srciext, const int& srcjext, const int& srckext, CCTK_REAL8* dst, @@ -397,77 +364,7 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); -} - -template<> -void data<CCTK_REAL8,3> -::interpolate_from_innerloop (const generic_data<3>* gsrc, const ibbox& box) -{ - const data* src = (const data*)gsrc; - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - assert (proc() == src->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const ibbox& sext = src->extent(); - const ibbox& dext = extent(); - - int srcshp[3], dstshp[3]; - int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; - - for (int d=0; d<3; ++d) { - srcshp[d] = (sext.shape() / sext.stride())[d]; - dstshp[d] = (dext.shape() / dext.stride())[d]; - - srcbbox[0][d] = sext.lower()[d]; - srcbbox[1][d] = sext.upper()[d]; - srcbbox[2][d] = sext.stride()[d]; - - dstbbox[0][d] = dext.lower()[d]; - dstbbox[1][d] = dext.upper()[d]; - dstbbox[2][d] = dext.stride()[d]; - - regbbox[0][d] = box.lower()[d]; - regbbox[1][d] = box.upper()[d]; - regbbox[2][d] = box.stride()[d]; - } - - assert (all(dext.stride() == box.stride())); - if (all(sext.stride() > dext.stride())) { - CCTK_FNAME(prolongate_3d_real8) ((const CCTK_REAL8*)src->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); - } else if (all(sext.stride() < dext.stride())) { - CCTK_FNAME(restrict_3d_real8) ((const CCTK_REAL8*)src->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); - } else { - abort(); - } -} - - - -extern "C" { void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl) (const CCTK_REAL8* src1, const int& t1, const CCTK_REAL8* src2, const int& t2, @@ -477,40 +374,67 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3) + (const CCTK_REAL8* src1, const int& t1, + const CCTK_REAL8* src2, const int& t2, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, const int& t, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); + + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl) + (const CCTK_REAL8* src1, const int& t1, + const CCTK_REAL8* src2, const int& t2, + const CCTK_REAL8* src3, const int& t3, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, const int& t, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3) + (const CCTK_REAL8* src1, const int& t1, + const CCTK_REAL8* src2, const int& t2, + const CCTK_REAL8* src3, const int& t3, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, const int& t, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); + } template<> void data<CCTK_REAL8,3> -::interpolate_from_innerloop (const generic_data<3>* gsrc1, const int t1, - const generic_data<3>* gsrc2, const int t2, - const ibbox& box, const int t) +::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 data* src1 = (const data*)gsrc1; - const data* src2 = (const data*)gsrc2; - assert (has_storage() && src1->has_storage() && src2->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src1->extent().lower() - && box.lower()>=src2->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src1->extent().upper() - && box.upper()<=src2->extent().upper())); + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src1->extent().lower())%box.stride() == 0 - && (box.lower()-src2->extent().lower())%box.stride() == 0)); - assert (src1->proc() == src2->proc()); + assert (all((box.lower()-extent().lower())%box.stride() == 0)); + vector<const data*> srcs(gsrcs.size()); + for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t]; + assert (srcs.size() == tls.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); + } - assert (proc() == src1->proc()); + assert (proc() == srcs[0]->proc()); int rank; MPI_Comm_rank (dist::comm, &rank); assert (rank == proc()); - assert (src1->extent() == src2->extent()); - - const ibbox& sext = src1->extent(); + const ibbox& sext = srcs[0]->extent(); const ibbox& dext = extent(); int srcshp[3], dstshp[3]; @@ -534,16 +458,102 @@ void data<CCTK_REAL8,3> } assert (all(dext.stride() == box.stride())); - if (all(sext.stride() > dext.stride())) { - CCTK_FNAME(prolongate_3d_real8_2tl) - ((const CCTK_REAL8*)src1->storage(), t1, - (const CCTK_REAL8*)src2->storage(), t2, + if (all(sext.stride() < dext.stride())) { + + assert (srcs.size() == 1); + CCTK_FNAME(restrict_3d_real8) + ((const CCTK_REAL8*)srcs[0]->storage(), srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), t, + (CCTK_REAL8*)storage(), dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); + srcbbox, dstbbox, regbbox); + + } else if (all(sext.stride() > dext.stride())) { + + switch (srcs.size()) { + + case 1: + // One timelevel + switch (order_space) { + case 1: + CCTK_FNAME(prolongate_3d_real8) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 3: + CCTK_FNAME(prolongate_3d_real8_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + abort(); + } + break; + + case 2: + // Two timelevels + switch (order_space) { + case 1: + CCTK_FNAME(prolongate_3d_real8_2tl) + ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], + (const CCTK_REAL8*)srcs[1]->storage(), tls[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), tl, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 3: + CCTK_FNAME(prolongate_3d_real8_2tl_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], + (const CCTK_REAL8*)srcs[1]->storage(), tls[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), tl, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + abort(); + } + break; + + case 3: + // Three timelevels + switch (order_space) { + case 1: + CCTK_FNAME(prolongate_3d_real8_3tl) + ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], + (const CCTK_REAL8*)srcs[1]->storage(), tls[1], + (const CCTK_REAL8*)srcs[2]->storage(), tls[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), tl, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 3: + CCTK_FNAME(prolongate_3d_real8_3tl_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), tls[0], + (const CCTK_REAL8*)srcs[1]->storage(), tls[1], + (const CCTK_REAL8*)srcs[2]->storage(), tls[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), tl, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + abort(); + } + break; + + default: + abort(); + } + } else { abort(); } @@ -564,7 +574,7 @@ void data<T,D>::write_ascii_output_element (ofstream& file, const ivect& index) // Output template<class T,int D> ostream& data<T,D>::out (ostream& os) const { - os << "data<" STR(T) "," << D << ">:" + os << "data<T," << D << ">:" << "extent=" << extent() << "," << "stride=" << stride() << ",size=" << size(); return os; diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index c4594c56d..5182d0379 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.4 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.5 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -99,13 +99,12 @@ public: } // Data manipulators - void copy_from_innerloop (const generic_data<D>* src, + void copy_from_innerloop (const generic_data<D>* gsrc, const ibbox& box); - void interpolate_from_innerloop (const generic_data<D>* src, - const ibbox& box); - void interpolate_from_innerloop (const generic_data<D>* src1, const int t1, - const generic_data<D>* src2, const int t2, - const ibbox& box, const int t); + 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); void write_ascii_output_element (ofstream& file, const ivect& index) const; // void write_ieee (const string name, const int time, diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index bf483c62f..aa25b17b8 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.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/defs.cc,v 1.3 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.4 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -79,6 +79,7 @@ ostream& operator<< (ostream& os, const vector<T>& v) { template ostream& operator<< (ostream& os, const list<bbox<int,3> >& l); template ostream& operator<< (ostream& os, const set<bbox<int,3> >& s); template ostream& operator<< (ostream& os, const set<bboxset<int,3> >& s); +template ostream& operator<< (ostream& os, const vector<int>& v); template ostream& operator<< (ostream& os, const vector<list<bbox<int,3> > >& v); template ostream& operator<< (ostream& os, const vector<vector<bbox<int,3> > >& v); template ostream& operator<< (ostream& os, const vector<vector<vector<bbox<int,3> > > >& v); diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index 5418db23a..6f2a2ecd3 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.9 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.10 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -88,78 +88,38 @@ void generic_data<D>::copy_from (const generic_data* src, const ibbox& box) template<int D> void generic_data<D> -::interpolate_from (const generic_data* src, const ibbox& box) +::interpolate_from (const vector<const generic_data*> srcs, + const vector<int> tls, + const ibbox& box, const int tl, + const int order_space) { - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); assert (all(box.stride()==extent().stride())); assert (all((box.lower()-extent().lower())%box.stride() == 0)); - - if (proc() == src->proc()) { - // interpolate on same processor - - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - interpolate_from_innerloop (src, box); - } - - } else { - // interpolate from other processor - - generic_data* const tmp = make_typed(); - tmp->allocate (box, src->proc()); - tmp->interpolate_from (src, box); - tmp->change_processor (proc()); - copy_from (tmp, box); - delete tmp; - + assert (srcs.size() == tls.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); } -} - - - -template<int D> -void generic_data<D> -::interpolate_from (const generic_data* src1, const int t1, - const generic_data* src2, const int t2, - const ibbox& box, const int t) -{ - assert (has_storage() && src1->has_storage() && src2->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.upper()<=extent().upper())); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src1->extent().lower() - && box.lower()>=src2->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src1->extent().upper() - && box.upper()<=src2->extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src1->extent().lower())%box.stride() == 0 - && (box.lower()-src2->extent().lower())%box.stride() == 0)); - assert (src1->proc() == src2->proc()); - if (proc() == src1->proc()) { + if (proc() == srcs[0]->proc()) { // interpolate on same processor int rank; MPI_Comm_rank (dist::comm, &rank); if (rank == proc()) { - interpolate_from_innerloop (src1, t1, src2, t2, box, t); + interpolate_from_innerloop (srcs, tls, box, tl, order_space); } } else { // interpolate from other processor generic_data* const tmp = make_typed(); - tmp->allocate (box, src1->proc()); - tmp->interpolate_from (src1, t1, src2, t2, box, t); + tmp->allocate (box, srcs[0]->proc()); + tmp->interpolate_from (srcs, tls, box, tl, order_space); 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 1a16f92c2..ddb1c0603 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.5 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.6 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -135,19 +135,18 @@ public: // Data manipulators void copy_from (const generic_data* src, const ibbox& box); - void interpolate_from (const generic_data* src, const ibbox& box); - void interpolate_from (const generic_data* src1, const int t1, - const generic_data* src2, const int t2, - const ibbox& box, const int t); + void interpolate_from (const vector<const generic_data*> srcs, + const vector<int> tls, + const ibbox& box, const int tl, + const int order_space); protected: virtual void copy_from_innerloop (const generic_data* src, const ibbox& box) = 0; virtual void - interpolate_from_innerloop (const generic_data* src, const ibbox& box) =0; - virtual void - interpolate_from_innerloop (const generic_data* src1, const int t1, - const generic_data* src2, const int t2, - const ibbox& box, const int t) = 0; + interpolate_from_innerloop (const vector<const generic_data*> srcs, + const vector<int> tls, + const ibbox& box, const int tl, + const int order_space) = 0; public: diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 09966ac69..2fbb0eae3 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.6 2001/03/22 18:42:05 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.7 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -20,6 +20,7 @@ ***************************************************************************/ #include <assert.h> +#include <stdlib.h> #include <iostream> #include <string> @@ -170,7 +171,7 @@ void generic_gf<D>::cycle (int rl, int c, int ml) { // Operations -// Copy region for a component (between time levels) +// Copy a region template<int D> void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, @@ -194,7 +195,7 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, (storage[tl2-tmin][rl2][c2][ml2], recv); } -// Copy regions for a component (between multigrid levels) +// Copy regions template<int D> void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, @@ -217,12 +218,12 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, r!=recv.end(); ++r, ++s) { // (use the send boxes for communication) // copy the content - storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + storage[tl1-tmin][rl1][c1][ml1]->copy_from (storage[tl2-tmin][rl2][c2][ml2], *r); } } -// Copy regions for a level (between refinement levels) +// Copy regions template<int D> void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, @@ -246,54 +247,75 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, r!=recv.end(); ++r, ++s) { // (use the send boxes for communication) // copy the content - storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + storage[tl1-tmin][rl1][c1][ml1]->copy_from (storage[tl2-tmin][rl2][c2][ml2], *r); } } } -// Interpolate a component (between time levels) +// Interpolate a region template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, - int tl2a, int tl2b, int rl2, int ml2, - const ibbox dh<D>::dboxes::* send_list) + const vector<int> tl2s, int rl2, int ml2, + const ibbox dh<D>::dboxes::* send_list, + int order_space) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); assert (c1>=0 && c1<h.components(rl1)); assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); - assert (tl2a>=tmin && tl2a<=tmax); - assert (tl2b>=tmin && tl2b<=tmax); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s[i]>=tmin && tl2s[i]<=tmax); + } assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; assert (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) { + gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; + tls[i] = tl2s[i] * t.get_delta(rl2,ml2); + } + const int tl = tl1 * t.get_delta(rl1,ml1); + const ibbox recv = d.boxes[rl1][c1][ml1].*recv_list; const ibbox send = d.boxes[rl2][c2][ml2].*send_list; assert (recv.size()==send.size()); // interpolate the content assert (recv==send); storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (storage[tl2a-tmin][rl2][c2][ml2], tl2a, - storage[tl2b-tmin][rl2][c2][ml2], tl2b, recv, tl1); + (gsrcs, tls, recv, tl, order_space); } -// Interpolate a component (between multigrid levels) +// Interpolate regions template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, - int tl2a, int tl2b, int rl2, int ml2, - const iblist dh<D>::dboxes::* send_list) + const vector<int> tl2s, int rl2, int ml2, + const iblist dh<D>::dboxes::* send_list, + int order_space) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); assert (c1>=0 && c1<h.components(rl1)); assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); - assert (tl2a>=tmin && tl2a<=tmax); - assert (tl2b>=tmin && tl2b<=tmax); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s[i]>=tmin && tl2s[i]<=tmax); + } assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; assert (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) { + gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; + tls[i] = tl2s[i] * t.get_delta(rl2,ml2); + } + const int tl = tl1 * t.get_delta(rl1,ml1); + const iblist recv = d.boxes[rl1][c1][ml1].*recv_list; const iblist send = d.boxes[rl2][c2][ml2].*send_list; assert (recv.size()==send.size()); @@ -303,28 +325,38 @@ 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 - (storage[tl2a-tmin][rl2][c2][ml2], tl2a, - storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1); + (gsrcs, tls, *r, tl, order_space); } } -// Interpolate a level (between refinement levels) +// Interpolate regions template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, - int tl2a, int tl2b, int rl2, int ml2, - const iblistvect dh<D>::dboxes::* send_listvect) + const vector<int> tl2s, int rl2, int ml2, + const iblistvect dh<D>::dboxes::* send_listvect, + int order_space) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1<h.reflevels()); assert (c1>=0 && c1<h.components(rl1)); assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); - assert (tl2a>=tmin && tl2a<=tmax); - assert (tl2b>=tmin && tl2b<=tmax); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s[i]>=tmin && tl2s[i]<=tmax); + } assert (rl2>=0 && rl2<h.reflevels()); // walk all components for (int c2=0; c2<h.components(rl2); ++c2) { assert (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) { + gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2]; + tls[i] = tl2s[i] * t.get_delta(rl2,ml2); + } + const int tl = tl1 * t.get_delta(rl1,ml1); + const iblist recv = (d.boxes[rl1][c1][ml1].*recv_listvect)[c2]; const iblist send = (d.boxes[rl2][c2][ml2].*send_listvect)[c1]; assert (recv.size()==send.size()); @@ -334,8 +366,7 @@ 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 - (storage[tl2a-tmin][rl2][c2][ml2], tl2a, - storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1); + (gsrcs, tls, *r, tl, order_space); } } } @@ -360,59 +391,75 @@ 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) { +void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml, + int order_space, int order_time) { // Interpolate assert (rl>=1); - double time = - (t.time(tl,rl,ml) - t.get_time(rl-1,ml)) / (double)t.get_delta(rl-1, ml); - const int tl2 = (int)floor(time + 1e-6); - assert (tl2>=tmin && tl2<=tmax); - if (time==tl2) { - copycat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, - tl2,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine); + 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 { - const int tl2a=tl2; - const int tl2b=tl2+1; - assert (tl2b>=tmin && tl2b<=tmax); - intercat (tl,rl,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, - tl2a, tl2b, rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine); + // 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; } + 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); } // Restrict a multigrid level template<int D> -void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml) { +void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml, + int order_space) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml-1)); - copycat (tl,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, - tl,rl, ml-1, &dh<D>::dboxes::send_mg_fine); + 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); } // Prolongate a multigrid level template<int D> -void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml) { +void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml, + int order_space) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml+1)); - copycat (tl,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, - tl,rl, ml+1, &dh<D>::dboxes::send_mg_fine); + 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); } // Restrict a refinement level template<int D> -void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml) { +void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml, + int order_space) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl+1,ml)); - copycat (tl,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine, - tl,rl+1, ml, &dh<D>::dboxes::send_ref_coarse); + 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); } // Prolongate a refinement level template<int D> -void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml) { +void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml, + int order_space) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl-1,ml)); - copycat (tl,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse, - tl,rl-1, ml, &dh<D>::dboxes::send_ref_fine); + const vector<int> tl2s(1,tl); + intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse, + tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine, + order_space); } diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 6e998c145..6cacc0724 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.3 2001/03/22 18:42:06 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.4 2001/03/24 22:38:48 eschnett Exp $ ***************************************************************************/ @@ -26,6 +26,7 @@ #include <iostream> #include <string> +#include <vector> #include "defs.hh" #include "dh.hh" @@ -110,41 +111,44 @@ protected: protected: - // Copy region for a component (between time levels) + // Copy a region void copycat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list); - // Copy regions for a component (between multigrid levels) + // Copy regions void copycat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list); - // Copy regions for a level (between refinement levels) + // Copy regions void copycat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, int tl2, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect); - // Interpolate a component (between time levels) + // Interpolate a region void intercat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, - int tl2a, int tl2b, int rl2, int ml2, - const ibbox dh<D>::dboxes::* send_list); + const vector<int> tl2s, int rl2, int ml2, + const ibbox dh<D>::dboxes::* send_list, + int order_space); - // Interpolate a component (between multigrid levels) + // Interpolate regions void intercat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, - int tl2a, int tl2b, int rl2, int ml2, - const iblist dh<D>::dboxes::* send_list); + const vector<int> tl2s, int rl2, int ml2, + const iblist dh<D>::dboxes::* send_list, + int order_space); - // Interpolate a level (between refinement levels) + // Interpolate regions void intercat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, - int tl2a, int tl2b, int rl2, int ml2, - const iblistvect dh<D>::dboxes::* send_listvect); + const vector<int> tl2s, int rl2, int ml2, + const iblistvect dh<D>::dboxes::* send_listvect, + int order_space); @@ -163,19 +167,20 @@ 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); + void ref_bnd_prolongate (int tl, int rl, int c, int ml, + int order_space=1, int order_time=1); // Restrict a multigrid level - void mg_restrict (int tl, int rl, int c, int ml); + void mg_restrict (int tl, int rl, int c, int ml, int order_space=1); // Prolongate a multigrid level - void mg_prolongate (int tl, int rl, int c, int ml); + void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1); // Restrict a refinement level - void ref_restrict (int tl, int rl, int c, int ml); + void ref_restrict (int tl, int rl, int c, int ml, int order_space=1); // Prolongate a refinement level - void ref_prolongate (int tl, int rl, int c, int ml); + void ref_prolongate (int tl, int rl, int c, int ml, int order_space=1); diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index dd1756780..b072b8c62 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -1,9 +1,16 @@ # Main make.code.defn file for thorn CarpetLib -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.3 2001/03/22 18:42:06 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.4 2001/03/24 22:38:48 eschnett Exp $ # Source files in this directory SRCS = bbox.cc bboxset.cc data.cc defs.cc dh.cc dist.cc gdata.cc gf.cc ggf.cc gh.cc th.cc vect.cc \ - copy_3d_real8.F77 prolongate_3d_real8.F77 prolongate_3d_real8_2tl.F77 restrict_3d_real8.F77 + copy_3d_real8.F77 \ + prolongate_3d_real8.F77 \ + prolongate_3d_real8_o3.F77 \ + prolongate_3d_real8_2tl.F77 \ + prolongate_3d_real8_2tl_o3.F77 \ + prolongate_3d_real8_3tl.F77 \ + prolongate_3d_real8_3tl_o3.F77 \ + restrict_3d_real8.F77 # Note: skipping io.cc # Subdirectories containing source files diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 index 04053c524..531a3a6a7 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.2 2001/03/22 18:42:06 eschnett Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $ #include "cctk.h" @@ -53,6 +53,12 @@ c bbox(:,3) is stride $ .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 + if (regbbox(d,1).lt.srcbbox(d,1) + $ .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") + end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 index f9b380c5d..0cb837b3d 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.2 2001/03/22 18:42:06 eschnett Exp $ +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 $ #include "cctk.h" @@ -62,6 +62,12 @@ c bbox(:,3) is stride $ .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 + if (regbbox(d,1).lt.srcbbox(d,1) + $ .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") + end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -93,6 +99,11 @@ c bbox(:,3) is stride +c Linear (first order) interpolation + if (t1.eq.t2) then + call CCTK_WARN (0, "Internal error: arrays have same time") + end if + s1fac = (t - t2) * one / (t1 - t2) s2fac = (t - t1) * one / (t2 - t1) diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 index 0f0509cff..554c4fc38 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.13 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.1 2001/03/24 22:38:48 eschnett Exp $ #include "cctk.h" - - - -#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ - if ((i).lt.1 .or. (i).gt.(imax) \ - .or. (j).lt.1 .or. (j).gt.(jmax) \ - .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ - write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ - (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ - call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ - end if - - subroutine prolongate_3d_real8_2tl_o3 ( $ src1, t1, src2, t2, srciext, srcjext, srckext, @@ -26,24 +13,19 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d CCTK_REAL8 one parameter (one = 1) - CCTK_REAL8 eps - parameter (eps = 1.0d-10) - integer srciext, srcjext, srckext CCTK_REAL8 src1(srciext,srcjext,srckext) - CCTK_REAL8 t1 + integer t1 CCTK_REAL8 src2(srciext,srcjext,srckext) - CCTK_REAL8 t2 + integer t2 integer dstiext, dstjext, dstkext CCTK_REAL8 dst(dstiext,dstjext,dstkext) - CCTK_REAL8 t + integer t c bbox(:,1) is lower boundary (inclusive) c bbox(:,2) is upper boundary (inclusive) c bbox(:,3) is stride integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3) - integer offsetlo, offsethi - integer regiext, regjext, regkext integer dstifac, dstjfac, dstkfac @@ -53,18 +35,14 @@ 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(4), jfac(4), kfac(4) integer ii, jj, kk - integer fac + CCTK_REAL8 fi, fj, fk + CCTK_REAL8 fac CCTK_REAL8 res integer d - character msg*1000 - do d=1,3 @@ -80,41 +58,15 @@ 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 - if (regbbox(d,1).gt.regbbox(d,2)) then -c This could be handled, but is likely to point to an error elsewhere - call CCTK_WARN (0, "Internal error: region extent is empty") - end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if - regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 - dstkfac = srcbbox(d,3) / dstbbox(d,3) - srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) - offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) then - offsetlo = 0 - if (regkext.gt.1) then - offsetlo = regbbox(d,3) - end if - end if - offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) then - offsethi = 0 - if (regkext.gt.1) then - offsethi = regbbox(d,3) - end if - end if - if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) - $ .or. regbbox(d,2)+offsethi.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 extent") + if (regbbox(d,1)-1.lt.srcbbox(d,1) + $ .or. regbbox(d,1)-1.lt.dstbbox(d,1) + $ .or. regbbox(d,2)+1.gt.srcbbox(d,2) + $ .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then + call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") end if end do @@ -151,65 +103,74 @@ c Linear (first order) interpolation if (t1.eq.t2) then call CCTK_WARN (0, "Internal error: arrays have same time") end if - if (t.lt.min(t1,t2)-eps .or. t.gt.max(t1,t2)+eps) then - call CCTK_WARN (0, "Internal error: extrapolation in time") - end if - s1fac = (t - t2) / (t1 - t2) - s2fac = (t - t1) / (t2 - t1) + s1fac = (t - t2) * one / (t1 - t2) + s2fac = (t - t1) * one / (t2 - t1) c Loop over fine region - dstdiv = one / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3) - 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 - 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 + + i0 = (srcioff + i) / dstifac + 1 + j0 = (srcjoff + j) / dstjfac + 1 + k0 = (srckoff + k) / dstkfac + 1 + + fi = mod(srcioff + i, dstifac) * one / dstifac + fj = mod(srcjoff + j, dstjfac) * one / dstjfac + fk = mod(srckoff + k, dstkfac) * one / dstkfac res = 0 - do kk=1,4 - do jj=1,4 - do ii=1,4 + do kk=-1,2 + do jj=-1,2 + do ii=-1,2 - fac = ifac(ii) * jfac(jj) * kfac(kk) + fac = 1 + + if (ii.eq.-1) then + fac = fac * (fi ) * (fi-1) * (fi-2) / (-6) + else if (ii.eq.0) then + fac = fac * (fi+1) * (fi-1) * (fi-2) / 2 + else if (ii.eq.1) then + fac = fac * (fi+1) * (fi ) * (fi-2) / (-2) + else + fac = fac * (fi+1) * (fi ) * (fi-1) / 6 + end if + + if (jj.eq.-1) then + fac = fac * (fj ) * (fj-1) * (fj-2) / (-6) + else if (jj.eq.0) then + fac = fac * (fj+1) * (fj-1) * (fj-2) / 2 + else if (jj.eq.1) then + fac = fac * (fj+1) * (fj ) * (fj-2) / (-2) + else + fac = fac * (fj+1) * (fj ) * (fj-1) / 6 + end if + + if (kk.eq.-1) then + fac = fac * (fk ) * (fk-1) * (fk-2) / (-6) + else if (kk.eq.0) then + fac = fac * (fk+1) * (fk-1) * (fk-2) / 2 + else if (kk.eq.1) then + fac = fac * (fk+1) * (fk ) * (fk-2) / (-2) + else + fac = fac * (fk+1) * (fk ) * (fk-1) / 6 + end if if (fac.ne.0) then - CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ - srciext,srcjext,srckext, "source") res = res - $ + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1) - $ + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1) + $ + fac * s1fac * src1(i0+ii,j0+jj,k0+kk) + $ + fac * s2fac * src2(i0+ii,j0+jj,k0+kk) end if end do end do end do - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res end do end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 index 9f3e3d944..17834ddd1 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.9 2004/03/11 12:03:09 schnetter Exp $ +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 $ #include "cctk.h" - - - -#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ - if ((i).lt.1 .or. (i).gt.(imax) \ - .or. (j).lt.1 .or. (j).gt.(jmax) \ - .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ - write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ - (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ - call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ - end if - - subroutine prolongate_3d_real8_3tl ( $ src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext, @@ -26,19 +13,16 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d CCTK_REAL8 one parameter (one = 1) - CCTK_REAL8 eps - parameter (eps = 1.0d-10) - integer srciext, srcjext, srckext CCTK_REAL8 src1(srciext,srcjext,srckext) - CCTK_REAL8 t1 + integer t1 CCTK_REAL8 src2(srciext,srcjext,srckext) - CCTK_REAL8 t2 + integer t2 CCTK_REAL8 src3(srciext,srcjext,srckext) - CCTK_REAL8 t3 + integer t3 integer dstiext, dstjext, dstkext CCTK_REAL8 dst(dstiext,dstjext,dstkext) - CCTK_REAL8 t + integer t c bbox(:,1) is lower boundary (inclusive) c bbox(:,2) is upper boundary (inclusive) c bbox(:,3) is stride @@ -53,18 +37,14 @@ c bbox(:,3) is stride CCTK_REAL8 s1fac, s2fac, s3fac - 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 integer d - character msg*1000 - do d=1,3 @@ -80,24 +60,15 @@ 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 - if (regbbox(d,1).gt.regbbox(d,2)) then -c This could be handled, but is likely to point to an error elsewhere - call CCTK_WARN (0, "Internal error: region extent is empty") - end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if if (regbbox(d,1).lt.srcbbox(d,1) $ .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 extent") + call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") end if end do @@ -134,61 +105,63 @@ c Quadratic (second order) interpolation if (t1.eq.t2 .or. t1.eq.t3 .or. t2.eq.t3) then call CCTK_WARN (0, "Internal error: arrays have same time") end if - if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then - call CCTK_WARN (0, "Internal error: extrapolation in time") - end if - s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3)) - s2fac = (t - t1) * (t - t3) / ((t2 - t1) * (t2 - t3)) - s3fac = (t - t1) * (t - t2) / ((t3 - t1) * (t3 - t2)) + s1fac = (t - t2) * (t - t3) * one / ((t1 - t2) * (t1 - t3)) + s2fac = (t - t1) * (t - t3) * one / ((t2 - t1) * (t2 - t3)) + s3fac = (t - t1) * (t - t2) * one / ((t3 - t1) * (t3 - t2)) 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 + + i0 = (srcioff + i) / dstifac + 1 + j0 = (srcjoff + j) / dstjfac + 1 + k0 = (srckoff + k) / dstkfac + 1 + fi = mod(srcioff + i, dstifac) - ifac(1) = (fi-dstifac) * (-1) - ifac(2) = (fi ) * 1 + fj = mod(srcjoff + j, dstjfac) + fk = mod(srckoff + k, dstkfac) res = 0 - do kk=1,2 - do jj=1,2 - do ii=1,2 + do kk=0,1 + do jj=0,1 + do ii=0,1 + + fac = 1 - fac = ifac(ii) * jfac(jj) * kfac(kk) + 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 if (fac.ne.0) then - CHKIDX (i0+ii, j0+jj, k0+kk, \ - srciext,srcjext,srckext, "source") res = res - $ + fac * s1fac * src1(i0+ii, j0+jj, k0+kk) - $ + fac * s2fac * src2(i0+ii, j0+jj, k0+kk) - $ + fac * s3fac * src3(i0+ii, j0+jj, k0+kk) + $ + fac * s1fac * src1(i0+ii,j0+jj,k0+kk) + $ + fac * s2fac * src2(i0+ii,j0+jj,k0+kk) + $ + fac * s3fac * src3(i0+ii,j0+jj,k0+kk) end if end do end do end do - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) + $ = res / (dstifac * dstjfac * dstkfac) end do 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 c44e1119f..ad360647c 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.13 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.1 2001/03/24 22:38:48 eschnett Exp $ #include "cctk.h" - - - -#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ - if ((i).lt.1 .or. (i).gt.(imax) \ - .or. (j).lt.1 .or. (j).gt.(jmax) \ - .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ - write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ - (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ - call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ - end if - - subroutine prolongate_3d_real8_3tl_o3 ( $ src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext, @@ -26,26 +13,21 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d CCTK_REAL8 one parameter (one = 1) - CCTK_REAL8 eps - parameter (eps = 1.0d-10) - integer srciext, srcjext, srckext CCTK_REAL8 src1(srciext,srcjext,srckext) - CCTK_REAL8 t1 + integer t1 CCTK_REAL8 src2(srciext,srcjext,srckext) - CCTK_REAL8 t2 + integer t2 CCTK_REAL8 src3(srciext,srcjext,srckext) - CCTK_REAL8 t3 + integer t3 integer dstiext, dstjext, dstkext CCTK_REAL8 dst(dstiext,dstjext,dstkext) - CCTK_REAL8 t + integer t c bbox(:,1) is lower boundary (inclusive) c bbox(:,2) is upper boundary (inclusive) c bbox(:,3) is stride integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3) - integer offsetlo, offsethi - integer regiext, regjext, regkext integer dstifac, dstjfac, dstkfac @@ -55,18 +37,14 @@ c bbox(:,3) is stride CCTK_REAL8 s1fac, s2fac, s3fac - CCTK_REAL8 dstdiv integer i, j, k integer i0, j0, k0 - integer fi, fj, fk - integer ifac(4), jfac(4), kfac(4) integer ii, jj, kk - integer fac + CCTK_REAL8 fi, fj, fk + CCTK_REAL8 fac CCTK_REAL8 res integer d - character msg*1000 - do d=1,3 @@ -82,41 +60,15 @@ 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 - if (regbbox(d,1).gt.regbbox(d,2)) then -c This could be handled, but is likely to point to an error elsewhere - call CCTK_WARN (0, "Internal error: region extent is empty") - end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if - regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 - dstkfac = srcbbox(d,3) / dstbbox(d,3) - srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) - offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) then - offsetlo = 0 - if (regkext.gt.1) then - offsetlo = regbbox(d,3) - end if - end if - offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) then - offsethi = 0 - if (regkext.gt.1) then - offsethi = regbbox(d,3) - end if - end if - if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) - $ .or. regbbox(d,2)+offsethi.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 extent") + if (regbbox(d,1)-1.lt.srcbbox(d,1) + $ .or. regbbox(d,1)-1.lt.dstbbox(d,1) + $ .or. regbbox(d,2)+1.gt.srcbbox(d,2) + $ .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then + call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") end if end do @@ -153,67 +105,76 @@ c Quadratic (second order) interpolation if (t1.eq.t2 .or. t1.eq.t3 .or. t2.eq.t3) then call CCTK_WARN (0, "Internal error: arrays have same time") end if - if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then - call CCTK_WARN (0, "Internal error: extrapolation in time") - end if - s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3)) - s2fac = (t - t1) * (t - t3) / ((t2 - t1) * (t2 - t3)) - s3fac = (t - t1) * (t - t2) / ((t3 - t1) * (t3 - t2)) + s1fac = (t - t2) * (t - t3) * one / ((t1 - t2) * (t1 - t3)) + s2fac = (t - t1) * (t - t3) * one / ((t2 - t1) * (t2 - t3)) + s3fac = (t - t1) * (t - t2) * one / ((t3 - t1) * (t3 - t2)) c Loop over fine region - dstdiv = one / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3) - 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 - 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 + + i0 = (srcioff + i) / dstifac + 1 + j0 = (srcjoff + j) / dstjfac + 1 + k0 = (srckoff + k) / dstkfac + 1 + + fi = mod(srcioff + i, dstifac) * one / dstifac + fj = mod(srcjoff + j, dstjfac) * one / dstjfac + fk = mod(srckoff + k, dstkfac) * one / dstkfac res = 0 - do kk=1,4 - do jj=1,4 - do ii=1,4 + do kk=-1,2 + do jj=-1,2 + do ii=-1,2 - fac = ifac(ii) * jfac(jj) * kfac(kk) + fac = 1 + + if (ii.eq.-1) then + fac = fac * (fi ) * (fi-1) * (fi-2) / (-6) + else if (ii.eq.0) then + fac = fac * (fi+1) * (fi-1) * (fi-2) / 2 + else if (ii.eq.1) then + fac = fac * (fi+1) * (fi ) * (fi-2) / (-2) + else + fac = fac * (fi+1) * (fi ) * (fi-1) / 6 + end if + + if (jj.eq.-1) then + fac = fac * (fj ) * (fj-1) * (fj-2) / (-6) + else if (jj.eq.0) then + fac = fac * (fj+1) * (fj-1) * (fj-2) / 2 + else if (jj.eq.1) then + fac = fac * (fj+1) * (fj ) * (fj-2) / (-2) + else + fac = fac * (fj+1) * (fj ) * (fj-1) / 6 + end if + + if (kk.eq.-1) then + fac = fac * (fk ) * (fk-1) * (fk-2) / (-6) + else if (kk.eq.0) then + fac = fac * (fk+1) * (fk-1) * (fk-2) / 2 + else if (kk.eq.1) then + fac = fac * (fk+1) * (fk ) * (fk-2) / (-2) + else + fac = fac * (fk+1) * (fk ) * (fk-1) / 6 + end if if (fac.ne.0) then - CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ - srciext,srcjext,srckext, "source") res = res - $ + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1) - $ + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1) - $ + fac * s3fac * src3(i0+ii-1, j0+jj-1, k0+kk-1) + $ + fac * s1fac * src1(i0+ii,j0+jj,k0+kk) + $ + fac * s2fac * src2(i0+ii,j0+jj,k0+kk) + $ + fac * s3fac * src3(i0+ii,j0+jj,k0+kk) end if end do end do end do - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res end do end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 index 21d82a733..c539ee7de 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.9 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.1 2001/03/24 22:38:48 eschnett Exp $ #include "cctk.h" - - - -#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ - if ((i).lt.1 .or. (i).gt.(imax) \ - .or. (j).lt.1 .or. (j).gt.(jmax) \ - .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ - write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ - (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ - call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ - end if - - subroutine prolongate_3d_real8_o3 ( $ src, srciext, srcjext, srckext, @@ -35,8 +22,6 @@ c bbox(:,2) is upper boundary (inclusive) c bbox(:,3) is stride integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3) - integer offsetlo, offsethi - integer regiext, regjext, regkext integer dstifac, dstjfac, dstkfac @@ -44,18 +29,14 @@ c bbox(:,3) is stride integer srcioff, srcjoff, srckoff integer dstioff, dstjoff, dstkoff - CCTK_REAL8 dstdiv integer i, j, k integer i0, j0, k0 - integer fi, fj, fk - integer ifac(4), jfac(4), kfac(4) integer ii, jj, kk - integer fac + CCTK_REAL8 fi, fj, fk + CCTK_REAL8 fac CCTK_REAL8 res integer d - character msg*1000 - do d=1,3 @@ -71,41 +52,15 @@ 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 - if (regbbox(d,1).gt.regbbox(d,2)) then -c This could be handled, but is likely to point to an error elsewhere - call CCTK_WARN (0, "Internal error: region extent is empty") - end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if - regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 - dstkfac = srcbbox(d,3) / dstbbox(d,3) - srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) - offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) then - offsetlo = 0 - if (regkext.gt.1) then - offsetlo = regbbox(d,3) - end if - end if - offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) then - offsethi = 0 - if (regkext.gt.1) then - offsethi = regbbox(d,3) - end if - end if - if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) - $ .or. regbbox(d,2)+offsethi.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 extent") + if (regbbox(d,1)-1.lt.srcbbox(d,1) + $ .or. regbbox(d,1)-1.lt.dstbbox(d,1) + $ .or. regbbox(d,2)+1.gt.srcbbox(d,2) + $ .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then + call CCTK_WARN (0, "Internal error: region extent is not contained in array extents") end if end do @@ -139,53 +94,65 @@ c This could be handled, but is likely to point to an error elsewhere c Loop over fine region - dstdiv = one / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3) - 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 - 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 + + i0 = (srcioff + i) / dstifac + 1 + j0 = (srcjoff + j) / dstjfac + 1 + k0 = (srckoff + k) / dstkfac + 1 + + fi = mod(srcioff + i, dstifac) * one / dstifac + fj = mod(srcjoff + j, dstjfac) * one / dstjfac + fk = mod(srckoff + k, dstkfac) * one / dstkfac res = 0 - do kk=1,4 - do jj=1,4 - do ii=1,4 + do kk=-1,2 + do jj=-1,2 + do ii=-1,2 + + fac = 1 + + if (ii.eq.-1) then + fac = fac * (fi ) * (fi-1) * (fi-2) / (-6) + else if (ii.eq.0) then + fac = fac * (fi+1) * (fi-1) * (fi-2) / 2 + else if (ii.eq.1) then + fac = fac * (fi+1) * (fi ) * (fi-2) / (-2) + else + fac = fac * (fi+1) * (fi ) * (fi-1) / 6 + end if + + if (jj.eq.-1) then + fac = fac * (fj ) * (fj-1) * (fj-2) / (-6) + else if (jj.eq.0) then + fac = fac * (fj+1) * (fj-1) * (fj-2) / 2 + else if (jj.eq.1) then + fac = fac * (fj+1) * (fj ) * (fj-2) / (-2) + else + fac = fac * (fj+1) * (fj ) * (fj-1) / 6 + end if - fac = ifac(ii) * jfac(jj) * kfac(kk) + if (kk.eq.-1) then + fac = fac * (fk ) * (fk-1) * (fk-2) / (-6) + else if (kk.eq.0) then + fac = fac * (fk+1) * (fk-1) * (fk-2) / 2 + else if (kk.eq.1) then + fac = fac * (fk+1) * (fk ) * (fk-2) / (-2) + else + fac = fac * (fk+1) * (fk ) * (fk-1) / 6 + end if if (fac.ne.0) then - CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ - srciext,srcjext,srckext, "source") - res = res + fac * src(i0+ii-1, j0+jj-1, k0+kk-1) + res = res + fac * src(i0+ii,j0+jj,k0+kk) end if end do end do end do - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res end do end do diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77 index 68ea98b11..d580f423c 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.2 2001/03/22 18:42:06 eschnett Exp $ +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 $ #include "cctk.h" @@ -48,6 +48,12 @@ c bbox(:,3) is stride $ .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 + if (regbbox(d,1).lt.srcbbox(d,1) + $ .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") + end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par index 94558c2de..db3d76f2c 100644 --- a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par +++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par @@ -7,7 +7,7 @@ # @enddesc # @@*/ # -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par,v 1.4 2001/03/17 00:36:14 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl2.par,v 1.5 2001/03/24 22:38:49 eschnett Exp $ ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetIOFlexIO CarpetLib CarpetSlab IDScalarWave WaveToyF77" @@ -20,6 +20,8 @@ driver::global_ny = 24 driver::global_nz = 32 Carpet::max_refinement_levels = 2 +Carpet::prolongation_order_space= 3 +Carpet::prolongation_order_time = 2 grid::type = byrange grid::xmin = -4.8 diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par index dd4f38ab7..7a3cfa981 100644 --- a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par +++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par @@ -7,7 +7,7 @@ # @enddesc # @@*/ # -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par,v 1.4 2001/03/17 00:36:14 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_coarse_rl3.par,v 1.5 2001/03/24 22:38:49 eschnett Exp $ ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetIOFlexIO CarpetLib CarpetSlab IDScalarWave WaveToyF77" @@ -20,6 +20,7 @@ driver::global_ny = 24 driver::global_nz = 32 Carpet::max_refinement_levels = 3 +Carpet::prolongation_order = 2 grid::type = byrange grid::xmin = -4.8 diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par index f8d868c5a..b8141e51c 100644 --- a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par +++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par @@ -7,7 +7,7 @@ # @enddesc # @@*/ # -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par,v 1.5 2001/03/17 00:36:14 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl2.par,v 1.6 2001/03/24 22:38:49 eschnett Exp $ ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetIOFlexIO CarpetLib CarpetSlab IDScalarWave WaveToyF77" @@ -20,6 +20,8 @@ driver::global_ny = 46 driver::global_nz = 62 Carpet::max_refinement_levels = 2 +Carpet::prolongation_order_space= 3 +Carpet::prolongation_order_time = 2 grid::type = byrange grid::xmin = -4.8 diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par index 85f145f60..eee5ba0d7 100644 --- a/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par +++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par @@ -7,7 +7,7 @@ # @enddesc # @@*/ # -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par,v 1.4 2001/03/17 00:36:14 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_periodic_rl3.par,v 1.5 2001/03/24 22:38:49 eschnett Exp $ ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetIOFlexIO CarpetLib CarpetSlab IDScalarWave WaveToyF77" @@ -20,6 +20,7 @@ driver::global_ny = 46 driver::global_nz = 62 Carpet::max_refinement_levels = 3 +Carpet::prolongation_order = 2 grid::type = byrange grid::xmin = -4.8 |