Skip to content

Commit

Permalink
Merged and refactored external fields to use new MultiFab register in…
Browse files Browse the repository at this point in the history 
… Cylindrical and Cartesian.
  • Loading branch information
@clarkse
clarkse committed Sep 25, 2024
1 parent ba9678a commit 344a3dd
Show file tree
Hide file tree
Showing 6 changed files with 117 additions and 192 deletions.
5 changes: 4 additions & 1 deletion Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,10 @@
/* Copyright 2020 Remi Lehe
/* Copyright 2020-2024 The WarpX Community
*
* This file is part of WarpX.
*
* Authors: Remi Lehe (LBNL)
* S. Eric Clark (Helion Energy)
*
* License: BSD-3-Clause-LBNL
*/

Expand Down
97 changes: 39 additions & 58 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
/* Copyright 2023 The WarpX Community
/* Copyright 2023-2024 The WarpX Community
*
* This file is part of WarpX.
*
* Authors: Roelof Groenewald (TAE Technologies)
* S. Eric Clark (Helion Energy)
*
* License: BSD-3-Clause-LBNL
*/
Expand All @@ -12,6 +13,8 @@

#include "HybridPICModel_fwd.H"

#include "Fields.H"

#include "Utils/WarpXAlgorithmSelection.H"

#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
Expand Down Expand Up @@ -42,8 +45,8 @@ public:
void ReadParameters ();

/** Allocate hybrid-PIC specific multifabs. Called in constructor. */
void AllocateLevelMFs (ablastr::fields::MultiFabRegister & fields,
void AllocateLevelMFs (
ablastr::fields::MultiFabRegister & fields,
int lev,
const amrex::BoxArray& ba,
const amrex::DistributionMapping& dm,
Expand Down Expand Up @@ -72,41 +75,19 @@ 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 (bool skip_check = false);

void GetCurrentExternal (
ablastr::fields::VectorField const& edge_lengths,
int lev
);

void GetFieldsExternal (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths
);

void GetFieldsExternal (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
amrex::Real t);

void GetFieldsExternal (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev
);
void GetFieldsExternal (amrex::Real t);

void GetExternalFieldFromExpression (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
warpx::fields::FieldType field_type,
std::array< amrex::ParserExecutor<4>, 3> const& expression,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev
);

int lev);

void GetExternalFieldFromExpression (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
warpx::fields::FieldType field_type,
std::array< amrex::ParserExecutor<4>, 3> const& expression,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev, amrex::Real t
);
int lev, amrex::Real t);

/**
* \brief
Expand Down Expand Up @@ -252,35 +233,35 @@ public:
std::array< amrex::ParserExecutor<4>, 3> m_E_external;

// Declare multifabs specifically needed for the hybrid-PIC model
amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_fp_temp;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_ampere;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_external;
amrex::Vector< std::unique_ptr<amrex::MultiFab> > electron_pressure_fp;

amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_hyb_external;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_hyb_external;
// amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_fp_temp;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_ampere;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_external;
// amrex::Vector< std::unique_ptr<amrex::MultiFab> > electron_pressure_fp;

// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_hyb_external;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_hyb_external;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_hyb_self;
// amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_hyb_self;

// Helper functions to retrieve hybrid-PIC multifabs
[[nodiscard]] amrex::MultiFab*
get_pointer_current_fp_ampere (int lev, int direction) const
{
return current_fp_ampere[lev][direction].get();
}

[[nodiscard]] amrex::MultiFab*
get_pointer_current_fp_external (int lev, int direction) const
{
return current_fp_external[lev][direction].get();
}

[[nodiscard]] amrex::MultiFab*
get_pointer_electron_pressure_fp (int lev) const
{
return electron_pressure_fp[lev].get();
}
// // Helper functions to retrieve hybrid-PIC multifabs
// [[nodiscard]] amrex::MultiFab*
// get_pointer_current_fp_ampere (int lev, int direction) const
// {
// return current_fp_ampere[lev][direction].get();
// }

// [[nodiscard]] amrex::MultiFab*
// get_pointer_current_fp_external (int lev, int direction) const
// {
// return current_fp_external[lev][direction].get();
// }

// [[nodiscard]] amrex::MultiFab*
// get_pointer_electron_pressure_fp (int lev) const
// {
// return electron_pressure_fp[lev].get();
// }

/** Gpu Vector with index type of the Jx multifab */
amrex::GpuArray<int, 3> Jx_IndexType;
Expand Down
156 changes: 45 additions & 111 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,17 @@
/* Copyright 2023 The WarpX Community
/* Copyright 2023-2024 The WarpX Community
*
* This file is part of WarpX.
*
* Authors: Roelof Groenewald (TAE Technologies)
* S. Eric Clark (Helion Energy)
*
* License: BSD-3-Clause-LBNL
*/

#include "HybridPICModel.H"

#include "EmbeddedBoundary/Enabled.H"
#include "FieldSolver/Fields.H"
#include "Fields.H"
#include "Particles/MultiParticleContainer.H"
#include "WarpX.H"

Expand Down Expand Up @@ -68,30 +69,8 @@ void HybridPICModel::ReadParameters ()
}
}

void HybridPICModel::AllocateLevelMFs (ablastr::fields::MultiFabRegister & fields,
int lev, const BoxArray& ba, const DistributionMapping& dm,
const int ncomps, const IntVect& ngJ, const IntVect& ngRho,
const IntVect& jx_nodal_flag,
const IntVect& jy_nodal_flag,
const IntVect& jz_nodal_flag,
const IntVect& rho_nodal_flag)
void HybridPICModel::AllocateMFs (int nlevs_max)
{
electron_pressure_fp.resize(nlevs_max);
rho_fp_temp.resize(nlevs_max);
current_fp_temp.resize(nlevs_max);
current_fp_ampere.resize(nlevs_max);
current_fp_external.resize(nlevs_max);

if (m_add_external_fields) {
Bfield_hyb_external.resize(nlevs_max);
Efield_hyb_external.resize(nlevs_max);
// Bfield_hyb_self.resize(nlevs_max);
// Efield_hyb_self.resize(nlevs_max);
}
}

void HybridPICModel::AllocateLevelMFs (
ablastr::fields::MultiFabRegister & fields,
int lev, const BoxArray& ba, const DistributionMapping& dm,
const int ncomps,
const IntVect& ngJ, const IntVect& ngRho,
Expand Down Expand Up @@ -158,18 +137,24 @@ void HybridPICModel::AllocateLevelMFs (

if (m_add_external_fields) {
// These are nodal to match when B-field is added in evaluation of Ohm's law
WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][0], amrex::convert(ba, Bx_nodal_flag),
dm, ncomps, ngB, lev, "Bfield_hyb_external[x]", 0.0_rt);
WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][1], amrex::convert(ba, By_nodal_flag),
dm, ncomps, ngB, lev, "Bfield_hyb_external[y]", 0.0_rt);
WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][2], amrex::convert(ba, Bz_nodal_flag),
dm, ncomps, ngB, lev, "Bfield_hyb_external[z]", 0.0_rt);
WarpX::AllocInitMultiFab(Efield_hyb_external[lev][0], amrex::convert(ba, Ex_nodal_flag),
dm, ncomps, ngE, lev, "Efield_hyb_external[x]", 0.0_rt);
WarpX::AllocInitMultiFab(Efield_hyb_external[lev][1], amrex::convert(ba, Ey_nodal_flag),
dm, ncomps, ngE, lev, "Efield_hyb_external[y]", 0.0_rt);
WarpX::AllocInitMultiFab(Efield_hyb_external[lev][2], amrex::convert(ba, Ez_nodal_flag),
dm, ncomps, ngE, lev, "Efield_hyb_external[z]", 0.0_rt);
fields.alloc_init(FieldType::hybrid_B_fp_external, Direction{0},
lev, amrex::convert(ba, Bx_nodal_flag),
dm, ncomps, ngB, 0.0_rt);
fields.alloc_init(FieldType::hybrid_B_fp_external, Direction{1},
lev, amrex::convert(ba, By_nodal_flag),
dm, ncomps, ngB, 0.0_rt);
fields.alloc_init(FieldType::hybrid_B_fp_external, Direction{2},
lev, amrex::convert(ba, Bz_nodal_flag),
dm, ncomps, ngB, 0.0_rt);
fields.alloc_init(FieldType::hybrid_E_fp_external, Direction{0},
lev, amrex::convert(ba, Ex_nodal_flag),
dm, ncomps, ngE, 0.0_rt);
fields.alloc_init(FieldType::hybrid_E_fp_external, Direction{1},
lev, amrex::convert(ba, Ey_nodal_flag),
dm, ncomps, ngE, 0.0_rt);
fields.alloc_init(FieldType::hybrid_E_fp_external, Direction{2},
lev, amrex::convert(ba, Ez_nodal_flag),
dm, ncomps, ngE, 0.0_rt);
}

#ifdef WARPX_DIM_RZ
Expand All @@ -179,21 +164,6 @@ void HybridPICModel::AllocateLevelMFs (
#endif
}

void HybridPICModel::ClearLevel (int lev)
{
electron_pressure_fp[lev].reset();
rho_fp_temp[lev].reset();
for (int i = 0; i < 3; ++i) {
current_fp_temp[lev][i].reset();
current_fp_ampere[lev][i].reset();
current_fp_external[lev][i].reset();
if (m_add_external_fields) {
Bfield_hyb_external[lev][i].reset();
Efield_hyb_external[lev][i].reset();
}
}
}

void HybridPICModel::InitData ()
{
m_resistivity_parser = std::make_unique<amrex::Parser>(
Expand Down Expand Up @@ -315,96 +285,60 @@ void HybridPICModel::InitData ()
// Initialize external current - note that this approach skips the check
// 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);
GetFieldsExternal(ablastr::fields::a2m(edge_lengths), lev);
}
GetCurrentExternal(true);
GetFieldsExternal(warpx.gett_new(0));
}

void HybridPICModel::GetCurrentExternal (
ablastr::fields::MultiLevelVectorField const& edge_lengths)
void HybridPICModel::GetCurrentExternal (bool skip_check /*=false*/)
{
if (!m_external_field_has_time_dependence) { return; }
if (!skip_check && !m_external_field_has_time_dependence) { return; }

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

void HybridPICModel::GetCurrentExternal (
ablastr::fields::VectorField const& edge_lengths,
int lev)
{
if (!m_external_field_has_time_dependence) { return; }
GetExternalFieldFromExpression(current_fp_external[lev], m_J_external, edge_lengths, lev);
}

void HybridPICModel::GetFieldsExternal (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
amrex::Real t)
void HybridPICModel::GetFieldsExternal (amrex::Real t)
{
auto& warpx = WarpX::GetInstance();
for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
{
GetExternalFieldFromExpression(Bfield_hyb_external[lev], m_B_external, edge_lengths[lev], lev, t);
GetExternalFieldFromExpression(Efield_hyb_external[lev], m_E_external, edge_lengths[lev], lev, t);
GetExternalFieldFromExpression(
FieldType::hybrid_B_fp_external,
m_B_external, lev, t);
GetExternalFieldFromExpression(
FieldType::hybrid_E_fp_external,
m_E_external, lev, t);
}
}

void HybridPICModel::GetFieldsExternal (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev)
{
GetExternalFieldFromExpression(Bfield_hyb_external[lev], m_B_external, edge_lengths, lev);
GetExternalFieldFromExpression(Efield_hyb_external[lev], m_E_external, edge_lengths, lev);
}

void HybridPICModel::GetExternalFieldFromExpression (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
FieldType field_type,
std::array< amrex::ParserExecutor<4>, 3> const& expression,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev)
{
auto & warpx = WarpX::GetInstance();
auto t = warpx.gett_new(lev);
GetExternalFieldFromExpression(field, expression, edge_lengths, lev, t);
GetExternalFieldFromExpression(field_type, expression, lev, warpx.gett_new(lev));
}

void HybridPICModel::GetExternalFieldFromExpression (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
FieldType field_type,
std::array< amrex::ParserExecutor<4>, 3> const& expression,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
int lev, amrex::Real t)
{
using ablastr::fields::Direction;

// This logic matches closely to WarpX::InitializeExternalFieldsOnGridUsingParser
// except that the parsers include time dependence.
auto & warpx = WarpX::GetInstance();
auto dx_lev = warpx.Geom(lev).CellSizeArray();
const RealBox& real_box = warpx.Geom(lev).ProbDomain();

auto& mfx = field[0];
auto& mfy = field[1];
auto& mfz = field[2];
auto const& mfx = warpx.m_fields.get(field_type, Direction{0}, lev);
auto const& mfy = warpx.m_fields.get(field_type, Direction{1}, lev);
auto const& mfz = warpx.m_fields.get(field_type, Direction{2}, lev);

const amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect();
const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect();
Expand All @@ -427,9 +361,9 @@ void HybridPICModel::GetExternalFieldFromExpression (

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);
lx = warpx.m_fields.get(FieldType::edge_lengths, Direction{0}, lev)->array(mfi);
ly = warpx.m_fields.get(FieldType::edge_lengths, Direction{1}, lev)->array(mfi);
lz = warpx.m_fields.get(FieldType::edge_lengths, Direction{2}, lev)->array(mfi);
}

amrex::ParallelFor (tbx, tby, tbz,
Expand Down
Loading

0 comments on commit 344a3dd

Please sign in to comment.