From 275e279967be343c0217336be792c746f0f851c9 Mon Sep 17 00:00:00 2001 From: Zachariah Etienne Date: Thu, 23 Aug 2018 16:20:50 -0400 Subject: [PATCH 1/2] Carpet/CarpetLib: Add Lagrange_third_order_prolong prolongation option. --- Carpet/src/Cycle.cc | 1 + Carpet/src/SetupGH.cc | 2 ++ CarpetLib/src/data.cc | 68 ++++++++++++++++++++++++++++++++++++-- CarpetLib/src/operators.hh | 6 ++++ 4 files changed, 74 insertions(+), 3 deletions(-) diff --git a/Carpet/src/Cycle.cc b/Carpet/src/Cycle.cc index 25fd2c790..4da6d44d8 100644 --- a/Carpet/src/Cycle.cc +++ b/Carpet/src/Cycle.cc @@ -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: diff --git a/Carpet/src/SetupGH.cc b/Carpet/src/SetupGH.cc index 2551271e1..9a3bb8e92 100644 --- a/Carpet/src/SetupGH.cc +++ b/Carpet/src/SetupGH.cc @@ -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")) { diff --git a/CarpetLib/src/data.cc b/CarpetLib/src/data.cc index 69c539a1b..3b3406e13 100644 --- a/CarpetLib/src/data.cc +++ b/CarpetLib/src/data.cc @@ -584,7 +584,7 @@ void data::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"); @@ -659,6 +659,64 @@ void data::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, + ibbox3 const &restrict srcregbbox, ibbox3 const &restrict dstregbbox, + void *const extraargs) = { + NULL, &prolongate_3d_rf2, NULL, &prolongate_3d_rf2, + NULL, &prolongate_3d_rf2, NULL, &prolongate_3d_rf2, + NULL, &prolongate_3d_rf2, NULL, &prolongate_3d_rf2, + }; + call_operator( + the_operators[3], static_cast(src->storage()), + src->padded_shape(), src->shape(), static_cast(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( + static_cast(src->storage()), src->padded_shape(), + src->shape(), static_cast(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, &prolongate_3d_cc_rf2, + &prolongate_3d_cc_rf2, &prolongate_3d_cc_rf2, + &prolongate_3d_cc_rf2, &prolongate_3d_cc_rf2}; + call_operator( + the_operators[3], static_cast(src->storage()), + src->padded_shape(), src->shape(), static_cast(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(); @@ -932,6 +990,7 @@ void data::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(); @@ -948,7 +1007,7 @@ void data::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; @@ -1029,6 +1088,7 @@ void data::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: @@ -1144,6 +1204,7 @@ void data::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: @@ -1196,7 +1257,8 @@ void data::time_interpolate(vector 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) { diff --git a/CarpetLib/src/operators.hh b/CarpetLib/src/operators.hh index aa07e1c4b..84bf477c9 100644 --- a/CarpetLib/src/operators.hh +++ b/CarpetLib/src/operators.hh @@ -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) From 15a4c4eed735e65b8636264d397cf33a630e1493 Mon Sep 17 00:00:00 2001 From: Zachariah Etienne Date: Thu, 23 Aug 2018 16:25:33 -0400 Subject: [PATCH 2/2] Carpet/src/SetupGH.cc: fix tiny typo --- Carpet/src/SetupGH.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Carpet/src/SetupGH.cc b/Carpet/src/SetupGH.cc index 9a3bb8e92..b201b7d5e 100644 --- a/Carpet/src/SetupGH.cc +++ b/Carpet/src/SetupGH.cc @@ -2197,7 +2197,7 @@ operator_type get_transport_operator(cGH const *const cctkGH, int const group, } 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:; + return op_Lagrange_third_order_prolong; } else if (CCTK_Equals(prolong_string, "ENO")) { return op_ENO; } else if (CCTK_Equals(prolong_string, "WENO")) {