From ccf5bcc2c3241584a74cea8ae2b49697c8abc220 Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Mon, 10 Mar 2025 10:45:29 -0700 Subject: [PATCH 1/3] added diagnostic_interval int for nonlinear solver diagnostic. --- Source/NonlinearSolvers/NewtonSolver.H | 4 +++- Source/NonlinearSolvers/NonlinearSolver.H | 1 + Source/NonlinearSolvers/PicardSolver.H | 4 +++- 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/Source/NonlinearSolvers/NewtonSolver.H b/Source/NonlinearSolvers/NewtonSolver.H index 9b87f2e8103..031a8acd805 100644 --- a/Source/NonlinearSolvers/NewtonSolver.H +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -245,6 +245,7 @@ void NewtonSolver::ParseParameters () pp_newton.query("max_iterations", m_maxits); pp_newton.query("require_convergence", m_require_convergence); pp_newton.query("diagnostic_file", this->m_diagnostic_file); + pp_newton.query("diagnostic_interval", this->m_diagnostic_interval); const amrex::ParmParse pp_gmres("gmres"); pp_gmres.query("verbose_int", m_gmres_verbose_int); @@ -358,7 +359,8 @@ void NewtonSolver::Solve ( Vec& a_U, } } - if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) { + if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor() && + a_step%this->m_diagnostic_interval==0) { std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::app}; diagnostic_file << std::setprecision(14); diagnostic_file << a_step; diff --git a/Source/NonlinearSolvers/NonlinearSolver.H b/Source/NonlinearSolvers/NonlinearSolver.H index d7175d8ea56..74f3fc19642 100644 --- a/Source/NonlinearSolvers/NonlinearSolver.H +++ b/Source/NonlinearSolvers/NonlinearSolver.H @@ -83,6 +83,7 @@ protected: bool m_is_defined = false; mutable bool m_verbose = true; std::string m_diagnostic_file; + int m_diagnostic_interval = 1; }; diff --git a/Source/NonlinearSolvers/PicardSolver.H b/Source/NonlinearSolvers/PicardSolver.H index ec53bf94587..12e359889db 100644 --- a/Source/NonlinearSolvers/PicardSolver.H +++ b/Source/NonlinearSolvers/PicardSolver.H @@ -156,6 +156,7 @@ void PicardSolver::ParseParameters () pp_picard.query("max_iterations", m_maxits); pp_picard.query("require_convergence", m_require_convergence); pp_picard.query("diagnostic_file", this->m_diagnostic_file); + pp_picard.query("diagnostic_interval", this->m_diagnostic_interval); } template @@ -244,7 +245,8 @@ void PicardSolver::Solve ( Vec& a_U, } } - if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) { + if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor() && + a_step%this->m_diagnostic_interval==0) { std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::app}; diagnostic_file << std::setprecision(14); diagnostic_file << a_step; From 649ee0d09e696ce731827f24bed611e1f53f2c27 Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Mon, 10 Mar 2025 14:27:21 -0700 Subject: [PATCH 2/3] 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. --- Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H | 7 +++++++ Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H | 2 +- Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp | 5 +++-- .../FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H | 2 +- .../ImplicitSolvers/StrangImplicitSpectralEM.cpp | 5 +++-- Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H | 5 ----- Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp | 2 +- Source/NonlinearSolvers/CurlCurlMLMGPC.H | 6 +++--- 8 files changed, 19 insertions(+), 15 deletions(-) diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H index 0d2083793ac..517252bff49 100644 --- a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H @@ -96,6 +96,8 @@ public: [[nodiscard]] amrex::Array GetLinOpBCLo () const; [[nodiscard]] amrex::Array 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 */ diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H index 62401d7b48f..d8039ef9623 100644 --- a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H +++ b/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; diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp index aae10b978ee..6bd37e35c50 100644 --- a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp +++ b/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 diff --git a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H b/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H index fcfcb9821a7..336ddcd0c53 100644 --- a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.H +++ b/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; diff --git a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp b/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp index debbff1b35e..13427f69d90 100644 --- a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp +++ b/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 diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H index c9de6b82a00..0196da4a00b 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H +++ b/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 diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp index 54644d357d2..418ed5721fc 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp +++ b/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 ); diff --git a/Source/NonlinearSolvers/CurlCurlMLMGPC.H b/Source/NonlinearSolvers/CurlCurlMLMGPC.H index b3fcc6fe38f..c525e03d4a4 100644 --- a/Source/NonlinearSolvers/CurlCurlMLMGPC.H +++ b/Source/NonlinearSolvers/CurlCurlMLMGPC.H @@ -272,8 +272,8 @@ void CurlCurlMLMGPC::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::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"; From aff2041701778050583cedcf91b6c0150712edea Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Mon, 10 Mar 2025 14:29:58 -0700 Subject: [PATCH 3/3] adjusting step and time for solver diagnostic to be consistent with other reduced diagnostics. --- Source/NonlinearSolvers/NewtonSolver.H | 6 +++--- Source/NonlinearSolvers/PicardSolver.H | 7 +++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/Source/NonlinearSolvers/NewtonSolver.H b/Source/NonlinearSolvers/NewtonSolver.H index 031a8acd805..f71785f99cf 100644 --- a/Source/NonlinearSolvers/NewtonSolver.H +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -360,12 +360,12 @@ void NewtonSolver::Solve ( Vec& a_U, } if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor() && - a_step%this->m_diagnostic_interval==0) { + ((a_step+1)%this->m_diagnostic_interval==0 || a_step==0)) { std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::app}; diagnostic_file << std::setprecision(14); - diagnostic_file << a_step; + diagnostic_file << a_step+1; diagnostic_file << " ";; - diagnostic_file << a_time; + diagnostic_file << a_time + a_dt; diagnostic_file << " ";; diagnostic_file << iter; diagnostic_file << " ";; diff --git a/Source/NonlinearSolvers/PicardSolver.H b/Source/NonlinearSolvers/PicardSolver.H index 12e359889db..fbbd6f46f5f 100644 --- a/Source/NonlinearSolvers/PicardSolver.H +++ b/Source/NonlinearSolvers/PicardSolver.H @@ -170,7 +170,6 @@ void PicardSolver::Solve ( Vec& a_U, WARPX_ALWAYS_ASSERT_WITH_MESSAGE( this->m_is_defined, "PicardSolver::Solve() called on undefined object"); - amrex::ignore_unused(a_dt); using namespace amrex::literals; // @@ -246,12 +245,12 @@ void PicardSolver::Solve ( Vec& a_U, } if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor() && - a_step%this->m_diagnostic_interval==0) { + (a_step%this->m_diagnostic_interval==0 || a_step==0)) { std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::app}; diagnostic_file << std::setprecision(14); - diagnostic_file << a_step; + diagnostic_file << a_step+1; diagnostic_file << " "; - diagnostic_file << a_time; + diagnostic_file << a_time + a_dt; diagnostic_file << " "; diagnostic_file << iter; diagnostic_file << " ";