Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lagrange3oprolong #2

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Carpet/src/Cycle.cc
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ void CycleTimeLevels(cGH *const cctkGH) {
if (activetimelevels > 1) {
switch (groupdata.AT(group).transport_operator) {
case op_Lagrange:
case op_Lagrange_third_order_prolong:
case op_ENO:
case op_WENO:
case op_Lagrange_monotone:
Expand Down
2 changes: 2 additions & 0 deletions Carpet/src/SetupGH.cc
Original file line number Diff line number Diff line change
Expand Up @@ -2196,6 +2196,8 @@ operator_type get_transport_operator(cGH const *const cctkGH, int const group,
return op_copy;
} else if (CCTK_Equals(prolong_string, "Lagrange")) {
return op_Lagrange;
} else if (CCTK_Equals(prolong_string, "Lagrange_third_order_prolong")) {
return op_Lagrange_third_order_prolong;
} else if (CCTK_Equals(prolong_string, "ENO")) {
return op_ENO;
} else if (CCTK_Equals(prolong_string, "WENO")) {
Expand Down
68 changes: 65 additions & 3 deletions CarpetLib/src/data.cc
Original file line number Diff line number Diff line change
Expand Up @@ -584,7 +584,7 @@ void data<T>::transfer_prolongate(data const *const src, ibbox const &dstbox,
#if CARPET_DIM == 3

switch (transport_operator) {

case op_copy:
case op_Lagrange: {
static Timer timer("prolongate_Lagrange");
Expand Down Expand Up @@ -659,6 +659,64 @@ void data<T>::transfer_prolongate(data const *const src, ibbox const &dstbox,
break;
}

case op_Lagrange_third_order_prolong: {
static Timer timer("prolongate_Lagrange_3o_prolong");
timer.start();
// enum centering { vertex_centered, cell_centered };
switch (cent) {
case vertex_centered: {
static void (*the_operators[])(
T const *restrict const src, ivect3 const &restrict srcpadext,
ivect3 const &restrict srcext, T *restrict const dst,
ivect3 const &restrict dstpadext, ivect3 const &restrict dstext,
ibbox3 const &restrict srcbbox, ibbox3 const &restrict dstbbox,
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eschnett says

This chooses third order interpolation even if the overall order is set to be first order. This will bypass the check for sufficient number of ghost zones. The operator’s check will still catch this, but with an error message that is difficult to interpret.

I suggest changing this, so that the the_operators array is set up to contain the actual desired operators, and then using the respective array element here.

Otherwise there is no point in setting up an array, and you can simple pass &prolongate_3d_rf2<T, 3> as argument here.

ibbox3 const &restrict srcregbbox, ibbox3 const &restrict dstregbbox,
void *const extraargs) = {
NULL, &prolongate_3d_rf2<T, 1>, NULL, &prolongate_3d_rf2<T, 3>,
NULL, &prolongate_3d_rf2<T, 5>, NULL, &prolongate_3d_rf2<T, 7>,
NULL, &prolongate_3d_rf2<T, 9>, NULL, &prolongate_3d_rf2<T, 11>,
};
call_operator<T>(
the_operators[3], static_cast<T const *>(src->storage()),
src->padded_shape(), src->shape(), static_cast<T *>(this->storage()),
this->padded_shape(), this->shape(), src->extent(), this->extent(),
srcbox, dstbox, NULL);
break;
}
case cell_centered: {
if (use_dgfe) {
// Don't use call_operator, because we parallelise ourselves
prolongate_3d_dgfe_rf2<T, 5>(
static_cast<T const *>(src->storage()), src->padded_shape(),
src->shape(), static_cast<T *>(this->storage()),
this->padded_shape(), this->shape(), src->extent(), this->extent(),
srcbox, dstbox, NULL);
break;
}
static void (*the_operators[])(
T const *restrict const src, ivect3 const &restrict srcpadext,
ivect3 const &restrict srcext, T *restrict const dst,
ivect3 const &restrict dstpadext, ivect3 const &restrict dstext,
ibbox3 const &restrict srcbbox, ibbox3 const &restrict dstbbox,
ibbox3 const &restrict srcregbbox, ibbox3 const &restrict dstregbbox,
void *const extraargs) = {
&prolongate_3d_cc_rf2<T, 0>, &prolongate_3d_cc_rf2<T, 1>,
&prolongate_3d_cc_rf2<T, 2>, &prolongate_3d_cc_rf2<T, 3>,
&prolongate_3d_cc_rf2<T, 4>, &prolongate_3d_cc_rf2<T, 5>};
call_operator<T>(
the_operators[3], static_cast<T const *>(src->storage()),
src->padded_shape(), src->shape(), static_cast<T *>(this->storage()),
this->padded_shape(), this->shape(), src->extent(), this->extent(),
srcbox, dstbox, NULL);
break;
}
default:
assert(0);
}
timer.stop(0);
break;
}

case op_ENO: {
static Timer timer("prolongate_ENO");
timer.start();
Expand Down Expand Up @@ -932,6 +990,7 @@ void data<T>::transfer_prolongate(data const *const src, ibbox const &dstbox,
switch (transport_operator) {

case op_copy:
case op_Lagrange_third_order_prolong:
case op_Lagrange: {
static Timer timer("prolongate_Lagrange");
timer.start();
Expand All @@ -948,7 +1007,7 @@ void data<T>::transfer_prolongate(data const *const src, ibbox const &dstbox,
break;
default:
CCTK_ERROR("There is no vertex-centred stencil for op=\"LAGRANGE\" "
"with order_space not in {1}");
"with order_space not in {1}, and CARPET_DIM==4.");
break;
}
break;
Expand Down Expand Up @@ -1029,6 +1088,7 @@ void data<T>::transfer_restrict(data const *const src, ibbox const &dstregbox,

case op_copy:
case op_Lagrange:
case op_Lagrange_third_order_prolong:
case op_ENO:
case op_WENO:
case op_TVD:
Expand Down Expand Up @@ -1144,6 +1204,7 @@ void data<T>::transfer_restrict(data const *const src, ibbox const &dstregbox,

case op_copy:
case op_Lagrange:
case op_Lagrange_third_order_prolong:
// enum centering { vertex_centered, cell_centered };
switch (cent) {
case vertex_centered:
Expand Down Expand Up @@ -1196,7 +1257,8 @@ void data<T>::time_interpolate(vector<data *> const &srcs, ibbox const &dstbox,
case op_STAGGER101:
case op_STAGGER110:
case op_STAGGER111:
case op_Lagrange: {
case op_Lagrange:
case op_Lagrange_third_order_prolong: {
static Timer timer("time_interpolate_Lagrange");
timer.start();
switch (order_time) {
Expand Down
6 changes: 6 additions & 0 deletions CarpetLib/src/operators.hh
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,12 @@ enum operator_type {
op_copy, // use simple copying for prolongation
// (needs only one time level)
op_Lagrange, // Lagrange interpolation (standard)
op_Lagrange_third_order_prolong, // Force third-order Lagrange interpolation:
// Useful for hydro/GRMHD evolved variables,
// where default, higher-order Lagrange
// (used for spacetime fields) may
// induce spurious oscillations in
// hydro/GRMHD fields.
op_ENO, // use ENO stencils (for hydro)
op_WENO, // use WENO stencils (for hydro)
op_TVD, // use TVD stencils (for hydro)
Expand Down