diff --git a/res/cmake/test.cmake b/res/cmake/test.cmake index 1f7331b1d..9b6e1653c 100644 --- a/res/cmake/test.cmake +++ b/res/cmake/test.cmake @@ -8,6 +8,7 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests add_subdirectory(tests/core/data/ndarray) + add_subdirectory(tests/core/data/field) add_subdirectory(tests/core/data/grid) add_subdirectory(tests/core/data/gridlayout) add_subdirectory(tests/core/data/vecfield) diff --git a/src/amr/data/field/refine/electric_field_refiner.hpp b/src/amr/data/field/refine/electric_field_refiner.hpp index 89debd3f4..1abac57e9 100644 --- a/src/amr/data/field/refine/electric_field_refiner.hpp +++ b/src/amr/data/field/refine/electric_field_refiner.hpp @@ -3,14 +3,14 @@ #include "core/def/phare_mpi.hpp" // IWYU pragma: keep - -#include - -#include "amr/resources_manager/amr_utils.hpp" #include "core/utilities/constants.hpp" #include "core/utilities/point/point.hpp" #include "core/data/grid/gridlayoutdefs.hpp" +#include "amr/resources_manager/amr_utils.hpp" + +#include + #include namespace PHARE::amr @@ -37,22 +37,23 @@ class ElectricFieldRefiner // electric field refinement strategy follows // fujimoto et al. 2011 : doi:10.1016/j.jcp.2011.08.002 - template - void operator()(FieldT const& coarseField, FieldT& fineField, - core::Point fineIndex) + void operator()(auto const& coarseField, auto& fineField, auto const& fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { TBOX_ASSERT(coarseField.physicalQuantity() == fineField.physicalQuantity()); - auto const locFineIdx = AMRToLocal(fineIndex, fineBox_); - auto const coarseIdx = toCoarseIndex(fineIndex); - auto const locCoarseIdx = AMRToLocal(coarseIdx, coarseBox_); + if (not std::isnan(fineVal)) + return; // KEEP! if constexpr (dimension == 1) - refine1D_(coarseField, fineField, locFineIdx, locCoarseIdx); + fineVal = coarseVal; else if constexpr (dimension == 2) - refine2D_(coarseField, fineField, fineIndex, locFineIdx, locCoarseIdx); + refine2D_(coarseField, fineField, fineIndex, coarseIndex, locFineIdx, locCoarseIdx, + fineVal, coarseVal); else if constexpr (dimension == 3) - refine3D_(coarseField, fineField, fineIndex, locFineIdx, locCoarseIdx); + refine3D_(coarseField, fineField, fineIndex, coarseIndex, locFineIdx, locCoarseIdx, + fineVal, coarseVal); } private: @@ -73,44 +74,13 @@ class ElectricFieldRefiner } - template - void refine1D_(FieldT const& coarseField, FieldT& fineField, - core::Point const& locFineIdx, - core::Point const& locCoarseIdx) - { - // if dual, then Ex - // if even fine index, we're on top of coarse, we take 100% coarse overlapped fieldValue - // e.g. fineIndex==100, we take coarse[100/2] - // e.g. fineIndex==101 we take coarse[101/2] = coarse[50] as well - // - // 49 50 51 - // o x o x o x o Ex on coarse - // o x o x o x o x o x o x o Ex on fine - // 98 99 100 101 102 103 - // - // - // if primal, then Ey,Ez we just copy the coarse value - // since they are on top of the fine one. - // - // therefore in all cases in 1D we just copy the coarse value - // - if (std::isnan(fineField(locFineIdx[dirX]))) - fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]); - } - - template - void refine2D_(FieldT const& coarseField, FieldT& fineField, - core::Point const& fineIndex, - core::Point const& locFineIdx, - core::Point const& locCoarseIdx) + void refine2D_(auto const& coarseField, auto& fineField, auto const& fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { // ilc: index local coarse - // ilf: index local fine auto const ilcx = locCoarseIdx[dirX]; auto const ilcy = locCoarseIdx[dirY]; - auto const ilfx = locFineIdx[dirX]; - auto const ilfy = locFineIdx[dirY]; - // this is Ex if (centerings_[dirX] == core::QtyCentering::dual @@ -120,16 +90,14 @@ class ElectricFieldRefiner { // we're on a fine edge shared with coarse mesh // take the coarse face value - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) = coarseField(ilcx, ilcy); + fineVal = coarseVal; } else { // we're on a fine edge in between two coarse edges // we take the average - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); + + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy + 1)); } } // Ey @@ -143,16 +111,15 @@ class ElectricFieldRefiner // both fine Ey e.g. at j=100 and 101 will take j=50 on coarse // so no need to look at whether jfine is even or odd // just take the value at the local coarse index - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) = coarseField(ilcx, ilcy); + + fineVal = coarseVal; } else { // we're on a fine edge in between two coarse ones // we take the average - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); + + fineVal = 0.5 * (coarseVal + coarseField(ilcx + 1, ilcy)); } } // and this is now Ez @@ -161,47 +128,34 @@ class ElectricFieldRefiner { if (onCoarseXFace_(fineIndex) and onCoarseYFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) = coarseField(ilcx, ilcy); + fineVal = coarseVal; } else if (onCoarseXFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy + 1)); } else if (onCoarseYFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx + 1, ilcy)); } else { - if (std::isnan(fineField(ilfx, ilfy))) - fineField(ilfx, ilfy) - = 0.25 - * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy) - + coarseField(ilcx, ilcy + 1) + coarseField(ilcx + 1, ilcy + 1)); + fineVal = 0.25 + * (coarseVal + coarseField(ilcx + 1, ilcy) + coarseField(ilcx, ilcy + 1) + + coarseField(ilcx + 1, ilcy + 1)); } } } - template - void refine3D_(FieldT const& coarseField, FieldT& fineField, - core::Point const& fineIndex, - core::Point const& locFineIdx, - core::Point const& locCoarseIdx) + void refine3D_(auto const& coarseField, auto& fineField, auto const& fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { // ilc: index local coarse - // ilf: index local fine auto const ilcx = locCoarseIdx[dirX]; auto const ilcy = locCoarseIdx[dirY]; auto const ilcz = locCoarseIdx[dirZ]; - auto const ilfx = locFineIdx[dirX]; - auto const ilfy = locFineIdx[dirY]; - auto const ilfz = locFineIdx[dirZ]; // Ex if (centerings_[dirX] == core::QtyCentering::dual @@ -212,34 +166,28 @@ class ElectricFieldRefiner // just copy the coarse value if (onCoarseYFace_(fineIndex) and onCoarseZFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz); + fineVal = coarseVal; } // we share the Y face but not the Z face // we must be one of the 2 X fine edges on a Y face // thus we take the average of the two surrounding edges at Z and Z+DZ else if (onCoarseYFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy, ilcz + 1)); } // we share a Z face but not the Y face // we must be one of the 2 X fine edges on a Z face // we thus take the average of the two X edges at y and y+dy else if (onCoarseZFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy + 1, ilcz)); } else { // we don't share any face thus we're on one of the 2 middle X edges // we take the average of the 4 surrounding X averages - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.25 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)) + + fineVal = 0.25 * (coarseVal + coarseField(ilcx, ilcy + 1, ilcz)) + 0.25 * (coarseField(ilcx, ilcy, ilcz + 1) + coarseField(ilcx, ilcy + 1, ilcz + 1)); @@ -254,8 +202,8 @@ class ElectricFieldRefiner if (onCoarseXFace_(fineIndex) and onCoarseZFace_(fineIndex)) { // we thus just copy the coarse value - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz); + + fineVal = coarseVal; } // now we only have same X face, but not (else) the Z face // so we're a new fine Y edge in between two coarse Y edges @@ -267,28 +215,23 @@ class ElectricFieldRefiner // this means we are on a Y edge that lies in between 2 coarse edges // at z and z+dz // take the average of these 2 coarse value - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); + + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy, ilcz + 1)); } // we're on a Z coarse face, but not on a X coarse face // we thus must be one of the 2 Y edges on a Z face // and thus we take the average of the 2 Y edges at X and X+dX else if (onCoarseZFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx + 1, ilcy, ilcz)); } // now we're not on any of the coarse faces // so we must be one of the two Y edge in the middle of the cell // we thus average over the 4 Y edges of the coarse cell else { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.25 - * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz) + fineVal = 0.25 + * (coarseVal + coarseField(ilcx + 1, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1) + coarseField(ilcx + 1, ilcy, ilcz + 1)); } @@ -302,36 +245,29 @@ class ElectricFieldRefiner // we thus copy the coarse value if (onCoarseXFace_(fineIndex) and onCoarseYFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz); + fineVal = coarseVal; } // here we're on a coarse X face, but not a Y face // we must be 1 of the 2 Z edges on a X face // thus we average the 2 surrounding Z coarse edges at Y and Y+dY else if (onCoarseXFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx, ilcy + 1, ilcz)); } // here we're on a coarse Y face, but not a X face // we must be 1 of the 2 Z edges on a Y face // thus we average the 2 surrounding Z coarse edges at X and X+dX else if (onCoarseYFace_(fineIndex)) { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); + fineVal = 0.5 * (coarseVal + coarseField(ilcx + 1, ilcy, ilcz)); } // we're not on any coarse face thus must be one of the 2 Z edges // in the middle of the coarse cell // we therefore take the average of the 4 surrounding Z edges else { - if (std::isnan(fineField(ilfx, ilfy, ilfz))) - fineField(ilfx, ilfy, ilfz) - = 0.25 - * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz) + fineVal = 0.25 + * (coarseVal + coarseField(ilcx + 1, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz) + coarseField(ilcx + 1, ilcy + 1, ilcz)); } diff --git a/src/amr/data/field/refine/field_moments_refiner.hpp b/src/amr/data/field/refine/field_moments_refiner.hpp index c5b04f442..66bf3dc19 100644 --- a/src/amr/data/field/refine/field_moments_refiner.hpp +++ b/src/amr/data/field/refine/field_moments_refiner.hpp @@ -4,8 +4,7 @@ #include "core/def/phare_mpi.hpp" // IWYU pragma: keep -#include "core/data/field/field.hpp" -#include "core/utilities/constants.hpp" + #include "core/utilities/point/point.hpp" #include "core/data/grid/gridlayoutdefs.hpp" @@ -14,7 +13,6 @@ #include #include -#include namespace PHARE::amr diff --git a/src/amr/data/field/refine/field_refine_operator.hpp b/src/amr/data/field/refine/field_refine_operator.hpp index 7e5df80b5..ab9fddbfb 100644 --- a/src/amr/data/field/refine/field_refine_operator.hpp +++ b/src/amr/data/field/refine/field_refine_operator.hpp @@ -1,17 +1,17 @@ #ifndef PHARE_FIELD_REFINE_OPERATOR_HPP #define PHARE_FIELD_REFINE_OPERATOR_HPP - - -#include "core/def/phare_mpi.hpp" // IWYU pragma: keep - #include "core/def.hpp" +#include "core/def/phare_mpi.hpp" // IWYU pragma: keep +#include "core/utilities/box/box_span.hpp" #include "amr/data/field/field_data.hpp" +#include "amr/resources_manager/amr_utils.hpp" #include "amr/data/tensorfield/tensor_field_data.hpp" -#include "field_linear_refine.hpp" + +// #include "field_linear_refine.hpp" #include #include @@ -27,13 +27,91 @@ using core::dirX; using core::dirY; using core::dirZ; +template +struct CrossLevelIndices +{ + auto constexpr static dim = Field_t::dimension; + using FieldRow_t = core::FieldBoxPointSpans; + using FieldRow_ct = core::FieldBoxPointSpans; + + core::FieldBox& dst; + core::FieldBox const& src; + core::FieldBoxSpan dst_lcl_span + = core::make_field_box_point_span(dst.lcl_box, dst.field); + core::FieldBoxSpan src_lcl_span + = core::make_field_box_point_span(src.lcl_box, src.field); + + core::BoxSpan dst_amr_span = core::make_box_span(dst.amr_box); + core::BoxSpan src_amr_span = core::make_box_span(src.amr_box); +}; -template -void refine_field(Dst& destinationField, auto& sourceField, auto& intersectionBox, auto& refiner) + + + +template +void refine_field(core::FieldBox& dst, core::FieldBox& src, auto& refiner) { - for (auto const bix : phare_box_from(intersectionBox)) - refiner(sourceField, destinationField, bix); + auto constexpr static IDX = Field_t::dimension - 1; + + CrossLevelIndices indices{dst, src}; + + auto d_f_slabs = indices.dst_lcl_span.begin(); + auto d_b_slabs = indices.dst_amr_span.begin(); + auto s_f_slabs = indices.src_lcl_span.begin(); + auto s_b_slabs = indices.src_amr_span.begin(); + + auto slab_idx = indices.dst_amr_span.slab_begin(); + + for (; d_f_slabs != indices.dst_lcl_span.end(); ++d_f_slabs, ++d_b_slabs, ++slab_idx) + { + auto d_f_spans = d_f_slabs.begin(); + auto s_f_spans = s_f_slabs.begin(); + auto d_b_spans = d_b_slabs.begin(); + auto s_b_spans = s_b_slabs.begin(); + + auto span_idx = d_b_slabs.span_begin(); + for (; d_f_spans != d_f_slabs.end(); ++d_f_spans, ++d_b_spans, ++span_idx) + { + auto const& [d_amr_point, d_size] = *d_b_spans; + auto const& [s_amr_point, s_size] = *s_b_spans; + + auto&& [d_span, d_lcl_point] = *d_f_spans; + auto&& [s_span, s_lcl_point] = *s_f_spans; + + std::size_t dst_idx = 0; + std::size_t src_idx = 0; + for (; dst_idx < d_span.size(); ++dst_idx) + { + assert(s_amr_point == toCoarseIndex(d_amr_point)); + + refiner(src.field, dst.field, d_amr_point, s_amr_point, d_lcl_point, s_lcl_point, + d_span[dst_idx], s_span[src_idx]); + + if (d_amr_point[IDX] % 2 != 0) + { + ++src_idx; + ++s_lcl_point[IDX]; + ++s_amr_point[IDX]; + } + + ++d_amr_point[IDX]; + ++d_lcl_point[IDX]; + } + + if (span_idx % 2 != 0) + { + ++s_f_spans; + ++s_b_spans; + } + } + + if (slab_idx % 2 != 0) + { + ++s_f_slabs; + ++s_b_slabs; + } + } } @@ -112,8 +190,10 @@ class FieldRefineOperator : public SAMRAI::hier::RefineOperator { // we compute the intersection with the destination, // and then we apply the refine operation on each fine index. - auto intersectionBox = destFieldBox * box; - refine_field(destinationField, sourceField, intersectionBox, refiner); + auto const dst_overlap = phare_box_from(destFieldBox * box); + core::FieldBox dst{destinationField, destLayout, dst_overlap}; + core::FieldBox src{sourceField, srcLayout, coarsen_box(dst_overlap)}; + refine_field(dst, src, refiner); } } }; @@ -178,6 +258,8 @@ class TensorFieldRefineOperator : public SAMRAI::hier::RefineOperator auto const& sourceFields = TensorFieldDataT::getFields(source, sourceId); auto const& srcLayout = TensorFieldDataT::getLayout(source, sourceId); + PHARE_LOG_SCOPE(2, "TensorFieldRefineOperator::refine::" + sourceFields[0].name()); + // We assume that quantity are all the same. // Note that an assertion will be raised in refineIt operator for (std::uint16_t c = 0; c < N; ++c) @@ -199,8 +281,13 @@ class TensorFieldRefineOperator : public SAMRAI::hier::RefineOperator { // we compute the intersection with the destination, // and then we apply the refine operation on each fine index. - auto const intersectionBox = destFieldBox * box; - refine_field(destinationFields[c], sourceFields[c], intersectionBox, refiner); + auto const dst_overlap = phare_box_from(destFieldBox * box); + core::FieldBox dst{destinationFields[c], destLayout, dst_overlap}; + + auto const src_overlap + = *(coarsen_box(dst_overlap) * srcLayout.AMRGhostBoxFor(sourceFields[c])); + core::FieldBox src{sourceFields[c], srcLayout, src_overlap}; + refine_field(dst, src, refiner); } } } diff --git a/src/amr/data/field/refine/field_refiner.hpp b/src/amr/data/field/refine/field_refiner.hpp index f4e74fc4a..229d232cf 100644 --- a/src/amr/data/field/refine/field_refiner.hpp +++ b/src/amr/data/field/refine/field_refiner.hpp @@ -43,8 +43,8 @@ namespace amr } - /** @brief Given a sourceField , a destinationField, and a fineIndex compute the - * interpolation from the coarseField(sourceField) to the fineFiled(destinationField) at the + /** @brief Given a coarseField , a fineField, and a fineIndex compute the + * interpolation from the coarseField(coarseField) to the fineFiled(fineField) at the * fineIndex index * * @@ -55,11 +55,15 @@ namespace amr * - we just have to know which one to use, depending on where the fineIndex is in the * coarse cell */ - template - void operator()(FieldT const& sourceField, FieldT& destinationField, - core::Point fineIndex) + void operator()(auto const& coarseField, auto& fineField, auto fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { - TBOX_ASSERT(sourceField.physicalQuantity() == destinationField.physicalQuantity()); + TBOX_ASSERT(coarseField.physicalQuantity() == fineField.physicalQuantity()); + + + if (not std::isnan(fineVal)) + return; // First we get the coarseStartIndex for a given fineIndex // then we get the index in weights table for a given fineIndex. @@ -75,29 +79,23 @@ namespace amr coarseStartIndex = AMRToLocal(coarseStartIndex, coarseBox_); fineIndex = AMRToLocal(fineIndex, fineBox_); - double fieldValue = 0.; - + assert(coarseField(coarseStartIndex) == coarseVal); + fineVal = 0.; if constexpr (dimension == 1) { - auto const& xStartIndex = coarseStartIndex[dirX]; - + auto const& xStartIndex = coarseStartIndex[dirX]; auto const& xWeights = indexesAndWeights_.weights(core::Direction::X); auto const& leftRightWeights = xWeights[iWeight[dirX]]; for (std::size_t iShiftX = 0; iShiftX < leftRightWeights.size(); ++iShiftX) - { - fieldValue += sourceField(xStartIndex + iShiftX) * leftRightWeights[iShiftX]; - } - if (std::isnan(destinationField(fineIndex[dirX]))) - destinationField(fineIndex[dirX]) = fieldValue; + fineVal += coarseField(xStartIndex + iShiftX) * leftRightWeights[iShiftX]; } - else if constexpr (dimension == 2) { auto const& xStartIndex = coarseStartIndex[dirX]; @@ -114,14 +112,11 @@ namespace amr double Yinterp = 0.; for (std::size_t iShiftY = 0; iShiftY < yLeftRightWeights.size(); ++iShiftY) { - Yinterp += sourceField(xStartIndex + iShiftX, yStartIndex + iShiftY) + Yinterp += coarseField(xStartIndex + iShiftX, yStartIndex + iShiftY) * yLeftRightWeights[iShiftY]; } - fieldValue += Yinterp * xLeftRightWeights[iShiftX]; + fineVal += Yinterp * xLeftRightWeights[iShiftX]; } - - if (std::isnan(destinationField(fineIndex[dirX], fineIndex[dirY]))) - destinationField(fineIndex[dirX], fineIndex[dirY]) = fieldValue; } @@ -150,18 +145,14 @@ namespace amr double Zinterp = 0.; for (std::size_t iShiftZ = 0; iShiftZ < zLeftRightWeights.size(); ++iShiftZ) { - Zinterp += sourceField(xStartIndex + iShiftX, yStartIndex + iShiftY, + Zinterp += coarseField(xStartIndex + iShiftX, yStartIndex + iShiftY, zStartIndex + iShiftZ) * zLeftRightWeights[iShiftZ]; } Yinterp += Zinterp * yLeftRightWeights[iShiftY]; } - fieldValue += Yinterp * xLeftRightWeights[iShiftX]; + fineVal += Yinterp * xLeftRightWeights[iShiftX]; } - - if (std::isnan(destinationField(fineIndex[dirX], fineIndex[dirY], fineIndex[dirZ]))) - destinationField(fineIndex[dirX], fineIndex[dirY], fineIndex[dirZ]) - = fieldValue; } } diff --git a/src/amr/data/field/refine/magnetic_field_init_refiner.hpp b/src/amr/data/field/refine/magnetic_field_init_refiner.hpp index c1fce019b..45c9fcb9c 100644 --- a/src/amr/data/field/refine/magnetic_field_init_refiner.hpp +++ b/src/amr/data/field/refine/magnetic_field_init_refiner.hpp @@ -48,16 +48,13 @@ class MagneticFieldInitRefiner // onto the 2 (1D), 4 (2/3D) colocated fine faces. This way the total flux on // these fine faces equals that on the overlaped coarse face. // see fujimoto et al. 2011 : doi:10.1016/j.jcp.2011.08.002 - template - void operator()(FieldT const& coarseField, FieldT& fineField, - core::Point fineIndex) + + void operator()(auto const& coarseField, auto& fineField, auto const& fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { TBOX_ASSERT(coarseField.physicalQuantity() == fineField.physicalQuantity()); - auto locFineIdx = AMRToLocal(fineIndex, fineBox_); - auto coarseIdx = toCoarseIndex(fineIndex); - auto locCoarseIdx = AMRToLocal(coarseIdx, coarseBox_); - if constexpr (dimension == 1) { @@ -76,10 +73,9 @@ class MagneticFieldInitRefiner if (centerings_[0] == core::QtyCentering::primal) { if (fineIndex[0] % 2 == 0) - { - fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]); - } + fineVal = coarseVal; } + // dual case, By, Bz // 49 50 51 // o + o + o + o Byz on coarse : + @@ -89,9 +85,7 @@ class MagneticFieldInitRefiner // 100 takes 50 = 100/2 // 101 takes 50 = 101/2 else - { - fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]); - } + fineVal = coarseVal; } @@ -107,8 +101,7 @@ class MagneticFieldInitRefiner { // we're on a coarse X face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -119,8 +112,7 @@ class MagneticFieldInitRefiner { // we're on a coarse Y face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -129,18 +121,13 @@ class MagneticFieldInitRefiner // Bz // we're always on a coarse Z face since there is no dual in z // all 4 fine Bz take the coarse Z value - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if constexpr (dimension == 3) { - auto ix = locCoarseIdx[dirX]; - auto iy = locCoarseIdx[dirY]; - auto iz = locCoarseIdx[dirZ]; - if (centerings_[dirX] == core::QtyCentering::primal and centerings_[dirY] == core::QtyCentering::dual and centerings_[dirZ] == core::QtyCentering::dual) @@ -150,8 +137,7 @@ class MagneticFieldInitRefiner { // we're on a coarse X face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -163,8 +149,7 @@ class MagneticFieldInitRefiner { // we're on a coarse Y face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -176,8 +161,7 @@ class MagneticFieldInitRefiner { // we're on a coarse X face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } } diff --git a/src/amr/data/field/refine/magnetic_field_refiner.hpp b/src/amr/data/field/refine/magnetic_field_refiner.hpp index 55e7cc970..ccf0b8c63 100644 --- a/src/amr/data/field/refine/magnetic_field_refiner.hpp +++ b/src/amr/data/field/refine/magnetic_field_refiner.hpp @@ -2,7 +2,8 @@ #define PHARE_MAGNETIC_FIELD_REFINER_HPP -#include "core/def/phare_mpi.hpp" + +#include "core/def/phare_mpi.hpp" // IWYU pragma: keep #include "core/utilities/constants.hpp" #include "core/utilities/point/point.hpp" #include "core/data/grid/gridlayoutdefs.hpp" @@ -45,20 +46,14 @@ class MagneticFieldRefiner // onto the 2 (1D), 4 (2/3D) colocated fine faces. This way the total flux on // these fine faces equals that on the overlaped coarse face. // see fujimoto et al. 2011 : doi:10.1016/j.jcp.2011.08.002 - template - void operator()(FieldT const& coarseField, FieldT& fineField, - core::Point fineIndex) + void operator()(auto const& coarseField, auto& fineField, auto const& fineIndex, + auto const& coarseIndex, auto const& locFineIdx, auto const& locCoarseIdx, + auto& fineVal, auto const& coarseVal) { TBOX_ASSERT(coarseField.physicalQuantity() == fineField.physicalQuantity()); - using core::dirX; - using core::dirY; - using core::dirZ; - - auto const locFineIdx = AMRToLocal(fineIndex, fineBox_); - auto const coarseIdx = toCoarseIndex(fineIndex); - auto const locCoarseIdx = AMRToLocal(coarseIdx, coarseBox_); - + if (not std::isnan(fineVal)) // KEEP + return; if constexpr (dimension == 1) { @@ -76,10 +71,8 @@ class MagneticFieldRefiner // if (centerings_[0] == core::QtyCentering::primal) { - if (fineIndex[0] % 2 == 0 && std::isnan(fineField(locFineIdx[dirX]))) - { - fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]); - } + if (fineIndex[0] % 2 == 0) + fineVal = coarseVal; } // dual case, By, Bz // 49 50 51 @@ -90,41 +83,33 @@ class MagneticFieldRefiner // 100 takes 50 = 100/2 // 101 takes 50 = 101/2 else - { - if (std::isnan(fineField(locFineIdx[dirX]))) - fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]); - } + fineVal = coarseVal; } - else if constexpr (dimension == 2) { if (centerings_[dirX] == core::QtyCentering::primal and centerings_[dirY] == core::QtyCentering::dual) { // Bx - if (fineIndex[dirX] % 2 == 0 - && std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY]))) + if (fineIndex[dirX] % 2 == 0) { // we're on a coarse X face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual and centerings_[dirY] == core::QtyCentering::primal) { // By - if (fineIndex[dirY] % 2 == 0 - && std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY]))) + if (fineIndex[dirY] % 2 == 0) { // we're on a coarse Y face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -133,31 +118,23 @@ class MagneticFieldRefiner // Bz // we're always on a coarse Z face since there is no dual in z // all 4 fine Bz take the coarse Z value - if (std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY]))) - fineField(locFineIdx[dirX], locFineIdx[dirY]) - = coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]); + fineVal = coarseVal; } } else if constexpr (dimension == 3) { - auto ix = locCoarseIdx[dirX]; - auto iy = locCoarseIdx[dirY]; - auto iz = locCoarseIdx[dirZ]; - if (centerings_[dirX] == core::QtyCentering::primal and centerings_[dirY] == core::QtyCentering::dual and centerings_[dirZ] == core::QtyCentering::dual) { // Bx - if (fineIndex[dirX] % 2 == 0 - && std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]))) + if (fineIndex[dirX] % 2 == 0) { // we're on a coarse X face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -165,13 +142,11 @@ class MagneticFieldRefiner and centerings_[dirZ] == core::QtyCentering::dual) { // By - if (fineIndex[dirY] % 2 == 0 - && std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]))) + if (fineIndex[dirY] % 2 == 0) { // we're on a coarse Y face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -179,13 +154,11 @@ class MagneticFieldRefiner and centerings_[dirZ] == core::QtyCentering::primal) { // Bz - if (fineIndex[dirZ] % 2 == 0 - && std::isnan(fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]))) + if (fineIndex[dirZ] % 2 == 0) { // we're on a coarse Z face // take the coarse face value - fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = coarseField(ix, iy, iz); + fineVal = coarseVal; } } } diff --git a/src/amr/data/field/refine/magnetic_refine_patch_strategy.hpp b/src/amr/data/field/refine/magnetic_refine_patch_strategy.hpp index 6fc0d3f60..0854dbb32 100644 --- a/src/amr/data/field/refine/magnetic_refine_patch_strategy.hpp +++ b/src/amr/data/field/refine/magnetic_refine_patch_strategy.hpp @@ -1,17 +1,22 @@ #ifndef PHARE_AMR_MAGNETIC_REFINE_PATCH_STRATEGY_HPP #define PHARE_AMR_MAGNETIC_REFINE_PATCH_STRATEGY_HPP -#include "core/utilities/constants.hpp" +#include "core/data/field/field_box_span.hpp" +#include "core/logger.hpp" #include "core/utilities/types.hpp" +#include "core/utilities/constants.hpp" #include "amr/data/field/field_geometry.hpp" #include "amr/utilities/box/amr_box.hpp" +#include "amr/data/field/field_geometry.hpp" #include "amr/resources_manager/amr_utils.hpp" + #include "SAMRAI/xfer/RefinePatchStrategy.h" #include #include +#include namespace PHARE::amr { @@ -19,6 +24,9 @@ using core::dirX; using core::dirY; using core::dirZ; + + + template class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy { @@ -54,18 +62,20 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy } - void preprocessRefine(SAMRAI::hier::Patch& fine, SAMRAI::hier::Patch const& coarse, - SAMRAI::hier::Box const& fine_box, - SAMRAI::hier::IntVector const& ratio) override + void preprocessRefine(SAMRAI::hier::Patch& /*fine*/, SAMRAI::hier::Patch const& /*coarse*/, + SAMRAI::hier::Box const& /*fine_box*/, + SAMRAI::hier::IntVector const& /*ratio*/) override { } // We compute the values of the new fine magnetic faces using what was already refined, ie // the values on the old coarse faces. - void postprocessRefine(SAMRAI::hier::Patch& fine, SAMRAI::hier::Patch const& coarse, + void postprocessRefine(SAMRAI::hier::Patch& fine, SAMRAI::hier::Patch const& /*coarse*/, SAMRAI::hier::Box const& fine_box, - SAMRAI::hier::IntVector const& ratio) override + SAMRAI::hier::IntVector const& /*ratio*/) override { + PHARE_LOG_SCOPE(2, "MagneticRefinePatchStrategy::postprocessRefine"); + assertIDsSet(); auto& fields = TensorFieldDataT::getFields(fine, b_id_); @@ -77,8 +87,9 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy auto const fine_field_box = core::for_N_make_array([&](auto i) { using PhysicalQuantity = std::decay_t; - return FieldGeometry::toFieldBox( - fine_box, fields[i].physicalQuantity(), fineBoxLayout); + return phare_box_from( + FieldGeometry::toFieldBox( + fine_box, fields[i].physicalQuantity(), fineBoxLayout)); }); if constexpr (dimension == 1) @@ -86,7 +97,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy // if we ever go to c++23 we could use std::views::zip to iterate both on the local and // global indices instead of passing the box to do an amr to local inside the function, // which is not obvious at call site - for (auto const& i : phare_box_from(fine_field_box[dirX])) + for (auto const& i : fine_field_box[dirX]) { postprocessBx1d(bx, layout, i); } @@ -94,12 +105,12 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy else if constexpr (dimension == 2) { - for (auto const& i : phare_box_from(fine_field_box[dirX])) + for (auto const& i : fine_field_box[dirX]) { postprocessBx2d(bx, by, layout, i); } - for (auto const& i : phare_box_from(fine_field_box[dirY])) + for (auto const& i : fine_field_box[dirY]) { postprocessBy2d(bx, by, layout, i); } @@ -109,20 +120,20 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy { auto meshSize = layout.meshSize(); - for (auto const& i : phare_box_from(fine_field_box[dirX])) - { - postprocessBx3d(bx, by, bz, meshSize, layout, i); - } + for (auto const& slab : core::make_box_span(fine_field_box[dirX])) + for (auto const& [point, size] : slab) + for (std::size_t i = 0; i < size; ++i, ++point[dimension - 1]) + postprocessBx3d(bx, by, bz, meshSize, layout, point); - for (auto const& i : phare_box_from(fine_field_box[dirY])) - { - postprocessBy3d(bx, by, bz, meshSize, layout, i); - } + for (auto const& slab : core::make_box_span(fine_field_box[dirY])) + for (auto const& [point, size] : slab) + for (std::size_t i = 0; i < size; ++i, ++point[dimension - 1]) + postprocessBy3d(bx, by, bz, meshSize, layout, point); - for (auto const& i : phare_box_from(fine_field_box[dirZ])) - { - postprocessBz3d(bx, by, bz, meshSize, layout, i); - } + for (auto const& slab : core::make_box_span(fine_field_box[dirZ])) + for (auto const& [point, size] : slab) + for (std::size_t i = 0; i < size; ++i, ++point[dimension - 1]) + postprocessBz3d(bx, by, bz, meshSize, layout, point); } } @@ -134,7 +145,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy return amrIdx[dir] % 2 != 0; } - static void postprocessBx1d(auto& bx, auto const& layout, core::Point idx) + static void postprocessBx1d(auto& bx, auto const& layout, auto const& idx) { auto const locIdx = layout.AMRToLocal(idx); auto const ix = locIdx[dirX]; @@ -142,8 +153,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy bx(ix) = 0.5 * (bx(ix - 1) + bx(ix + 1)); } - static void postprocessBx2d(auto& bx, auto& by, auto const& layout, - core::Point idx) + static void postprocessBx2d(auto& bx, auto& by, auto const& layout, auto const& idx) { auto const locIdx = layout.AMRToLocal(idx); auto const ix = locIdx[dirX]; @@ -157,8 +167,8 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy // modifying, but dual for the field we are indexing to compute // second and third order terms, then the formula reduces to offset // = 1 - int xoffset = 1; - int yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; + int const xoffset = 1; + int const yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; bx(ix, iy) = 0.5 * (bx(ix - 1, iy) + bx(ix + 1, iy)) + 0.25 @@ -169,8 +179,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy } } - static void postprocessBy2d(auto& bx, auto& by, auto const& layout, - core::Point idx) + static void postprocessBy2d(auto& bx, auto& by, auto const& layout, auto const& idx) { auto const locIdx = layout.AMRToLocal(idx); auto const ix = locIdx[dirX]; @@ -180,8 +189,8 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy // | if (isNewFineFace(idx, dirY)) { - int xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; - int yoffset = 1; + int const xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; + int const yoffset = 1; by(ix, iy) = 0.5 * (by(ix, iy - 1) + by(ix, iy + 1)) + 0.25 @@ -193,7 +202,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy } static void postprocessBx3d(auto& bx, auto& by, auto& bz, auto const& meshSize, - auto const& layout, core::Point idx) + auto const& layout, auto const& idx) { auto const Dx = meshSize[dirX]; auto const Dy = meshSize[dirY]; @@ -206,9 +215,9 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy if (isNewFineFace(idx, dirX)) { - int xoffset = 1; - int yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; - int zoffset = (idx[dirZ] % 2 == 0) ? 0 : 1; + int const xoffset = 1; + int const yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; + int const zoffset = (idx[dirZ] % 2 == 0) ? 0 : 1; bx(ix, iy, iz) = 0.5 * (bx(ix - 1, iy, iz) + bx(ix + 1, iy, iz)) @@ -252,7 +261,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy }; static void postprocessBy3d(auto& bx, auto& by, auto& bz, auto const& meshSize, - auto const& layout, core::Point idx) + auto const& layout, auto const& idx) { auto const Dx = meshSize[dirX]; auto const Dy = meshSize[dirY]; @@ -265,9 +274,9 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy if (isNewFineFace(idx, dirY)) { - int xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; - int yoffset = 1; - int zoffset = (idx[dirZ] % 2 == 0) ? 0 : 1; + int const xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; + int const yoffset = 1; + int const zoffset = (idx[dirZ] % 2 == 0) ? 0 : 1; by(ix, iy, iz) = 0.5 * (by(ix, iy - 1, iz) + by(ix, iy + 1, iz)) @@ -311,7 +320,7 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy }; static void postprocessBz3d(auto& bx, auto& by, auto& bz, auto const& meshSize, - auto const& layout, core::Point idx) + auto const& layout, auto const& idx) { auto const Dx = meshSize[dirX]; auto const Dy = meshSize[dirY]; @@ -324,9 +333,9 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy if (isNewFineFace(idx, dirZ)) { - int xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; - int yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; - int zoffset = 1; + int const xoffset = (idx[dirX] % 2 == 0) ? 0 : 1; + int const yoffset = (idx[dirY] % 2 == 0) ? 0 : 1; + int const zoffset = 1; bz(ix, iy, iz) = 0.5 * (bz(ix, iy, iz - 1) + bz(ix, iy, iz + 1)) diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp index feaebf92d..122ebc616 100644 --- a/src/amr/resources_manager/amr_utils.hpp +++ b/src/amr/resources_manager/amr_utils.hpp @@ -345,4 +345,16 @@ namespace amr } // namespace amr } // namespace PHARE +namespace PHARE::amr +{ + +template +NO_DISCARD Box_t coarsen_box(Box_t const& box) +{ + return Box_t{toCoarseIndex(box.lower), toCoarseIndex(box.upper)}; +} + +} // namespace PHARE::amr + + #endif // UTILS_HPP diff --git a/src/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index 8b5b22f3e..7df0b236a 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -286,24 +286,22 @@ void SolverPPC::accumulateFluxSum(IPhysicalModel_t& mode PHARE_LOG_SCOPE(1, "SolverPPC::accumulateFluxSum"); auto& hybridModel = dynamic_cast(model); + auto& Eavg = electromagAvg_.E; + auto& rm = *hybridModel.resourcesManager; - for (auto& patch : level) + for (auto& patch : rm.enumerate(level, fluxSumE_, Eavg)) { - auto& Eavg = electromagAvg_.E; - auto const& layout = amr::layoutFromPatch(*patch); - auto _ = hybridModel.resourcesManager->setOnPatch(*patch, fluxSumE_, Eavg); + auto const&& [Fx, Fy, Fz] = fluxSumE_(); + auto const&& [Ex, Ey, Ez] = Eavg(); - layout.evalOnGhostBox(fluxSumE_(core::Component::X), [&](auto const&... args) mutable { - fluxSumE_(core::Component::X)(args...) += Eavg(core::Component::X)(args...) * coef; - }); + for (std::size_t i = 0; i < Fx.size(); ++i) + Fx.data()[i] += Ex.data()[i] * coef; - layout.evalOnGhostBox(fluxSumE_(core::Component::Y), [&](auto const&... args) mutable { - fluxSumE_(core::Component::Y)(args...) += Eavg(core::Component::Y)(args...) * coef; - }); + for (std::size_t i = 0; i < Fy.size(); ++i) + Fy.data()[i] += Ey.data()[i] * coef; - layout.evalOnGhostBox(fluxSumE_(core::Component::Z), [&](auto const&... args) mutable { - fluxSumE_(core::Component::Z)(args...) += Eavg(core::Component::Z)(args...) * coef; - }); + for (std::size_t i = 0; i < Fz.size(); ++i) + Fz.data()[i] += Ez.data()[i] * coef; } } diff --git a/src/amr/solvers/solver_ppc_model_view.hpp b/src/amr/solvers/solver_ppc_model_view.hpp index c6c25e2ae..09f423b87 100644 --- a/src/amr/solvers/solver_ppc_model_view.hpp +++ b/src/amr/solvers/solver_ppc_model_view.hpp @@ -1,11 +1,12 @@ #ifndef PHARE_SOLVER_SOLVER_PPC_MODEL_VIEW_HPP #define PHARE_SOLVER_SOLVER_PPC_MODEL_VIEW_HPP +#include "core/numerics/ohm/ohm.hpp" #include "core/numerics/ampere/ampere.hpp" #include "core/numerics/faraday/faraday.hpp" -#include "core/numerics/ohm/ohm.hpp" #include "amr/solvers/solver.hpp" +#include "amr/resources_manager/amr_utils.hpp" namespace PHARE::solver { diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp index eb4eb64ee..2ab169a1f 100644 --- a/src/core/data/electrons/electrons.hpp +++ b/src/core/data/electrons/electrons.hpp @@ -119,9 +119,7 @@ class StandardHybridElectronFluxComputer auto& Vez = Ve_(Component::Z); // from Ni because all components defined on primal - layout.evalOnBox(Ne, [&](auto const&... args) { - auto arr = std::array{args...}; - + layout.evalOnBox(Ne, [&](auto const& arr) { auto const JxOnVx = GridLayout::template project(Jx, arr); auto const JyOnVy = GridLayout::template project(Jy, arr); auto const JzOnVz = GridLayout::template project(Jz, arr); diff --git a/src/core/data/field/field_box.hpp b/src/core/data/field/field_box.hpp index 6868b3414..cf88d2153 100644 --- a/src/core/data/field/field_box.hpp +++ b/src/core/data/field/field_box.hpp @@ -1,34 +1,38 @@ #ifndef PHARE_CORE_DATA_FIELD_FIELD_BOX_HPP #define PHARE_CORE_DATA_FIELD_FIELD_BOX_HPP -#include "core/def.hpp" -#include "core/logger.hpp" #include "core/utilities/types.hpp" #include "core/utilities/box/box.hpp" +#include "core/data/field/field_box_span.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" + #include #include +#include +#include #include namespace PHARE::core { -template +template class FieldBox { +public: + using Field_t = Field_t_; using value_type = std::decay_t; -public: auto constexpr static dimension = Field_t::dimension; Field_t& field; - Box amr_ghost_box; + Box amr_box; Box lcl_box; template FieldBox(Field_t& field_, GridLayout_t const& layout) : field{field_} - , amr_ghost_box{layout.AMRGhostBoxFor(field.physicalQuantity())} + , amr_box{layout.AMRGhostBoxFor(field.physicalQuantity())} , lcl_box{layout.ghostBoxFor(field)} { } @@ -37,7 +41,7 @@ class FieldBox FieldBox(Field_t& field_, GridLayout_t const& layout, Box const& selection) : field{field_} - , amr_ghost_box{layout.AMRGhostBoxFor(field.physicalQuantity())} + , amr_box{layout.localToAMR(selection)} , lcl_box{selection} { } @@ -45,7 +49,7 @@ class FieldBox template FieldBox(Field_t& field_, GridLayout_t const& layout, Box const& selection) : field{field_} - , amr_ghost_box{layout.AMRGhostBoxFor(field.physicalQuantity())} + , amr_box{selection} , lcl_box{layout.AMRToLocal(selection)} { } @@ -58,36 +62,83 @@ class FieldBox }; + +template +void operate_on_span(auto& dst, T const* src_data) + requires(std::is_same_v>) +{ + std::memcpy(dst.data(), src_data, dst.size() * sizeof(T)); +} + +template +void operate_on_span(auto& dst, T const* src_data) +{ + for (std::size_t i = 0; i < dst.size(); ++i, ++src_data) + Operator{dst[i]}(*src_data); +} + + + template void operate_on_fields(auto& dst, auto const& src) { assert(dst.lcl_box.size() == src.lcl_box.size()); - auto src_it = src.lcl_box.begin(); - auto dst_it = dst.lcl_box.begin(); - for (; dst_it != dst.lcl_box.end(); ++src_it, ++dst_it) - Operator{dst.field(*dst_it)}(src.field(*src_it)); + auto d_span = make_field_box_span(dst.lcl_box, dst.field); + auto const s_span = make_field_box_span(src.lcl_box, src.field); + + auto d_slabs = d_span.begin(); + auto s_slabs = s_span.begin(); + for (; s_slabs != s_span.end(); ++s_slabs, ++d_slabs) + { + auto d_spans = d_slabs.begin(); + auto s_spans = s_slabs.begin(); + for (; s_spans != s_slabs.end(); ++s_spans, ++d_spans) + operate_on_span(*d_spans, (*s_spans).data()); + } } + template template void FieldBox::set_from(std::vector const& vec, std::size_t seek) { - auto dst_it = lcl_box.begin(); - for (; dst_it != lcl_box.end(); ++seek, ++dst_it) - Operator{field(*dst_it)}(vec[seek]); + for (auto& slab : make_field_box_span(lcl_box, field)) + for (auto& span : slab) + { + operate_on_span(span, vec.data() + seek); + seek += span.size(); + } } + template void FieldBox::append_to(std::vector& vec) { // reserve vec before use! - auto src_it = lcl_box.begin(); - for (; src_it != lcl_box.end(); ++src_it) - vec.push_back(field(*src_it)); + std::size_t seek = vec.size(); + vec.resize(vec.size() + lcl_box.size()); + + for (auto const& slab : make_field_box_span(lcl_box, field)) + for (auto const& span : slab) + { + std::memcpy(vec.data() + seek, span.data(), span.size() * sizeof(value_type)); + seek += span.size(); + } +} + + + +template +void fill_ghost(Field_t& field, auto const& layout, auto const v) +{ + std::size_t const nghosts = layout.nbrGhosts(); + std::size_t const max = nghosts - 1; + field[NdArrayMask{0, max}] = v; } + } // namespace PHARE::core diff --git a/src/core/data/field/field_box_span.hpp b/src/core/data/field/field_box_span.hpp new file mode 100644 index 000000000..bb9404d4c --- /dev/null +++ b/src/core/data/field/field_box_span.hpp @@ -0,0 +1,179 @@ +#ifndef PHARE_CORE_DATA_FIELD_FIELD_BOX_SPAN_HPP +#define PHARE_CORE_DATA_FIELD_FIELD_BOX_SPAN_HPP + + +#include "core/utilities/span.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/box/box_span.hpp" + +#include + + +namespace PHARE::core +{ + +template +class FieldBoxSpans : public BoxSpans +{ + constexpr auto static dim = Array_t::dimension; + using Super = BoxSpans; + using raw_value_type = Array_t::value_type; + using Super::span_size; + +protected: + using Super::point; + +public: + using value_type + = std::conditional_t, raw_value_type const, raw_value_type>; + + FieldBoxSpans(Array_t& arr, auto&&... args) + : Super{args...} + , arr{arr} + { + } + + + auto& operator*() { return (span = Span{&arr(point()), span_size}); } + +private: + Array_t& arr; + Span span{0, 0}; +}; + + +template +class FieldBoxPointSpans : public FieldBoxSpans +{ + auto constexpr static dim = Array_t::dimension; + using Super = FieldBoxSpans; + +public: + using value_type = Super::value_type; + + FieldBoxPointSpans(auto&&... args) + : Super{args...} + , tup{std::forward_as_tuple(*super(), Super::_point)} + { + } + + Super& super() { return *this; } + auto& operator*() + { + std::get<0>(tup) = *super(); + std::get<1>(tup) = this->point(); + return tup; + } + +private: + std::tuple&, Point&> tup; +}; + + + +template> +class FieldBoxSlab : public BoxSlab +{ + auto constexpr static dim = Array_t::dimension; + using FieldBoxSpans_t = Spans_t; + using BaseBoxSpans_t = BoxSpans; + using Super = BoxSlab; + using Super::box; + using Super::slab_idx; + using Super::span_begin; + using Super::span_end; + + +public: + FieldBoxSlab(Array_t& ar, auto&&... args) + : Super{args...} + , arr{ar} + { + } + + FieldBoxSpans_t begin() { return {arr, box, slab_idx, span_begin()}; } + FieldBoxSpans_t begin() const { return {arr, box, slab_idx, span_begin()}; } + BaseBoxSpans_t end() { return {box, slab_idx, span_end()}; } + BaseBoxSpans_t end() const { return {box, slab_idx, span_end()}; } + + auto& operator*() { return *this; } + auto& operator*() const { return *this; } + +private: + Array_t& arr; +}; + + + +template> +class FieldBoxSpan : public BoxSpan +{ + constexpr auto static dim = Array_t::dimension; + using FieldBoxSlab_t = FieldBoxSlab; + using Super = BoxSpan; + using Super::box; + using Super::slab_begin; + using Super::slab_end; + +public: + FieldBoxSpan(Array_t& ar, auto&&... args) + : Super{args...} + , arr{ar} + { + } + + FieldBoxSlab_t begin() { return {arr, box, slabs_begin()}; } + FieldBoxSlab_t begin() const { return {arr, box, slabs_begin()}; } + FieldBoxSlab_t end() { return {arr, box, slabs_end()}; } + FieldBoxSlab_t end() const { return {arr, box, slabs_end()}; } + + std::uint32_t slabs_begin() const + { + if constexpr (dim > 2) + return box.lower[0]; + return 0; + } + std::uint32_t slabs_end() const + { + if constexpr (dim > 2) + return box.upper[0] + 1; + return 1; + } + +private: + Array_t& arr; +}; + + + +template +auto make_field_box_span(Box const box, Array_t& arr) +{ + return FieldBoxSpan{arr, box}; +} + +template +auto make_field_box_span(Box const box, Array_t const& arr) +{ + return FieldBoxSpan{arr, box}; +} + + +template +auto make_field_box_point_span(Box const box, Array_t& arr) +{ + using Row_t = FieldBoxPointSpans; + return FieldBoxSpan{arr, box}; +} + +template +auto make_field_box_point_span(Box const box, Array_t const& arr) +{ + using Row_t = FieldBoxPointSpans; + return FieldBoxSpan{arr, box}; +} + +} // namespace PHARE::core + + +#endif // PHARE_CORE_DATA_FIELD_FIELD_BOX_SPAN_HPP diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index e53ae379e..f6c7cfa0c 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -6,6 +6,7 @@ #include "core/utilities/types.hpp" #include "core/utilities/box/box.hpp" #include "core/utilities/constants.hpp" +#include "core/data/field/field_box.hpp" #include "core/utilities/index/index.hpp" #include "core/utilities/point/point.hpp" #include "core/hybrid/hybrid_quantities.hpp" @@ -1142,22 +1143,27 @@ namespace core template Box ghostBoxFor(Field const& field) const { - return BoxFor(field, [&](auto const& centering, auto const direction) { - return this->ghostStartToEnd(centering, direction); - }); + return boxFor(field, [&](auto&&... args) { return this->ghostStartToEnd(args...); }); } + template + Box domainBoxFor(Field const& field) const + { + return boxFor(field, [&](auto&&... args) { return this->physicalStartToEnd(args...); }); + } - auto AMRBoxFor(auto const& field) const + auto FieldBoxFor(auto const& field, Box box) const { auto const centerings = centering(field); - auto box = AMRBox_; + for (std::uint8_t i = 0; i < dimension; ++i) box.upper[i] += (centerings[i] == QtyCentering::primal) ? 1 : 0; return box; } + auto AMRBoxFor(auto const& field) const { return FieldBoxFor(field, AMRBox_); } + auto AMRGhostBoxFor(auto const& field) const { auto const centerings = centering(field); @@ -1166,24 +1172,16 @@ namespace core } - template - void evalOnBox(Field& field, Fn&& fn) const + template + void evalOnBox(auto const& field, auto&& fn, Args&&... args) const { - auto indices = [&](auto const& centering, auto const direction) { - return this->physicalStartToEnd(centering, direction); - }; - - evalOnBox_(field, fn, indices); + evalOnAnyBox(domainBoxFor(field), fn, std::forward(args)...); } - template - void evalOnGhostBox(Field& field, Fn&& fn) const + template + void evalOnGhostBox(auto const& field, auto&& fn, Args&&... args) const { - auto indices = [&](auto const& centering, auto const direction) { - return this->ghostStartToEnd(centering, direction); - }; - - evalOnBox_(field, fn, indices); + evalOnAnyBox(ghostBoxFor(field), fn, std::forward(args)...); } auto levelNumber() const { return levelNumber_; } @@ -1198,41 +1196,20 @@ namespace core } private: - template - static void evalOnBox_(Field& field, Fn& fn, IndicesFn& startToEnd) + template + void evalOnAnyBox(auto const& box, auto&& fn, Args&&... args) const { - auto const [ix0, ix1] = startToEnd(field, Direction::X); - for (auto ix = ix0; ix <= ix1; ++ix) - { - if constexpr (dimension == 1) - { - fn(ix); - } - else - { - auto const [iy0, iy1] = startToEnd(field, Direction::Y); - - for (auto iy = iy0; iy <= iy1; ++iy) - { - if constexpr (dimension == 2) - { - fn(ix, iy); - } - else - { - auto const [iz0, iz1] = startToEnd(field, Direction::Z); - - for (auto iz = iz0; iz <= iz1; ++iz) - fn(ix, iy, iz); - } - } - } - } + auto constexpr static IDX = dimension - 1; // Fast index + + for (auto const& slab : make_box_span(box)) + for (auto const& [point, _] : slab) + for (; point[IDX] <= box.upper[IDX]; ++point[IDX]) + fn(point, std::forward(args)...); } template - auto BoxFor(Field const& field, Fn startToEnd) const + auto boxFor(Field const& field, Fn const startToEnd) const { std::array lower, upper; diff --git a/src/core/numerics/ampere/ampere.hpp b/src/core/numerics/ampere/ampere.hpp index 1f973d011..777b33e37 100644 --- a/src/core/numerics/ampere/ampere.hpp +++ b/src/core/numerics/ampere/ampere.hpp @@ -1,17 +1,15 @@ #ifndef PHARE_CORE_NUMERICS_AMPERE_AMPERE_HPP #define PHARE_CORE_NUMERICS_AMPERE_AMPERE_HPP -#include -#include - #include "core/data/grid/gridlayoutdefs.hpp" #include "core/data/grid/gridlayout_utils.hpp" #include "core/data/vecfield/vecfield_component.hpp" -#include "core/utilities/index/index.hpp" namespace PHARE::core { + + template class Ampere : public LayoutHolder { @@ -32,53 +30,51 @@ class Ampere : public LayoutHolder auto& Jy = J(Component::Y); auto& Jz = J(Component::Z); - layout_->evalOnBox(Jx, [&](auto&... args) mutable { JxEq_(Jx, B, args...); }); - layout_->evalOnBox(Jy, [&](auto&... args) mutable { JyEq_(Jy, B, args...); }); - layout_->evalOnBox(Jz, [&](auto&... args) mutable { JzEq_(Jz, B, args...); }); + auto& layout = *layout_; + layout.evalOnBox(Jx, [](auto&&... args) { JxEq_(args...); }, Jx, B, layout); + layout.evalOnBox(Jy, [](auto&&... args) { JyEq_(args...); }, Jy, B, layout); + layout.evalOnBox(Jz, [](auto&&... args) { JzEq_(args...); }, Jz, B, layout); } private: - template - void JxEq_(Field& Jx, VecField const& B, Indexes const&... ijk) const + static void JxEq_(auto const& ijk, auto& Jx, auto const& B, auto const& layout) { auto const& [_, By, Bz] = B(); if constexpr (dimension == 1) - Jx(ijk...) = 0.0; + Jx(ijk) = 0.0; if constexpr (dimension == 2) - Jx(ijk...) = layout_->template deriv(Bz, {ijk...}); + Jx(ijk) = layout.template deriv(Bz, ijk); if constexpr (dimension == 3) - Jx(ijk...) = layout_->template deriv(Bz, {ijk...}) - - layout_->template deriv(By, {ijk...}); + Jx(ijk) = layout.template deriv(Bz, ijk) + - layout.template deriv(By, ijk); } - template - void JyEq_(Field& Jy, VecField const& B, Indexes const&... ijk) const + static void JyEq_(auto const& ijk, auto& Jy, auto const& B, auto const& layout) { auto const& [Bx, By, Bz] = B(); if constexpr (dimension == 1 || dimension == 2) - Jy(ijk...) = -layout_->template deriv(Bz, {ijk...}); + Jy(ijk) = -layout.template deriv(Bz, ijk); if constexpr (dimension == 3) - Jy(ijk...) = layout_->template deriv(Bx, {ijk...}) - - layout_->template deriv(Bz, {ijk...}); + Jy(ijk) = layout.template deriv(Bx, ijk) + - layout.template deriv(Bz, ijk); } - template - void JzEq_(Field& Jz, VecField const& B, Indexes const&... ijk) const + static void JzEq_(auto const& ijk, auto& Jz, auto const& B, auto const& layout) { auto const& [Bx, By, Bz] = B(); if constexpr (dimension == 1) - Jz(ijk...) = layout_->template deriv(By, {ijk...}); + Jz(ijk) = layout.template deriv(By, ijk); else - Jz(ijk...) = layout_->template deriv(By, {ijk...}) - - layout_->template deriv(Bx, {ijk...}); + Jz(ijk) = layout.template deriv(By, ijk) + - layout.template deriv(Bx, ijk); } }; diff --git a/src/core/numerics/faraday/faraday.hpp b/src/core/numerics/faraday/faraday.hpp index 0fe28cc77..0ca75bb40 100644 --- a/src/core/numerics/faraday/faraday.hpp +++ b/src/core/numerics/faraday/faraday.hpp @@ -39,59 +39,66 @@ class Faraday : public LayoutHolder auto& Bynew = Bnew(Component::Y); auto& Bznew = Bnew(Component::Z); - layout_->evalOnBox(Bxnew, [&](auto&... args) mutable { BxEq_(Bx, E, Bxnew, args...); }); - layout_->evalOnBox(Bynew, [&](auto&... args) mutable { ByEq_(By, E, Bynew, args...); }); - layout_->evalOnBox(Bznew, [&](auto&... args) mutable { BzEq_(Bz, E, Bznew, args...); }); + auto const& layout = *layout_; + + layout.evalOnGhostBox( + Bxnew, [](auto&&... args) { BxEq_(args...); }, Bx, E, Bxnew, layout, dt_); + layout.evalOnGhostBox( + Bynew, [](auto&&... args) { ByEq_(args...); }, By, E, Bynew, layout, dt_); + layout.evalOnGhostBox( + Bznew, [](auto&&... args) { BzEq_(args...); }, Bz, E, Bznew, layout, dt_); } - private: double dt_; - template - void BxEq_(Field const& Bx, VecField const& E, Field& Bxnew, Indexes const&... ijk) const + static void BxEq_(auto const& ijk, auto&&... args) { - auto const& [_, Ey, Ez] = E(); + auto const& [Bx, E, Bxnew, layout, dt] = std::forward_as_tuple(args...); + auto const& [_, Ey, Ez] = E(); if constexpr (dimension == 1) - Bxnew(ijk...) = Bx(ijk...); + Bxnew(ijk) = Bx(ijk); if constexpr (dimension == 2) - Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_->template deriv(Ez, {ijk...}); + Bxnew(ijk) = Bx(ijk) - dt * layout.template deriv(Ez, ijk); if constexpr (dimension == 3) - Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_->template deriv(Ez, {ijk...}) - + dt_ * layout_->template deriv(Ey, {ijk...}); + Bxnew(ijk) = Bx(ijk) - dt * layout.template deriv(Ez, ijk) + + dt * layout.template deriv(Ey, ijk); } - template - void ByEq_(Field const& By, VecField const& E, Field& Bynew, Indexes const&... ijk) const + + static void ByEq_(auto const& ijk, auto&&... args) { - auto const& [Ex, _, Ez] = E(); + auto const& [By, E, Bynew, layout, dt] = std::forward_as_tuple(args...); + auto const& [Ex, _, Ez] = E(); if constexpr (dimension == 1 || dimension == 2) - Bynew(ijk...) = By(ijk...) + dt_ * layout_->template deriv(Ez, {ijk...}); + Bynew(ijk) = By(ijk) + dt * layout.template deriv(Ez, ijk); if constexpr (dimension == 3) - Bynew(ijk...) = By(ijk...) - dt_ * layout_->template deriv(Ex, {ijk...}) - + dt_ * layout_->template deriv(Ez, {ijk...}); + Bynew(ijk) = By(ijk) - dt * layout.template deriv(Ex, ijk) + + dt * layout.template deriv(Ez, ijk); } - template - void BzEq_(Field const& Bz, VecField const& E, Field& Bznew, Indexes const&... ijk) const + + static void BzEq_(auto const& ijk, auto&&... args) { - auto const& [Ex, Ey, _] = E(); + auto const& [Bz, E, Bznew, layout, dt] = std::forward_as_tuple(args...); + auto const& [Ex, Ey, _] = E(); if constexpr (dimension == 1) - Bznew(ijk...) = Bz(ijk...) - dt_ * layout_->template deriv(Ey, {ijk...}); + Bznew(ijk) = Bz(ijk) - dt * layout.template deriv(Ey, ijk); else - Bznew(ijk...) = Bz(ijk...) - dt_ * layout_->template deriv(Ey, {ijk...}) - + dt_ * layout_->template deriv(Ex, {ijk...}); + Bznew(ijk) = Bz(ijk) - dt * layout.template deriv(Ey, ijk) + + dt * layout.template deriv(Ex, ijk); } }; + } // namespace PHARE::core diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index c67cc391a..790068410 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -35,23 +35,21 @@ class Ohm : public LayoutHolder void operator()(Field const& n, VecField const& Ve, Field const& Pe, VecField const& B, VecField const& J, VecField& Enew) { - using Pack = OhmPack; - if (!this->hasLayout()) throw std::runtime_error( "Error - Ohm - GridLayout not set, cannot proceed to calculate ohm()"); auto const& [Exnew, Eynew, Eznew] = Enew(); - layout_->evalOnBox(Exnew, [&](auto&... args) mutable { - this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); - }); - layout_->evalOnBox(Eynew, [&](auto&... args) mutable { - this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); - }); - layout_->evalOnBox(Eznew, [&](auto&... args) mutable { - this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); - }); + layout_->evalOnBox( + Exnew, [](auto&&... args) { E_Eq_(args...); }, n, Ve, Pe, B, J, Enew, + *this); + layout_->evalOnBox( + Eynew, [](auto&&... args) { E_Eq_(args...); }, n, Ve, Pe, B, J, Enew, + *this); + layout_->evalOnBox( + Eznew, [](auto&&... args) { E_Eq_(args...); }, n, Ve, Pe, B, J, Enew, + *this); } @@ -62,33 +60,25 @@ class Ohm : public LayoutHolder double const nu_; HyperMode const hyper_mode; - template - struct OhmPack - { - VecField& Exyz; - Field const &n, &Pe; - VecField const &Ve, &B, &J; - }; - - template - void E_Eq_(OhmPack&& pack, IDXs const&... ijk) const + template + static void E_Eq_(auto const& ijk, auto&&... args) { - auto const& [E, n, Pe, Ve, B, J] = pack; - auto& Exyz = E(Tag); + auto const& [n, Ve, Pe, B, J, E, self] = std::forward_as_tuple(args...); + auto& Exyz = E(Tag); static_assert(Components::check()); - Exyz(ijk...) = ideal_(Ve, B, {ijk...}) // - + pressure_(n, Pe, {ijk...}) // - + resistive_(J, {ijk...}) // - + hyperresistive_(J, B, n, {ijk...}); + Exyz(ijk) = self.template ideal_(Ve, B, ijk) // + + self.template pressure_(n, Pe, ijk) // + + self.template resistive_(J, ijk) // + + self.template hyperresistive_(J, B, n, ijk); } template - auto ideal_(VecField const& Ve, VecField const& B, MeshIndex index) const + auto ideal_(VecField const& Ve, VecField const& B, auto const index) const { if constexpr (component == Component::X) { @@ -139,7 +129,7 @@ class Ohm : public LayoutHolder template - auto pressure_(Field const& n, Field const& Pe, MeshIndex index) const + auto pressure_(Field const& n, Field const& Pe, auto const index) const { if constexpr (component == Component::X) { @@ -189,7 +179,7 @@ class Ohm : public LayoutHolder template - auto resistive_(VecField const& J, MeshIndex index) const + auto resistive_(VecField const& J, auto const index) const { auto const& Jxyx = J(component); @@ -214,7 +204,7 @@ class Ohm : public LayoutHolder template auto hyperresistive_(VecField const& J, VecField const& B, Field const& n, - MeshIndex index) const + auto const index) const { if (hyper_mode == HyperMode::constant) return constant_hyperresistive_(J, index); @@ -226,7 +216,7 @@ class Ohm : public LayoutHolder template - auto constant_hyperresistive_(VecField const& J, MeshIndex index) const + auto constant_hyperresistive_(VecField const& J, auto const index) const { // TODO : https://github.com/PHAREHUB/PHARE/issues/3 return -nu_ * layout_->laplacian(J(component), index); } @@ -234,7 +224,7 @@ class Ohm : public LayoutHolder template auto spatial_hyperresistive_(VecField const& J, VecField const& B, Field const& n, - MeshIndex index) const + auto const index) const { auto const lvlCoeff = 1. / std::pow(4, layout_->levelNumber()); auto constexpr min_density = 0.1; @@ -263,7 +253,7 @@ class Ohm : public LayoutHolder GridLayout::BzToEz, GridLayout::momentsToEz>(); } } -}; +}; // namespace PHARE::core } // namespace PHARE::core diff --git a/src/core/utilities/box/box_span.hpp b/src/core/utilities/box/box_span.hpp new file mode 100644 index 000000000..035950b99 --- /dev/null +++ b/src/core/utilities/box/box_span.hpp @@ -0,0 +1,209 @@ +#ifndef PHARE_CORE_UTILITIES_BOX_BOX_SPAN_HPP +#define PHARE_CORE_UTILITIES_BOX_BOX_SPAN_HPP + + +// #include "core/def.hpp" +// #include "core/logger.hpp" +// #include "core/utilities/span.hpp" +#include "core/utilities/box/box.hpp" + +#include + + +namespace PHARE::core +{ + +template +class BoxSpans +{ + auto constexpr static dim = dim_; + +public: + using value_type = T; + + BoxSpans(Box const& box, T const slab_idx, T const span_idx) + : box{box} + , slab_idx{slab_idx} + , span_idx{span_idx} + , span_size{_span_size()} + { + } + + BoxSpans& operator++() + { + ++span_idx; + return *this; + } + BoxSpans operator++(int) // postfix increment + { + auto copy = *this; + ++(*this); + return copy; + } + + BoxSpans& operator--() + { + --span_idx; + return *this; + } + BoxSpans operator--(int) // postfix decrement + { + auto copy = *this; + --(*this); + return copy; + } + + auto& point() + { + if constexpr (dim == 1) + return (_point = {box.lower[0]}); + if constexpr (dim == 2) + return (_point = {span_idx, box.lower[1]}); + if constexpr (dim == 3) + return (_point = {slab_idx, span_idx, box.lower[2]}); + } + + auto operator*() { return std::forward_as_tuple(point(), span_size); } + + bool operator==(BoxSpans const& that) const { return span_idx == that.span_idx; } + bool operator!=(BoxSpans const& that) const { return span_idx != that.span_idx; } + bool operator<(BoxSpans const& that) const { return span_idx < that.span_idx; } + bool operator>(BoxSpans const& that) const { return span_idx > that.span_idx; } + + +protected: + std::uint32_t _span_size() const { return box.upper[dim - 1] - box.lower[dim - 1] + 1; } + + Box const& box; + T const slab_idx; + T span_idx; + std::uint32_t const span_size; + Point _point{}; +}; + + +template +class BoxSlab +{ + auto constexpr static dim = dim_; + using BoxSpans_t = BoxSpans; + +public: + using value_type = T; + + BoxSlab(Box const& _box, T const _slab_idx) + : box{_box} + , slab_idx{_slab_idx} + { + } + + BoxSpans_t begin() { return {box, slab_idx, span_begin()}; } + BoxSpans_t begin() const { return {box, slab_idx, span_begin()}; } + BoxSpans_t end() { return {box, slab_idx, span_end()}; } + BoxSpans_t end() const { return {box, slab_idx, span_end()}; } + + auto& operator*() { return *this; } + auto& operator*() const { return *this; } + + bool operator==(BoxSlab const& that) const { return slab_idx == that.slab_idx; } + bool operator!=(BoxSlab const& that) const { return slab_idx != that.slab_idx; } + bool operator<(BoxSlab const& that) const { return slab_idx < that.slab_idx; } + bool operator>(BoxSlab const& that) const { return slab_idx > that.slab_idx; } + + BoxSlab& operator++() + { + ++slab_idx; + return *this; + } + BoxSlab operator++(int) // postfix increment + { + auto copy = *this; + ++(*this); + return copy; + } + + BoxSlab& operator--() + { + --slab_idx; + return *this; + } + BoxSlab operator--(int) // postfix decrement + { + auto copy = *this; + --(*this); + return copy; + } + + T span_begin() const + { + if constexpr (dim > 1) + return box.lower[dim - 2]; + return 0; + } + T span_end() const + { + if constexpr (dim > 1) + return box.upper[dim - 2] + 1; + return 1; + } + +protected: + Box const& box; + T slab_idx; +}; + +template +class BoxSpan +{ + auto constexpr static dim = dim_; + using BoxSlab_t = BoxSlab; + +public: + using value_type = T; + + BoxSpan(Box const& _box) + : box{_box} + { + } + + BoxSlab_t begin() { return {box, slab_begin()}; } + BoxSlab_t begin() const { return {box, slab_begin()}; } + BoxSlab_t end() { return {box, slab_end()}; } + BoxSlab_t end() const { return {box, slab_end()}; } + + T slab_begin() const + { + if constexpr (dim > 2) + return box.lower[0]; + return 0; + } + T slab_end() const + { + if constexpr (dim > 2) + return box.upper[0] + 1; + return 1; + } + + std::uint32_t size() const + { + auto s = slab_end() - slab_begin(); + assert(s > 0); + return s; + } + +protected: + Box const box; +}; + + +template +auto make_box_span(Box const box) +{ + return BoxSpan{box}; +} + + +} // namespace PHARE::core + + +#endif // PHARE_CORE_UTILITIES_BOX_BOX_SPAN_HPP diff --git a/src/core/utilities/span.hpp b/src/core/utilities/span.hpp index dab49c849..1f1906d36 100644 --- a/src/core/utilities/span.hpp +++ b/src/core/utilities/span.hpp @@ -21,20 +21,33 @@ concept Spannable = requires(T t) { template - struct Span { - using value_type = T; + using value_type = std::decay_t; + + Span(T* ptr_ = nullptr, SIZE s_ = 0) + : ptr{ptr_} + , s{s_} + { + } + + Span(Span&&) = default; + Span(Span const&) = default; + Span& operator=(Span&&) = default; + Span& operator=(Span const&) = default; NO_DISCARD auto& operator[](SIZE i) { return ptr[i]; } - NO_DISCARD auto& operator[](SIZE i) const { return ptr[i]; } - NO_DISCARD T const* const& data() const { return ptr; } - NO_DISCARD T const* const& begin() const { return ptr; } - NO_DISCARD T* end() const { return ptr + s; } + NO_DISCARD value_type const& operator[](SIZE i) const { return ptr[i]; } + NO_DISCARD value_type const* data() const { return ptr; } + NO_DISCARD auto data() { return ptr; } + NO_DISCARD auto begin() { return ptr; } + NO_DISCARD auto begin() const { return ptr; } + NO_DISCARD auto end() { return ptr + s; } + NO_DISCARD auto end() const { return ptr + s; } NO_DISCARD SIZE const& size() const { return s; } - T const* ptr = nullptr; - SIZE s = 0; + T* ptr = nullptr; + SIZE s = 0; }; @@ -88,7 +101,11 @@ struct SpanSet { } - NO_DISCARD Span operator[](SIZE i) const + NO_DISCARD Span operator[](SIZE i) + { + return {this->vec.data() + displs[i], this->sizes[i]}; + } + NO_DISCARD Span operator[](SIZE i) const { return {this->vec.data() + displs[i], this->sizes[i]}; } diff --git a/tests/core/data/field/CMakeLists.txt b/tests/core/data/field/CMakeLists.txt new file mode 100644 index 000000000..beeeb9629 --- /dev/null +++ b/tests/core/data/field/CMakeLists.txt @@ -0,0 +1,22 @@ +cmake_minimum_required (VERSION 3.20.1) + +project(test-field) + +function(_field_test test_name) + + add_executable(${test_name} ${test_name}.cpp) + + target_include_directories(${test_name} PRIVATE + ${GTEST_INCLUDE_DIRS} + ) + + target_link_libraries(${test_name} PRIVATE + phare_initializer + ${GTEST_LIBS}) + + add_no_mpi_phare_test(${test_name} ${CMAKE_CURRENT_BINARY_DIR}) + +endfunction(_field_test) + +_field_test(test_field_box_span) + diff --git a/tests/core/data/field/test_field_box_span.cpp b/tests/core/data/field/test_field_box_span.cpp new file mode 100644 index 000000000..c2a5433fb --- /dev/null +++ b/tests/core/data/field/test_field_box_span.cpp @@ -0,0 +1,94 @@ + +#include "core/data/field/field_box_span.hpp" + +#include "phare_core.hpp" +#include "phare_simulator_options.hpp" + +#include "tests/core/data/gridlayout/test_gridlayout.hpp" + +#include "gtest/gtest.h" + +#include + + +namespace PHARE::core +{ + +template +struct TestParam +{ + auto constexpr static dim = _dim; +}; + + +template +class FieldBoxSpanTest : public ::testing::Test +{ + static constexpr std::size_t dim = Param::dim; + static constexpr std::size_t interp = 1; + + static constexpr PHARE::SimOpts opts{dim, interp}; + using PHARE_Types = PHARE::core::PHARE_Types; + using GridLayout_t = TestGridLayout; + using Grid_t = PHARE_Types::Grid_t; + +public: + GridLayout_t layout{9}; + Grid_t dst{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 0}; + Grid_t src{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 0}; +}; + +using Permutations_t = testing::Types, TestParam<2>, TestParam<3>>; + +TYPED_TEST_SUITE(FieldBoxSpanTest, Permutations_t, ); + +TYPED_TEST(FieldBoxSpanTest, test_field_box_span_nd) +{ + auto& layout = this->layout; + auto& dst = this->dst; + auto& src = this->src; + auto const domain_box = layout.domainBoxFor(dst); + + for (auto const& bix : domain_box) + src(bix) = product(bix); + + auto d_span = make_field_box_span(domain_box, dst); + auto const s_span = make_field_box_span(domain_box, src); + + auto d_slabs = d_span.begin(); + auto s_slabs = s_span.begin(); + auto const send = s_span.end(); + for (; s_slabs != send; ++s_slabs, ++d_slabs) + { + auto d_rows = d_slabs.begin(); + auto s_rows = s_slabs.begin(); + auto const srend = s_slabs.end(); + + for (; s_rows != srend; ++s_rows, ++d_rows) + { + auto& s_row = *s_rows; + auto& d_row = *d_rows; + + for (std::size_t i = 0; i < s_row.size(); ++i) + d_row[i] += s_row[i]; + } + } + + + for (auto const& bix : domain_box) + { + EXPECT_EQ(dst(bix), product(bix)); + } + + EXPECT_EQ(sum(dst), sum(src)); +} + + +} // namespace PHARE::core + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/core/utilities/box/CMakeLists.txt b/tests/core/utilities/box/CMakeLists.txt index c0eb43f56..ff512c54a 100644 --- a/tests/core/utilities/box/CMakeLists.txt +++ b/tests/core/utilities/box/CMakeLists.txt @@ -2,18 +2,22 @@ cmake_minimum_required (VERSION 3.20.1) project(test-box) -set(SOURCES test_box.cpp) -add_executable(${PROJECT_NAME} ${SOURCES}) +function(_box_test test_name) -target_include_directories(${PROJECT_NAME} PRIVATE - ${GTEST_INCLUDE_DIRS} + add_executable(${test_name} ${test_name}.cpp) + + target_include_directories(${test_name} PRIVATE + ${GTEST_INCLUDE_DIRS} ) -target_link_libraries(${PROJECT_NAME} PRIVATE - phare_core - ${GTEST_LIBS}) + target_link_libraries(${test_name} PRIVATE + phare_initializer + ${GTEST_LIBS}) -add_no_mpi_phare_test(${PROJECT_NAME} ${CMAKE_CURRENT_BINARY_DIR}) + add_no_mpi_phare_test(${test_name} ${CMAKE_CURRENT_BINARY_DIR}) +endfunction(_box_test) +_box_test(test_box) +_box_test(test_box_span) diff --git a/tests/core/utilities/box/test_box_span.cpp b/tests/core/utilities/box/test_box_span.cpp new file mode 100644 index 000000000..6ccd6624c --- /dev/null +++ b/tests/core/utilities/box/test_box_span.cpp @@ -0,0 +1,193 @@ + + +#include "core/utilities/box/box.hpp" +#include "core/utilities/box/box_span.hpp" +#include "core/data/field/field_box_span.hpp" + +#include "phare_core.hpp" +#include "phare_simulator_options.hpp" + +#include "tests/core/data/gridlayout/test_gridlayout.hpp" + +#include "gtest/gtest.h" +#include + + +namespace PHARE::core +{ + + + +TEST(BoxSpanTest, test_reverse_iterator) +{ + static constexpr std::size_t dim = 3; + static constexpr std::size_t IDX = dim - 1; + static constexpr PHARE::SimOpts opts{dim, 1}; + using PHARE_Types = PHARE::core::PHARE_Types; + using GridLayout_t = TestGridLayout; + + GridLayout_t layout{3}; + + std::vector> points; + + auto const slabs = make_box_span(layout.AMRBox()); + + for (auto slabit = slabs.end(); slabit-- > slabs.begin();) + { + auto const& slab = *slabit; + for (auto spans = slab.end(); spans-- > slab.begin();) + { + auto const&& [start, size] = *spans; + auto point = start; + point[IDX] += size; + for (; point[IDX]-- > start[IDX];) + points.emplace_back(point); + } + } + + std::vector> const expected{ + {2, 2, 2}, {2, 2, 1}, {2, 2, 0}, // + {2, 1, 2}, {2, 1, 1}, {2, 1, 0}, // + {2, 0, 2}, {2, 0, 1}, {2, 0, 0}, // + {1, 2, 2}, {1, 2, 1}, {1, 2, 0}, // + {1, 1, 2}, {1, 1, 1}, {1, 1, 0}, // + {1, 0, 2}, {1, 0, 1}, {1, 0, 0}, // + {0, 2, 2}, {0, 2, 1}, {0, 2, 0}, // + {0, 1, 2}, {0, 1, 1}, {0, 1, 0}, // + {0, 0, 2}, {0, 0, 1}, {0, 0, 0}, + }; + + EXPECT_EQ(expected, points); +} + + +TEST(BoxSpanTest, test_reverse_iterator_unsigned_box) +{ + static constexpr std::size_t dim = 3; + static constexpr std::size_t IDX = dim - 1; + + Box const box{{0, 0, 0}, {2, 2, 2}}; + + std::vector> points; + + auto const slabs = make_box_span(box); + + for (auto slabit = slabs.end(); slabit-- > slabs.begin();) + { + auto const& slab = *slabit; + for (auto spans = slab.end(); spans-- > slab.begin();) + { + auto const&& [start, size] = *spans; + auto point = start; + point[IDX] += size; + for (; point[IDX]-- > start[IDX];) + points.emplace_back(point); + } + } + + std::vector> const expected{ + {2, 2, 2}, {2, 2, 1}, {2, 2, 0}, // + {2, 1, 2}, {2, 1, 1}, {2, 1, 0}, // + {2, 0, 2}, {2, 0, 1}, {2, 0, 0}, // + {1, 2, 2}, {1, 2, 1}, {1, 2, 0}, // + {1, 1, 2}, {1, 1, 1}, {1, 1, 0}, // + {1, 0, 2}, {1, 0, 1}, {1, 0, 0}, // + {0, 2, 2}, {0, 2, 1}, {0, 2, 0}, // + {0, 1, 2}, {0, 1, 1}, {0, 1, 0}, // + {0, 0, 2}, {0, 0, 1}, {0, 0, 0}, + }; + + EXPECT_EQ(expected, points); +} + + +TEST(BoxSpanTest, test_range_loop) +{ + std::size_t static constexpr dim = 3; + Box box{{0, 0, 0}, {9, 9, 9}}; + std::size_t elements = 0; + + for (auto const& slab : make_box_span(box)) + for (auto const& [start, size] : slab) + elements += size; + + EXPECT_EQ(elements, 10 * 10 * 10); +} + + + +TEST(BoxSpanTest, test_iter_loop_dim1) +{ + std::size_t static constexpr dim = 1; + Box box{{0}, {9}}; + std::size_t elements = 0; + + auto const& slabs = make_box_span(box); + + for (auto slabit = slabs.begin(); slabit != slabs.end(); ++slabit) + { + auto const& slab = *slabit; + + for (auto rowit = slab.begin(); rowit != slab.end(); ++rowit) + { + auto const& [start, size] = *rowit; + + elements += size; + } + } + + EXPECT_EQ(elements, 10); +} + +TEST(BoxSpanTest, test_iter_loop_dim2) +{ + std::size_t static constexpr dim = 2; + Box box{{0, 0}, {9, 9}}; + std::size_t elements = 0; + + auto const& slabs = make_box_span(box); + + for (auto slabit = slabs.begin(); slabit != slabs.end(); ++slabit) + { + auto const& slab = *slabit; + + for (auto rowit = slab.begin(); rowit != slab.end(); ++rowit) + { + auto const& [start, size] = *rowit; + elements += size; + } + } + + EXPECT_EQ(elements, 10 * 10); +} + +TEST(BoxSpanTest, test_iter_loop_dim3) +{ + std::size_t static constexpr dim = 3; + Box box{{0, 0, 0}, {9, 9, 9}}; + std::size_t elements = 0; + + auto const& slabs = make_box_span(box); + + for (auto slabit = slabs.begin(); slabit != slabs.end(); ++slabit) + { + auto const& slab = *slabit; + + for (auto rowit = slab.begin(); rowit != slab.end(); ++rowit) + { + auto const& [start, size] = *rowit; + elements += size; + } + } + + EXPECT_EQ(elements, 10 * 10 * 10); +} + +} // namespace PHARE::core + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +}