From d1439aff72130d0e0626ff82f7b8c9e45b56063e Mon Sep 17 00:00:00 2001
From: Justin Angus <angus1@llnl.gov>
Date: Mon, 4 Nov 2024 13:32:01 -0800
Subject: [PATCH] dt passed to m_nlsolver->Solve() and eventually used by the
 PC is now the appropriate fractional time step. Removed theta() function from
 ImplicitSolver.H.

---
 Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H    | 4 +---
 Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp  | 2 +-
 Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H   | 2 --
 Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp | 2 +-
 Source/NonlinearSolvers/CurlCurlMLMGPC.H               | 5 +++--
 5 files changed, 6 insertions(+), 9 deletions(-)

diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
index 94b6733d3c8..f8f0390e17a 100644
--- a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
+++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H
@@ -86,8 +86,6 @@ public:
                               int              a_nl_iter,
                               bool             a_from_jacobian ) = 0;
 
-    [[nodiscard]] virtual amrex::Real theta () const { return 1.0; }
-
     [[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; }
 
     [[nodiscard]] const amrex::Geometry& GetGeometry (int) const;
@@ -113,7 +111,7 @@ protected:
     /**
      * \brief Time step
      */
-    mutable amrex::Real m_dt = DBL_MAX;
+    mutable amrex::Real m_dt = 0.0;
 
     /**
      * \brief Nonlinear solver type and object
diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
index 516de694eda..117c3baecaa 100644
--- a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
+++ b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp
@@ -79,7 +79,7 @@ void SemiImplicitEM::OneStep ( amrex::Real  a_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, m_dt );
+    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 );
diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
index dae942bb22d..7461b77fb51 100644
--- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
+++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H
@@ -77,8 +77,6 @@ public:
                       int              a_nl_iter,
                       bool             a_from_jacobian ) override;
 
-    [[nodiscard]] amrex::Real theta () const override { return m_theta; }
-
 private:
 
     /**
diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
index 40bc6dfbb6d..8ca592517ac 100644
--- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
+++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp
@@ -107,7 +107,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real  a_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, theta_time, m_dt );
+    m_nlsolver->Solve( m_E, m_Eold, theta_time, m_theta*m_dt );
 
     // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
     UpdateWarpXFields( m_E, theta_time );
diff --git a/Source/NonlinearSolvers/CurlCurlMLMGPC.H b/Source/NonlinearSolvers/CurlCurlMLMGPC.H
index 47d7310995c..b3fcc6fe38f 100644
--- a/Source/NonlinearSolvers/CurlCurlMLMGPC.H
+++ b/Source/NonlinearSolvers/CurlCurlMLMGPC.H
@@ -272,7 +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
-    const RT alpha = (m_ops->theta()*this->m_dt*PhysConst::c) * (m_ops->theta()*this->m_dt*PhysConst::c);
+    // (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 beta = RT(1.0);
 
 // currently not implemented in 1D
@@ -282,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)
-                       << ": dt = " << this->m_dt << ", "
+                       << ": theta*dt = " << this->m_dt << ", "
                        << " coefficients: "
                        << "alpha = " << alpha << ", "
                        << "beta = " << beta << "\n";