Skip to content
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 res/cmake/test.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
164 changes: 50 additions & 114 deletions src/amr/data/field/refine/electric_field_refiner.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,14 @@


#include "core/def/phare_mpi.hpp" // IWYU pragma: keep

#include <SAMRAI/hier/Box.h>

#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 <SAMRAI/hier/Box.h>

#include <cstddef>

namespace PHARE::amr
Expand All @@ -37,22 +37,23 @@ class ElectricFieldRefiner

// electric field refinement strategy follows
// fujimoto et al. 2011 : doi:10.1016/j.jcp.2011.08.002
template<typename FieldT>
void operator()(FieldT const& coarseField, FieldT& fineField,
core::Point<int, dimension> 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:
Expand All @@ -73,44 +74,13 @@ class ElectricFieldRefiner
}


template<typename FieldT>
void refine1D_(FieldT const& coarseField, FieldT& fineField,
core::Point<int, dimension> const& locFineIdx,
core::Point<int, dimension> 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<typename FieldT>
void refine2D_(FieldT const& coarseField, FieldT& fineField,
core::Point<int, dimension> const& fineIndex,
core::Point<int, dimension> const& locFineIdx,
core::Point<int, dimension> 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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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<typename FieldT>
void refine3D_(FieldT const& coarseField, FieldT& fineField,
core::Point<int, dimension> const& fineIndex,
core::Point<int, dimension> const& locFineIdx,
core::Point<int, dimension> 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
Expand All @@ -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));
Expand All @@ -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
Expand All @@ -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));
}
Expand All @@ -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));
}
Expand Down
4 changes: 1 addition & 3 deletions src/amr/data/field/refine/field_moments_refiner.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand All @@ -14,7 +13,6 @@
#include <SAMRAI/hier/Box.h>

#include <array>
#include <vector>


namespace PHARE::amr
Expand Down
Loading
Loading