Skip to content

Commit

Permalink
second major refactor - create proper magnetostatic solver class
Browse files Browse the repository at this point in the history 
Signed-off-by: roelof-groenewald <[email protected]>
  • Loading branch information
@roelof-groenewald
roelof-groenewald committed Aug 30, 2024
1 parent cee7a95 commit ee2bd89
Show file tree
Hide file tree
Showing 10 changed files with 326 additions and 152 deletions.
6 changes: 3 additions & 3 deletions Python/pywarpx/fields.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -737,23 +737,23 @@ def GFPWrapper(level=0, include_ghosts=False):

def AxFPWrapper(level=0, include_ghosts=False):
return _MultiFABWrapper(
mf_name="Afield_fp_nodal[x]",
mf_name="Afield_fp[x]",
level=level,
include_ghosts=include_ghosts,
)


def AyFPWrapper(level=0, include_ghosts=False):
return _MultiFABWrapper(
mf_name="Afield_fp_nodal[y]",
mf_name="Afield_fp[y]",
level=level,
include_ghosts=include_ghosts,
)


def AzFPWrapper(level=0, include_ghosts=False):
return _MultiFABWrapper(
mf_name="Afield_fp_nodal[z]",
mf_name="Afield_fp[z]",
level=level,
include_ghosts=include_ghosts,
)
Expand Down
10 changes: 5 additions & 5 deletions Source/Diagnostics/FullDiagnostics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -671,7 +671,7 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
} else if ( m_varnames[comp] == "jz_displacement" ) {
m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(2, lev, m_crse_ratio, true);
} else if ( m_varnames[comp] == "Az" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 2), lev, m_crse_ratio);
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 2), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "rho" ){
// Initialize rho functor to dump total rho
m_all_field_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio, true);
Expand Down Expand Up @@ -720,9 +720,9 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
} else if (m_varnames[comp] == "jt_displacement" ){
m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(1, lev, m_crse_ratio, true);
} else if ( m_varnames[comp] == "Ar" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 0), lev, m_crse_ratio);
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 0), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "At" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 1), lev, m_crse_ratio);
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 1), lev, m_crse_ratio);
} else {
WARPX_ABORT_WITH_MESSAGE(m_varnames[comp] + " is not a known field output type for RZ geometry");
}
Expand All @@ -747,9 +747,9 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
} else if ( m_varnames[comp] == "jy_displacement" ){
m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(1, lev, m_crse_ratio);
} else if ( m_varnames[comp] == "Ax" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 0), lev, m_crse_ratio);
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 0), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "Ay" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 1), lev, m_crse_ratio);
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 1), lev, m_crse_ratio);
} else {
std::cout << "Error on component " << m_varnames[comp] << std::endl;
WARPX_ABORT_WITH_MESSAGE(m_varnames[comp] + " is not a known field output type for this geometry");
Expand Down
2 changes: 1 addition & 1 deletion Source/FieldSolver/Fields.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ namespace warpx::fields
F_fp,
G_fp,
phi_fp,
vector_potential_fp,
Afield_fp,
Efield_cp,
Bfield_cp,
current_cp,
Expand Down
149 changes: 142 additions & 7 deletions Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,34 +2,169 @@
*
* This file is part of WarpX.
*
* Authors: S. Eric Clark (LLNL)
* Authors: S. Eric Clark (LLNL), Roelof Groenewald (TAE)
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_MAGNETOSTATICSOLVER_H_
#define WARPX_MAGNETOSTATICSOLVER_H_

#include "MagnetostaticSolver_fwd.H"

#include "Utils/Parser/ParserUtils.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXAlgorithmSelection.H"

#include <AMReX_Array.H>
#include <AMReX_MultiFab.H>
#include <AMReX_MLMG.H>
#include <AMReX_REAL.H>

using namespace amrex;


/**
* \brief This class contains functionality involved in solving the self
* magnetic field of a plasma due to internal currents, i.e., the magnetostatic
* solver. The solver uses the plasma current to calculate the vector potential
* (in the Coulomb gauge) and evaluates B = ∇ x A.
*/
class MagnetostaticSolver
{
public:
MagnetostaticSolver (int nlevs_max); // constructor

/** Read user-defined model parameters. Called in constructor. */
void ReadParameters ();

/** Allocate member multifabs. Called in constructor. */
void AllocateMFs (int nlevs_max);
void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho,
const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag,
const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag);

/** Helper function to clear values from member multifabs. */
void ClearLevel (int lev);

namespace MagnetostaticSolver {
void InitData ();

/** Boundary Handler for the Vector Potential Poisson Solver
* This only will handle homogeneous Dirichlet boundary conditions on
* embedded boundaries, and homogeneous dirichlet/Neumann or periodic boundary conditions
* embedded boundaries, and homogeneous dirichlet/Neumann or periodic
* boundary conditions on domain boundaries.
*/
class VectorPoissonBoundaryHandler {
public:
struct BoundaryHandler {
amrex::Array<amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM>, 3> lobc, hibc;
bool bcs_set = false;
std::array<std::array<bool, AMREX_SPACEDIM * 2>, 3> dirichlet_flag;
bool has_non_periodic = false;

void defineVectorPotentialBCs ();
void defineVectorPotentialBCs (
Vector<FieldBoundaryType> field_boundary_lo,
Vector<FieldBoundaryType> field_boundary_hi,
[[maybe_unused]] Real prob_lo
) {
for (int adim = 0; adim < 3; adim++) {
#ifdef WARPX_DIM_RZ
int dim_start = 0;
if (prob_lo == 0){
lobc[adim][0] = LinOpBCType::Neumann;
dirichlet_flag[adim][0] = false;
dim_start = 1;

// handle the r_max boundary explicitly
if (field_boundary_hi[0] == FieldBoundaryType::PEC) {
if (adim == 0) {
hibc[adim][0] = LinOpBCType::Neumann;
dirichlet_flag[adim][1] = false;
} else{
hibc[adim][0] = LinOpBCType::Dirichlet;
dirichlet_flag[adim][1] = true;
}
}
else if (field_boundary_hi[0] == FieldBoundaryType::Neumann) {
hibc[adim][0] = LinOpBCType::Neumann;
dirichlet_flag[adim][1] = false;
}
}
#else
const int dim_start = 0;
#endif
bool ndotA = false;
for (int idim=dim_start; idim<AMREX_SPACEDIM; idim++){
ndotA = (adim == idim);

#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
if (idim == 1) { ndotA = (adim == 2); }
#endif

if ( field_boundary_lo[idim] == FieldBoundaryType::Periodic
&& field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
lobc[adim][idim] = LinOpBCType::Periodic;
hibc[adim][idim] = LinOpBCType::Periodic;
dirichlet_flag[adim][idim*2] = false;
dirichlet_flag[adim][idim*2+1] = false;
}
else {
has_non_periodic = true;
if ( field_boundary_lo[idim] == FieldBoundaryType::PEC ) {
if (ndotA) {
lobc[adim][idim] = LinOpBCType::Neumann;
dirichlet_flag[adim][idim*2] = false;
} else {
lobc[adim][idim] = LinOpBCType::Dirichlet;
dirichlet_flag[adim][idim*2] = true;
}
}

else if ( field_boundary_lo[idim] == FieldBoundaryType::Neumann ) {
lobc[adim][idim] = LinOpBCType::Neumann;
dirichlet_flag[adim][idim*2] = false;
}
else {
WARPX_ABORT_WITH_MESSAGE(
"Field boundary conditions have to be either periodic, PEC or neumann "
"when using the magnetostatic solver."
);
}

if ( field_boundary_hi[idim] == FieldBoundaryType::PEC ) {
if (ndotA) {
hibc[adim][idim] = LinOpBCType::Neumann;
dirichlet_flag[adim][idim*2+1] = false;
} else {
hibc[adim][idim] = LinOpBCType::Dirichlet;
dirichlet_flag[adim][idim*2+1] = true;
}
}
else if ( field_boundary_hi[idim] == FieldBoundaryType::Neumann ) {
hibc[adim][idim] = LinOpBCType::Neumann;
dirichlet_flag[adim][idim*2+1] = false;
}
else {
WARPX_ABORT_WITH_MESSAGE(
"Field boundary conditions have to be either periodic, PEC or neumann "
"when using the magnetostatic solver."
);
}
}
}
}
bcs_set = true;
}
};
} // namespace MagnetostaticSolver

BoundaryHandler m_boundary_handler = BoundaryHandler();

// numerical parameters
amrex::Real required_precision = 1.e-11_rt;
amrex::Real absolute_tolerance = 0.0_rt;
int max_iters = 200;
int verbosity = 2;

// Declare multifabs specifically needed for the magnetostatic solver
Vector<std::array< std::unique_ptr<MultiFab>, 3 > > Afield_fp_nodal;
Vector<std::array< std::unique_ptr<MultiFab>, 3 > > current_fp_temp;
};
#endif //WARPX_MAGNETOSTATICSOLVER_H_
Loading

0 comments on commit ee2bd89

Please sign in to comment.