Skip to content

Commit

Permalink
General parser function for external fields (#5349)
Browse files Browse the repository at this point in the history 
This feature replaces the InitializeExternalFieldsOnGridUsingParser with
a more general function ComputeExternalFieldOnGridUsingParser. The new
function takes parsed functions with four dimensions (x,y,z,t), so that
it can be used to evaluate functions throughout the simulation, for
example time-dependent external fields.

The function ComputeExternalFieldOnGridUsingParser has optional edge
length and face areas.

---------

Co-authored-by: Kristoffer Lindvall <kristoffer.lindvall@novatronfusion.com>
  • Loading branch information
@kli-jfp
kli-jfp and Kristoffer Lindvall authored Oct 2, 2024
1 parent 6757b4c commit 841d7ee
Show file tree
Hide file tree
Showing 8 changed files with 144 additions and 298 deletions.
8 changes: 1 addition & 7 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,7 @@ public:
* external current multifab. Note the external current can be a function
* of time and therefore this should be re-evaluated at every step.
*/
void GetCurrentExternal (
ablastr::fields::MultiLevelVectorField const& edge_lengths
);
void GetCurrentExternal (
ablastr::fields::VectorField const& edge_lengths,
int lev
);
void GetCurrentExternal ();

/**
* \brief
Expand Down
169 changes: 17 additions & 152 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,167 +221,32 @@ void HybridPICModel::InitData ()
// if the current is time dependent which is what needs to be done to
// write time independent fields on the first step.
for (int lev = 0; lev <= warpx.finestLevel(); ++lev) {
auto edge_lengths = std::array<std::unique_ptr<amrex::MultiFab>, 3>();
#ifdef AMREX_USE_EB
if (EB::enabled()) {
using ablastr::fields::Direction;
auto const & edge_lengths_x = *warpx.m_fields.get(FieldType::edge_lengths, Direction{0}, lev);
auto const & edge_lengths_y = *warpx.m_fields.get(FieldType::edge_lengths, Direction{1}, lev);
auto const & edge_lengths_z = *warpx.m_fields.get(FieldType::edge_lengths, Direction{2}, lev);

edge_lengths = std::array< std::unique_ptr<amrex::MultiFab>, 3 >{
std::make_unique<amrex::MultiFab>(
edge_lengths_x, amrex::make_alias, 0, edge_lengths_x.nComp()),
std::make_unique<amrex::MultiFab>(
edge_lengths_y, amrex::make_alias, 0, edge_lengths_y.nComp()),
std::make_unique<amrex::MultiFab>(
edge_lengths_z, amrex::make_alias, 0, edge_lengths_z.nComp())
};
}
#endif
GetCurrentExternal(ablastr::fields::a2m(edge_lengths), lev);
warpx.ComputeExternalFieldOnGridUsingParser(
FieldType::hybrid_current_fp_external,
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));
}
}

void HybridPICModel::GetCurrentExternal (
ablastr::fields::MultiLevelVectorField const& edge_lengths)
void HybridPICModel::GetCurrentExternal ()
{
if (!m_external_field_has_time_dependence) { return; }

auto& warpx = WarpX::GetInstance();
for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
{
GetCurrentExternal(edge_lengths[lev], lev);
}
}


void HybridPICModel::GetCurrentExternal (
ablastr::fields::VectorField const& edge_lengths,
int lev)
{
// This logic matches closely to WarpX::InitializeExternalFieldsOnGridUsingParser
// except that the parsers include time dependence.
auto & warpx = WarpX::GetInstance();

auto t = warpx.gett_new(lev);

auto dx_lev = warpx.Geom(lev).CellSizeArray();
const RealBox& real_box = warpx.Geom(lev).ProbDomain();

using ablastr::fields::Direction;
amrex::MultiFab * mfx = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{0}, lev);
amrex::MultiFab * mfy = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{1}, lev);
amrex::MultiFab * mfz = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{2}, lev);

const amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect();
const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect();
const amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect();

// avoid implicit lambda capture
auto Jx_external = m_J_external[0];
auto Jy_external = m_J_external[1];
auto Jz_external = m_J_external[2];

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() );
const amrex::Box& tbz = mfi.tilebox( z_nodal_flag, mfz->nGrowVect() );

auto const& mfxfab = mfx->array(mfi);
auto const& mfyfab = mfy->array(mfi);
auto const& mfzfab = mfz->array(mfi);

amrex::Array4<amrex::Real> lx, ly, lz;
if (EB::enabled()) {
lx = edge_lengths[0]->array(mfi);
ly = edge_lengths[1]->array(mfi);
lz = edge_lengths[2]->array(mfi);
}

amrex::ParallelFor (tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// skip if node is covered by an embedded boundary
if (lx && lx(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)
const amrex::Real x = 0._rt;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
const amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
#else
const amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real fac_y = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
const amrex::Real fac_z = (1._rt - x_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
#endif
// Initialize the x-component of the field.
mfxfab(i,j,k) = Jx_external(x,y,z,t);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// skip if node is covered by an embedded boundary
if (ly && ly(i, j, k) <= 0) { return; }

#if defined(WARPX_DIM_1D_Z)
const amrex::Real x = 0._rt;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
const amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
#elif defined(WARPX_DIM_3D)
const amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real fac_y = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
const amrex::Real fac_z = (1._rt - y_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
#endif
// Initialize the y-component of the field.
mfyfab(i,j,k) = Jy_external(x,y,z,t);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// skip if node is covered by an embedded boundary
if (lz && lz(i, j, k) <= 0) { return; }

#if defined(WARPX_DIM_1D_Z)
const amrex::Real x = 0._rt;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
const amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real y = 0._rt;
const amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
#elif defined(WARPX_DIM_3D)
const amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
const amrex::Real fac_y = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
const amrex::Real fac_z = (1._rt - z_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
#endif
// Initialize the z-component of the field.
mfzfab(i,j,k) = Jz_external(x,y,z,t);
}
);
warpx.ComputeExternalFieldOnGridUsingParser(
FieldType::hybrid_current_fp_external,
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));
}
}

Expand Down
3 changes: 1 addition & 2 deletions Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,8 +68,7 @@ void WarpX::HybridPICEvolveFields ()
const int sub_steps = m_hybrid_pic_model->m_substeps;

// Get the external current
m_hybrid_pic_model->GetCurrentExternal(
m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level));
m_hybrid_pic_model->GetCurrentExternal();

// Reference hybrid-PIC multifabs
ablastr::fields::MultiLevelScalarField rho_fp_temp = m_fields.get_mr_levels(FieldType::hybrid_rho_fp_temp, finest_level);
Expand Down
17 changes: 8 additions & 9 deletions Source/Fluids/WarpXFluidContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1010,24 +1010,23 @@ void WarpXFluidContainer::GatherAndPush (
// External field parsers
external_e_fields = (m_E_ext_s == "parse_e_ext_function");
external_b_fields = (m_B_ext_s == "parse_b_ext_function");

amrex::ParserExecutor<4> Exfield_parser;
amrex::ParserExecutor<4> Eyfield_parser;
amrex::ParserExecutor<4> Ezfield_parser;
amrex::ParserExecutor<4> Bxfield_parser;
amrex::ParserExecutor<4> Byfield_parser;
amrex::ParserExecutor<4> Bzfield_parser;

if (external_e_fields){
constexpr int num_arguments = 4; //x,y,z,t
Exfield_parser = m_Ex_parser->compile<num_arguments>();
Eyfield_parser = m_Ey_parser->compile<num_arguments>();
Ezfield_parser = m_Ez_parser->compile<num_arguments>();
Exfield_parser = m_Ex_parser->compile<4>();
Eyfield_parser = m_Ey_parser->compile<4>();
Ezfield_parser = m_Ez_parser->compile<4>();
}

if (external_b_fields){
constexpr int num_arguments = 4; //x,y,z,t
Bxfield_parser = m_Bx_parser->compile<num_arguments>();
Byfield_parser = m_By_parser->compile<num_arguments>();
Bzfield_parser = m_Bz_parser->compile<num_arguments>();
Bxfield_parser = m_Bx_parser->compile<4>();
Byfield_parser = m_By_parser->compile<4>();
Bzfield_parser = m_Bz_parser->compile<4>();
}


Expand Down
12 changes: 6 additions & 6 deletions Source/Initialization/ExternalField.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,11 +127,11 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx)
str_Bz_ext_grid_function);

Bxfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_Bx_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_Bx_ext_grid_function,{"x","y","z","t"}));
Byfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_By_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_By_ext_grid_function,{"x","y","z","t"}));
Bzfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_Bz_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_Bz_ext_grid_function,{"x","y","z","t"}));
}
//___________________________________________________________________________

Expand Down Expand Up @@ -163,11 +163,11 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx)
str_Ez_ext_grid_function);

Exfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_Ex_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_Ex_ext_grid_function,{"x","y","z","t"}));
Eyfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_Ey_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_Ey_ext_grid_function,{"x","y","z","t"}));
Ezfield_parser = std::make_unique<amrex::Parser>(
utils::parser::makeParser(str_Ez_ext_grid_function,{"x","y","z"}));
utils::parser::makeParser(str_Ez_ext_grid_function,{"x","y","z","t"}));
}
//___________________________________________________________________________

Expand Down
Loading

0 comments on commit 841d7ee

Please sign in to comment.