Skip to content

Commit f9ad5a4

Browse files
lucafedeli88lucafedeli88
committed
move ApplyBCKCorrection out of WarpX class
1 parent eaaaef5 commit f9ad5a4

File tree

2 files changed

+151
-97
lines changed

2 files changed

+151
-97
lines changed

‎Source/EmbeddedBoundary/WarpXFaceExtensions.cpp

+151-88
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,144 @@
2121
using namespace ablastr::fields;
2222
using warpx::fields::FieldType;
2323

24+
namespace
25+
{
26+
27+
/**
28+
* \brief Compute the minimal area for stability for the face i, j, k with normal 'dim'.
29+
*
30+
* \tparam dim normal direction to the plane in consideration (0 for x, 1 for y, 2 for z)
31+
*
32+
* \param[in] i, j, k the indices of the cell
33+
* \param[in] lx, ly, lz arrays containing the edge lengths
34+
* \param[in] dx, dy, dz the mesh with in each direction
35+
*/
36+
template <int dim>
37+
[[nodiscard]]
38+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
39+
amrex::Real
40+
ComputeSStab(const int i, const int j, const int k,
41+
const amrex::Array4<const amrex::Real> lx,
42+
const amrex::Array4<const amrex::Real> ly,
43+
const amrex::Array4<const amrex::Real> lz,
44+
const amrex::Real dx, const amrex::Real dy, const amrex::Real dz)
45+
{
46+
47+
using namespace amrex::literals;
48+
49+
if constexpr (dim == 0) {
50+
return 0.5_rt * std::max({ly(i, j, k) * dz, ly(i, j, k + 1) * dz,
51+
lz(i, j, k) * dy, lz(i, j + 1, k) * dy});
52+
}
53+
else if constexpr (dim == 1)
54+
{
55+
#ifdef WARPX_DIM_XZ
56+
return 0.5_rt * std::max({lx(i, j, k) * dz, lx(i, j + 1, k) * dz,
57+
lz(i, j, k) * dx, lz(i + 1, j, k) * dx});
58+
#elif defined(WARPX_DIM_3D)
59+
return 0.5_rt * std::max({lx(i, j, k) * dz, lx(i, j, k + 1) * dz,
60+
lz(i, j, k) * dx, lz(i + 1, j, k) * dx});
61+
#else
62+
amrex::Abort("ComputeSStab: Only implemented in 2D3V and 3D3V");
63+
#endif
64+
}
65+
else if constexpr(dim == 2){
66+
return 0.5_rt * std::max({lx(i, j, k) * dy, lx(i, j + 1, k) * dy,
67+
ly(i, j, k) * dx, ly(i + 1, j, k) * dx});
68+
}
69+
70+
amrex::Abort("ComputeSStab: dim must be 0, 1 or 2");
71+
72+
return -1;
73+
}
74+
75+
/**
76+
* \brief Compute the minimal area for stability for the face i, j, k with normal 'dim'.
77+
* (wrapper to allow using ComputeSStab as a non-templated function, by passing 'dim' as an argument)
78+
*
79+
* \param[in] i, j, k the indices of the cell
80+
* \param[in] lx, ly, lz arrays containing the edge lengths
81+
* \param[in] dx, dy, dz the mesh with in each direction
82+
* \param[in] dim normal direction to the plane in consideration (0 for x, 1 for y, 2 for z)
83+
*/
84+
[[nodiscard]]
85+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
86+
amrex::Real
87+
ComputeSStab(const int i, const int j, const int k,
88+
const amrex::Array4<const amrex::Real> lx,
89+
const amrex::Array4<const amrex::Real> ly,
90+
const amrex::Array4<const amrex::Real> lz,
91+
const amrex::Real dx, const amrex::Real dy, const amrex::Real dz,
92+
const int dim)
93+
{
94+
if (dim == 0) {
95+
return ::ComputeSStab<0>(i,j,k,lx,ly,lz,dx,dy,dz);
96+
}
97+
else if (dim == 1) {
98+
return ::ComputeSStab<1>(i,j,k,lx,ly,lz,dx,dy,dz);
99+
}
100+
else if (dim == 2) {
101+
return ::ComputeSStab<2>(i,j,k,lx,ly,lz,dx,dy,dz);
102+
}
103+
return -1;
104+
}
105+
106+
107+
/**
108+
* \brief Whenever an unstable cell cannot be extended we increase its area to be the minimal for stability.
109+
* This is the method Benkler-Chavannes-Kuster method and it is less accurate than the regular ECT but it
110+
* still works better than staircasing. (see https://ieeexplore.ieee.org/document/1638381)
111+
*
112+
* @tparam idim Integer indicating the dimension (x->0, y->1, z->2) for which the BCK correction is done
113+
*
114+
* @param[in] cell_size_max_lev The cell size at the maximum refinement level
115+
* @param[in, out] all_fields The field manager
116+
* @param[in] flag_ext_face The extension flag used by the ECT solver
117+
* @param[out] flag_info_face The info flag used by the ECT solver
118+
*/
119+
template <int idim>
120+
void ApplyBCKCorrection(
121+
const std::array<amrex::Real,3> &cell_size_max_lev,
122+
ablastr::fields::MultiFabRegister& all_fields,
123+
const int max_level,
124+
const amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>, 3 > >& flag_ext_face,
125+
const amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>, 3 > >& flag_info_face)
126+
{
127+
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EB::enabled(),
128+
"ApplyBCKCorrection only works when EBs are enabled at runtime");
129+
130+
#if defined(AMREX_USE_EB) and !defined(WARPX_DIM_RZ)
131+
132+
const amrex::Real dx = cell_size_max_lev[0];
133+
const amrex::Real dy = cell_size_max_lev[1];
134+
const amrex::Real dz = cell_size_max_lev[2];
135+
136+
for (amrex::MFIter mfi(*all_fields.get(FieldType::Bfield_fp, Direction{idim}, max_level), amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
137+
138+
const amrex::Box &box = mfi.tilebox();
139+
const amrex::Array4<int> &flag_ext_face_max_lev_idim = flag_ext_face[max_level][idim]->array(mfi);
140+
const amrex::Array4<int> &flag_info_face_max_lev_idim = flag_info_face[max_level][idim]->array(mfi);
141+
const amrex::Array4<amrex::Real> &S = all_fields.get(FieldType::face_areas, Direction{idim}, max_level)->array(mfi);
142+
const amrex::Array4<amrex::Real> &lx = all_fields.get(FieldType::face_areas, Direction{0}, max_level)->array(mfi);
143+
const amrex::Array4<amrex::Real> &ly = all_fields.get(FieldType::face_areas, Direction{1}, max_level)->array(mfi);
144+
const amrex::Array4<amrex::Real> &lz = all_fields.get(FieldType::face_areas, Direction{2}, max_level)->array(mfi);
145+
146+
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
147+
if (flag_ext_face_max_lev_idim(i, j, k)) {
148+
// Modify the area according to the BCK algorithm
149+
S(i, j, k) = ::ComputeSStab<idim>(i, j, k, lx, ly, lz, dx, dy, dz);
150+
// Update the face info so that the solver doesn't think that this face is being extended
151+
flag_info_face_max_lev_idim(i, j, k) = -1;
152+
}
153+
});
154+
}
155+
#else
156+
amrex::ignore_unused(idim, cell_size_max_lev, all_fields, max_level, flag_ext_face, flag_info_face);
157+
#endif
158+
}
159+
160+
}
161+
24162
/**
25163
* \brief Get the value of arr in the neighbor (i_n, j_n) on the plane with normal 'dim'.
26164
*
@@ -124,50 +262,6 @@ SetNeigh(const amrex::Array4<T>& arr, const T val,
124262
amrex::Abort("SetNeigh: dim must be 0, 1 or 2");
125263
}
126264

127-
128-
/**
129-
* \brief Compute the minimal area for stability for the face i, j, k with normal 'dim'.
130-
*
131-
* \param[in] i, j, k the indices of the cell
132-
* \param[in] lx, ly, lz arrays containing the edge lengths
133-
* \param[in] dx, dy, dz the mesh with in each direction
134-
* \param[in] dim normal direction to the plane in consideration (0 for x, 1 for y, 2 for z)
135-
*/
136-
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
137-
amrex::Real
138-
ComputeSStab(const int i, const int j, const int k,
139-
const amrex::Array4<const amrex::Real> lx,
140-
const amrex::Array4<const amrex::Real> ly,
141-
const amrex::Array4<const amrex::Real> lz,
142-
const amrex::Real dx, const amrex::Real dy, const amrex::Real dz,
143-
const int dim){
144-
145-
using namespace amrex::literals;
146-
147-
if(dim == 0) {
148-
return 0.5_rt * std::max({ly(i, j, k) * dz, ly(i, j, k + 1) * dz,
149-
lz(i, j, k) * dy, lz(i, j + 1, k) * dy});
150-
}else if(dim == 1){
151-
#ifdef WARPX_DIM_XZ
152-
return 0.5_rt * std::max({lx(i, j, k) * dz, lx(i, j + 1, k) * dz,
153-
lz(i, j, k) * dx, lz(i + 1, j, k) * dx});
154-
#elif defined(WARPX_DIM_3D)
155-
return 0.5_rt * std::max({lx(i, j, k) * dz, lx(i, j, k + 1) * dz,
156-
lz(i, j, k) * dx, lz(i + 1, j, k) * dx});
157-
#else
158-
amrex::Abort("ComputeSStab: Only implemented in 2D3V and 3D3V");
159-
#endif
160-
}else if(dim == 2){
161-
return 0.5_rt * std::max({lx(i, j, k) * dy, lx(i, j + 1, k) * dy,
162-
ly(i, j, k) * dx, ly(i + 1, j, k) * dx});
163-
}
164-
165-
amrex::Abort("ComputeSStab: dim must be 0, 1 or 2");
166-
167-
return -1;
168-
}
169-
170-
171265
amrex::Array1D<int, 0, 2>
172266
WarpX::CountExtFaces () {
173267
amrex::Array1D<int, 0, 2> sums{0, 0, 0};
@@ -254,19 +348,25 @@ WarpX::ComputeFaceExtensions ()
254348
// If any cell could not be extended we use the BCK method to stabilize them
255349
#if !defined(WARPX_DIM_XZ) && !defined(WARPX_DIM_RZ)
256350
if (N_ext_faces_after_eight_ways(0) > 0) {
257-
ApplyBCKCorrection(0);
351+
::ApplyBCKCorrection<0>(
352+
CellSize(maxLevel()), m_fields, maxLevel(),
353+
m_flag_ext_face, m_flag_info_face);
258354
using_bck = true;
259355
}
260356
#endif
261357

262358
if (N_ext_faces_after_eight_ways(1) > 0) {
263-
ApplyBCKCorrection(1);
359+
::ApplyBCKCorrection<1>(
360+
CellSize(maxLevel()), m_fields, maxLevel(),
361+
m_flag_ext_face, m_flag_info_face);
264362
using_bck = true;
265363
}
266364

267365
#if !defined(WARPX_DIM_XZ) && !defined(WARPX_DIM_RZ)
268366
if (N_ext_faces_after_eight_ways(2) > 0) {
269-
ApplyBCKCorrection(2);
367+
::ApplyBCKCorrection<2>(
368+
CellSize(maxLevel()), m_fields, maxLevel(),
369+
m_flag_ext_face, m_flag_info_face);
270370
using_bck = true;
271371
}
272372
#endif
@@ -490,7 +590,7 @@ WarpX::ComputeOneWayExtensions ()
490590
return 0;
491591
}
492592

493-
const amrex::Real S_stab = ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
593+
const amrex::Real S_stab = ::ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
494594

495595
const amrex::Real S_ext = S_stab - S(i, j, k);
496596
const int n_borrow =
@@ -512,7 +612,7 @@ WarpX::ComputeOneWayExtensions ()
512612
} else{
513613
borrowing_inds_pointer(i, j, k) = borrowing_inds + ps;
514614

515-
const amrex::Real S_stab = ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
615+
const amrex::Real S_stab = ::ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
516616

517617
const amrex::Real S_ext = S_stab - S(i, j, k);
518618
for (int i_n = -1; i_n < 2; i_n++) {
@@ -617,7 +717,7 @@ WarpX::ComputeEightWaysExtensions ()
617717
if (!flag_ext_face(i, j, k)) {
618718
return 0;
619719
}
620-
const amrex::Real S_stab = ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
720+
const amrex::Real S_stab = ::ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
621721

622722
const amrex::Real S_ext = S_stab - S(i, j, k);
623723
const int n_borrow = ComputeNBorrowEightFacesExtension(cell, S_ext, S_mod, S,
@@ -646,7 +746,7 @@ WarpX::ComputeEightWaysExtensions ()
646746
borrowing_inds_pointer(i, j, k) = borrowing_inds + ps;
647747

648748
S_mod(i, j, k) = S(i, j, k);
649-
const amrex::Real S_stab = ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
749+
const amrex::Real S_stab = ::ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
650750

651751
const amrex::Real S_ext = S_stab - S(i, j, k);
652752
amrex::Array2D<amrex::Real, 0, 2, 0, 2> local_avail{};
@@ -724,43 +824,6 @@ WarpX::ComputeEightWaysExtensions ()
724824
#endif
725825
}
726826

727-
void
728-
WarpX::ApplyBCKCorrection (const int idim)
729-
{
730-
if (!EB::enabled()) {
731-
throw std::runtime_error("ApplyBCKCorrection only works when EBs are enabled at runtime");
732-
}
733-
#if defined(AMREX_USE_EB) and !defined(WARPX_DIM_RZ)
734-
const std::array<amrex::Real,3> &cell_size = CellSize(maxLevel());
735-
736-
const amrex::Real dx = cell_size[0];
737-
const amrex::Real dy = cell_size[1];
738-
const amrex::Real dz = cell_size[2];
739-
740-
for (amrex::MFIter mfi(*m_fields.get(FieldType::Bfield_fp, Direction{idim}, maxLevel()), amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
741-
742-
const amrex::Box &box = mfi.tilebox();
743-
const amrex::Array4<int> &flag_ext_face = m_flag_ext_face[maxLevel()][idim]->array(mfi);
744-
const amrex::Array4<int> &flag_info_face = m_flag_info_face[maxLevel()][idim]->array(mfi);
745-
const amrex::Array4<amrex::Real> &S = m_fields.get(FieldType::face_areas, Direction{idim}, maxLevel())->array(mfi);
746-
const amrex::Array4<amrex::Real> &lx = m_fields.get(FieldType::face_areas, Direction{0}, maxLevel())->array(mfi);;
747-
const amrex::Array4<amrex::Real> &ly = m_fields.get(FieldType::face_areas, Direction{1}, maxLevel())->array(mfi);;
748-
const amrex::Array4<amrex::Real> &lz = m_fields.get(FieldType::face_areas, Direction{2}, maxLevel())->array(mfi);;
749-
750-
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
751-
if (flag_ext_face(i, j, k)) {
752-
// Modify the area according to the BCK algorithm
753-
S(i, j, k) = ComputeSStab(i, j, k, lx, ly, lz, dx, dy, dz, idim);
754-
// Update the face info so that the solver doesn't think that this face is being extended
755-
flag_info_face(i, j, k) = -1;
756-
}
757-
});
758-
}
759-
#else
760-
amrex::ignore_unused(idim);
761-
#endif
762-
}
763-
764827
void
765828
WarpX::ShrinkBorrowing ()
766829
{

‎Source/WarpX.H

-9
Original file line numberDiff line numberDiff line change
@@ -1122,15 +1122,6 @@ public:
11221122
* \brief Do the eight-ways extension
11231123
*/
11241124
void ComputeEightWaysExtensions();
1125-
/**
1126-
* \brief Whenever an unstable cell cannot be extended we increase its area to be the minimal for stability.
1127-
* This is the method Benkler-Chavannes-Kuster method and it is less accurate than the regular ECT but it
1128-
* still works better than staircasing. (see https://ieeexplore.ieee.org/document/1638381)
1129-
*
1130-
* @param idim Integer indicating the dimension (x->0, y->1, z->2) for which the BCK correction is done
1131-
*
1132-
*/
1133-
void ApplyBCKCorrection(int idim);
11341125

11351126
/**
11361127
* \brief Subtract the average of the cumulative sums of the preliminary current D

0 commit comments

Comments
 (0)