moved theta to implicit solve base class. time and dt passed to nonlinear solver is time at beginning of time step and the full time step.

Author: JustinRayAngus

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

@@ -96,6 +96,8 @@ public:
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCHi () const;
 
+    [[nodiscard]] amrex::Real  GetTheta () const { return m_theta; }
+
 protected:
 
     /**
@@ -115,6 +117,11 @@ protected:
      */
     mutable amrex::Real m_dt = 0.0;
 
+    /**
+     * \brief Time-biasing parameter for fields used on RHS to advance system.
+     */
+    amrex::Real m_theta = 0.5;
+
     /**
      * \brief Nonlinear solver type and object
      */

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H

@@ -63,7 +63,7 @@ public:
 
     void ComputeRHS ( WarpXSolverVec&  a_RHS,
                 const WarpXSolverVec&  a_E,
-                      amrex::Real      half_time,
+                      amrex::Real      start_time,
                       int              a_nl_iter,
                       bool             a_from_jacobian ) override;
 

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp

@@ -80,7 +80,7 @@ void SemiImplicitEM::OneStep ( amrex::Real  start_time,
     // 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, a_step );
+    m_nlsolver->Solve( m_E, m_Eold, start_time, m_dt, a_step );
 
     // Update WarpX owned Efield_fp to t_{n+1/2}
     m_WarpX->SetElectricFieldAndApplyBCs( m_E, half_time );
@@ -103,12 +103,13 @@ void SemiImplicitEM::OneStep ( amrex::Real  start_time,
 
 void SemiImplicitEM::ComputeRHS ( WarpXSolverVec&  a_RHS,
                             const WarpXSolverVec&  a_E,
-                                  amrex::Real      half_time,
+                                  amrex::Real      start_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
+    const amrex::Real half_time = start_time + 0.5_rt*m_dt;
     m_WarpX->SetElectricFieldAndApplyBCs( a_E, half_time );
 
     // Update particle positions and velocities using the current state

Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H

@@ -64,7 +64,7 @@ public:
 
     void ComputeRHS ( WarpXSolverVec&  a_RHS,
                 const WarpXSolverVec&  a_E,
-                      amrex::Real      half_time,
+                      amrex::Real      start_time,
                       int              a_nl_iter,
                       bool             a_from_jacobian ) override;
 

Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp

@@ -81,7 +81,7 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real start_time,
 
     // Solve nonlinear system for E at t_{n+1/2}
     // Particles will be advanced to t_{n+1/2}
-    m_nlsolver->Solve( m_E, m_Eold, half_time, 0.5_rt*m_dt, a_step );
+    m_nlsolver->Solve( m_E, m_Eold, start_time, m_dt, a_step );
 
     // Update WarpX owned Efield_fp and Bfield_fp to t_{n+1/2}
     UpdateWarpXFields( m_E, half_time );
@@ -101,12 +101,13 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real start_time,
 
 void StrangImplicitSpectralEM::ComputeRHS ( WarpXSolverVec& a_RHS,
                                             WarpXSolverVec const & a_E,
-                                            amrex::Real half_time,
+                                            amrex::Real start_time,
                                             int a_nl_iter,
                                             bool a_from_jacobian )
 {
     // Update WarpX-owned Efield_fp and Bfield_fp using current state of
     // E from the nonlinear solver at time n+1/2
+    const amrex::Real half_time = start_time + 0.5_rt*m_dt;
     UpdateWarpXFields( a_E, half_time );
 
     // Self consistently update particle positions and velocities using the

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H

@@ -86,11 +86,6 @@ public:
 
 private:
 
-    /**
-     * \brief Time-biasing parameter for fields used on RHS to advance system
-     */
-    amrex::Real m_theta = 0.5;
-
     /**
      * \brief Solver vectors to be used in the nonlinear solver to solve for the
      * electric field E. The main logic for determining which variables should be

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp

@@ -106,7 +106,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real  start_time,
     // Solve nonlinear system for Eg at t_{n+theta}
     // Particles will be advanced to t_{n+1/2}
     m_E.Copy(m_Eold); // initial guess for Eg^{n+theta}
-    m_nlsolver->Solve( m_E, m_Eold, start_time, m_theta*m_dt, a_step );
+    m_nlsolver->Solve( m_E, m_Eold, start_time, m_dt, a_step );
 
     // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
     UpdateWarpXFields( m_E, start_time );

Source/NonlinearSolvers/CurlCurlMLMGPC.H

@@ -272,8 +272,8 @@ void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)
     amrex::ignore_unused(a_U);
 
     // set the coefficients alpha and beta for curl-curl op
-    // (m_dt here is actually theta<=0.5 times simulation dt)
-    const RT alpha = (this->m_dt*PhysConst::c) * (this->m_dt*PhysConst::c);
+    const RT thetaDt = m_ops->GetTheta()*this->m_dt;
+    const RT alpha = (thetaDt*PhysConst::c) * (thetaDt*PhysConst::c);
     const RT beta = RT(1.0);
 
 // currently not implemented in 1D
@@ -283,7 +283,7 @@ void CurlCurlMLMGPC<T,Ops>::Update (const T& a_U)
 
     if (m_verbose) {
         amrex::Print() << "Updating " << amrex::getEnumNameString(PreconditionerType::pc_curl_curl_mlmg)
-                       << ": theta*dt = " << this->m_dt << ", "
+                       << ": theta*dt = " << thetaDt << ", "
                        << " coefficients: "
                        << "alpha = " << alpha << ", "
                        << "beta = " << beta << "\n";