diff options
Diffstat (limited to 'Carpet/CarpetLib')
30 files changed, 950 insertions, 523 deletions
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc index de281f582..f90c14d4e 100644 --- a/Carpet/CarpetLib/src/bbox.cc +++ b/Carpet/CarpetLib/src/bbox.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/bbox.cc,v 1.4 2001/03/12 16:54:25 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.5 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include "defs.hh" @@ -28,6 +29,8 @@ # include "bbox.hh" #endif +using namespace std; + // Constructors @@ -223,19 +226,37 @@ bbox<T,D>::iterator bbox<T,D>::end () const { // Output +#ifndef SGI +// This doesn't work on SGIs. Is this legal C++? template<class T,int D> ostream& operator<< (ostream& os, const bbox<T,D>& b) { os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")"; return os; } +#else +ostream& operator<< (ostream& os, const bbox<int,1>& b) { + os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")"; + return os; +} +ostream& operator<< (ostream& os, const bbox<int,2>& b) { + os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")"; + return os; +} +ostream& operator<< (ostream& os, const bbox<int,3>& b) { + os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")"; + return os; +} +#endif #if defined(TMPL_EXPLICIT) template class bbox<int,1>; template ostream& operator<< (ostream& os, const bbox<int,1>& b); + template class bbox<int,2>; template ostream& operator<< (ostream& os, const bbox<int,2>& b); + template class bbox<int,3>; template ostream& operator<< (ostream& os, const bbox<int,3>& b); #endif diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index db43519b4..396a9c5ca 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.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/bbox.hh,v 1.6 2001/03/14 11:00:26 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.7 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -26,13 +26,15 @@ #include "defs.hh" #include "vect.hh" +using namespace std; + // Forward definition template<class T, int D> class bbox; // Output -template<class T,int D> +template<class T, int D> ostream& operator<< (ostream& os, const bbox<T,D>& b); @@ -102,7 +104,7 @@ public: // Iterators class iterator { protected: - bbox box; + const bbox& box; vect<T,D> pos; public: iterator (const bbox& box, const vect<T,D>& pos); @@ -115,7 +117,7 @@ public: iterator end () const; // Output - friend ostream& operator<< <>(ostream& os, const bbox& b); + friend ostream& operator<< <>(ostream& os, const bbox& s); }; diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 489856298..f18385c1b 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.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/bboxset.cc,v 1.4 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.5 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include <set> @@ -28,6 +29,8 @@ # include "bboxset.hh" #endif +using namespace std; + // Constructors @@ -210,6 +213,8 @@ bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s) { // Difference +#ifndef SGI +// This doesn't work on SGIs. Is this legal C++? template<class T, int D> bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) { assert (b1.is_aligned_with(b2)); @@ -243,6 +248,104 @@ bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) { assert (r.invariant()); return r; } +#else +bboxset<int,1> operator- (const bbox<int,1>& b1, const bbox<int,1>& b2) { + assert (b1.is_aligned_with(b2)); + if (b1.empty()) return bboxset<int,1>(); + if (b2.empty()) return bboxset<int,1>(b1); + const vect<int,1> str = b1.stride(); + bboxset<int,1> r; + for (int d=0; d<1; ++d) { + // make resulting bboxes as large as possible in x-direction (for + // better consumption by Fortranly ordered arrays) + vect<int,1> lb, ub; + bbox<int,1> b; + for (int dd=0; dd<1; ++dd) { + if (dd<d) { + lb[dd] = b2.lower()[dd]; + ub[dd] = b2.upper()[dd]; + } else if (dd>d) { + lb[dd] = b1.lower()[dd]; + ub[dd] = b1.upper()[dd]; + } + } + lb[d] = b1.lower()[d]; + ub[d] = b2.lower()[d] - str[d]; + b = bbox<int,1>(lb,ub,str) & b1; + r += b; + lb[d] = b2.upper()[d] + str[d]; + ub[d] = b1.upper()[d]; + b = bbox<int,1>(lb,ub,str) & b1; + r += b; + } + assert (r.invariant()); + return r; +} +bboxset<int,2> operator- (const bbox<int,2>& b1, const bbox<int,2>& b2) { + assert (b1.is_aligned_with(b2)); + if (b1.empty()) return bboxset<int,2>(); + if (b2.empty()) return bboxset<int,2>(b1); + const vect<int,2> str = b1.stride(); + bboxset<int,2> r; + for (int d=0; d<2; ++d) { + // make resulting bboxes as large as possible in x-direction (for + // better consumption by Fortranly ordered arrays) + vect<int,2> lb, ub; + bbox<int,2> b; + for (int dd=0; dd<2; ++dd) { + if (dd<d) { + lb[dd] = b2.lower()[dd]; + ub[dd] = b2.upper()[dd]; + } else if (dd>d) { + lb[dd] = b1.lower()[dd]; + ub[dd] = b1.upper()[dd]; + } + } + lb[d] = b1.lower()[d]; + ub[d] = b2.lower()[d] - str[d]; + b = bbox<int,2>(lb,ub,str) & b1; + r += b; + lb[d] = b2.upper()[d] + str[d]; + ub[d] = b1.upper()[d]; + b = bbox<int,2>(lb,ub,str) & b1; + r += b; + } + assert (r.invariant()); + return r; +} +bboxset<int,3> operator- (const bbox<int,3>& b1, const bbox<int,3>& b2) { + assert (b1.is_aligned_with(b2)); + if (b1.empty()) return bboxset<int,3>(); + if (b2.empty()) return bboxset<int,3>(b1); + const vect<int,3> str = b1.stride(); + bboxset<int,3> r; + for (int d=0; d<3; ++d) { + // make resulting bboxes as large as possible in x-direction (for + // better consumption by Fortranly ordered arrays) + vect<int,3> lb, ub; + bbox<int,3> b; + for (int dd=0; dd<3; ++dd) { + if (dd<d) { + lb[dd] = b2.lower()[dd]; + ub[dd] = b2.upper()[dd]; + } else if (dd>d) { + lb[dd] = b1.lower()[dd]; + ub[dd] = b1.upper()[dd]; + } + } + lb[d] = b1.lower()[d]; + ub[d] = b2.lower()[d] - str[d]; + b = bbox<int,3>(lb,ub,str) & b1; + r += b; + lb[d] = b2.upper()[d] + str[d]; + ub[d] = b1.upper()[d]; + b = bbox<int,3>(lb,ub,str) & b1; + r += b; + } + assert (r.invariant()); + return r; +} +#endif template<class T, int D> bboxset<T,D> bboxset<T,D>::operator- (const box& b) const { @@ -281,21 +384,38 @@ bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const { return r; } +#ifndef SGI +// This doesn't work on SGIs. Is this legal C++? template<class T, int D> bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) { bboxset<T,D> r = bboxset<T,D>(b) - s; assert (r.invariant()); return r; } +#else +bboxset<int,1> operator- (const bbox<int,1>& b, const bboxset<int,1>& s) { + bboxset<int,1> r = bboxset<int,1>(b) - s; + assert (r.invariant()); + return r; +} +bboxset<int,2> operator- (const bbox<int,2>& b, const bboxset<int,2>& s) { + bboxset<int,2> r = bboxset<int,2>(b) - s; + assert (r.invariant()); + return r; +} +bboxset<int,3> operator- (const bbox<int,3>& b, const bboxset<int,3>& s) { + bboxset<int,3> r = bboxset<int,3>(b) - s; + assert (r.invariant()); + return r; +} +#endif // Output template<class T,int D> ostream& operator<< (ostream& os, const bboxset<T,D>& s) { -// os << "bboxset<" STR(T) "," << D << ">:size=" << s.size() << "," -// << "set=" << s.bs; - os << s.bs; + os << "bboxset<T," << D << ">:size=" << s.size() << "," << "set=" << s.bs; return os; } @@ -303,7 +423,11 @@ ostream& operator<< (ostream& os, const bboxset<T,D>& s) { #if defined(TMPL_EXPLICIT) template class bboxset<int,3>; -template bboxset<int,3> operator- (const bbox<int,3>& b1, const bbox<int,3>& b3); -template bboxset<int,3> operator- (const bbox<int,3>& b, const bboxset<int,3>& s); -template ostream& operator<< (ostream& os, const bboxset<int,3>& b); + +template +bboxset<int,3> operator- (const bbox<int,3>& b1, const bbox<int,3>& b3); +template +bboxset<int,3> operator- (const bbox<int,3>& b, const bboxset<int,3>& s); +template +ostream& operator<< (ostream& os, const bboxset<int,3>& b); #endif diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index dae3d3f3f..4883cc6e9 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.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/bboxset.hh,v 1.3 2001/03/12 16:54:25 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -21,7 +21,8 @@ #ifndef BBOXSET_HH #define BBOXSET_HH -#include <cassert> +#include <assert.h> + #include <iostream> #include <set> @@ -29,6 +30,8 @@ #include "defs.hh" #include "vect.hh" +using namespace std; + // Forward definition @@ -105,8 +108,8 @@ public: // friend bboxset operator- <T,D>(const box& b, const bboxset& s); // Iterators - typedef bset::const_iterator const_iterator; - typedef bset::iterator iterator; + typedef typename bset::const_iterator const_iterator; + typedef typename bset::iterator iterator; iterator begin () const { return bs.begin(); } iterator end () const { return bs.end(); } diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77 index 716ffd135..ae5d1c34f 100644 --- a/Carpet/CarpetLib/src/copy_3d_real8.F77 +++ b/Carpet/CarpetLib/src/copy_3d_real8.F77 @@ -1,10 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.8 2004/03/11 12:03:09 schnetter Exp $ +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 $ #include "cctk.h" -#include "cctk_Parameters.h" - - subroutine copy_3d_real8 ( $ src, srciext, srcjext, srckext, @@ -13,8 +10,6 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8 implicit none - DECLARE_CCTK_PARAMETERS - integer srciext, srcjext, srckext CCTK_REAL8 src(srciext,srcjext,srckext) integer dstiext, dstjext, dstkext @@ -32,8 +27,6 @@ c bbox(:,3) is stride integer i, j, k integer d - character msg*1000 - do d=1,3 @@ -46,25 +39,10 @@ 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 - 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") - end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -93,19 +71,12 @@ c This could be handled, but is likely to point to an error elsewhere c Loop over region - do k = 1, regkext - do j = 1, regjext - do i = 1, regiext - - if (check_array_accesses.ne.0) then - call checkindex (srcioff+i, srcjoff+j+1, srckoff+k+1, 1,1,1, - $ "source") - call checkindex (dstioff+i, dstjoff+j+1, dstkoff+k+1, 1,1,1, - $ "destination") - end if + do k = 0, regkext-1 + do j = 0, regjext-1 + do i = 0, regiext-1 - dst (dstioff+i, dstjoff+j, dstkoff+k) - $ = src (srcioff+i, srcjoff+j, srckoff+k) + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) + $ = src (srcioff+i+1, srcjoff+j+1, srckoff+k+1) end do end do diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 932b692e6..23cd28b57 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.7 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.8 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,12 +18,15 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <fstream> #include <string> #include <mpi.h> +#include "cctk.h" + #include "bbox.hh" #include "defs.hh" #include "dist.hh" @@ -33,6 +36,8 @@ # include "data.hh" #endif +using namespace std; + // Constructors @@ -65,7 +70,7 @@ data<T,D>* data<T,D>::make_typed () const { // Storage management template<class T, int D> void data<T,D>::allocate (const ibbox& extent, const int proc, - void* const mem=0) { + void* const mem) { assert (!_has_storage); _has_storage = true; // data @@ -110,7 +115,7 @@ void data<T,D>::transfer_from (generic_data<D>* gsrc) { // Processor management template<class T, int D> -void data<T,D>::change_processor (const int newproc, void* const mem=0) { +void data<T,D>::change_processor (const int newproc, void* const mem) { if (newproc == _proc) { assert (!mem); return; @@ -162,7 +167,9 @@ void data<T,D>::change_processor (const int newproc, void* const mem=0) { // Data manipulators template<class T, int D> -void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) { +void data<T,D> +::copy_from_innerloop (const generic_data<D>* gsrc, const ibbox& box) +{ const data* src = (const data*)gsrc; assert (has_storage() && src->has_storage()); assert (all(box.lower()>=extent().lower() @@ -174,164 +181,371 @@ void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) { assert (all((box.lower()-extent().lower())%box.stride() == 0 && (box.lower()-src->extent().lower())%box.stride() == 0)); - if (proc() == src->proc()) { - // copy on same processor + assert (proc() == src->proc()); + + int rank; + MPI_Comm_rank (dist::comm, &rank); + assert (rank == proc()); + + for (ibbox::iterator it=box.begin(); it!=box.end(); ++it) { + const ivect index = *it; + (*this)[index] = (*src)[index]; + } + +} + + + +template<class T, int D> +void data<T,D> +::interpolate_from_innerloop (const generic_data<D>* 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()); + + for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { + const ivect& pos = *posi; - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - - for (ibbox::iterator it=box.begin(); it!=box.end(); ++it) { - const ivect index = *it; - (*this)[index] = (*src)[index]; + // 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; - } else { + } // 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 (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)); + + for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { + const ivect& pos = *posi; - // copy to different processor - data* const tmp = new data(box, src->proc()); - tmp->copy_from (src, box); - tmp->change_processor (proc()); - copy_from (tmp, box); - delete tmp; + // get box around current position + const ibbox frombox + = ibbox(pos,pos,extent().stride()).expanded_for(src1->extent()); - } + // interpolate from box to position + T src1_sum = 0; + T src2_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]; + } + (*this)[pos] = src1_fac * src1_sum + src2_fac * src2_sum; + + } // for pos + } -template<class T, int D> -void data<T,D>::interpolate_from (const generic_data<D>* gsrc, - const ibbox& box) + + +extern "C" { + void CCTK_FCALL CCTK_FNAME(copy_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]); +} + +template<> +void data<CCTK_REAL8,3> +::copy_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() - /* && box.stride()<=src->extent().stride() */ )); + && box.stride()==src->extent().stride())); assert (all((box.lower()-extent().lower())%box.stride() == 0 - /* && (box.lower()-src->extent().lower())%box.stride() == 0 */ )); + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + assert (proc() == src->proc()); - if (proc() == src->proc()) { - // interpolate on same processor + 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]; - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (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 - - } + srcbbox[0][d] = sext.lower()[d]; + srcbbox[1][d] = sext.upper()[d]; + srcbbox[2][d] = sext.stride()[d]; - } else { - // interpolate from other processor + dstbbox[0][d] = dext.lower()[d]; + dstbbox[1][d] = dext.upper()[d]; + dstbbox[2][d] = dext.stride()[d]; - data* const tmp = new data(box, src->proc()); - tmp->interpolate_from (src, box); - tmp->change_processor (proc()); - copy_from (tmp, box); - delete tmp; + 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(copy_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(); } } -template<class T, int D> -void data<T,D>::interpolate_from (const generic_data<D>* gsrc, - const double sfact, - const generic_data<D>* gtrc, - const double tfact, - const ibbox& box) + + +extern "C" { + void CCTK_FCALL CCTK_FNAME(prolongate_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(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]); +} + +template<> +void data<CCTK_REAL8,3> +::interpolate_from_innerloop (const generic_data<3>* gsrc, const ibbox& box) { const data* src = (const data*)gsrc; - const data* trc = (const data*)gtrc; - assert (has_storage() && src->has_storage() && trc->has_storage()); + 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() - && box.lower()>=trc->extent().lower())); + && box.lower()>=src->extent().lower())); assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper() - && box.upper()<=trc->extent().upper())); - assert (all(box.stride()==extent().stride() - /* && box.stride()<=src->extent().stride() */ - /* && trc->extent().stride()==src->extent().stride() */ )); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0 - && (box.lower()-trc->extent().lower())%box.stride() == 0)); - assert (src->proc() == trc->proc()); + && box.upper()<=src->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 + 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]; - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (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 src_sum = 0; - T trc_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)); - src_sum += f * (*src)[frompos]; - trc_sum += f * (*trc)[frompos]; - } - (*this)[pos] = (T)sfact * src_sum + (T)tfact * trc_sum; - - } // for pos - - } + 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 { - // interpolate from other processor + 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, + 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) +{ + 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 (proc() == src1->proc()); + + int rank; + MPI_Comm_rank (dist::comm, &rank); + assert (rank == proc()); + + assert (src1->extent() == src2->extent()); + + const ibbox& sext = src1->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]; - data* const tmp = new data(box, src->proc()); - tmp->interpolate_from (src, sfact, trc, tfact, box); - tmp->change_processor (proc()); - copy_from (tmp, box); - delete tmp; + 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_2tl) + ((const CCTK_REAL8*)src1->storage(), t1, + (const CCTK_REAL8*)src2->storage(), t2, + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), t, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, + dstbbox, + regbbox); + } else { + abort(); } } @@ -361,7 +575,7 @@ ostream& data<T,D>::out (ostream& os) const { #if defined(TMPL_EXPLICIT) #define INSTANTIATE(T) \ -template data<T,3>; +template class data<T,3>; #include "instantiate" diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index 84e4c92b7..c4594c56d 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.3 2001/03/10 20:55:06 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -21,7 +21,8 @@ #ifndef DATA_HH #define DATA_HH -#include <cassert> +#include <assert.h> + #include <string> #include "defs.hh" @@ -30,6 +31,8 @@ #include "gdata.hh" #include "vect.hh" +using namespace std; + // Forward definition @@ -50,7 +53,7 @@ class data: public generic_data<D> { typedef bbox<int,D> ibbox; // Fields - T* restrict _storage; // the data (if located on this processor) + T* _storage; // the data (if located on this processor) public: @@ -74,12 +77,12 @@ public: virtual void change_processor (const int newproc, void* const mem=0); // Accessors - virtual const T* storage () const { + virtual const void* storage () const { assert (_has_storage); return _storage; } - virtual T* storage () { + virtual void* storage () { assert (_has_storage); return _storage; } @@ -96,14 +99,13 @@ public: } // Data manipulators - virtual void copy_from (const generic_data<D>* gsrc, const ibbox& b); - virtual void interpolate_from (const generic_data<D>* gsrc, - const ibbox& box); - virtual void interpolate_from (const generic_data<D>* gsrc, - const double sfact, - const generic_data<D>* gtrc, - const double tfact, - const ibbox& box); + void copy_from_innerloop (const generic_data<D>* src, + 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 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 bbfdf990a..bf483c62f 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.2 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.3 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include <list> #include <set> @@ -28,6 +29,8 @@ # include "defs.hh" #endif +using namespace std; + // List output diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh index 9071bad2b..ec59237f8 100644 --- a/Carpet/CarpetLib/src/defs.hh +++ b/Carpet/CarpetLib/src/defs.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/defs.hh,v 1.3 2001/03/18 05:20:24 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -31,6 +31,8 @@ #include <set> #include <vector> +using namespace std; + // Stringification #define STR(s) #s diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 09962b96d..c9e03753d 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.9 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.10 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -19,7 +19,7 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> #include "defs.hh" #include "dist.hh" @@ -32,6 +32,8 @@ #undef DEBUG_OUTPUT +using namespace std; + // Constructors diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index 91634d1ca..aa9e331da 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.3 2001/03/10 20:55:06 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -22,7 +22,8 @@ #ifndef DH_HH #define DH_HH -#include <cassert> +#include <assert.h> + #include <iostream> #include <list> #include <string> @@ -34,6 +35,8 @@ #include "gh.hh" #include "vect.hh" +using namespace std; + // Forward declaration @@ -57,6 +60,7 @@ class dh { typedef list<ibbox> iblist; typedef vector<iblist> iblistvect; // vector of lists +public: struct dboxes { ibbox exterior; // whole region (including boundaries) @@ -77,6 +81,7 @@ class dh { iblistvect recv_ref_bnd_coarse; // received from coarser grids ibset sync_not; // not received while syncing (outer boundary) }; +private: struct dbases { ibbox exterior; // whole region (including boundaries) diff --git a/Carpet/CarpetLib/src/dist.cc b/Carpet/CarpetLib/src/dist.cc index 9aa613764..99cc2e83a 100644 --- a/Carpet/CarpetLib/src/dist.cc +++ b/Carpet/CarpetLib/src/dist.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/dist.cc,v 1.3 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.cc,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,7 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> #include <mpi.h> @@ -31,6 +31,8 @@ # include "dist.hh" #endif +using namespace std; + namespace dist { diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh index dac83076f..9c488e0fd 100644 --- a/Carpet/CarpetLib/src/dist.hh +++ b/Carpet/CarpetLib/src/dist.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/dist.hh,v 1.3 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -21,9 +21,9 @@ #ifndef DIST_HH #define DIST_HH -#include <cassert> -#include <cstdio> -#include <cstdlib> +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> #include <complex> @@ -31,6 +31,8 @@ #include "defs.hh" +using namespace std; + namespace dist { diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index eb19701ce..5418db23a 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.8 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.9 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <fstream> #include <iomanip> @@ -31,6 +32,8 @@ # include "gdata.hh" #endif +using namespace std; + // Constructors @@ -45,6 +48,127 @@ generic_data<D>::~generic_data () { } +// Data manipulators +template<int D> +void generic_data<D>::copy_from (const generic_data* src, const ibbox& box) +{ + assert (has_storage() && src->has_storage()); + 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() + && box.stride()==src->extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0 + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + if (proc() == src->proc()) { + // copy on same processor + + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == proc()) { + copy_from_innerloop (src, box); + } + + } else { + + // copy to different processor + generic_data* const tmp = make_typed(); + tmp->allocate (box, src->proc()); + tmp->copy_from (src, box); + tmp->change_processor (proc()); + copy_from (tmp, box); + delete tmp; + + } +} + + + +template<int D> +void generic_data<D> +::interpolate_from (const generic_data* src, const ibbox& box) +{ + 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)); + + 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; + + } +} + + + +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()) { + // interpolate on same processor + + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == proc()) { + interpolate_from_innerloop (src1, t1, src2, t2, box, t); + } + + } 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->change_processor (proc()); + copy_from (tmp, box); + delete tmp; + + } +} + + + // Output template<int D> template<int DD> @@ -75,7 +199,8 @@ void generic_data<D>::write_ascii (const string name, const int time, << " component " << c << " multigrid level " << ml << endl << "# column format: it tl rl c ml"; assert (D>=1 && D<=3); - for (int d=0; d<D; ++d) file << " " << "xyz"[d]; + const char* const coords = "xyz"; + for (int d=0; d<D; ++d) file << " " << coords[d]; file << " data" << endl; const vect<int,DD> lo = extent().lower()[dirs]; diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh index e765cc7c6..1a16f92c2 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.4 2001/03/14 11:00:26 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.5 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -21,8 +21,9 @@ #ifndef GDATA_HH #define GDATA_HH -#include <cassert> -#include <cstdlib> +#include <assert.h> +#include <stdlib.h> + #include <fstream> #include <iostream> #include <string> @@ -32,6 +33,8 @@ #include "bbox.hh" #include "vect.hh" +using namespace std; + // Forward declaration @@ -131,13 +134,22 @@ public: } // Data manipulators - virtual void copy_from (const generic_data* src, - const ibbox& b) = 0; - virtual void interpolate_from (const generic_data* src, - const ibbox& box) = 0; - virtual void interpolate_from (const generic_data* src, const double sfact, - const generic_data* trc, const double tfact, - const ibbox& box) = 0; + 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); +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; + +public: // Output template<int DD> diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc index 86a2d518d..6847d8506 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.2 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.3 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -19,7 +19,7 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> #include "defs.hh" @@ -27,6 +27,8 @@ # include "gf.hh" #endif +using namespace std; + // Constructors @@ -81,6 +83,7 @@ ostream& gf<T,D>::out (ostream& os) const { template class gf<T,3>; #include "instantiate" + #undef INSTANTIATE #endif diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh index 298c26589..bf916bb82 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.1 2001/03/01 13:40:10 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.2 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -22,8 +22,9 @@ #ifndef GF_HH #define GF_HH -#include <cassert> -#include <cmath> +#include <assert.h> +#include <math.h> + #include <iostream> #include <string> @@ -36,6 +37,8 @@ #include "th.hh" #include "vect.hh" +using namespace std; + // A real grid function diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 53aa15179..09966ac69 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.5 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.6 2001/03/22 18:42:05 eschnett Exp $ ***************************************************************************/ @@ -19,7 +19,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include <string> @@ -31,6 +32,8 @@ # include "ggf.hh" #endif +using namespace std; + // Constructors @@ -150,6 +153,19 @@ void generic_gf<D>::recompose () { } // for tl } +// Cycle the time levels by rotating the data sets +template<int D> +void generic_gf<D>::cycle (int rl, int c, int ml) { + assert (rl>=0 && rl<h.reflevels()); + assert (c>=0 && c<h.components(rl)); + assert (ml>=0 && ml<h.mglevels(rl,c)); + generic_data<D>* tmpdata = storage[tmin-tmin][rl][c][ml]; + for (int tl=tmin; tl<=tmax-1; ++tl) { + storage[tl-tmin][rl][c][ml] = storage[tl+1-tmin][rl][c][ml]; + } + storage[tmax-tmin][rl][c][ml] = tmpdata; +} + // Operations @@ -240,17 +256,15 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1, template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, + int tl2a, int tl2b, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list) { 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 (tl2>=tmin && tl2<=tmax); - assert (tl3>=tmin && tl3<=tmax); + assert (tl2a>=tmin && tl2a<=tmax); + assert (tl2b>=tmin && tl2b<=tmax); assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; assert (ml2<h.mglevels(rl2,c2)); @@ -260,25 +274,23 @@ 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 - (storage[tl2-tmin][rl2][c2][ml2], fact2, - storage[tl3-tmin][rl2][c2][ml2], fact3, recv); + (storage[tl2a-tmin][rl2][c2][ml2], tl2a, + storage[tl2b-tmin][rl2][c2][ml2], tl2b, recv, tl1); } // Interpolate a component (between multigrid levels) template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, + int tl2a, int tl2b, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list) { 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 (tl2>=tmin && tl2<=tmax); - assert (tl3>=tmin && tl3<=tmax); + assert (tl2a>=tmin && tl2a<=tmax); + assert (tl2b>=tmin && tl2b<=tmax); assert (rl2>=0 && rl2<h.reflevels()); const int c2=c1; assert (ml2<h.mglevels(rl2,c2)); @@ -291,8 +303,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 - (storage[tl2-tmin][rl2][c2][ml2], fact2, - storage[tl3-tmin][rl2][c2][ml2], fact3, *r); + (storage[tl2a-tmin][rl2][c2][ml2], tl2a, + storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1); } } @@ -300,18 +312,16 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, template<int D> void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, + int tl2a, int tl2b, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect) { 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 (tl2>=tmin && tl2<=tmax); + assert (tl2a>=tmin && tl2a<=tmax); + assert (tl2b>=tmin && tl2b<=tmax); assert (rl2>=0 && rl2<h.reflevels()); - assert (tl3>=tmin && tl3<=tmax); // walk all components for (int c2=0; c2<h.components(rl2); ++c2) { assert (ml2<h.mglevels(rl2,c2)); @@ -324,8 +334,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 - (storage[tl2-tmin][rl2][c2][ml2], fact2, - storage[tl3-tmin][rl2][c2][ml2], fact3, *r); + (storage[tl2a-tmin][rl2][c2][ml2], tl2a, + storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1); } } } @@ -361,13 +371,11 @@ void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml) { copycat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, tl2,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine); } else { - int tl3=tl2+1; - assert (tl3>=tmin && tl3<=tmax); - const double fact2 = 1 - (time - tl2); - const double fact3 = 1 - fact2; + 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, - tl2,fact2, tl3,fact3, - rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine); + tl2a, tl2b, rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine); } } diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index f59d4e415..6e998c145 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.2 2001/03/15 09:59:43 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -22,7 +22,8 @@ #ifndef GGF_HH #define GGF_HH -#include <cassert> +#include <assert.h> + #include <iostream> #include <string> @@ -32,6 +33,8 @@ #include "gh.hh" #include "th.hh" +using namespace std; + // Forward declaration @@ -83,12 +86,15 @@ public: virtual ~generic_gf (); // Comparison - virtual bool operator== (const generic_gf<D>& f) const; + bool operator== (const generic_gf<D>& f) const; // Modifiers - virtual void recompose (); + void recompose (); + + // Cycle the time levels by rotating the data sets + void cycle (int rl, int c, int ml); @@ -105,46 +111,40 @@ protected: protected: // Copy region for a component (between time levels) - virtual 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); + 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) - virtual 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); + 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) - virtual 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); + 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) - virtual void intercat (int tl1, int rl1, int c1, int ml1, - const ibbox dh<D>::dboxes::* recv_list, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, - const ibbox dh<D>::dboxes::* send_list); + 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); // Interpolate a component (between multigrid levels) - virtual void intercat (int tl1, int rl1, int c1, int ml1, - const iblist dh<D>::dboxes::* recv_list, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, - const iblist dh<D>::dboxes::* send_list); + 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); // Interpolate a level (between refinement levels) - virtual void intercat (int tl1, int rl1, int c1, int ml1, - const iblistvect dh<D>::dboxes::* recv_listvect, - int tl2, const double fact2, - int tl3, const double fact3, - int rl2, int ml2, - const iblistvect dh<D>::dboxes::* send_listvect); + 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); @@ -157,25 +157,25 @@ public: // synchronised. They don't need to be prolongated. // Copy a component from the next time level - virtual void copy (int tl, int rl, int c, int ml); + void copy (int tl, int rl, int c, int ml); // Synchronise the boundaries of a component - virtual void sync (int tl, int rl, int c, int ml); + void sync (int tl, int rl, int c, int ml); // Prolongate the boundaries of a component - virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml); + void ref_bnd_prolongate (int tl, int rl, int c, int ml); // Restrict a multigrid level - virtual void mg_restrict (int tl, int rl, int c, int ml); + void mg_restrict (int tl, int rl, int c, int ml); // Prolongate a multigrid level - virtual void mg_prolongate (int tl, int rl, int c, int ml); + void mg_prolongate (int tl, int rl, int c, int ml); // Restrict a refinement level - virtual void ref_restrict (int tl, int rl, int c, int ml); + void ref_restrict (int tl, int rl, int c, int ml); // Prolongate a refinement level - virtual void ref_prolongate (int tl, int rl, int c, int ml); + void ref_prolongate (int tl, int rl, int c, int ml); diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index b1608b6b7..ee10251b1 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -7,7 +7,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.4 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.5 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -20,8 +20,8 @@ * * ***************************************************************************/ -#include <cassert> -#include <cstdlib> +#include <assert.h> +#include <stdlib.h> #include <iostream> #include "defs.hh" @@ -32,6 +32,8 @@ # include "gh.hh" #endif +using namespace std; + // Constructors diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh index ff90a9476..841abbdfb 100644 --- a/Carpet/CarpetLib/src/gh.hh +++ b/Carpet/CarpetLib/src/gh.hh @@ -7,7 +7,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.2 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -23,7 +23,8 @@ #ifndef GH_HH #define GH_HH -#include <cassert> +#include <assert.h> + #include <iostream> #include <list> #include <vector> @@ -33,6 +34,8 @@ #include "dist.hh" #include "vect.hh" +using namespace std; + // Forward declaration diff --git a/Carpet/CarpetLib/src/instantiate b/Carpet/CarpetLib/src/instantiate index aafc645d2..f85d27c58 100644 --- a/Carpet/CarpetLib/src/instantiate +++ b/Carpet/CarpetLib/src/instantiate @@ -1,7 +1,7 @@ // Instantiate templates for all available types -*-C++-*- // (C) 2001 Erik Schnetter <schnetter@uni-tuebingen.de> -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.2 2001/03/19 21:30:19 eschnett Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.3 2001/03/22 18:42:06 eschnett Exp $ // Usage: // Define the macro INSTANTIATE(T) to instantiate for the type T, @@ -55,7 +55,9 @@ #endif #if !defined(CARPET_BYTE) && !defined(CARPET_INT2) && !defined(CARPET_INT4) && !defined(CARPET_INT8) && !defined(CARPET_REAL4) && !defined(CARPET_REAL8) && !defined(CARPET_REAL16) && !defined(CARPET_COMPLEX8) && !defined(CARPET_COMPLEX16) && !defined(CARPET_COMPLEX32) -// Assume the user just wants REAL +// Assume the user just wants INT and REAL +# undef CARPET_INT +# define CARPET_INT # undef CARPET_REAL # define CARPET_REAL #endif diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 121d29247..dd1756780 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -1,8 +1,9 @@ # Main make.code.defn file for thorn CarpetLib -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.2 2001/03/17 22:37:56 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.3 2001/03/22 18:42:06 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 +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 # 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 8eae50332..04053c524 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 @@ -1,21 +1,8 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.11 2004/03/11 12:03:09 schnetter Exp $ +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 $ #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 ( $ src, srciext, srcjext, srckext, $ dst, dstiext, dstjext, dstkext, @@ -23,9 +10,6 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d implicit none - CCTK_REAL8 one - parameter (one = 1) - integer srciext, srcjext, srckext CCTK_REAL8 src(srciext,srcjext,srckext) integer dstiext, dstjext, dstkext @@ -42,18 +26,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(2), jfac(2), kfac(2) integer ii, jj, kk integer fac CCTK_REAL8 res integer d - character msg*1000 - do d=1,3 @@ -69,25 +49,10 @@ 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") - end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -120,71 +85,52 @@ c This could be handled, but is likely to point to an error elsewhere 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 = ifac(ii) * jfac(jj) * kfac(kk) + fac = 1 + + 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 * src(i0+ii, j0+jj, k0+kk) + res = res + fac * src(i0+ii,j0+jj,k0+kk) end if end do end do end do -c$$$ fac = ifac(1) * jfac(1) * kfac(1) -c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1) -c$$$ -c$$$ fac = ifac(2) * jfac(1) * kfac(1) -c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1) -c$$$ -c$$$ fac = ifac(1) * jfac(2) * kfac(1) -c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1) -c$$$ -c$$$ fac = ifac(2) * jfac(2) * kfac(1) -c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1) -c$$$ -c$$$ fac = ifac(1) * jfac(1) * kfac(2) -c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2) -c$$$ -c$$$ fac = ifac(2) * jfac(1) * kfac(2) -c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2) -c$$$ -c$$$ fac = ifac(1) * jfac(2) * kfac(2) -c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2) -c$$$ -c$$$ fac = ifac(2) * jfac(2) * kfac(2) -c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2) - - 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_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 index 8fd69178f..f9b380c5d 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.12 2004/03/11 12:03:09 schnetter Exp $ +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 $ #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 ( $ src1, t1, src2, t2, srciext, srcjext, srckext, @@ -26,17 +13,14 @@ 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 @@ -51,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(2), jfac(2), kfac(2) integer ii, jj, kk integer fac CCTK_REAL8 res integer d - character msg*1000 - do d=1,3 @@ -78,25 +58,10 @@ 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") - end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -128,63 +93,60 @@ c This could be handled, but is likely to point to an error elsewhere -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 / (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 = ifac(ii) * jfac(jj) * kfac(kk) + fac = 1 + + 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 * 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 / (dstifac * dstjfac * dstkfac) end do end do diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77 index 81f4cfd0a..68ea98b11 100644 --- a/Carpet/CarpetLib/src/restrict_3d_real8.F77 +++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77 @@ -1,20 +1,7 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.7 2004/03/11 12:03:09 schnetter Exp $ +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 $ #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 restrict_3d_real8 ( $ src, srciext, srcjext, srckext, @@ -42,8 +29,6 @@ c bbox(:,3) is stride integer i, j, k integer d - character msg*1000 - do d=1,3 @@ -59,25 +44,10 @@ 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 - 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") - end if end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -114,10 +84,6 @@ c Loop over coarse region do j = 0, regjext-1 do i = 0, regiext-1 - CHKIDX (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1, \ - srciext,srcjext,srckext, "source") - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) $ = src (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1) diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc index aab64bc78..17aaff924 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.2 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.3 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -19,7 +19,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include "defs.hh" @@ -29,6 +30,8 @@ # include "th.hh" #endif +using namespace std; + // Constructors diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh index 6089e37cd..7c39a7d18 100644 --- a/Carpet/CarpetLib/src/th.hh +++ b/Carpet/CarpetLib/src/th.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/th.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.2 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -22,13 +22,16 @@ #ifndef TH_HH #define TH_HH -#include <cassert> +#include <assert.h> + #include <iostream> #include <vector> #include "defs.hh" #include "gh.hh" +using namespace std; + // Forward declaration diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc index 6c8d012a2..5a622fd1c 100644 --- a/Carpet/CarpetLib/src/vect.cc +++ b/Carpet/CarpetLib/src/vect.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/vect.cc,v 1.2 2001/03/19 21:30:19 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.3 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -18,7 +18,8 @@ * * ***************************************************************************/ -#include <cassert> +#include <assert.h> + #include <iostream> #include "defs.hh" @@ -27,9 +28,13 @@ # include "vect.hh" #endif +using namespace std; + // Output +#ifndef SGI +// This doesn't work on SGIs. Is this legal C++? template<class T,int D> ostream& operator<< (ostream& os, const vect<T,D>& a) { os << "["; @@ -40,6 +45,35 @@ ostream& operator<< (ostream& os, const vect<T,D>& a) { os << "]"; return os; } +#else +ostream& operator<< (ostream& os, const vect<int,1>& a) { + os << "["; + for (int d=0; d<1; ++d) { + if (d>0) os << ","; + os << a[d]; + } + os << "]"; + return os; +} +ostream& operator<< (ostream& os, const vect<int,2>& a) { + os << "["; + for (int d=0; d<2; ++d) { + if (d>0) os << ","; + os << a[d]; + } + os << "]"; + return os; +} +ostream& operator<< (ostream& os, const vect<int,3>& a) { + os << "["; + for (int d=0; d<3; ++d) { + if (d>0) os << ","; + os << a[d]; + } + os << "]"; + return os; +} +#endif diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index b4b2a2bba..d5ad5cbab 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.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/vect.hh,v 1.2 2001/03/17 22:37:56 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $ ***************************************************************************/ @@ -21,10 +21,13 @@ #ifndef VECT_HH #define VECT_HH -#include <cassert> -#include <cmath> +#include <assert.h> +#include <math.h> + #include <iostream> +using namespace std; + // Forward definition @@ -201,25 +204,25 @@ public: // Non-modifying operators vect operator+ () const { - vect r(*this); + vect r; for (int d=0; d<D; ++d) r[d]=+r[d]; return r; } vect operator- () const { - vect r(*this); + vect r; for (int d=0; d<D; ++d) r[d]=-r[d]; return r; } vect<bool,D> operator! () const { - vect r(*this); + vect<bool,D> r; for (int d=0; d<D; ++d) r[d]=!r[d]; return r; } vect operator~ () const { - vect r(*this); + vect r; for (int d=0; d<D; ++d) r[d]=~r[d]; return r; } @@ -273,14 +276,14 @@ public: } vect<bool,D> operator&& (const T x) const { - vect r(*this); - for (int d=0; d<D; ++d) r[d]=r[d]&&x; + vect<bool,D> r; + for (int d=0; d<D; ++d) r[d]=(*this)[d]&&x; return r; } vect<bool,D> operator|| (const T x) const { - vect r(*this); - for (int d=0; d<D; ++d) r[d]=r[d]||x; + vect<bool,D> r; + for (int d=0; d<D; ++d) r[d]=(*this)[d]||x; return r; } @@ -402,8 +405,6 @@ public: }; #endif - // Output - friend ostream& operator<< <>(ostream& os, const vect& a); }; |