aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetLib/src/data.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetLib/src/data.cc')
-rw-r--r--Carpet/CarpetLib/src/data.cc229
1 files changed, 183 insertions, 46 deletions
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index f6158bae8..6f65158b8 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,10 +1,15 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.38 2004/01/21 16:32:04 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.39 2004/01/25 14:57:29 schnetter Exp $
#include <assert.h>
#include <limits.h>
+#include <stdlib.h>
+#include <math.h>
+#include <algorithm>
#include <iostream>
+#include <sstream>
#include <string>
+#include <vector>
#include <mpi.h>
@@ -34,16 +39,17 @@ static int nexttag ()
// Constructors
template<class T, int D>
-data<T,D>::data (const int varindex_)
- : gdata<D>(varindex_),
+data<T,D>::data (const int varindex_, const operator_type transport_operator_)
+ : gdata<D>(varindex_, transport_operator_),
_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_),
+data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+ const ibbox& extent_, const int proc_)
+ : gdata<D>(varindex_, transport_operator_),
_storage(0),
comm_active(false),
tag(nexttag())
@@ -59,8 +65,11 @@ data<T,D>::~data () {
// Pseudo constructors
template<class T, int D>
-data<T,D>* data<T,D>::make_typed (const int varindex_) const {
- return new data(varindex_);
+data<T,D>* data<T,D>::make_typed (const int varindex_,
+ const operator_type transport_operator_)
+ const
+{
+ return new data(varindex_, transport_operator_);
}
@@ -107,7 +116,7 @@ void data<T,D>::transfer_from (gdata<D>* gsrc) {
data* src = (data*)gsrc;
assert (!_storage);
*this = *src;
- *src = data(this->varindex);
+ *src = data(this->varindex, this->transport_operator);
}
@@ -159,9 +168,12 @@ void data<T,D>::change_processor_recv (const int newproc, void* const mem)
_storage = (T*)mem;
}
+ const double wtime1 = MPI_Wtime();
T dummy;
MPI_Irecv (_storage, this->_size, dist::datatype(dummy), this->_proc,
this->tag, dist::comm, &request);
+ const double wtime2 = MPI_Wtime();
+ this->wtime_irecv += wtime2 - wtime1;
} else if (rank == this->_proc) {
// copy to other processor
@@ -197,9 +209,12 @@ void data<T,D>::change_processor_send (const int newproc, void* const mem)
assert (!mem);
assert (_storage);
+ const double wtime1 = MPI_Wtime();
T dummy;
MPI_Isend (_storage, this->_size, dist::datatype(dummy), newproc,
this->tag, dist::comm, &request);
+ const double wtime2 = MPI_Wtime();
+ this->wtime_isend += wtime2 - wtime1;
} else {
assert (!mem);
@@ -227,8 +242,11 @@ void data<T,D>::change_processor_wait (const int newproc, void* const mem)
if (rank == newproc) {
// copy from other processor
+ const double wtime1 = MPI_Wtime();
MPI_Status status;
MPI_Wait (&request, &status);
+ const double wtime2 = MPI_Wtime();
+ this->wtime_irecvwait += wtime2 - wtime1;
} else if (rank == this->_proc) {
// copy to other processor
@@ -236,8 +254,11 @@ void data<T,D>::change_processor_wait (const int newproc, void* const mem)
assert (!mem);
assert (_storage);
+ const double wtime1 = MPI_Wtime();
MPI_Status status;
MPI_Wait (&request, &status);
+ const double wtime2 = MPI_Wtime();
+ this->wtime_isendwait += wtime2 - wtime1;
if (this->_owns_storage) {
delete [] _storage;
@@ -553,6 +574,15 @@ extern "C" {
const int srcbbox[3][3],
const int dstbbox[3][3],
const int regbbox[3][3]);
+ void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_rf2)
+ (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+ const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+ 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_2tl_o3)
(const CCTK_REAL8* src1, const CCTK_REAL8& t1,
const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -562,6 +592,15 @@ extern "C" {
const int srcbbox[3][3],
const int dstbbox[3][3],
const int regbbox[3][3]);
+ void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2)
+ (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+ const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+ 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_2tl_minmod)
(const CCTK_REAL8* src1, const CCTK_REAL8& t1,
const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -651,6 +690,8 @@ void data<CCTK_REAL8,3>
const int order_space,
const int order_time)
{
+ const CCTK_REAL eps = 1.0e-10;
+
assert (has_storage());
assert (all(box.lower()>=extent().lower()));
assert (all(box.upper()<=extent().upper()));
@@ -694,8 +735,46 @@ void data<CCTK_REAL8,3>
regbbox[2][d] = box.stride()[d];
}
+ // Check that the times are consistent
+ assert (times.size() > 0);
+ CCTK_REAL min_time = times[0];
+ CCTK_REAL max_time = times[0];
+ for (size_t tl=1; tl<times.size(); ++tl) {
+ min_time = min(min_time, times[tl]);
+ max_time = max(max_time, times[tl]);
+ }
+ if (time < min_time - eps || time > max_time + eps) {
+ ostringstream buf;
+ buf << "Internal error: extrapolation in time."
+ << " time=" << time
+ << " times=" << times;
+ CCTK_WARN (0, buf.str().c_str());
+ }
+
+ // Is it necessary to interpolate in time?
+ if (times.size() > 1) {
+ for (size_t tl=0; tl<times.size(); ++tl) {
+ assert (abs(1.5) > 1.4);
+ if (abs(times[tl] - time) < eps) {
+ // It is not.
+ vector<const gdata<3>*> my_gsrcs(1);
+ vector<CCTK_REAL> my_times(1);
+ my_gsrcs[0] = gsrcs[tl];
+ my_times[0] = times[tl];
+ const int my_order_time = 0;
+ this->interpolate_from_innerloop
+ (my_gsrcs, my_times, box, time, order_space, my_order_time);
+ return;
+ }
+ }
+ }
+
assert (all(dext.stride() == box.stride()));
if (all(sext.stride() < dext.stride())) {
+ // Restrict
+
+ assert (times.size() == 1);
+ assert (abs(times[0] - time) < eps);
switch (transport_operator) {
@@ -722,9 +801,10 @@ void data<CCTK_REAL8,3>
default:
assert (0);
- } // switch (prolong_method)
+ } // switch (transport_operator)
} else if (all(sext.stride() > dext.stride())) {
+ // Prolongate
switch (transport_operator) {
@@ -732,6 +812,8 @@ void data<CCTK_REAL8,3>
switch (order_time) {
case 0:
+ assert (times.size() == 1);
+ assert (abs(times[0] - time) < eps);
assert (srcs.size()>=1);
switch (order_space) {
case 0:
@@ -789,23 +871,43 @@ void data<CCTK_REAL8,3>
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);
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_2tl_rf2)
+ ((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);
+ } else {
+ 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);
+ if (all (sext.stride() == dext.stride() * 2)) {
+ CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2)
+ ((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);
+ } else {
+ 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:
@@ -893,40 +995,75 @@ void data<CCTK_REAL8,3>
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);
+ assert (times.size() == 1);
+ assert (abs(times[0] - time) < eps);
+ switch (order_space) {
+ case 0:
+ case 1:
+ CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+ break;
+ case 2:
+ case 3:
+ 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;
+ default:
+ assert (0);
+ }
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);
+ switch (order_space) {
+ case 0:
+ case 1:
+ CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+ break;
+ case 2:
+ case 3:
+ 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);
+ break;
+ default:
+ assert (0);
+ }
break;
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);
+ switch (order_space) {
+ case 0:
+ case 1:
+ CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+ break;
+ case 2:
+ case 3:
+ 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 (order_time)
break;
default:
assert(0);
- } // switch (prolong_method)
+ } // switch (transport_operator)
} else {
assert (0);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-02 01:22:45 +0000