Skip to content

Commit

Permalink
Add SemiImplicitPicard scheme
Browse files Browse the repository at this point in the history
  • Loading branch information
@dpgrote
dpgrote committed Nov 13, 2023
1 parent a78579b commit dc7b37a
Show file tree
Hide file tree
Showing 7 changed files with 78 additions and 50 deletions.
3 changes: 3 additions & 0 deletions Docs/source/usage/parameters.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,9 @@ Overall simulation parameters
:math:`\omega_{pe} \delta t` is close to or greater than one.
The method is described in `Angus et al., On numerical energy conservation for an implicit particle-in-cell method coupled with a binary Monte-Carlo algorithm for Coulomb collisions <https://doi.org/10.1016/j.jcp.2022.111030>`__.

* ``semi_implicit_picard``: Use an energy conserving semi-implicit solver that uses a Picard iteration to solve the system.
Note that this method has the CFL limitation :math:`\omega_{pe} \delta t < 1`.
The method is described in `Chen et al., A semi-implicit, energy- and charge-conserving particle-in-cell algorithm for the relativistic Vlasov-Maxwell equations <https://doi.org/10.1016/j.jcp.2020.109228>`__.

* ``warpx.do_electrostatic`` (`string`) optional (default `none`)
Specifies the electrostatic mode. When turned on, instead of updating
Expand Down
3 changes: 2 additions & 1 deletion Source/Evolve/WarpXEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ WarpX::Evolve (int numsteps)
{
if (evolve_scheme == EvolveScheme::Explicit) {
EvolveExplicit(numsteps);
} else if (evolve_scheme == EvolveScheme::ImplicitPicard) {
} else if (evolve_scheme == EvolveScheme::ImplicitPicard ||
evolve_scheme == EvolveScheme::SemiImplicitPicard) {
EvolveImplicitPicard(numsteps);
} else {
amrex::Abort(Utils::TextMsg::Err("Unknown evolve scheme"));
Expand Down
104 changes: 62 additions & 42 deletions Source/Evolve/WarpXEvolveImplicitPicard.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,26 +77,30 @@ WarpX::EvolveImplicitPicardInit (const int lev)
// Only one refinement level is supported
const int nlevs_max = maxLevel() + 1;
Efield_n.resize(nlevs_max);
Bfield_n.resize(nlevs_max);
Efield_save.resize(nlevs_max);
Bfield_save.resize(nlevs_max);
if (evolve_scheme == EvolveScheme::ImplicitPicard) {
Bfield_n.resize(nlevs_max);
Bfield_save.resize(nlevs_max);
}

// Strange, the WarpX::DistributionMap(0) is not consistent with Ex_fp.DistributionMap()???

// Note that the *_fp will be the n+theta and n+1 time level
AllocInitMultiFabFromModel(Efield_n[lev][0], *Efield_fp[0][0], lev, "Efield_n[0]");
AllocInitMultiFabFromModel(Efield_n[lev][1], *Efield_fp[0][1], lev, "Efield_n[1]");
AllocInitMultiFabFromModel(Efield_n[lev][2], *Efield_fp[0][2], lev, "Efield_n[2]");
AllocInitMultiFabFromModel(Bfield_n[lev][0], *Bfield_fp[0][0], lev, "Bfield_n[0]");
AllocInitMultiFabFromModel(Bfield_n[lev][1], *Bfield_fp[0][1], lev, "Bfield_n[1]");
AllocInitMultiFabFromModel(Bfield_n[lev][2], *Bfield_fp[0][2], lev, "Bfield_n[2]");

AllocInitMultiFabFromModel(Efield_save[lev][0], *Efield_fp[0][0], lev, "Efield_save[0]");
AllocInitMultiFabFromModel(Efield_save[lev][1], *Efield_fp[0][1], lev, "Efield_save[1]");
AllocInitMultiFabFromModel(Efield_save[lev][2], *Efield_fp[0][2], lev, "Efield_save[2]");
AllocInitMultiFabFromModel(Bfield_save[lev][0], *Bfield_fp[0][0], lev, "Bfield_save[0]");
AllocInitMultiFabFromModel(Bfield_save[lev][1], *Bfield_fp[0][1], lev, "Bfield_save[1]");
AllocInitMultiFabFromModel(Bfield_save[lev][2], *Bfield_fp[0][2], lev, "Bfield_save[2]");

if (evolve_scheme == EvolveScheme::ImplicitPicard) {
AllocInitMultiFabFromModel(Bfield_n[lev][0], *Bfield_fp[0][0], lev, "Bfield_n[0]");
AllocInitMultiFabFromModel(Bfield_n[lev][1], *Bfield_fp[0][1], lev, "Bfield_n[1]");
AllocInitMultiFabFromModel(Bfield_n[lev][2], *Bfield_fp[0][2], lev, "Bfield_n[2]");
AllocInitMultiFabFromModel(Bfield_save[lev][0], *Bfield_fp[0][0], lev, "Bfield_save[0]");
AllocInitMultiFabFromModel(Bfield_save[lev][1], *Bfield_fp[0][1], lev, "Bfield_save[1]");
AllocInitMultiFabFromModel(Bfield_save[lev][2], *Bfield_fp[0][2], lev, "Bfield_save[2]");
}

}

Expand Down Expand Up @@ -178,9 +182,17 @@ WarpX::EvolveImplicitPicard (int numsteps)
amrex::MultiFab::Copy(*Efield_n[0][0], *Efield_fp[0][0], 0, 0, ncomps, Efield_fp[0][0]->nGrowVect());
amrex::MultiFab::Copy(*Efield_n[0][1], *Efield_fp[0][1], 0, 0, ncomps, Efield_fp[0][1]->nGrowVect());
amrex::MultiFab::Copy(*Efield_n[0][2], *Efield_fp[0][2], 0, 0, ncomps, Efield_fp[0][2]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_n[0][0], *Bfield_fp[0][0], 0, 0, ncomps, Bfield_fp[0][0]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_n[0][1], *Bfield_fp[0][1], 0, 0, ncomps, Bfield_fp[0][1]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_n[0][2], *Bfield_fp[0][2], 0, 0, ncomps, Bfield_fp[0][2]->nGrowVect());
if (evolve_scheme == EvolveScheme::ImplicitPicard) {
amrex::MultiFab::Copy(*Bfield_n[0][0], *Bfield_fp[0][0], 0, 0, ncomps, Bfield_fp[0][0]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_n[0][1], *Bfield_fp[0][1], 0, 0, ncomps, Bfield_fp[0][1]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_n[0][2], *Bfield_fp[0][2], 0, 0, ncomps, Bfield_fp[0][2]->nGrowVect());
} else if (evolve_scheme == EvolveScheme::SemiImplicitPicard) {
// This updates Bfield_fp so it holds the new B at n+1/2
EvolveB(dt[0], DtType::Full);
// WarpX::sync_nodal_points is used to avoid instability
FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
ApplyBfieldBoundary(0, PatchType::fine, DtType::Full);
}

// Start the iterations
amrex::Real deltaE = 1._rt;
Expand Down Expand Up @@ -219,24 +231,26 @@ WarpX::EvolveImplicitPicard (int numsteps)
FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
ApplyEfieldBoundary(0, PatchType::fine);

if (picard_iteration_tolerance > 0. || iteration_count == max_picard_iterations) {
// Save the B at n+1/2 from the previous iteration so that the change
// in this iteration can be calculated
amrex::MultiFab::Copy(*Bfield_save[0][0], *Bfield_fp[0][0], 0, 0, ncomps, 0);
amrex::MultiFab::Copy(*Bfield_save[0][1], *Bfield_fp[0][1], 0, 0, ncomps, 0);
amrex::MultiFab::Copy(*Bfield_save[0][2], *Bfield_fp[0][2], 0, 0, ncomps, 0);
}
if (evolve_scheme == EvolveScheme::ImplicitPicard) {
if (picard_iteration_tolerance > 0. || iteration_count == max_picard_iterations) {
// Save the B at n+1/2 from the previous iteration so that the change
// in this iteration can be calculated
amrex::MultiFab::Copy(*Bfield_save[0][0], *Bfield_fp[0][0], 0, 0, ncomps, 0);
amrex::MultiFab::Copy(*Bfield_save[0][1], *Bfield_fp[0][1], 0, 0, ncomps, 0);
amrex::MultiFab::Copy(*Bfield_save[0][2], *Bfield_fp[0][2], 0, 0, ncomps, 0);
}

// Copy Bfield_n into Bfield_fp since EvolveB updates Bfield_fp in place
amrex::MultiFab::Copy(*Bfield_fp[0][0], *Bfield_n[0][0], 0, 0, ncomps, Bfield_n[0][0]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_fp[0][1], *Bfield_n[0][1], 0, 0, ncomps, Bfield_n[0][1]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_fp[0][2], *Bfield_n[0][2], 0, 0, ncomps, Bfield_n[0][2]->nGrowVect());
// Copy Bfield_n into Bfield_fp since EvolveB updates Bfield_fp in place
amrex::MultiFab::Copy(*Bfield_fp[0][0], *Bfield_n[0][0], 0, 0, ncomps, Bfield_n[0][0]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_fp[0][1], *Bfield_n[0][1], 0, 0, ncomps, Bfield_n[0][1]->nGrowVect());
amrex::MultiFab::Copy(*Bfield_fp[0][2], *Bfield_n[0][2], 0, 0, ncomps, Bfield_n[0][2]->nGrowVect());

// This updates Bfield_fp so it holds the new B at n+1/2
EvolveB(0.5_rt*dt[0], DtType::Full);
// WarpX::sync_nodal_points is used to avoid instability
FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
ApplyBfieldBoundary(0, PatchType::fine, DtType::Full);
// This updates Bfield_fp so it holds the new B at n+1/2
EvolveB(0.5_rt*dt[0], DtType::Full);
// WarpX::sync_nodal_points is used to avoid instability
FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
ApplyBfieldBoundary(0, PatchType::fine, DtType::Full);
}

// The B field update needs
if (num_mirrors>0){
Expand All @@ -259,17 +273,21 @@ WarpX::EvolveImplicitPicard (int numsteps)
amrex::Real deltaE2 = Efield_save[0][2]->norm0(0, 0)/maxE2;
/* amrex::Print() << "deltaE " << iteration_count << " " << deltaE0 << " " << deltaE1 << " " << deltaE2 << "\n"; */
deltaE = std::max(std::max(deltaE0, deltaE1), deltaE2);
Bfield_save[0][0]->minus(*Bfield_fp[0][0], 0, ncomps, 0);
Bfield_save[0][1]->minus(*Bfield_fp[0][1], 0, ncomps, 0);
Bfield_save[0][2]->minus(*Bfield_fp[0][2], 0, ncomps, 0);
amrex::Real maxB0 = std::max(1._rt, Bfield_fp[0][0]->norm0(0, 0));
amrex::Real maxB1 = std::max(1._rt, Bfield_fp[0][1]->norm0(0, 0));
amrex::Real maxB2 = std::max(1._rt, Bfield_fp[0][2]->norm0(0, 0));
amrex::Real deltaB0 = Bfield_save[0][0]->norm0(0, 0)/maxB0;
amrex::Real deltaB1 = Bfield_save[0][1]->norm0(0, 0)/maxB1;
amrex::Real deltaB2 = Bfield_save[0][2]->norm0(0, 0)/maxB2;
/* amrex::Print() << "deltaB " << iteration_count << " " << deltaB0 << " " << deltaB1 << " " << deltaB2 << "\n"; */
deltaB = std::max(std::max(deltaB0, deltaB1), deltaB2);
if (evolve_scheme == EvolveScheme::ImplicitPicard) {
Bfield_save[0][0]->minus(*Bfield_fp[0][0], 0, ncomps, 0);
Bfield_save[0][1]->minus(*Bfield_fp[0][1], 0, ncomps, 0);
Bfield_save[0][2]->minus(*Bfield_fp[0][2], 0, ncomps, 0);
amrex::Real maxB0 = std::max(1._rt, Bfield_fp[0][0]->norm0(0, 0));
amrex::Real maxB1 = std::max(1._rt, Bfield_fp[0][1]->norm0(0, 0));
amrex::Real maxB2 = std::max(1._rt, Bfield_fp[0][2]->norm0(0, 0));
amrex::Real deltaB0 = Bfield_save[0][0]->norm0(0, 0)/maxB0;
amrex::Real deltaB1 = Bfield_save[0][1]->norm0(0, 0)/maxB1;
amrex::Real deltaB2 = Bfield_save[0][2]->norm0(0, 0)/maxB2;
/* amrex::Print() << "deltaB " << iteration_count << " " << deltaB0 << " " << deltaB1 << " " << deltaB2 << "\n"; */
deltaB = std::max(std::max(deltaB0, deltaB1), deltaB2);
} else {
deltaB = 0.;
}
amrex::Print() << "Max delta " << iteration_count << " " << deltaE << " " << deltaB << "\n";
}

Expand All @@ -285,11 +303,13 @@ WarpX::EvolveImplicitPicard (int numsteps)
}

// Advance fields to step n+1
FinishImplicitFieldUpdate(Efield_fp, Efield_n);
FinishImplicitFieldUpdate(Bfield_fp, Bfield_n);
// WarpX::sync_nodal_points is used to avoid instability
FinishImplicitFieldUpdate(Efield_fp, Efield_n);
FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
if (evolve_scheme == EvolveScheme::ImplicitPicard) {
FinishImplicitFieldUpdate(Bfield_fp, Bfield_n);
FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
}

// Run multi-physics modules:
// ionization, Coulomb collisions, QED
Expand Down
3 changes: 2 additions & 1 deletion Source/Utils/WarpXAlgorithmSelection.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@ struct MediumForEM {
struct EvolveScheme {
enum {
Explicit = 0,
ImplicitPicard = 1
ImplicitPicard = 1,
SemiImplicitPicard = 2
};
};

Expand Down
7 changes: 4 additions & 3 deletions Source/Utils/WarpXAlgorithmSelection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,10 @@
// and corresponding integer for use inside the code

const std::map<std::string, int> evolve_scheme_to_int = {
{"explicit", EvolveScheme::Explicit },
{"implicit_picard", EvolveScheme::ImplicitPicard },
{"default", EvolveScheme::Explicit }
{"explicit", EvolveScheme::Explicit },
{"implicit_picard", EvolveScheme::ImplicitPicard },
{"semi_implicit_picard", EvolveScheme::SemiImplicitPicard },
{"default", EvolveScheme::Explicit }
};

const std::map<std::string, int> grid_to_int = {
Expand Down
2 changes: 1 addition & 1 deletion Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ public:
static short particle_pusher_algo;
//! Integer that corresponds to the type of Maxwell solver (Yee, CKC, PSATD, ECT)
static short electromagnetic_solver_id;
//! Integer that corresponds to the evolve scheme (explicit, implicit_picard)
//! Integer that corresponds to the evolve scheme (explicit, implicit_picard, semi_implicit_picard)
static short evolve_scheme;
//! The maximum number of Picard iterations to do each time step
static int max_picard_iterations;
Expand Down
6 changes: 4 additions & 2 deletions Source/WarpX.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2013,7 +2013,8 @@ void
WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm)
{
const bool aux_is_nodal = ((field_gathering_algo == GatheringAlgo::MomentumConserving) &&
(evolve_scheme != EvolveScheme::ImplicitPicard));
(evolve_scheme != EvolveScheme::ImplicitPicard &&
evolve_scheme != EvolveScheme::SemiImplicitPicard));

#if defined(WARPX_DIM_1D_Z)
const amrex::RealVect dx(WarpX::CellSize(lev)[2]);
Expand Down Expand Up @@ -2085,7 +2086,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
AllocLevelMFs(lev, ba, dm, guard_cells.ng_alloc_EB, guard_cells.ng_alloc_J,
guard_cells.ng_alloc_Rho, guard_cells.ng_alloc_F, guard_cells.ng_alloc_G, aux_is_nodal);

if (evolve_scheme == EvolveScheme::ImplicitPicard) {
if (evolve_scheme == EvolveScheme::ImplicitPicard ||
evolve_scheme == EvolveScheme::SemiImplicitPicard) {
EvolveImplicitPicardInit(lev);
}

Expand Down

0 comments on commit dc7b37a

Please sign in to comment.