Skip to content

Commit

Permalink
Complete BfieldEvolveRK (#2)
Browse files Browse the repository at this point in the history 
* continued transition to RK4 in B-field time integration in hybrid-PIC algorithm

* RK-integration over all levels
  • Loading branch information
@aveksler1
aveksler1 authored Dec 1, 2023
1 parent 1185283 commit b6f6d06
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 47 deletions.
9 changes: 9 additions & 0 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,15 @@ public:
amrex::Real dt, DtType a_dt_type,
amrex::IntVect ng, std::optional<bool> nodal_sync);

void BfieldEvolveRK (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
amrex::Real const dt, int lev, DtType dt_type,
amrex::IntVect ng, std::optional<bool> nodal_sync);

/**
* \brief
* Function to calculate the electron pressure at a given timestep type
Expand Down
110 changes: 91 additions & 19 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -543,42 +543,61 @@ void HybridPICModel::BfieldEvolveRK (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
amrex::Real const 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,
ng, nodal_sync
);
}
}

void HybridPICModel::BfieldEvolveRK (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
amrex::Real const dt, int lev, DtType dt_type,
IntVect ng, std::optional<bool> nodal_sync )
{
// make copies of the B-field multifabs at t = n
std::array< MultiFab, 3 > B_old;
B_old[0] = MultiFab(
Bfield[0][0]->boxArray(), Bfield[0][0]->DistributionMap(), 1,
Bfield[0][0]->nGrowVect()
Bfield[lev][0]->boxArray(), Bfield[lev][0]->DistributionMap(), 1,
Bfield[lev][0]->nGrowVect()
);
MultiFab::Copy(B_old[0], *Bfield[0][0], 0, 0, 1, ng);
MultiFab::Copy(B_old[0], *Bfield[lev][0], 0, 0, 1, ng);
B_old[1] = MultiFab(
Bfield[0][1]->boxArray(), Bfield[0][1]->DistributionMap(), 1,
Bfield[0][1]->nGrowVect()
Bfield[lev][1]->boxArray(), Bfield[lev][1]->DistributionMap(), 1,
Bfield[lev][1]->nGrowVect()
);
MultiFab::Copy(B_old[1], *Bfield[0][1], 0, 0, 1, ng);
MultiFab::Copy(B_old[1], *Bfield[lev][1], 0, 0, 1, ng);
B_old[2] = MultiFab(
Bfield[0][2]->boxArray(), Bfield[0][2]->DistributionMap(), 1,
Bfield[0][2]->nGrowVect()
Bfield[lev][2]->boxArray(), Bfield[lev][2]->DistributionMap(), 1,
Bfield[lev][2]->nGrowVect()
);
MultiFab::Copy(B_old[2], *Bfield[0][2], 0, 0, 1, ng);
MultiFab::Copy(B_old[2], *Bfield[lev][2], 0, 0, 1, ng);

// Create multifabs for each direction to store the Runge-Kutta
// intermediate terms. Each multifab has 2 components for the different
// terms that need to be stored.
std::array< amrex::MultiFab, 3 > K;
K[0] = amrex::MultiFab(
Bfield[0][0]->boxArray(), Bfield[0][0]->DistributionMap(), 2,
Bfield[0][0]->nGrowVect()
Bfield[lev][0]->boxArray(), Bfield[lev][0]->DistributionMap(), 2,
Bfield[lev][0]->nGrowVect()
);
K[0].setVal(0.0);
K[1] = amrex::MultiFab(
Bfield[0][1]->boxArray(), Bfield[0][1]->DistributionMap(), 2,
Bfield[0][1]->nGrowVect()
Bfield[lev][1]->boxArray(), Bfield[lev][1]->DistributionMap(), 2,
Bfield[lev][1]->nGrowVect()
);
K[1].setVal(0.0);
K[2] = amrex::MultiFab(
Bfield[0][2]->boxArray(), Bfield[0][2]->DistributionMap(), 2,
Bfield[0][2]->nGrowVect()
Bfield[lev][2]->boxArray(), Bfield[lev][2]->DistributionMap(), 2,
Bfield[lev][2]->nGrowVect()
);
K[2].setVal(0.0);

Expand All @@ -601,7 +620,7 @@ void HybridPICModel::BfieldEvolveRK (
{
// Extract 0.5 * dt * K0 for each direction into index 0 of K.
MultiFab::LinComb(
K[ii], 1._rt, *Bfield[0][ii], 0, -1._rt, B_old[ii], 0, 0, 1, ng
K[ii], 1._rt, *Bfield[lev][ii], 0, -1._rt, B_old[ii], 0, 0, 1, ng
);
}

Expand All @@ -622,12 +641,65 @@ void HybridPICModel::BfieldEvolveRK (
{
// Subtract 0.5 * dt * K0 from the Bfield for each direction, to get
// B_new = B_old + 0.5 * dt * K1.
MultiFab::Subtract(*Bfield[0][ii], K[ii], 0, 0, 1, ng);
MultiFab::Subtract(*Bfield[lev][ii], K[ii], 0, 0, 1, ng);
// Extract 0.5 * dt * K1 for each direction into index 1 of K.
MultiFab::LinComb(
K[ii], 1._rt, *Bfield[0][ii], 0, -1._rt, B_old[ii], 0, 1, 1, ng
K[ii], 1._rt, *Bfield[lev][ii], 0, -1._rt, B_old[ii], 0, 1, 1, ng
);
}

// TODO INSERT STEPS 3 & 4
// Step 3:
// Calculate J = curl x B / mu0
CalculateCurrentAmpere(Bfield, edge_lengths);
// Calculate the E-field form Ohm's law
HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
warpx.FillBoundaryE(ng, nodal_sync);
// Push forward the B-field using Faraday's law.
warpx.EvolveB(dt, dt_type);
warpx.FillBoundaryB(ng, nodal_sync);

// The Bfield is now given by:
// B_new = B_old + 0.5 * dt * K1
// + dt * [-curl x E(B_old + 0.5 * dt * K1)]
// = B_old + 0.5 * dt * K1 + dt * K2
for (int ii = 0; ii < 3; ii++)
{
// Subtract 0.5 * dt * K1 from the Bfield for each direction to get
// B_new = B_old + dt * K2.
MultiFab::Subtract(*Bfield[lev][ii], K[ii], 1, 0, 1, ng);
}

// Step 4:
// Calculate J = curl x B / mu0
CalculateCurrentAmpere(Bfield, edge_lengths);
// Calculate the E-field form Ohm's law
HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
warpx.FillBoundaryE(ng, nodal_sync);
// Push forward the B-field using Faraday's law.
warpx.EvolveB(0.5_rt*dt, dt_type);
warpx.FillBoundaryB(ng, nodal_sync);

// The Bfield is now given by:
// B_new = B_old + dt * K2 + 0.5 * dt * [-curl x E(B_old + dt * K2)]
// = B_old + dt * K2 + 0.5 * dt * K3
for (int ii = 0; ii < 3; ii++)
{
// Subtract B_old from the Bfield for each direction, to get
// B = dt * K2 + 0.5 * dt * K3.
MultiFab::Subtract(*Bfield[lev][ii], B_old[ii], 0, 0, 1, ng);

// Add dt * K2 + 0.5 * dt * K3 to index 0 of K (= 0.5 * dt * K0).
MultiFab::Add(K[ii], *Bfield[lev][ii], 0, 0, 1, ng);

// Add 2 * 0.5 * dt * K1 to index 0 of K.
MultiFab::LinComb(
K[ii], 1.0, K[ii], 0, 2.0, K[ii], 1, 0, 1, ng
);

// Overwrite the Bfield with the Runge-Kutta sum:
// B_new = B_old + 1/3 * dt * (0.5 * K0 + K1 + K2 + 0.5 * K3).
MultiFab::LinComb(
*Bfield[lev][ii], 1.0, B_old[ii], 0, 1.0/3.0, K[ii], 0, 0, 1, ng
);
}
}
56 changes: 28 additions & 28 deletions Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,21 +94,21 @@ void WarpX::HybridPICEvolveFields ()
// momentum equation
for (int sub_step = 0; sub_step < sub_steps; sub_step++)
{
m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
m_hybrid_pic_model->HybridPICSolveE(
Efield_fp, current_fp_temp, Bfield_fp, rho_fp_temp, m_edge_lengths,
true
);
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
EvolveB(0.5_rt / sub_steps * dt[0], DtType::FirstHalf);
FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

// m_hybrid_pic_model->BfieldEvolveRK(
// Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp,
// m_edge_lengths, 0.5_rt/sub_steps*dt[0],
// DtType::FirstHalf, guard_cells.ng_FieldSolver,
// WarpX::sync_nodal_points
// m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
// m_hybrid_pic_model->HybridPICSolveE(
// Efield_fp, current_fp_temp, Bfield_fp, rho_fp_temp, m_edge_lengths,
// true
// );
// FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
// EvolveB(0.5_rt / sub_steps * dt[0], DtType::FirstHalf);
// FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

m_hybrid_pic_model->BfieldEvolveRK(
Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp,
m_edge_lengths, 0.5_rt/sub_steps*dt[0],
DtType::FirstHalf, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points
);
}

// Average rho^{n} and rho^{n+1} to get rho^{n+1/2} in rho_fp_temp
Expand All @@ -129,21 +129,21 @@ void WarpX::HybridPICEvolveFields ()
// Now push the B field from t=n+1/2 to t=n+1 using the n+1/2 quantities
for (int sub_step = 0; sub_step < sub_steps; sub_step++)
{
m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
m_hybrid_pic_model->HybridPICSolveE(
Efield_fp, current_fp, Bfield_fp, rho_fp_temp, m_edge_lengths,
true
);
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
EvolveB(0.5_rt / sub_steps * dt[0], DtType::SecondHalf);
FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

// m_hybrid_pic_model->BfieldEvolveRK(
// Bfield_fp, Efield_fp, current_fp, rho_fp_temp,
// m_edge_lengths, 0.5_rt/sub_steps*dt[0],
// DtType::SecondHalf, guard_cells.ng_FieldSolver,
// WarpX::sync_nodal_points
// m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
// m_hybrid_pic_model->HybridPICSolveE(
// Efield_fp, current_fp, Bfield_fp, rho_fp_temp, m_edge_lengths,
// true
// );
// FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
// EvolveB(0.5_rt / sub_steps * dt[0], DtType::SecondHalf);
// FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

m_hybrid_pic_model->BfieldEvolveRK(
Bfield_fp, Efield_fp, current_fp, rho_fp_temp,
m_edge_lengths, 0.5_rt/sub_steps*dt[0],
DtType::SecondHalf, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points
);
}

// Extrapolate the ion current density to t=n+1 using
Expand Down

0 comments on commit b6f6d06

Please sign in to comment.