Skip to content

Commit

Permalink
Update initialization of external fields
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiLehe
RemiLehe committed Jan 11, 2025
1 parent ab87809 commit 78bb8c1
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 89 deletions.
10 changes: 4 additions & 6 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,9 +226,8 @@ void HybridPICModel::InitData ()
m_J_external[0],
m_J_external[1],
m_J_external[2],
lev, PatchType::fine, 'e',
warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev),
warpx.m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine,
warpx.m_eb_update_E);
}
}

Expand All @@ -244,9 +243,8 @@ void HybridPICModel::GetCurrentExternal ()
m_J_external[0],
m_J_external[1],
m_J_external[2],
lev, PatchType::fine, 'e',
warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev),
warpx.m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine,
warpx.m_eb_update_E);
}
}

Expand Down
105 changes: 27 additions & 78 deletions Source/Initialization/WarpXInitData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -976,23 +976,21 @@ WarpX::InitLevelData (int lev, Real /*time*/)
// The default maxlevel_extEMfield_init value is the total number of levels in the simulation
if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
&& (lev > 0) && (lev <= maxlevel_extEMfield_init)) {

// TODO: raise error when EB is activated
ComputeExternalFieldOnGridUsingParser(
FieldType::Bfield_aux,
m_p_ext_field_params->Bxfield_parser->compile<4>(),
m_p_ext_field_params->Byfield_parser->compile<4>(),
m_p_ext_field_params->Bzfield_parser->compile<4>(),
lev, PatchType::fine, 'f',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine, m_eb_update_B);

ComputeExternalFieldOnGridUsingParser(
FieldType::Bfield_cp,
m_p_ext_field_params->Bxfield_parser->compile<4>(),
m_p_ext_field_params->Byfield_parser->compile<4>(),
m_p_ext_field_params->Bzfield_parser->compile<4>(),
lev, PatchType::coarse, 'f',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_mr_levels_alldirs(FieldType::face_areas, max_level)[lev]);
lev, PatchType::coarse, m_eb_update_E);
}

// if the input string for the E-field is "parse_e_ext_grid_function",
Expand Down Expand Up @@ -1023,18 +1021,14 @@ WarpX::InitLevelData (int lev, Real /*time*/)
m_p_ext_field_params->Exfield_parser->compile<4>(),
m_p_ext_field_params->Eyfield_parser->compile<4>(),
m_p_ext_field_params->Ezfield_parser->compile<4>(),
lev, PatchType::fine, 'e',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine, m_eb_update_E);

ComputeExternalFieldOnGridUsingParser(
FieldType::Efield_cp,
m_p_ext_field_params->Exfield_parser->compile<4>(),
m_p_ext_field_params->Eyfield_parser->compile<4>(),
m_p_ext_field_params->Ezfield_parser->compile<4>(),
lev, PatchType::coarse, 'e',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::coarse, m_eb_update_E);
#ifdef AMREX_USE_EB
if (eb_enabled) {
if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) {
Expand Down Expand Up @@ -1068,9 +1062,8 @@ void WarpX::ComputeExternalFieldOnGridUsingParser (
amrex::ParserExecutor<4> const& fx_parser,
amrex::ParserExecutor<4> const& fy_parser,
amrex::ParserExecutor<4> const& fz_parser,
int lev, PatchType patch_type, [[maybe_unused]] const char topology,
std::optional<ablastr::fields::VectorField> const& edge_lengths,
std::optional<ablastr::fields::VectorField> const& face_areas)
int lev, PatchType patch_type,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > const& eb_update_field)
{
auto t = gett_new(lev);

Expand All @@ -1093,8 +1086,6 @@ void WarpX::ComputeExternalFieldOnGridUsingParser (
const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect();
const amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect();

const bool eb_enabled = EB::enabled();

for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
const amrex::Box& tbx = mfi.tilebox( x_nodal_flag, mfx->nGrowVect() );
const amrex::Box& tby = mfi.tilebox( y_nodal_flag, mfy->nGrowVect() );
Expand All @@ -1104,44 +1095,19 @@ void WarpX::ComputeExternalFieldOnGridUsingParser (
auto const& mfyfab = mfy->array(mfi);
auto const& mfzfab = mfz->array(mfi);

amrex::Array4<amrex::Real> lx, ly, lz, Sx, Sy, Sz;
if (eb_enabled) {
if (edge_lengths.has_value()) {
const auto& edge_lengths_array = edge_lengths.value();
lx = edge_lengths_array[0]->array(mfi);
ly = edge_lengths_array[1]->array(mfi);
lz = edge_lengths_array[2]->array(mfi);
}
if (face_areas.has_value()) {
const auto& face_areas_array = face_areas.value();
Sx = face_areas_array[0]->array(mfi);
Sy = face_areas_array[1]->array(mfi);
Sz = face_areas_array[2]->array(mfi);
}
}

#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
amrex::Dim3 lx_lo, lx_hi, lz_lo, lz_hi;
#endif
if (eb_enabled) {
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
lx_lo = amrex::lbound(lx);
lx_hi = amrex::ubound(lx);
lz_lo = amrex::lbound(lz);
lz_hi = amrex::ubound(lz);
#endif
amrex::Array4<int> update_fx_arr, update_fy_arr, update_fz_arr;
if (EB::enabled()) {
update_fx_arr = eb_update_field[lev][0]->array(mfi);
update_fy_arr = eb_update_field[lev][1]->array(mfi);
update_fz_arr = eb_update_field[lev][2]->array(mfi);
}

amrex::ParallelFor (tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef AMREX_USE_EB
#ifdef WARPX_DIM_3D
if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and Sx(i, j, k)<=0))) { return; }
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
//In XZ and RZ Ex is associated with a x-edge, while Bx is associated with a z-edge
if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and lz(i, j, k)<=0))) { return; }
#endif
#endif

// Do not set fields inside the embedded boundary
if (update_fx_arr && update_fx_arr(i,j,k) == 0) { return; }

// Shift required in the x-, y-, or z- position
// depending on the index type of the multifab
#if defined(WARPX_DIM_1D_Z)
Expand All @@ -1167,19 +1133,10 @@ void WarpX::ComputeExternalFieldOnGridUsingParser (
mfxfab(i,j,k) = fx_parser(x,y,z,t);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef AMREX_USE_EB
#ifdef WARPX_DIM_3D
if(ly && ((topology=='e' and ly(i, j, k)<=0) or (topology=='f' and Sy(i, j, k)<=0))) { return; }
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
//In XZ and RZ Ey is associated with a mesh node, so we need to check if the mesh node is covered
if(lx &&
((topology=='e' and (lx(std::min(i , lx_hi.x), std::min(j , lx_hi.y), k)<=0
|| lx(std::max(i-1, lx_lo.x), std::min(j , lx_hi.y), k)<=0
|| lz(std::min(i , lz_hi.x), std::min(j , lz_hi.y), k)<=0
|| lz(std::min(i , lz_hi.x), std::max(j-1, lz_lo.y), k)<=0)) or
(topology=='f' and Sy(i,j,k)<=0))) { return; }
#endif
#endif

// Do not set fields inside the embedded boundary
if (update_fy_arr && update_fy_arr(i,j,k) == 0) { return; }

#if defined(WARPX_DIM_1D_Z)
const amrex::Real x = 0._rt;
const amrex::Real y = 0._rt;
Expand All @@ -1203,14 +1160,10 @@ void WarpX::ComputeExternalFieldOnGridUsingParser (
mfyfab(i,j,k) = fy_parser(x,y,z,t);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef AMREX_USE_EB
#ifdef WARPX_DIM_3D
if(lz && ((topology=='e' and lz(i, j, k)<=0) or (topology=='f' and Sz(i, j, k)<=0))) { return; }
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
//In XZ and RZ Ez is associated with a z-edge, while Bz is associated with a x-edge
if(lz && ((topology=='e' and lz(i, j, k)<=0) or (topology=='f' and lx(i, j, k)<=0))) { return; }
#endif
#endif

// Do not set fields inside the embedded boundary
if (update_fz_arr && update_fz_arr(i,j,k) == 0) { return; }

#if defined(WARPX_DIM_1D_Z)
const amrex::Real x = 0._rt;
const amrex::Real y = 0._rt;
Expand Down Expand Up @@ -1407,9 +1360,7 @@ WarpX::LoadExternalFields (int const lev)
m_p_ext_field_params->Bxfield_parser->compile<4>(),
m_p_ext_field_params->Byfield_parser->compile<4>(),
m_p_ext_field_params->Bzfield_parser->compile<4>(),
lev, PatchType::fine, 'f',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine, m_eb_update_B);
}
else if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) {
#if defined(WARPX_DIM_RZ)
Expand All @@ -1432,9 +1383,7 @@ WarpX::LoadExternalFields (int const lev)
m_p_ext_field_params->Exfield_parser->compile<4>(),
m_p_ext_field_params->Eyfield_parser->compile<4>(),
m_p_ext_field_params->Ezfield_parser->compile<4>(),
lev, PatchType::fine, 'e',
m_fields.get_alldirs(FieldType::edge_lengths, lev),
m_fields.get_alldirs(FieldType::face_areas, lev));
lev, PatchType::fine, m_eb_update_E );
}
else if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) {
#if defined(WARPX_DIM_RZ)
Expand Down
12 changes: 7 additions & 5 deletions Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -937,9 +937,8 @@ public:
amrex::ParserExecutor<4> const& fx_parser,
amrex::ParserExecutor<4> const& fy_parser,
amrex::ParserExecutor<4> const& fz_parser,
int lev, PatchType patch_type, [[maybe_unused]] char topology,
std::optional<ablastr::fields::VectorField> const& edge_lengths = std::nullopt,
std::optional<ablastr::fields::VectorField> const& face_areas = std::nullopt);
int lev, PatchType patch_type,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > const& eb_update_field);

/**
* \brief Load field values from a user-specified openPMD file,
Expand Down Expand Up @@ -1416,8 +1415,11 @@ private:
/** EB: Flag to indicate whether a gridpoint is inside the embedded boundary and therefore
* whether the E or B should not be updated. (One array per level and per direction, due to staggering)
*/
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > m_eb_update_E;
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > m_eb_update_B;
// TODO: leave as private, add setters, getters
public:
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > m_eb_update_E;
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > m_eb_update_B;
private:

/** EB: for every mesh face flag_info_face contains a:
* * 0 if the face needs to be extended
Expand Down

0 comments on commit 78bb8c1

Please sign in to comment.