From de8c64ae25d8e65f7d521ac19efb7302e0e854a8 Mon Sep 17 00:00:00 2001 From: schnetter <> Date: Wed, 21 Jan 2004 13:25:00 +0000 Subject: Check the stencil size before calling the minmod prolongation Check the stencil size before calling the minmod prolongation operators. darcs-hash:20040121132509-07bb3-7eb951d88ae3b66a3ad2338fecc632c8c77cb312.gz --- Carpet/CarpetLib/src/data.cc | 229 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 183 insertions(+), 46 deletions(-) diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 468602265..61aa92f20 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.36 2003/11/21 13:55:46 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.37 2004/01/21 14:25:09 schnetter Exp $ #include #include +#include +#include +#include #include +#include #include +#include #include @@ -34,16 +39,17 @@ static int nexttag () // Constructors template -data::data (const int varindex_) - : gdata(varindex_), +data::data (const int varindex_, const operator_type transport_operator_) + : gdata(varindex_, transport_operator_), _storage(0), comm_active(false), tag(nexttag()) { } template -data::data (const int varindex_, const ibbox& extent_, const int proc_) - : gdata(varindex_), +data::data (const int varindex_, const operator_type transport_operator_, + const ibbox& extent_, const int proc_) + : gdata(varindex_, transport_operator_), _storage(0), comm_active(false), tag(nexttag()) @@ -59,8 +65,11 @@ data::~data () { // Pseudo constructors template -data* data::make_typed (const int varindex_) const { - return new data(varindex_); +data* data::make_typed (const int varindex_, + const operator_type transport_operator_) + const +{ + return new data(varindex_, transport_operator_); } @@ -107,7 +116,7 @@ void data::transfer_from (gdata* 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::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::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::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::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 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 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 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 1.4); + if (abs(times[tl] - time) < eps) { + // It is not. + vector*> my_gsrcs(1); + vector 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 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 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 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 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); -- cgit v1.2.3