Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix the time for the insulator boundary field evaluation #5488

Merged
merged 12 commits into from
Dec 7, 2024
Merged

Fix the time for the insulator boundary field evaluation #5488

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
e7db4ac
Pass time into ApplyEfieldBoundary and ApplyBfieldBoundary
dpgrote Nov 26, 2024
7214cd8
Fix use of m_dt for spectral solver
dpgrote Nov 26, 2024
a129cc3
Fix unused variable
dpgrote Nov 26, 2024
8761d94
Fix initialization order in SpectralSolverRZ
dpgrote Dec 2, 2024
8ac2710
Uodate test_2d_pec_field_insulator.json
dpgrote Dec 2, 2024
08c1b13
Fix const
dpgrote Dec 2, 2024
d9afee2
Fix routine parameter name
dpgrote Dec 2, 2024
2c45815
Fix times in ThetaImplicitEM.H
dpgrote Dec 6, 2024
793059b
Multiple fixes
dpgrote Dec 6, 2024
dbf3af8
Merge branch 'development' into fix_insulator_BC_time
dpgrote Dec 6, 2024
518bf97
Fix argument name of PushPSATD
dpgrote Dec 7, 2024
baf5d03
Fix name of argument for ThetaImplicitEM::FinishFieldUpdate
dpgrote Dec 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json
View file
Open in desktop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you confirm that these changes in the checksums are expected? In particular, do you know which call(s) specifically were using the wrong time value and were thus fixed in this PR, leading to these changes in the checksums?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes were expected. Previously, in each call to EvolveB the time used when setting the boundary conditions was the start of the time step. With this PR, the time is the correct time, at the half time step and end time step for the split calls to EvolveB in the explicit solver.

Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
{
"lev=0": {
"Bx": 0.0,
"By": 0.33065279639752304,
"By": 0.3605690508580849,
"Bz": 0.0,
"Ex": 31873416.396984838,
"Ex": 37101418.32484252,
"Ey": 0.0,
"Ez": 99285542.27022335,
"Ez": 108228802.9434361,
"jx": 0.0,
"jy": 0.0,
"jz": 0.0
Expand Down
18 changes: 8 additions & 10 deletions Source/BoundaryConditions/WarpXFieldBoundaries.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ namespace

}

void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type, amrex::Real time)
{
using ablastr::fields::Direction;

Expand Down Expand Up @@ -94,23 +94,22 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
}

if (::isAnyBoundary<FieldBoundaryType::PECInsulator>(field_boundary_lo, field_boundary_hi)) {
amrex::Real const tnew = gett_new(lev);
if (patch_type == PatchType::fine) {
pec_insulator_boundary->ApplyPEC_InsulatortoEfield(
{m_fields.get(FieldType::Efield_fp,Direction{0},lev),
m_fields.get(FieldType::Efield_fp,Direction{1},lev),
m_fields.get(FieldType::Efield_fp,Direction{2},lev)},
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew);
lev, patch_type, ref_ratio, time);
if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
// apply pec on split E-fields in PML region
const bool split_pml_field = true;
pec_insulator_boundary->ApplyPEC_InsulatortoEfield(
m_fields.get_alldirs(FieldType::pml_E_fp, lev),
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew,
lev, patch_type, ref_ratio, time,
split_pml_field);
}
} else {
Expand All @@ -120,15 +119,15 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
m_fields.get(FieldType::Efield_cp,Direction{2},lev)},
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew);
lev, patch_type, ref_ratio, time);
if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
// apply pec on split E-fields in PML region
const bool split_pml_field = true;
pec_insulator_boundary->ApplyPEC_InsulatortoEfield(
m_fields.get_alldirs(FieldType::pml_E_cp, lev),
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew,
lev, patch_type, ref_ratio, time,
split_pml_field);
}
}
Expand All @@ -147,7 +146,7 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
#endif
}

void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_dt_type)
void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_dt_type, amrex::Real time)
{
using ablastr::fields::Direction;

Expand All @@ -172,23 +171,22 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d
}

if (::isAnyBoundary<FieldBoundaryType::PECInsulator>(field_boundary_lo, field_boundary_hi)) {
amrex::Real const tnew = gett_new(lev);
if (patch_type == PatchType::fine) {
pec_insulator_boundary->ApplyPEC_InsulatortoBfield(
{m_fields.get(FieldType::Bfield_fp,Direction{0},lev),
m_fields.get(FieldType::Bfield_fp,Direction{1},lev),
m_fields.get(FieldType::Bfield_fp,Direction{2},lev)},
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew);
lev, patch_type, ref_ratio, time);
} else {
pec_insulator_boundary->ApplyPEC_InsulatortoBfield(
{m_fields.get(FieldType::Bfield_cp,Direction{0},lev),
m_fields.get(FieldType::Bfield_cp,Direction{1},lev),
m_fields.get(FieldType::Bfield_cp,Direction{2},lev)},
field_boundary_lo, field_boundary_hi,
get_ng_fieldgather(), Geom(lev),
lev, patch_type, ref_ratio, tnew);
lev, patch_type, ref_ratio, time);
}
}

Expand Down
46 changes: 23 additions & 23 deletions Source/Evolve/WarpXEvolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -390,7 +390,7 @@ WarpX::OneStep_nosub (Real cur_time)
WarpX::Hybrid_QED_Push(dt);
FillBoundaryE(guard_cells.ng_alloc_EB);
}
PushPSATD();
PushPSATD(cur_time);

if (do_pml) {
DampPML();
Expand Down Expand Up @@ -418,23 +418,23 @@ WarpX::OneStep_nosub (Real cur_time)
FillBoundaryF(guard_cells.ng_FieldSolverF);
FillBoundaryG(guard_cells.ng_FieldSolverG);

EvolveB(0.5_rt * dt[0], DtType::FirstHalf); // We now have B^{n+1/2}
EvolveB(0.5_rt * dt[0], DtType::FirstHalf, cur_time); // We now have B^{n+1/2}
FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

if (WarpX::em_solver_medium == MediumForEM::Vacuum) {
// vacuum medium
EvolveE(dt[0]); // We now have E^{n+1}
EvolveE(dt[0], cur_time); // We now have E^{n+1}
} else if (WarpX::em_solver_medium == MediumForEM::Macroscopic) {
// macroscopic medium
MacroscopicEvolveE(dt[0]); // We now have E^{n+1}
MacroscopicEvolveE(dt[0], cur_time); // We now have E^{n+1}
} else {
WARPX_ABORT_WITH_MESSAGE("Medium for EM is unknown");
}
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

EvolveF(0.5_rt * dt[0], DtType::SecondHalf);
EvolveG(0.5_rt * dt[0], DtType::SecondHalf);
EvolveB(0.5_rt * dt[0], DtType::SecondHalf); // We now have B^{n+1}
EvolveB(0.5_rt * dt[0], DtType::SecondHalf, cur_time + 0.5_rt * dt[0]); // We now have B^{n+1}

if (do_pml) {
DampPML();
Expand Down Expand Up @@ -808,10 +808,10 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
{
pml[lev]->PushPSATD(m_fields, lev);
}
ApplyEfieldBoundary(lev, PatchType::fine);
if (lev > 0) { ApplyEfieldBoundary(lev, PatchType::coarse); }
ApplyBfieldBoundary(lev, PatchType::fine, DtType::FirstHalf);
if (lev > 0) { ApplyBfieldBoundary(lev, PatchType::coarse, DtType::FirstHalf); }
ApplyEfieldBoundary(lev, PatchType::fine, cur_time + dt[0]);
if (lev > 0) { ApplyEfieldBoundary(lev, PatchType::coarse, cur_time + dt[0]); }
ApplyBfieldBoundary(lev, PatchType::fine, DtType::FirstHalf, cur_time + dt[0]);
if (lev > 0) { ApplyBfieldBoundary(lev, PatchType::coarse, DtType::FirstHalf, cur_time + dt[0]); }
}

// Damp fields in PML before exchanging guard cells
Expand Down Expand Up @@ -886,17 +886,17 @@ WarpX::OneStep_sub1 (Real cur_time)
m_fields.get_mr_levels(FieldType::rho_cp, finest_level),
fine_lev, PatchType::fine, 0, 2*ncomps);

EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf);
EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf, cur_time);
EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf);
FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points);
FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_alloc_F,
WarpX::sync_nodal_points);

EvolveE(fine_lev, PatchType::fine, dt[fine_lev]);
EvolveE(fine_lev, PatchType::fine, dt[fine_lev], cur_time);
FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldGather);

EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf);
EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf, cur_time + 0.5_rt * dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf);

if (do_pml) {
Expand All @@ -922,22 +922,22 @@ WarpX::OneStep_sub1 (Real cur_time)
m_fields.get_mr_levels(FieldType::rho_buf, finest_level),
coarse_lev, 0, ncomps);

EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf);
EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf, cur_time);
EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf);
FillBoundaryB(fine_lev, PatchType::coarse, guard_cells.ng_FieldGather);
FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolverF);

EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]);
EvolveE(fine_lev, PatchType::coarse, dt[fine_lev], cur_time);
FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ng_FieldGather);

EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::FirstHalf);
EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::FirstHalf, cur_time);
EvolveF(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::FirstHalf);
FillBoundaryB(coarse_lev, PatchType::fine, guard_cells.ng_FieldGather,
WarpX::sync_nodal_points);
FillBoundaryF(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolverF,
WarpX::sync_nodal_points);

EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev]);
EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], cur_time);
FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ng_FieldGather);

// TODO Remove call to FillBoundaryAux before UpdateAuxilaryData?
Expand All @@ -961,16 +961,16 @@ WarpX::OneStep_sub1 (Real cur_time)
m_fields.get_mr_levels(FieldType::rho_cp, finest_level),
fine_lev, PatchType::fine, 0, ncomps);

EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf);
EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf, cur_time + dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf);
FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver);
FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_FieldSolverF);

EvolveE(fine_lev, PatchType::fine, dt[fine_lev]);
EvolveE(fine_lev, PatchType::fine, dt[fine_lev], cur_time + dt[fine_lev]);
FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points);

EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf);
EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf, cur_time + 1.5_rt*dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::SecondHalf);

if (do_pml) {
Expand All @@ -997,11 +997,11 @@ WarpX::OneStep_sub1 (Real cur_time)
m_fields.get_mr_levels(FieldType::rho_buf, finest_level),
coarse_lev, ncomps, ncomps);

EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]);
EvolveE(fine_lev, PatchType::coarse, dt[fine_lev], cur_time + 0.5_rt * dt[fine_lev]);
FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points);

EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf);
EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf, cur_time + 0.5_rt * dt[fine_lev]);
EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf);

if (do_pml) {
Expand All @@ -1016,11 +1016,11 @@ WarpX::OneStep_sub1 (Real cur_time)
FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolverF,
WarpX::sync_nodal_points);

EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev]);
EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], cur_time + 0.5_rt*dt[coarse_lev]);
FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolver,
WarpX::sync_nodal_points);

EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::SecondHalf);
EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::SecondHalf, cur_time + 0.5_rt*dt[coarse_lev]);
EvolveF(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], DtType::SecondHalf);

if (do_pml) {
Expand Down
6 changes: 4 additions & 2 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -348,7 +348,8 @@ void HybridPICModel::HybridPICSolveE (
Efield, current_fp_plasma, Jfield, Bfield, rhofield,
*electron_pressure_fp, edge_lengths, lev, this, solve_for_Faraday
);
warpx.ApplyEfieldBoundary(lev, patch_type);
amrex::Real const time = warpx.gett_old(0) + warpx.getdt(0);
warpx.ApplyEfieldBoundary(lev, patch_type, time);
}

void HybridPICModel::CalculateElectronPressure() const
Expand Down Expand Up @@ -556,6 +557,7 @@ void HybridPICModel::FieldPush (
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);
amrex::Real const t_old = warpx.gett_old(0);
warpx.EvolveB(dt, dt_type, t_old);
warpx.FillBoundaryB(ng, nodal_sync);
}
4 changes: 2 additions & 2 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,13 +57,13 @@ public:

void PrintParameters () const override;

void OneStep ( amrex::Real a_time,
void OneStep ( amrex::Real start_time,
amrex::Real a_dt,
int a_step ) override;

void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real half_time,
int a_nl_iter,
bool a_from_jacobian ) override;

Expand Down
19 changes: 10 additions & 9 deletions Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ void SemiImplicitEM::PrintParameters () const
amrex::Print() << "-----------------------------------------------------------\n\n";
}

void SemiImplicitEM::OneStep ( amrex::Real a_time,
void SemiImplicitEM::OneStep ( amrex::Real start_time,
amrex::Real a_dt,
int a_step )
{
Expand All @@ -71,46 +71,47 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time,
m_Eold.Copy( FieldType::Efield_fp );

// Advance WarpX owned Bfield_fp from t_{n} to t_{n+1/2}
m_WarpX->EvolveB(0.5_rt*m_dt, DtType::FirstHalf);
m_WarpX->EvolveB(0.5_rt*m_dt, DtType::FirstHalf, start_time);
m_WarpX->FillBoundaryB(m_WarpX->getngEB(), true);

const amrex::Real half_time = a_time + 0.5_rt*m_dt;
const amrex::Real half_time = start_time + 0.5_rt*m_dt;

// Solve nonlinear system for Eg at t_{n+1/2}
// Particles will be advanced to t_{n+1/2}
m_E.Copy(m_Eold); // initial guess for Eg^{n+1/2}
m_nlsolver->Solve( m_E, m_Eold, half_time, 0.5_rt*m_dt );

// Update WarpX owned Efield_fp to t_{n+1/2}
m_WarpX->SetElectricFieldAndApplyBCs( m_E );
m_WarpX->SetElectricFieldAndApplyBCs( m_E, half_time );

// Advance particles from time n+1/2 to time n+1
m_WarpX->FinishImplicitParticleUpdate();

// Advance Eg from time n+1/2 to time n+1
// Eg^{n+1} = 2.0*Eg^{n+1/2} - Eg^n
m_E.linComb( 2._rt, m_E, -1._rt, m_Eold );
m_WarpX->SetElectricFieldAndApplyBCs( m_E );
const amrex::Real new_time = start_time + m_dt;
m_WarpX->SetElectricFieldAndApplyBCs( m_E, new_time );

// Advance WarpX owned Bfield_fp from t_{n+1/2} to t_{n+1}
m_WarpX->EvolveB(0.5_rt*m_dt, DtType::SecondHalf);
m_WarpX->EvolveB(0.5_rt*m_dt, DtType::SecondHalf, half_time);
m_WarpX->FillBoundaryB(m_WarpX->getngEB(), true);

}

void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real half_time,
int a_nl_iter,
bool a_from_jacobian )
{
// Update WarpX-owned Efield_fp using current state of Eg from
// the nonlinear solver at time n+theta
m_WarpX->SetElectricFieldAndApplyBCs( a_E );
m_WarpX->SetElectricFieldAndApplyBCs( a_E, half_time );

// Update particle positions and velocities using the current state
// of Eg and Bg. Deposit current density at time n+1/2
m_WarpX->ImplicitPreRHSOp( a_time, m_dt, a_nl_iter, a_from_jacobian );
m_WarpX->ImplicitPreRHSOp( half_time, m_dt, a_nl_iter, a_from_jacobian );

// RHS = cvac^2*0.5*dt*( curl(Bg^{n+1/2}) - mu0*Jg^{n+1/2} )
m_WarpX->ImplicitComputeRHSE(0.5_rt*m_dt, a_RHS);
Expand Down
6 changes: 3 additions & 3 deletions Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,13 +58,13 @@ public:

void PrintParameters () const override;

void OneStep ( amrex::Real a_time,
void OneStep ( amrex::Real start_time,
amrex::Real a_dt,
int a_step ) override;

void ComputeRHS ( WarpXSolverVec& a_RHS,
const WarpXSolverVec& a_E,
amrex::Real a_time,
amrex::Real half_time,
int a_nl_iter,
bool a_from_jacobian ) override;

Expand Down Expand Up @@ -92,7 +92,7 @@ private:
* \brief Update the E and B fields owned by WarpX
*/
void UpdateWarpXFields ( WarpXSolverVec const& a_E,
amrex::Real a_time );
amrex::Real half_time );

/**
* \brief Nonlinear solver is for the time-centered values of E. After
Expand Down
Loading
Loading