aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetLib/src
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetLib/src')
-rw-r--r--Carpet/CarpetLib/src/copy_3d_int4.F777
-rw-r--r--Carpet/CarpetLib/src/copy_3d_real8.F7737
-rw-r--r--Carpet/CarpetLib/src/data.cc679
-rw-r--r--Carpet/CarpetLib/src/data.hh17
-rw-r--r--Carpet/CarpetLib/src/dh.cc8
-rw-r--r--Carpet/CarpetLib/src/dh.hh4
-rw-r--r--Carpet/CarpetLib/src/dist.hh4
-rw-r--r--Carpet/CarpetLib/src/gdata.cc384
-rw-r--r--Carpet/CarpetLib/src/gdata.hh97
-rw-r--r--Carpet/CarpetLib/src/gf.cc4
-rw-r--r--Carpet/CarpetLib/src/ggf.cc141
-rw-r--r--Carpet/CarpetLib/src/ggf.hh39
-rw-r--r--Carpet/CarpetLib/src/gh.cc10
-rw-r--r--Carpet/CarpetLib/src/gh.hh5
-rw-r--r--Carpet/CarpetLib/src/make.code.defn11
-rw-r--r--Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F779
-rw-r--r--Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F779
-rw-r--r--Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F777
-rw-r--r--Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F777
-rw-r--r--Carpet/CarpetLib/src/restrict_3d_real8_rf2.F779
-rw-r--r--Carpet/CarpetLib/src/vect.hh6
21 files changed, 1063 insertions, 431 deletions
diff --git a/Carpet/CarpetLib/src/copy_3d_int4.F77 b/Carpet/CarpetLib/src/copy_3d_int4.F77
index b440a32b2..4cd070277 100644
--- a/Carpet/CarpetLib/src/copy_3d_int4.F77
+++ b/Carpet/CarpetLib/src/copy_3d_int4.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_int4.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_int4.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -54,11 +54,6 @@ c bbox(:,3) is stride
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)
diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77
index d115418b8..43507cdd3 100644
--- a/Carpet/CarpetLib/src/copy_3d_real8.F77
+++ b/Carpet/CarpetLib/src/copy_3d_real8.F77
@@ -1,18 +1,8 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.6 2003/02/25 22:57:00 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.7 2003/11/05 16:18:39 schnetter 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
+#include "cctk_Parameters.h"
@@ -23,6 +13,8 @@ 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
@@ -96,16 +88,19 @@ c This could be handled, but is likely to point to an error elsewhere
c Loop over region
- do k = 0, regkext-1
- do j = 0, regjext-1
- do i = 0, regiext-1
+ 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
- CHKIDX (srcioff+i+1, srcjoff+j+1, srckoff+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+i+1, srcjoff+j+1, srckoff+k+1)
+ dst (dstioff+i, dstjoff+j, dstkoff+k)
+ $ = src (srcioff+i, srcjoff+j, srckoff+k)
end do
end do
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index d2177364a..46ff5df79 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.32 2003/10/17 08:14:41 cvs_anon Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.33 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
#include <limits.h>
@@ -17,24 +17,36 @@
#include "data.hh"
-#include "util_ErrorCodes.h"
-#include "util_Table.h"
-
using namespace std;
+// Hand out the next MPI tag
+static int nexttag ()
+{
+ static int last = 100;
+ ++last;
+ if (last > 30000) last = 100;
+ return last;
+}
+
+
+
// Constructors
template<class T, int D>
data<T,D>::data (const int varindex_)
: gdata<D>(varindex_),
- _storage(0)
+ _storage(0),
+ comm_active(false),
+ tag(nexttag())
{ }
template<class T, int D>
data<T,D>::data (const int varindex_, const ibbox& extent_, const int proc_)
: gdata<D>(varindex_),
- _storage(0)
+ _storage(0),
+ comm_active(false),
+ tag(nexttag())
{
allocate(extent_, proc_);
}
@@ -102,7 +114,32 @@ void data<T,D>::transfer_from (gdata<D>* gsrc) {
// Processor management
template<class T, int D>
-void data<T,D>::change_processor (const int newproc, void* const mem) {
+void data<T,D>::change_processor (comm_state<D>& state,
+ const int newproc, void* const mem)
+{
+ switch (state.thestate) {
+ case state_recv:
+ change_processor_recv (newproc, mem);
+ break;
+ case state_send:
+ change_processor_send (newproc, mem);
+ break;
+ case state_wait:
+ change_processor_wait (newproc, mem);
+ break;
+ default:
+ assert(0);
+ }
+}
+
+
+
+template<class T, int D>
+void data<T,D>::change_processor_recv (const int newproc, void* const mem)
+{
+ assert (!comm_active);
+ comm_active = true;
+
if (newproc == this->_proc) {
assert (!mem);
return;
@@ -121,21 +158,87 @@ void data<T,D>::change_processor (const int newproc, void* const mem) {
} else {
_storage = (T*)mem;
}
-
+
T dummy;
- MPI_Status status;
- MPI_Recv (_storage, this->_size, dist::datatype(dummy), this->_proc,
- dist::tag, dist::comm, &status);
+ MPI_Irecv (_storage, this->_size, dist::datatype(dummy), this->_proc,
+ this->tag, dist::comm, &request);
+
+ } else if (rank == this->_proc) {
+ // copy to other processor
+
+ } else {
+ assert (!mem);
+ assert (!_storage);
+ }
+ }
+}
+
+
+template<class T, int D>
+void data<T,D>::change_processor_send (const int newproc, void* const mem)
+{
+ assert (comm_active);
+
+ if (newproc == this->_proc) {
+ assert (!mem);
+ return;
+ }
+
+ if (this->_has_storage) {
+ int rank;
+ MPI_Comm_rank (dist::comm, &rank);
+ if (rank == newproc) {
+ // copy from other processor
+
} else if (rank == this->_proc) {
// copy to other processor
assert (!mem);
assert (_storage);
+
T dummy;
- MPI_Send (_storage, this->_size, dist::datatype(dummy), newproc,
- dist::tag, dist::comm);
+ MPI_Isend (_storage, this->_size, dist::datatype(dummy), newproc,
+ this->tag, dist::comm, &request);
+
+ } else {
+ assert (!mem);
+ assert (!_storage);
+ }
+ }
+}
+
+
+template<class T, int D>
+void data<T,D>::change_processor_wait (const int newproc, void* const mem)
+{
+ assert (comm_active);
+ comm_active = false;
+
+ if (newproc == this->_proc) {
+ assert (!mem);
+ return;
+ }
+
+ if (this->_has_storage) {
+ int rank;
+ MPI_Comm_rank (dist::comm, &rank);
+ if (rank == newproc) {
+ // copy from other processor
+
+ MPI_Status status;
+ MPI_Wait (&request, &status);
+
+ } else if (rank == this->_proc) {
+ // copy to other processor
+
+ assert (!mem);
+ assert (_storage);
+
+ MPI_Status status;
+ MPI_Wait (&request, &status);
+
if (this->_owns_storage) {
delete [] _storage;
}
@@ -146,7 +249,7 @@ void data<T,D>::change_processor (const int newproc, void* const mem) {
assert (!_storage);
}
}
-
+
this->_proc = newproc;
}
@@ -169,6 +272,13 @@ void data<T,D>
&& (box.lower()-src->extent().lower())%box.stride() == 0));
assert (this->proc() == src->proc());
+
+ const int groupindex = CCTK_GroupIndexFromVarI(varindex);
+ const int group_tags_table = CCTK_GroupTagsTableI(groupindex);
+ assert (group_tags_table >= 0);
+
+ // Disallow this.
+ assert (0);
int rank;
MPI_Comm_rank (dist::comm, &rank);
@@ -212,15 +322,27 @@ void data<T,D>
MPI_Comm_rank (dist::comm, &rank);
assert (rank == this->proc());
- T Tdummy;
- CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
- "There is no interpolator available for variable type %s, dimension %d, spatial interpolation order %d, temporal interpolation order %d. The interpolation will not be done.",
- typestring(Tdummy), D, order_space, order_time);
+ assert (varindex >= 0);
+ const int groupindex = CCTK_GroupIndexFromVarI (varindex);
+ assert (groupindex >= 0);
+ char* groupname = CCTK_GroupName(groupindex);
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "There is no interpolator available for the group \"%s\" with variable type %s, dimension %d, spatial interpolation order %d, temporal interpolation order %d.",
+ groupname, D, order_space, order_time);
+ ::free (groupname);
}
extern "C" {
+ void CCTK_FCALL CCTK_FNAME(copy_3d_int4)
+ (const CCTK_INT4* src,
+ const int& srciext, const int& srcjext, const int& srckext,
+ CCTK_INT4* 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(copy_3d_real8)
(const CCTK_REAL8* src,
const int& srciext, const int& srcjext, const int& srckext,
@@ -232,6 +354,65 @@ extern "C" {
}
template<>
+void data<CCTK_INT4,3>
+::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box)
+{
+ const data* src = (const data*)gsrc;
+ 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));
+
+ 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(copy_3d_int4) ((const CCTK_INT4*)src->storage(),
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_INT4*)storage(),
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox,
+ dstbbox,
+ regbbox);
+
+ } else {
+ assert (0);
+ }
+}
+
+template<>
void data<CCTK_REAL8,3>
::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box)
{
@@ -302,6 +483,14 @@ 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_rf2)
+ (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]);
@@ -313,6 +502,14 @@ 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_rf2)
+ (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_o3)
(const CCTK_REAL8* src,
const int& srciext, const int& srcjext, const int& srckext,
@@ -321,6 +518,14 @@ 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_o3_rf2)
+ (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_minmod)
(const CCTK_REAL8* src,
const int& srciext, const int& srcjext, const int& srckext,
@@ -385,6 +590,16 @@ 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_3tl_rf2)
+ (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+ const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+ const CCTK_REAL8* src3, const CCTK_REAL8& t3,
+ const int& srciext, const int& srcjext, const int& srckext,
+ CCTK_REAL8* dst, const CCTK_REAL8& 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 CCTK_REAL8& t1,
const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -395,6 +610,16 @@ 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_3tl_o3_rf2)
+ (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+ const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+ const CCTK_REAL8* src3, const CCTK_REAL8& t3,
+ const int& srciext, const int& srcjext, const int& srckext,
+ CCTK_REAL8* dst, const CCTK_REAL8& 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_minmod)
(const CCTK_REAL8* src1, const CCTK_REAL8& t1,
const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -415,7 +640,6 @@ extern "C" {
const int srcbbox[3][3],
const int dstbbox[3][3],
const int regbbox[3][3]);
-
}
template<>
@@ -468,230 +692,151 @@ void data<CCTK_REAL8,3>
regbbox[1][d] = box.upper()[d];
regbbox[2][d] = box.stride()[d];
}
-
- const int varindex = (gsrcs[0])->var_index();
- const int groupindex = CCTK_GroupIndexFromVarI(varindex);
- const int group_tags_table = CCTK_GroupTagsTableI(groupindex);
- assert(group_tags_table >= 0);
-
- int typecode = -1;
- int keysize = -1;
- if (! Util_TableQueryValueInfo(group_tags_table,
- &typecode,
- &keysize,
- "Prolongation"))
- {
- CCTK_VWarn(4, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Tags table for group %s does not contain 'Prolongation'"
- " tag", CCTK_GroupName(groupindex));
- }
- else
- {
- assert(typecode == CCTK_VARIABLE_CHAR);
- }
-
-#define PROLONG_NONE 0
-#define PROLONG_LAGRANGE 1
-#define PROLONG_TVD 2
-
- int prolong_method;
- if (keysize < 0) { // No key in table - default to Lagrange.
- prolong_method = 1;
- }
- else {
- char prolong_string[keysize+10];
- const int error = Util_TableGetString(group_tags_table,
- keysize+10,
- prolong_string,
- "Prolongation");
- if (error < 0) {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Error code %d getting prolongation method"
- " from tags table for group %s.",
- error, CCTK_GroupName(groupindex));
- }
- if (CCTK_Equals(prolong_string, "None")){
- prolong_method = PROLONG_NONE;
- }
- else if (CCTK_Equals(prolong_string, "Lagrange")){
- prolong_method = PROLONG_LAGRANGE;
- }
- else if (CCTK_Equals(prolong_string, "TVD")){
- prolong_method = PROLONG_TVD;
- }
- else {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Prolongation method %s for group %s not recognized.\n"
- "Is this is a typo? Known values are:\n"
- "None\nLagrange\nTVD",
- prolong_method, CCTK_GroupName(groupindex));
- }
-
- }
-
-// CCTK_VInfo(CCTK_THORNSTRING,
-// "Variable %d is %s with method %d",
-// varindex, CCTK_VarName(varindex), prolong_method);
assert (all(dext.stride() == box.stride()));
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(),
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
+ switch (transport_operator) {
+
+ case op_Lagrange:
+ case op_TVD:
+ assert (srcs.size() == 1);
+ if (all (dext.stride() == sext.stride() * 2)) {
+ CCTK_FNAME(restrict_3d_real8_rf2)
+ ((const CCTK_REAL8*)srcs[0]->storage(),
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(),
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ } else {
+ CCTK_FNAME(restrict_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;
+
+ default:
+ assert (0);
+
+ } // switch (prolong_method)
} else if (all(sext.stride() > dext.stride())) {
-
- switch (prolong_method) {
-
- case PROLONG_NONE: // PROLONG_NONE: NOTHING
- break;
- case PROLONG_LAGRANGE: // PROLONG_LAGRANGE: DEFAULT
- switch (order_time) {
-
- case 0:
- assert (srcs.size()>=1);
- switch (order_space) {
- case 0:
- 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 2:
- 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;
- case 4:
- case 5:
- CCTK_FNAME(prolongate_3d_real8_o5)
- ((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:
- assert (0);
- }
- break;
-
- case 1:
- assert (srcs.size()>=2);
- switch (order_space) {
- case 0:
- case 1:
- CCTK_FNAME(prolongate_3d_real8_2tl)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- case 2:
- case 3:
- CCTK_FNAME(prolongate_3d_real8_2tl_o3)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- case 4:
- case 5:
- CCTK_FNAME(prolongate_3d_real8_2tl_o5)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- default:
- assert (0);
- }
- break;
-
- case 2:
- assert (srcs.size()>=3);
- switch (order_space) {
- case 0:
- case 1:
- CCTK_FNAME(prolongate_3d_real8_3tl)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- (const CCTK_REAL8*)srcs[2]->storage(), times[2],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- case 2:
- case 3:
- CCTK_FNAME(prolongate_3d_real8_3tl_o3)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- (const CCTK_REAL8*)srcs[2]->storage(), times[2],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- case 4:
- case 5:
- CCTK_FNAME(prolongate_3d_real8_3tl_o5)
- ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
- (const CCTK_REAL8*)srcs[1]->storage(), times[1],
- (const CCTK_REAL8*)srcs[2]->storage(), times[2],
- srcshp[0], srcshp[1], srcshp[2],
- (CCTK_REAL8*)storage(), time,
- dstshp[0], dstshp[1], dstshp[2],
- srcbbox, dstbbox, regbbox);
- break;
- default:
- assert (0);
- }
- break;
-
- default:
- assert (0);
- } // switch (order_time)
- break;
- case PROLONG_TVD: // PROLONG_TVD
- switch (order_time) {
- case 0:
- CCTK_FNAME(prolongate_3d_real8_minmod)
+
+ switch (transport_operator) {
+
+ case op_Lagrange:
+ switch (order_time) {
+
+ case 0:
+ assert (srcs.size()>=1);
+ switch (order_space) {
+ case 0:
+ case 1:
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_rf2)
+ ((const CCTK_REAL8*)srcs[0]->storage(),
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(),
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ } else {
+ 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 2:
+ case 3:
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_o3_rf2)
((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 1:
- CCTK_FNAME(prolongate_3d_real8_2tl_minmod)
+ } else {
+ 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;
+ case 4:
+ case 5:
+ CCTK_FNAME(prolongate_3d_real8_o5)
+ ((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:
+ assert (0);
+ }
+ break;
+
+ case 1:
+ assert (srcs.size()>=2);
+ switch (order_space) {
+ case 0:
+ case 1:
+ CCTK_FNAME(prolongate_3d_real8_2tl)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ break;
+ case 2:
+ case 3:
+ CCTK_FNAME(prolongate_3d_real8_2tl_o3)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ break;
+ case 4:
+ case 5:
+ CCTK_FNAME(prolongate_3d_real8_2tl_o5)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ break;
+ default:
+ assert (0);
+ }
+ break;
+
+ case 2:
+ assert (srcs.size()>=3);
+ switch (order_space) {
+ case 0:
+ case 1:
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_3tl_rf2)
((const CCTK_REAL8*)srcs[0]->storage(), times[0],
(const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ (const CCTK_REAL8*)srcs[2]->storage(), times[2],
srcshp[0], srcshp[1], srcshp[2],
(CCTK_REAL8*)storage(), time,
dstshp[0], dstshp[1], dstshp[2],
srcbbox, dstbbox, regbbox);
- break;
- case 2:
- CCTK_FNAME(prolongate_3d_real8_3tl_minmod)
+ } else {
+ CCTK_FNAME(prolongate_3d_real8_3tl)
((const CCTK_REAL8*)srcs[0]->storage(), times[0],
(const CCTK_REAL8*)srcs[1]->storage(), times[1],
(const CCTK_REAL8*)srcs[2]->storage(), times[2],
@@ -699,16 +844,88 @@ void data<CCTK_REAL8,3>
(CCTK_REAL8*)storage(), time,
dstshp[0], dstshp[1], dstshp[2],
srcbbox, dstbbox, regbbox);
- break;
- default:
- assert (0);
+ }
+ break;
+ case 2:
+ case 3:
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_3tl_o3_rf2)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ } else {
+ CCTK_FNAME(prolongate_3d_real8_3tl_o3)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ }
+ break;
+ case 4:
+ case 5:
+ CCTK_FNAME(prolongate_3d_real8_3tl_o5)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ break;
+ default:
+ assert (0);
}
break;
+
default:
- assert(0);
+ assert (0);
+ } // switch (order_time)
+ break;
+
+ case op_TVD:
+ switch (order_time) {
+ case 0:
+ CCTK_FNAME(prolongate_3d_real8_minmod)
+ ((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 1:
+ CCTK_FNAME(prolongate_3d_real8_2tl_minmod)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ case 2:
+ CCTK_FNAME(prolongate_3d_real8_3tl_minmod)
+ ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+ (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+ (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+ srcshp[0], srcshp[1], srcshp[2],
+ (CCTK_REAL8*)storage(), time,
+ dstshp[0], dstshp[1], dstshp[2],
+ srcbbox, dstbbox, regbbox);
+ break;
+ default:
+ assert (0);
+ }
+ break;
+
+ default:
+ assert(0);
} // switch (prolong_method)
-
} else {
assert (0);
}
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index f126c6c47..bb8e5c497 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.14 2003/10/14 16:39:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.15 2003/11/05 16:18:39 schnetter Exp $
#ifndef DATA_HH
#define DATA_HH
@@ -30,7 +30,12 @@ class data: public gdata<D> {
// Fields
T* _storage; // the data (if located on this processor)
-
+
+ bool comm_active;
+ MPI_Request request;
+
+ int tag; // MPI tag for this object
+
public:
// Constructors
@@ -50,7 +55,13 @@ public:
virtual void transfer_from (gdata<D>* gsrc);
// Processor management
- virtual void change_processor (const int newproc, void* const mem=0);
+ virtual void change_processor (comm_state<D>& state,
+ const int newproc, void* const mem=0);
+ private:
+ virtual void change_processor_recv (const int newproc, void* const mem=0);
+ virtual void change_processor_send (const int newproc, void* const mem=0);
+ virtual void change_processor_wait (const int newproc, void* const mem=0);
+ public:
// Accessors
virtual const void* storage () const {
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 83d43932b..4cdc1ae28 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.43 2003/10/16 17:00:04 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.44 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
@@ -31,7 +31,7 @@ dh<D>::dh (gh<D>& h,
assert (buffer_width>=0);
h.add(this);
CHECKPOINT;
- recompose(false);
+ recompose (0, true);
}
// Destructors
@@ -51,7 +51,7 @@ int dh<D>::prolongation_stencil_size () const {
// Modifiers
template<int D>
-void dh<D>::recompose (const int initialise_upto) {
+void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
DECLARE_CCTK_PARAMETERS;
CHECKPOINT;
@@ -562,7 +562,7 @@ void dh<D>::recompose (const int initialise_upto) {
for (typename list<ggf<D>*>::iterator f=gfs.begin();
f!=gfs.end(); ++f) {
- (*f)->recompose(initialise_upto);
+ (*f)->recompose (initialise_from, do_prolongate);
}
}
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 4a23d5cba..0f9cbaaa5 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.15 2003/07/20 21:03:43 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.16 2003/11/05 16:18:39 schnetter Exp $
#ifndef DH_HH
#define DH_HH
@@ -115,7 +115,7 @@ public:
int prolongation_stencil_size () const;
// Modifiers
- void recompose (const int initialise_upto);
+ void recompose (const int initialise_from, const bool do_prolongate);
// Grid function management
void add (ggf<D>* f);
diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh
index 10d4a3cfc..cb3d422ca 100644
--- a/Carpet/CarpetLib/src/dist.hh
+++ b/Carpet/CarpetLib/src/dist.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.8 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.9 2003/11/05 16:18:39 schnetter Exp $
#ifndef DIST_HH
#define DIST_HH
@@ -19,8 +19,6 @@ using namespace std;
namespace dist {
- const int tag = 1;
-
extern MPI_Comm comm;
extern MPI_Datatype mpi_complex_float;
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index e438a3147..4236b4b41 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -1,11 +1,15 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.23 2003/10/14 16:39:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.24 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
+#include <stdlib.h>
#include <iostream>
#include "cctk.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
#include "bbox.hh"
#include "defs.hh"
#include "dist.hh"
@@ -17,11 +21,43 @@ using namespace std;
+// Communication state control
+template<int D>
+comm_state<D>::comm_state ()
+ : thestate(state_recv),
+ current(0)
+{
+}
+
+template<int D>
+void comm_state<D>::step ()
+{
+ assert (thestate!=state_done);
+ assert (current==tmps.size());
+ thestate=astate(size_t(thestate)+1);
+ current=0;
+}
+
+template<int D>
+bool comm_state<D>::done ()
+{
+ return thestate==state_done;
+}
+
+template<int D>
+comm_state<D>::~comm_state ()
+{
+ assert (thestate==state_recv || thestate==state_done);
+}
+
+
+
// Constructors
template<int D>
gdata<D>::gdata (const int varindex_)
: varindex(varindex_),
- _has_storage(false)
+ _has_storage(false),
+ transport_operator (find_transport_operator(varindex_))
{ }
// Destructors
@@ -30,9 +66,93 @@ gdata<D>::~gdata () { }
+// Transport operator types
+template<int D>
+gdata<D>::operator_type gdata<D>::find_transport_operator (const int varindex) {
+ const operator_type default_operator = op_Lagrange;
+ if (varindex == -1) return op_error;
+ assert (varindex >= 0);
+ const int groupindex = CCTK_GroupIndexFromVarI (varindex);
+ assert (groupindex >= 0);
+ const int group_tags_table = CCTK_GroupTagsTableI(groupindex);
+ assert (group_tags_table >= 0);
+ char prolong_string[100];
+ const int ierr = Util_TableGetString
+ (group_tags_table, sizeof prolong_string, prolong_string, "Prolongation");
+ if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ char* groupname = CCTK_GroupName(groupindex);
+ CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Tags table for group \"%s\" does not contain \"Prolongation\" tag", groupname);
+ ::free (groupname);
+ return default_operator;
+ }
+ assert (ierr >= 0);
+ if (CCTK_Equals(prolong_string, "None")) {
+ return op_none;
+ } else if (CCTK_Equals(prolong_string, "Lagrange")) {
+ return op_Lagrange;
+ } else if (CCTK_Equals(prolong_string, "TVD")) {
+ return op_TVD;
+ } else {
+ assert (0);
+ }
+ return op_error;
+}
+
+
+
// Data manipulators
template<int D>
-void gdata<D>::copy_from (const gdata* src, const ibbox& box)
+void gdata<D>::copy_from (comm_state<D>& state,
+ const gdata* src, const ibbox& box)
+{
+ switch (state.thestate) {
+ case state_recv:
+ copy_from_recv (state, src, box);
+ break;
+ case state_send:
+ copy_from_send (state, src, box);
+ break;
+ case state_wait:
+ copy_from_wait (state, src, box);
+ break;
+ default:
+ assert(0);
+ }
+}
+
+
+
+template<int D>
+void gdata<D>::copy_from_nocomm (const gdata* 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 (box.empty()) return;
+
+ assert (proc() == src->proc());
+
+ // copy on same processor
+ int rank;
+ MPI_Comm_rank (dist::comm, &rank);
+ if (rank == proc()) {
+ copy_from_innerloop (src, box);
+ }
+}
+
+
+
+template<int D>
+void gdata<D>::copy_from_recv (comm_state<D>& state,
+ const gdata* src, const ibbox& box)
{
assert (has_storage() && src->has_storage());
assert (all(box.lower()>=extent().lower()
@@ -49,20 +169,81 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box)
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
- gdata* const tmp = make_typed(varindex);
+ gdata<D>* const tmp = make_typed(varindex);
+ // TODO: is this efficient?
+ state.tmps.push_back (tmp);
+ ++state.current;
tmp->allocate (box, src->proc());
- tmp->copy_from (src, box);
- tmp->change_processor (proc());
- copy_from (tmp, box);
+ tmp->change_processor_recv (proc());
+
+ }
+}
+
+
+
+template<int D>
+void gdata<D>::copy_from_send (comm_state<D>& state,
+ const gdata* 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 (box.empty()) return;
+
+ if (proc() == src->proc()) {
+ // copy on same processor
+
+ copy_from_nocomm (src, box);
+
+ } else {
+
+ // copy to different processor
+ gdata<D>* const tmp = state.tmps.at(state.current++);
+ assert (tmp);
+ tmp->copy_from_nocomm (src, box);
+ tmp->change_processor_send (proc());
+
+ }
+}
+
+
+
+template<int D>
+void gdata<D>::copy_from_wait (comm_state<D>& state,
+ const gdata* 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 (box.empty()) return;
+
+ if (proc() == src->proc()) {
+ // copy on same processor
+
+ } else {
+
+ // copy to different processor
+ gdata<D>* const tmp = state.tmps.at(state.current++);
+ assert (tmp);
+ tmp->change_processor_wait (proc());
+ copy_from_nocomm (tmp, box);
delete tmp;
}
@@ -72,11 +253,79 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box)
template<int D>
void gdata<D>
-::interpolate_from (const vector<const gdata*> srcs,
- const vector<CCTK_REAL> times,
- const ibbox& box, const CCTK_REAL time,
- const int order_space,
- const int order_time)
+::interpolate_from (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time)
+{
+ assert (transport_operator != op_error);
+ if (transport_operator == op_none) return;
+ switch (state.thestate) {
+ case state_recv:
+ interpolate_from_recv (state, srcs, times, box, time, order_space, order_time);
+ break;
+ case state_send:
+ interpolate_from_send (state, srcs, times, box, time, order_space, order_time);
+ break;
+ case state_wait:
+ interpolate_from_wait (state, srcs, times, box, time, order_space, order_time);
+ break;
+ default:
+ assert(0);
+ }
+}
+
+
+
+template<int D>
+void gdata<D>
+::interpolate_from_nocomm (const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time)
+{
+ 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 (srcs.size() == times.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 (! box.empty());
+ if (box.empty()) return;
+
+ assert (proc() == srcs[0]->proc());
+
+ assert (transport_operator != op_error);
+ assert (transport_operator != op_none);
+
+ // interpolate on same processor
+ int rank;
+ MPI_Comm_rank (dist::comm, &rank);
+ if (rank == proc()) {
+ interpolate_from_innerloop
+ (srcs, times, box, time, order_space, order_time);
+ }
+}
+
+
+
+template<int D>
+void gdata<D>
+::interpolate_from_recv (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time)
{
assert (has_storage());
assert (all(box.lower()>=extent().lower()));
@@ -96,21 +345,97 @@ void gdata<D>
if (proc() == srcs[0]->proc()) {
// interpolate on same processor
- int rank;
- MPI_Comm_rank (dist::comm, &rank);
- if (rank == proc()) {
- interpolate_from_innerloop
- (srcs, times, box, time, order_space, order_time);
- }
-
} else {
// interpolate from other processor
- gdata* const tmp = make_typed(varindex);
+ gdata<D>* const tmp = make_typed(varindex);
+ // TODO: is this efficient?
+ state.tmps.push_back (tmp);
+ ++state.current;
tmp->allocate (box, srcs[0]->proc());
- tmp->interpolate_from (srcs, times, box, time, order_space, order_time);
- tmp->change_processor (proc());
- copy_from (tmp, box);
+ tmp->change_processor_recv (proc());
+
+ }
+}
+
+
+
+template<int D>
+void gdata<D>
+::interpolate_from_send (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time)
+{
+ 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 (srcs.size() == times.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 (! box.empty());
+ if (box.empty()) return;
+
+ if (proc() == srcs[0]->proc()) {
+ // interpolate on same processor
+
+ interpolate_from_nocomm (srcs, times, box, time, order_space, order_time);
+
+ } else {
+ // interpolate from other processor
+
+ gdata<D>* const tmp = state.tmps.at(state.current++);
+ assert (tmp);
+ tmp->interpolate_from_nocomm (srcs, times, box, time, order_space, order_time);
+ tmp->change_processor_send (proc());
+
+ }
+}
+
+
+
+template<int D>
+void gdata<D>
+::interpolate_from_wait (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time)
+{
+ 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 (srcs.size() == times.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 (! box.empty());
+ if (box.empty()) return;
+
+ if (proc() == srcs[0]->proc()) {
+ // interpolate on same processor
+
+ } else {
+ // interpolate from other processor
+
+ gdata<D>* const tmp = state.tmps.at(state.current++);
+ assert (tmp);
+ tmp->change_processor_wait (proc());
+ copy_from_nocomm (tmp, box);
delete tmp;
}
@@ -118,4 +443,5 @@ void gdata<D>
+template class comm_state<3>;
template class gdata<3>;
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index edc08c8e7..808a8eb51 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.19 2003/10/15 07:14:01 hawke Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.20 2003/11/05 16:18:39 schnetter Exp $
#ifndef GDATA_HH
#define GDATA_HH
@@ -8,6 +8,7 @@
#include <iostream>
#include <string>
+#include <vector>
#include "cctk.h"
@@ -19,6 +20,29 @@
using namespace std;
+
+
+template<int D>
+class gdata;
+
+
+
+// State information for communications
+enum astate { state_recv, state_send, state_wait, state_done };
+
+template<int D>
+struct comm_state {
+ astate thestate;
+ comm_state ();
+ void step ();
+ bool done ();
+ ~comm_state ();
+
+ vector<gdata<D>*> tmps;
+ size_t current;
+};
+
+
// A generic data storage without type information
template<int D>
@@ -57,7 +81,13 @@ public:
virtual gdata<D>* make_typed (const int varindex) const = 0;
// Processor management
- virtual void change_processor (const int newproc, void* const mem=0) = 0;
+ virtual void change_processor (comm_state<D>& state,
+ const int newproc, void* const mem=0) = 0;
+ protected:
+ virtual void change_processor_recv (const int newproc, void* const mem=0) = 0;
+ virtual void change_processor_send (const int newproc, void* const mem=0) = 0;
+ virtual void change_processor_wait (const int newproc, void* const mem=0) = 0;
+ public:
// Storage management
virtual void transfer_from (gdata<D>* src) = 0;
@@ -68,10 +98,6 @@ public:
// Accessors
- int var_index () const {
- return varindex;
- }
-
bool has_storage () const {
return _has_storage;
}
@@ -117,14 +143,59 @@ public:
assert (all(ind>=0 && ind<=shape()));
return dot(ind, stride());
}
-
+
+ // Transport operator types
+ protected:
+ enum operator_type { op_error, op_none, op_Lagrange, op_TVD };
+ // readonly
+ operator_type transport_operator;
+ private:
+ static operator_type find_transport_operator (const int varindex_);
+
// Data manipulators
- void copy_from (const gdata* src, const ibbox& box);
- void interpolate_from (const vector<const gdata*> srcs,
- const vector<CCTK_REAL> times,
- const ibbox& box, const CCTK_REAL time,
- const int order_space,
- const int order_time);
+ public:
+ void copy_from (comm_state<D>& state, const gdata* src, const ibbox& box);
+ private:
+ void copy_from_nocomm (const gdata* src, const ibbox& box);
+ void copy_from_recv (comm_state<D>& state,
+ const gdata* src, const ibbox& box);
+ void copy_from_send (comm_state<D>& state,
+ const gdata* src, const ibbox& box);
+ void copy_from_wait (comm_state<D>& state,
+ const gdata* src, const ibbox& box);
+ public:
+ void interpolate_from (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time);
+ private:
+ void interpolate_from_nocomm (const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time);
+ void interpolate_from_recv (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time);
+ void interpolate_from_send (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time);
+ void interpolate_from_wait (comm_state<D>& state,
+ const vector<const gdata*> srcs,
+ const vector<CCTK_REAL> times,
+ const ibbox& box, const CCTK_REAL time,
+ const int order_space,
+ const int order_time);
+ public:
+
protected:
virtual void
copy_from_innerloop (const gdata* src, const ibbox& box) = 0;
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 8de652672..4a913d199 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.13 2003/10/14 16:39:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.14 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
@@ -20,7 +20,7 @@ gf<T,D>::gf (const int varindex, th<D>& t, dh<D>& d,
: ggf<D>(varindex, t, d, tmin, tmax, prolongation_order_time)
{
// VGF
- this->recompose();
+ this->recompose (0, true);
}
// Destructors
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index a096ce0fb..6786e5eb7 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.28 2003/10/14 16:39:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.29 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
#include <stdlib.h>
@@ -51,7 +51,9 @@ bool ggf<D>::operator== (const ggf<D>& f) const {
// Modifiers
// VGF
template<int D>
-void ggf<D>::recompose (const int initialise_upto) {
+void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
+
+ // TODO: restructure storage only when needed
// Retain storage that might be needed
fdata oldstorage = storage;
@@ -60,19 +62,19 @@ void ggf<D>::recompose (const int initialise_upto) {
storage.resize(tmax-tmin+1);
for (int tl=tmin; tl<=tmax; ++tl) {
storage[tl-tmin].resize(h.reflevels());
- for (int rl=0; rl<h.reflevels(); ++rl) {
+ for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
storage[tl-tmin][rl].resize(h.components(rl));
for (int c=0; c<h.components(rl); ++c) {
storage[tl-tmin][rl][c].resize(h.mglevels(rl,c));
for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
- storage[tl-tmin][rl][c][ml] = 0;
+ storage[tl-tmin][rl][c][ml] = 0;
} // for ml
} // for c
} // for rl
} // for tl
// Initialise the new storage
- for (int rl=0; rl<h.reflevels(); ++rl) {
+ for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
for (int tl=tmin; tl<=tmax; ++tl) {
for (int c=0; c<h.components(rl); ++c) {
for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
@@ -84,30 +86,32 @@ void ggf<D>::recompose (const int initialise_upto) {
storage[tl-tmin][rl][c][ml]->allocate
(d.boxes[rl][c][ml].exterior, h.proc(rl,c));
- // Initialise only if desired
- if (initialise_upto >= 0 && rl <= initialise_upto) {
-
+ if (do_prolongate) {
// Initialise from coarser level, if possible
// TODO: init only un-copied regions
if (rl>0) {
- const CCTK_REAL time = t.time(tl,rl,ml);
- ref_prolongate (tl,rl,c,ml,time);
+ for (comm_state<D> state; !state.done(); state.step()) {
+ const CCTK_REAL time = t.time(tl,rl,ml);
+ ref_prolongate (state,tl,rl,c,ml,time);
+ }
} // if rl
-
- // Copy from old storage, if possible
- if (rl<(int)oldstorage[tl-tmin].size()) {
+ } // if do_prolongate
+
+ // Copy from old storage, if possible
+ // todo: copy only from interior regions?
+ if (rl<(int)oldstorage[tl-tmin].size()) {
+ for (comm_state<D> state; !state.done(); state.step()) {
for (int cc=0; cc<(int)oldstorage[tl-tmin][rl].size(); ++cc) {
if (ml<(int)oldstorage[tl-tmin][rl][cc].size()) {
const ibbox ovlp =
(d.boxes[rl][c][ml].exterior
& oldstorage[tl-tmin][rl][cc][ml]->extent());
storage[tl-tmin][rl][c][ml]->copy_from
- (oldstorage[tl-tmin][rl][cc][ml], ovlp);
+ (state, oldstorage[tl-tmin][rl][cc][ml], ovlp);
} // if ml
} // for cc
- } // if rl
-
- } // if initialise
+ } // for step
+ } // if rl
} // for ml
} // for c
@@ -117,7 +121,7 @@ void ggf<D>::recompose (const int initialise_upto) {
// Delete old storage
for (int tl=tmin; tl<=tmax; ++tl) {
- for (int rl=0; rl<(int)oldstorage[tl-tmin].size(); ++rl) {
+ for (int rl=initialise_from; rl<(int)oldstorage[tl-tmin].size(); ++rl) {
for (int c=0; c<(int)oldstorage[tl-tmin][rl].size(); ++c) {
for (int ml=0; ml<(int)oldstorage[tl-tmin][rl][c].size(); ++ml) {
delete oldstorage[tl-tmin][rl][c][ml];
@@ -126,30 +130,29 @@ void ggf<D>::recompose (const int initialise_upto) {
} // for rl
} // for tl
- for (int rl=0; rl<h.reflevels(); ++rl) {
- for (int tl=tmin; tl<=tmax; ++tl) {
+ for (int tl=tmin; tl<=tmax; ++tl) {
+ for (int rl=initialise_from; rl<h.reflevels(); ++rl) {
// Set boundaries
for (int c=0; c<h.components(rl); ++c) {
for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
- // Initialise only if desired
- if (initialise_upto >= 0 && rl <= initialise_upto) {
-
- sync (tl,rl,c,ml);
- // TODO: assert that reflevel 0 boundaries are copied
- if (rl>0) {
+ for (comm_state<D> state; !state.done(); state.step()) {
+ sync (state,tl,rl,c,ml);
+ }
+ // TODO: assert that reflevel 0 boundaries are copied
+ if (rl>0) {
+ for (comm_state<D> state; !state.done(); state.step()) {
const CCTK_REAL time = t.time(tl,rl,ml);
- ref_bnd_prolongate (tl,rl,c,ml,time);
- } // if rl
-
- } // if initialise
-
+ ref_bnd_prolongate (state,tl,rl,c,ml,time);
+ }
+ } // if rl
+
} // for ml
} // for c
- } // for tl
- } // for rl
+ } // for rl
+ } // for tl
}
// Cycle the time levels by rotating the data sets
@@ -181,6 +184,7 @@ void ggf<D>::flip (int rl, int c, int ml) {
}
}
+#if 0
// Copy data from current time level to all previous levels
template<int D>
void ggf<D>::copytoprevs (int rl, int c, int ml) {
@@ -192,6 +196,7 @@ void ggf<D>::copytoprevs (int rl, int c, int ml) {
(storage[tmax-tmin][rl][c][ml], storage[tmax-tmin][rl][c][ml]->extent());
}
}
+#endif
@@ -201,7 +206,8 @@ void ggf<D>::copytoprevs (int rl, int c, int ml) {
// Copy a region
template<int D>
-void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::copycat (comm_state<D>& state,
+ 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)
@@ -220,12 +226,13 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
// copy the content
assert (recv==send);
storage[tl1-tmin][rl1][c1][ml1]->copy_from
- (storage[tl2-tmin][rl2][c2][ml2], recv);
+ (state, storage[tl2-tmin][rl2][c2][ml2], recv);
}
// Copy regions
template<int D>
-void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::copycat (comm_state<D>& state,
+ 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)
@@ -247,13 +254,14 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
// (use the send boxes for communication)
// copy the content
storage[tl1-tmin][rl1][c1][ml1]->copy_from
- (storage[tl2-tmin][rl2][c2][ml2], *r);
+ (state, storage[tl2-tmin][rl2][c2][ml2], *r);
}
}
// Copy regions
template<int D>
-void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::copycat (comm_state<D>& state,
+ 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)
@@ -276,14 +284,15 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1,
// (use the send boxes for communication)
// copy the content
storage[tl1-tmin][rl1][c1][ml1]->copy_from
- (storage[tl2-tmin][rl2][c2][ml2], *r);
+ (state, storage[tl2-tmin][rl2][c2][ml2], *r);
}
}
}
// Interpolate a region
template<int D>
-void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const ibbox dh<D>::dboxes::* recv_list,
const vector<int> tl2s, int rl2, int ml2,
const ibbox dh<D>::dboxes::* send_list,
@@ -317,13 +326,14 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
// interpolate the content
assert (recv==send);
storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
- (gsrcs, times, recv, time,
+ (state, gsrcs, times, recv, time,
d.prolongation_order_space, prolongation_order_time);
}
// Interpolate regions
template<int D>
-void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const iblist dh<D>::dboxes::* recv_list,
const vector<int> tl2s, int rl2, int ml2,
const iblist dh<D>::dboxes::* send_list,
@@ -360,14 +370,15 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
// (use the send boxes for communication)
// interpolate the content
storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
- (gsrcs, times, *r, time,
+ (state, gsrcs, times, *r, time,
d.prolongation_order_space, prolongation_order_time);
}
}
// Interpolate regions
template<int D>
-void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
+void ggf<D>::intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const iblistvect dh<D>::dboxes::* recv_listvect,
const vector<int> tl2s, int rl2, int ml2,
const iblistvect dh<D>::dboxes::* send_listvect,
@@ -405,7 +416,7 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
// (use the send boxes for communication)
// interpolate the content
storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
- (gsrcs, times, *r, time,
+ (state, gsrcs, times, *r, time,
d.prolongation_order_space, prolongation_order_time);
}
}
@@ -415,25 +426,28 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1,
// Copy a component from the next time level
template<int D>
-void ggf<D>::copy (int tl, int rl, int c, int ml)
+void ggf<D>::copy (comm_state<D>& state, int tl, int rl, int c, int ml)
{
// Copy
- copycat (tl ,rl,c,ml, &dh<D>::dboxes::exterior,
+ copycat (state,
+ tl ,rl,c,ml, &dh<D>::dboxes::exterior,
tl+1,rl, ml, &dh<D>::dboxes::exterior);
}
// Synchronise the boundaries a component
template<int D>
-void ggf<D>::sync (int tl, int rl, int c, int ml)
+void ggf<D>::sync (comm_state<D>& state, int tl, int rl, int c, int ml)
{
// Copy
- copycat (tl,rl,c,ml, &dh<D>::dboxes::recv_sync,
+ copycat (state,
+ tl,rl,c,ml, &dh<D>::dboxes::recv_sync,
tl,rl, ml, &dh<D>::dboxes::send_sync);
}
// Prolongate the boundaries of a component
template<int D>
-void ggf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml,
+void ggf<D>::ref_bnd_prolongate (comm_state<D>& state,
+ int tl, int rl, int c, int ml,
CCTK_REAL time)
{
// Interpolate
@@ -443,53 +457,61 @@ void ggf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml,
assert (tmax-tmin+1 >= prolongation_order_time+1);
tl2s.resize(prolongation_order_time+1);
for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
- intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
+ intercat (state,
+ tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine,
time);
}
// Restrict a multigrid level
template<int D>
-void ggf<D>::mg_restrict (int tl, int rl, int c, int ml,
+void ggf<D>::mg_restrict (comm_state<D>& state,
+ int tl, int rl, int c, int ml,
CCTK_REAL time)
{
// Require same times
assert (t.get_time(rl,ml) == t.get_time(rl,ml-1));
const vector<int> tl2s(1,tl);
- intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse,
+ intercat (state,
+ tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse,
tl2s,rl, ml-1, &dh<D>::dboxes::send_mg_fine,
time);
}
// Prolongate a multigrid level
template<int D>
-void ggf<D>::mg_prolongate (int tl, int rl, int c, int ml,
+void ggf<D>::mg_prolongate (comm_state<D>& state,
+ int tl, int rl, int c, int ml,
CCTK_REAL time)
{
// Require same times
assert (t.get_time(rl,ml) == t.get_time(rl,ml+1));
const vector<int> tl2s(1,tl);
- intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse,
+ intercat (state,
+ tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse,
tl2s,rl, ml+1, &dh<D>::dboxes::send_mg_fine,
time);
}
// Restrict a refinement level
template<int D>
-void ggf<D>::ref_restrict (int tl, int rl, int c, int ml,
+void ggf<D>::ref_restrict (comm_state<D>& state,
+ int tl, int rl, int c, int ml,
CCTK_REAL time)
{
// Require same times
assert (t.get_time(rl,ml) == t.get_time(rl+1,ml));
const vector<int> tl2s(1,tl);
- intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine,
+ intercat (state,
+ tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine,
tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse,
time);
}
// Prolongate a refinement level
template<int D>
-void ggf<D>::ref_prolongate (int tl, int rl, int c, int ml,
+void ggf<D>::ref_prolongate (comm_state<D>& state,
+ int tl, int rl, int c, int ml,
CCTK_REAL time)
{
assert (rl>=1);
@@ -498,7 +520,8 @@ void ggf<D>::ref_prolongate (int tl, int rl, int c, int ml,
assert (tmax-tmin+1 >= prolongation_order_time+1);
tl2s.resize(prolongation_order_time+1);
for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
- intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
+ intercat (state,
+ tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine,
time);
}
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 98c556800..b99c7aa10 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.16 2003/10/14 16:39:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.17 2003/11/05 16:18:39 schnetter Exp $
#ifndef GGF_HH
#define GGF_HH
@@ -80,7 +80,7 @@ public:
// Modifiers
// VGF
- void recompose (const int initialise_upto = -1);
+ void recompose (const int initialise_from, const bool do_prolongate);
// Cycle the time levels by rotating the data sets
void cycle (int rl, int c, int ml);
@@ -88,8 +88,11 @@ public:
// Flip the time levels by exchanging the data sets
void flip (int rl, int c, int ml);
+#if 0
// Copy data from current time level to all previous time levels
void copytoprevs (int rl, int c, int ml);
+#endif
+
// Helpers
@@ -105,39 +108,45 @@ protected:
protected:
// Copy a region
- void copycat (int tl1, int rl1, int c1, int ml1,
+ void copycat (comm_state<D>& state,
+ 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
- void copycat (int tl1, int rl1, int c1, int ml1,
+ void copycat (comm_state<D>& state,
+ 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
- void copycat (int tl1, int rl1, int c1, int ml1,
+ void copycat (comm_state<D>& state,
+ 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 region
- void intercat (int tl1, int rl1, int c1, int ml1,
+ void intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const ibbox dh<D>::dboxes::* recv_list,
const vector<int> tl2s, int rl2, int ml2,
const ibbox dh<D>::dboxes::* send_list,
CCTK_REAL time);
// Interpolate regions
- void intercat (int tl1, int rl1, int c1, int ml1,
+ void intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const iblist dh<D>::dboxes::* recv_list,
const vector<int> tl2s, int rl2, int ml2,
const iblist dh<D>::dboxes::* send_list,
CCTK_REAL time);
// Interpolate regions
- void intercat (int tl1, int rl1, int c1, int ml1,
+ void intercat (comm_state<D>& state,
+ int tl1, int rl1, int c1, int ml1,
const iblistvect dh<D>::dboxes::* recv_listvect,
const vector<int> tl2s, int rl2, int ml2,
const iblistvect dh<D>::dboxes::* send_listvect,
@@ -154,25 +163,25 @@ public:
// synchronised. They don't need to be prolongated.
// Copy a component from the next time level
- void copy (int tl, int rl, int c, int ml);
+ void copy (comm_state<D>& state, int tl, int rl, int c, int ml);
// Synchronise the boundaries of a component
- void sync (int tl, int rl, int c, int ml);
+ void sync (comm_state<D>& state, 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, CCTK_REAL time);
+ void ref_bnd_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time);
// Restrict a multigrid level
- void mg_restrict (int tl, int rl, int c, int ml, CCTK_REAL time);
+ void mg_restrict (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time);
// Prolongate a multigrid level
- void mg_prolongate (int tl, int rl, int c, int ml, CCTK_REAL time);
+ void mg_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time);
// Restrict a refinement level
- void ref_restrict (int tl, int rl, int c, int ml, CCTK_REAL time);
+ void ref_restrict (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time);
// Prolongate a refinement level
- void ref_prolongate (int tl, int rl, int c, int ml, CCTK_REAL time);
+ void ref_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index f8eeea4c0..b37fb3f15 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.22 2003/09/19 16:06:41 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.23 2003/11/05 16:18:39 schnetter Exp $
#include <assert.h>
#include <stdlib.h>
@@ -37,7 +37,8 @@ template<int D>
void gh<D>::recompose (const rexts& exts,
const rbnds& outer_bounds,
const rprocs& procs,
- const int initialise_upto)
+ const int initialise_from,
+ const bool do_prolongate)
{
DECLARE_CCTK_PARAMETERS;
@@ -81,6 +82,9 @@ void gh<D>::recompose (const rexts& exts,
assert (all(extents[rl][c][ml].stride()
== extents[rl][0][ml].stride()));
assert (extents[rl][c][ml].is_aligned_with(extents[rl][0][ml]));
+ for (int cc=c+1; cc<components(rl); ++cc) {
+ assert ((extents[rl][c][ml] & extents[rl][cc][ml]).empty());
+ }
}
}
}
@@ -161,7 +165,7 @@ void gh<D>::recompose (const rexts& exts,
}
for (typename list<dh<D>*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) {
- (*d)->recompose(initialise_upto);
+ (*d)->recompose (initialise_from, do_prolongate);
}
}
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index c748e93ee..db046185a 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.13 2003/05/02 14:23:12 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.14 2003/11/05 16:18:39 schnetter Exp $
#ifndef GH_HH
#define GH_HH
@@ -86,7 +86,8 @@ public:
// Modifiers
void recompose (const rexts& exts, const rbnds& outer_bounds,
const rprocs& procs,
- const int initialise_upto = -1);
+ const int initialise_from,
+ const bool do_prolongate);
// Helpers
cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index ab3874e22..0b718747a 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,5 +1,5 @@
# Main make.code.defn file for thorn CarpetLib -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.9 2003/10/14 14:53:46 hawke Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.10 2003/11/05 16:18:39 schnetter Exp $
# Source files in this directory
SRCS = bbox.cc \
@@ -14,20 +14,27 @@ SRCS = bbox.cc \
gh.cc \
th.cc \
vect.cc \
+ checkindex.F77 \
+ copy_3d_int4.F77 \
copy_3d_real8.F77 \
prolongate_3d_real8.F77 \
+ prolongate_3d_real8_rf2.F77 \
prolongate_3d_real8_o3.F77 \
+ prolongate_3d_real8_o3_rf2.F77 \
prolongate_3d_real8_o5.F77 \
prolongate_3d_real8_2tl.F77 \
prolongate_3d_real8_2tl_o3.F77 \
prolongate_3d_real8_2tl_o5.F77 \
prolongate_3d_real8_3tl.F77 \
+ prolongate_3d_real8_3tl_rf2.F77 \
prolongate_3d_real8_3tl_o3.F77 \
+ prolongate_3d_real8_3tl_o3_rf2.F77 \
prolongate_3d_real8_3tl_o5.F77 \
prolongate_3d_real8_minmod.F77 \
prolongate_3d_real8_2tl_minmod.F77 \
prolongate_3d_real8_3tl_minmod.F77 \
- restrict_3d_real8.F77
+ restrict_3d_real8.F77 \
+ restrict_3d_real8_rf2.F77
# Subdirectories containing source files
SUBDIRS =
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77
index daa130251..6c3cec79f 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -87,11 +87,6 @@ c bbox(:,3) is stride
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
srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
offsetlo = regbbox(d,3)
@@ -146,7 +141,7 @@ c Quadratic (second order) time interpolation
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 in time")
+ call CCTK_WARN (0, "Internal error: extrapolation")
end if
s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3))
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77
index a77498fef..2395b0821 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -77,11 +77,6 @@ c bbox(:,3) is stride
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)
@@ -120,7 +115,7 @@ c Quadratic (second order) time interpolation
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")
+ call CCTK_WARN (0, "Internal error: extrapolation")
end if
s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3))
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77
index 8a3b2629d..ad101b6c6 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -74,11 +74,6 @@ c bbox(:,3) is stride
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
srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
offsetlo = regbbox(d,3)
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77
index 556f4d092..6c59cc533 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -66,11 +66,6 @@ c bbox(:,3) is stride
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)
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77 b/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77
index b7f6a0454..71a8edfba 100644
--- a/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77
+++ b/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77
@@ -1,5 +1,5 @@
c -*-Fortran-*-
-c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -43,7 +43,7 @@ c bbox(:,3) is stride
call CCTK_WARN (0, "Internal error: strides disagree")
end if
if (dstbbox(d,3).ne.srcbbox(d,3)*2) then
- call CCTK_WARN (0, "Internal error: destination strides are not twice the source strides")
+ call CCTK_WARN (0, "Internal error: destination strides are not twice the sourcd strides")
end if
if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
$ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
@@ -54,11 +54,6 @@ c bbox(:,3) is stride
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)
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index 22d8ef054..dfb724694 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.19 2003/08/28 21:07:27 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.20 2003/11/05 16:18:39 schnetter Exp $
#ifndef VECT_HH
#define VECT_HH
@@ -58,7 +58,7 @@ public:
elt[0]=x; elt[1]=y;
}
- /** Constructor for 3-element vectors from 4 elements. */
+ /** Constructor for 3-element vectors from 3 elements. */
vect (const T x, const T y, const T z) {
assert (D==3);
elt[0]=x; elt[1]=y; elt[2]=z;
@@ -149,7 +149,7 @@ public:
/** Return a non-writable element of a vector. */
// Don't return a reference; *this might be a temporary
- const T operator[] (const int d) const {
+ T operator[] (const int d) const {
assert(d>=0 && d<D);
return elt[d];
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-01 05:40:17 +0000