Fix the time for the insulator boundary field evaluation (#5488)

In the application of the insulator boundary conditions, the field in
the insulator is specified using an expression that can be a function of
time. This PR updates the calls to that routine to pass in the correct
time, for example the time at the half step or the end of the step. This
requires passing the time down the call chain to the routines that apply
the insulator boundaries which makes this a fairly sizable set of
changes.

Note that in some places the variable names were changed to be more
precise about which time is passed in. For example, using `start_time`
whenever the time at the start of the step is expected.

Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json

@@ -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

Source/BoundaryConditions/WarpXFieldBoundaries.cpp

@@ -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;
 
@@ -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 {
@@ -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);
             }
         }
@@ -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;
 
@@ -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);
         }
     }
 

Source/Evolve/WarpXEvolve.cpp

@@ -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();
@@ -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();
@@ -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
@@ -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) {
@@ -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?
@@ -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) {
@@ -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) {
@@ -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) {

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp

@@ -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
@@ -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);
 }

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H

@@ -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;
 

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp

@@ -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 )
 {
@@ -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);

Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H

@@ -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;
 
@@ -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