Skip to content

Commit

Permalink
Use new stair-case approximation in hybrid solver (#5558)
Browse files Browse the repository at this point in the history 
Merge #5534 first.

This extends the changes of #5534 and #5574 to the hybrid solver.
Note that this effectively changes the definition of the stair-cased
embedded boundary for the hybrid solver.

---------

Co-authored-by: Roelof Groenewald <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
  • Loading branch information
@RemiLehe @roelof-groenewald @pre-commit-ci
3 people authored Jan 22, 2025
1 parent 88412f3 commit eaaaef5
Show file tree
Hide file tree
Showing 8 changed files with 128 additions and 113 deletions.
18 changes: 9 additions & 9 deletions Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ class FiniteDifferenceSolver
* \param[out] Efield vector of electric field MultiFabs updated at a given level
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] Jfield vector of current density MultiFabs at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] eb_update_E indicate in which cell E should be updated (related to embedded boundaries)
* \param[in] dt timestep of the simulation
* \param[in] macroscopic_properties contains user-defined properties of the medium.
*/
Expand Down Expand Up @@ -147,7 +147,7 @@ class FiniteDifferenceSolver
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] rhofield scalar ion charge density Multifab at a given level
* \param[in] Pefield scalar electron pressure MultiFab at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] eb_update_E indicate in which cell E should be updated (related to embedded boundaries)
* \param[in] lev level number for the calculation
* \param[in] hybrid_model instance of the hybrid-PIC model
* \param[in] solve_for_Faraday boolean flag for whether the E-field is solved to be used in Faraday's equation
Expand All @@ -158,7 +158,7 @@ class FiniteDifferenceSolver
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
amrex::MultiFab const& Pefield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev, HybridPICModel const* hybrid_model,
bool solve_for_Faraday );

Expand All @@ -168,13 +168,13 @@ class FiniteDifferenceSolver
*
* \param[out] Jfield vector of current MultiFabs at a given level
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] eb_update_E indicate in which cell E should be updated (related to embedded boundaries)
* \param[in] lev level number for the calculation
*/
void CalculateCurrentAmpere (
ablastr::fields::VectorField& Jfield,
ablastr::fields::VectorField const& Bfield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev );

private:
Expand Down Expand Up @@ -243,15 +243,15 @@ class FiniteDifferenceSolver
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
amrex::MultiFab const& Pefield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev, HybridPICModel const* hybrid_model,
bool solve_for_Faraday );

template<typename T_Algo>
void CalculateCurrentAmpereCylindrical (
ablastr::fields::VectorField& Jfield,
ablastr::fields::VectorField const& Bfield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev
);

Expand Down Expand Up @@ -347,15 +347,15 @@ class FiniteDifferenceSolver
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
amrex::MultiFab const& Pefield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev, HybridPICModel const* hybrid_model,
bool solve_for_Faraday );

template<typename T_Algo>
void CalculateCurrentAmpereCartesian (
ablastr::fields::VectorField& Jfield,
ablastr::fields::VectorField const& Bfield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
int lev
);
#endif
Expand Down
18 changes: 9 additions & 9 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,15 +62,15 @@ public:
* subtracted as well. Used in the Ohm's law solver (kinetic-fluid hybrid model).
*
* \param[in] Bfield Magnetic field from which the current is calculated.
* \param[in] edge_lengths Length of cell edges taking embedded boundaries into account
* \param[in] eb_update_E Indicate in which cell J should be calculated (related to embedded boundaries).
*/
void CalculatePlasmaCurrent (
ablastr::fields::MultiLevelVectorField const& Bfield,
ablastr::fields::MultiLevelVectorField const& edge_lengths
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E
);
void CalculatePlasmaCurrent (
ablastr::fields::VectorField const& Bfield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
int lev
);

Expand All @@ -83,31 +83,31 @@ public:
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelVectorField const& Bfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
bool solve_for_Faraday) const;

void HybridPICSolveE (
ablastr::fields::VectorField const& Efield,
ablastr::fields::VectorField const& Jfield,
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
int lev, bool solve_for_Faraday) const;

void HybridPICSolveE (
ablastr::fields::VectorField const& Efield,
ablastr::fields::VectorField const& Jfield,
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
int lev, PatchType patch_type, bool solve_for_Faraday) const;

void BfieldEvolveRK (
ablastr::fields::MultiLevelVectorField const& Bfield,
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, DtType a_dt_type,
amrex::IntVect ng, std::optional<bool> nodal_sync);

Expand All @@ -116,7 +116,7 @@ public:
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, int lev, DtType dt_type,
amrex::IntVect ng, std::optional<bool> nodal_sync);

Expand All @@ -125,7 +125,7 @@ public:
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, DtType dt_type,
amrex::IntVect ng, std::optional<bool> nodal_sync);

Expand Down
40 changes: 20 additions & 20 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,26 +250,26 @@ void HybridPICModel::GetCurrentExternal ()

void HybridPICModel::CalculatePlasmaCurrent (
ablastr::fields::MultiLevelVectorField const& Bfield,
ablastr::fields::MultiLevelVectorField const& edge_lengths)
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E)
{
auto& warpx = WarpX::GetInstance();
for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
{
CalculatePlasmaCurrent(Bfield[lev], edge_lengths[lev], lev);
CalculatePlasmaCurrent(Bfield[lev], eb_update_E[lev], lev);
}
}

void HybridPICModel::CalculatePlasmaCurrent (
ablastr::fields::VectorField const& Bfield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
const int lev)
{
WARPX_PROFILE("HybridPICModel::CalculatePlasmaCurrent()");

auto& warpx = WarpX::GetInstance();
ablastr::fields::VectorField current_fp_plasma = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_plasma, lev);
warpx.get_pointer_fdtd_solver_fp(lev)->CalculateCurrentAmpere(
current_fp_plasma, Bfield, edge_lengths, lev
current_fp_plasma, Bfield, eb_update_E, lev
);

// we shouldn't apply the boundary condition to J since J = J_i - J_e but
Expand All @@ -293,15 +293,15 @@ void HybridPICModel::HybridPICSolveE (
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelVectorField const& Bfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
const bool solve_for_Faraday) const
{
auto& warpx = WarpX::GetInstance();
for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
{
HybridPICSolveE(
Efield[lev], Jfield[lev], Bfield[lev], *rhofield[lev],
edge_lengths[lev], lev, solve_for_Faraday
eb_update_E[lev], lev, solve_for_Faraday
);
}
}
Expand All @@ -311,13 +311,13 @@ void HybridPICModel::HybridPICSolveE (
ablastr::fields::VectorField const& Jfield,
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
const int lev, const bool solve_for_Faraday) const
{
WARPX_PROFILE("WarpX::HybridPICSolveE()");

HybridPICSolveE(
Efield, Jfield, Bfield, rhofield, edge_lengths, lev,
Efield, Jfield, Bfield, rhofield, eb_update_E, lev,
PatchType::fine, solve_for_Faraday
);
if (lev > 0)
Expand All @@ -332,7 +332,7 @@ void HybridPICModel::HybridPICSolveE (
ablastr::fields::VectorField const& Jfield,
ablastr::fields::VectorField const& Bfield,
amrex::MultiFab const& rhofield,
ablastr::fields::VectorField const& edge_lengths,
std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
const int lev, PatchType patch_type,
const bool solve_for_Faraday) const
{
Expand All @@ -344,7 +344,7 @@ void HybridPICModel::HybridPICSolveE (
// Solve E field in regular cells
warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE(
Efield, current_fp_plasma, Jfield, Bfield, rhofield,
*electron_pressure_fp, edge_lengths, lev, this, solve_for_Faraday
*electron_pressure_fp, eb_update_E, lev, this, solve_for_Faraday
);
amrex::Real const time = warpx.gett_old(0) + warpx.getdt(0);
warpx.ApplyEfieldBoundary(lev, patch_type, time);
Expand Down Expand Up @@ -411,15 +411,15 @@ void HybridPICModel::BfieldEvolveRK (
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, DtType dt_type,
IntVect ng, std::optional<bool> nodal_sync )
{
auto& warpx = WarpX::GetInstance();
for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
{
BfieldEvolveRK(
Bfield, Efield, Jfield, rhofield, edge_lengths, dt, lev, dt_type,
Bfield, Efield, Jfield, rhofield, eb_update_E, dt, lev, dt_type,
ng, nodal_sync
);
}
Expand All @@ -430,7 +430,7 @@ void HybridPICModel::BfieldEvolveRK (
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, int lev, DtType dt_type,
IntVect ng, std::optional<bool> nodal_sync )
{
Expand All @@ -457,7 +457,7 @@ void HybridPICModel::BfieldEvolveRK (
// The Runge-Kutta scheme begins here.
// Step 1:
FieldPush(
Bfield, Efield, Jfield, rhofield, edge_lengths,
Bfield, Efield, Jfield, rhofield, eb_update_E,
0.5_rt*dt, dt_type, ng, nodal_sync
);

Expand All @@ -473,7 +473,7 @@ void HybridPICModel::BfieldEvolveRK (

// Step 2:
FieldPush(
Bfield, Efield, Jfield, rhofield, edge_lengths,
Bfield, Efield, Jfield, rhofield, eb_update_E,
0.5_rt*dt, dt_type, ng, nodal_sync
);

Expand All @@ -493,7 +493,7 @@ void HybridPICModel::BfieldEvolveRK (

// Step 3:
FieldPush(
Bfield, Efield, Jfield, rhofield, edge_lengths,
Bfield, Efield, Jfield, rhofield, eb_update_E,
dt, dt_type, ng, nodal_sync
);

Expand All @@ -509,7 +509,7 @@ void HybridPICModel::BfieldEvolveRK (

// Step 4:
FieldPush(
Bfield, Efield, Jfield, rhofield, edge_lengths,
Bfield, Efield, Jfield, rhofield, eb_update_E,
0.5_rt*dt, dt_type, ng, nodal_sync
);

Expand Down Expand Up @@ -543,16 +543,16 @@ void HybridPICModel::FieldPush (
ablastr::fields::MultiLevelVectorField const& Efield,
ablastr::fields::MultiLevelVectorField const& Jfield,
ablastr::fields::MultiLevelScalarField const& rhofield,
ablastr::fields::MultiLevelVectorField const& edge_lengths,
amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
amrex::Real dt, DtType dt_type,
IntVect ng, std::optional<bool> nodal_sync )
{
auto& warpx = WarpX::GetInstance();

// Calculate J = curl x B / mu0 - J_ext
CalculatePlasmaCurrent(Bfield, edge_lengths);
CalculatePlasmaCurrent(Bfield, eb_update_E);
// Calculate the E-field from Ohm's law
HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
HybridPICSolveE(Efield, Jfield, Bfield, rhofield, eb_update_E, true);
warpx.FillBoundaryE(ng, nodal_sync);
// Push forward the B-field using Faraday's law
amrex::Real const t_old = warpx.gett_old(0);
Expand Down
Loading

0 comments on commit eaaaef5

Please sign in to comment.