[WIP]Add diagnostic file output for the implicit solvers

Docs/source/install/dependencies.rst

@@ -6,7 +6,7 @@ Dependencies
 WarpX depends on the following popular third party software.
 Please see installation instructions below.
 
-- a mature `C++17 <https://en.wikipedia.org/wiki/C%2B%2B17>`__ compiler, e.g., GCC 8.4+, Clang 7, NVCC 11.0, MSVC 19.15 or newer
+- a mature `C++17 <https://en.wikipedia.org/wiki/C%2B%2B17>`__ compiler, e.g., GCC 9.1+, Clang 7, NVCC 11.0, MSVC 19.15 or newer
 - `CMake 3.24.0+ <https://cmake.org>`__
 - `Git 2.18+ <https://git-scm.com>`__
 - `AMReX <https://amrex-codes.github.io>`__: we automatically download and compile a copy of AMReX

Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb

@@ -52,13 +52,15 @@ implicit_evolve.particle_tolerance = 1.0e-12
 #picard.relative_tolerance = 0.0 #1.0e-12
 #picard.absolute_tolerance = 0.0 #1.0e-24
 #picard.require_convergence = false
+#picard.diagnostic_file = "diags/picard_solver.txt"
 
 implicit_evolve.nonlinear_solver = "newton"
 newton.verbose = true
 newton.max_iterations = 20
 newton.relative_tolerance = 1.0e-12
 newton.absolute_tolerance = 0.0
 newton.require_convergence = false
+newton.diagnostic_file = "diags/newton_solver.txt"
 
 gmres.verbose_int = 2
 gmres.max_iterations = 1000

Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp

@@ -76,7 +76,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, a_dt );
+    m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt, a_step );
 
     // Update WarpX owned Efield_fp to t_{n+1/2}
     m_WarpX->SetElectricFieldAndApplyBCs( m_E );

Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp

@@ -77,7 +77,7 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real a_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, a_dt );
+    m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt, a_step );
 
     // Update WarpX owned Efield_fp and Bfield_fp to t_{n+1/2}
     UpdateWarpXFields( m_E, half_time, a_dt );

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp

@@ -104,7 +104,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, a_dt );
+    m_nlsolver->Solve( m_E, m_Eold, theta_time, a_dt, a_step );
 
     // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta}
     UpdateWarpXFields( m_E, theta_time, a_dt );

Source/NonlinearSolvers/NewtonSolver.H

@@ -16,6 +16,8 @@
 #include <AMReX_ParmParse.H>
 
 #include <vector>
+#include <istream>
+#include <filesystem>
 
 /**
  * \brief Newton method to solve nonlinear equation of form:
@@ -45,7 +47,8 @@ public:
     void Solve ( Vec&   a_U,
            const Vec&   a_b,
            amrex::Real  a_time,
-           amrex::Real  a_dt ) const override;
+           amrex::Real  a_dt,
+           int          a_step) const override;
 
     void GetSolverParams ( amrex::Real&  a_rtol,
                            amrex::Real&  a_atol,
@@ -201,6 +204,34 @@ void NewtonSolver<Vec,Ops>::Define ( const Vec&  a_U,
 
     this->m_is_defined = true;
 
+    // Create diagnostic file and write header
+    if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) {
+
+        std::filesystem::path diagnostic_path(this->m_diagnostic_file);
+        std::filesystem::path diagnostic_dir = diagnostic_path.parent_path();
+        if (!diagnostic_dir.empty()) {
+            std::filesystem::create_directories(diagnostic_dir);
+        }
+
+        std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::trunc};
+        int c = 0;
+        diagnostic_file << "#";
+        diagnostic_file << "[" << c++ << "]step()";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]time(s)";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]iters";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]norm_abs";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]norm_rel";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]gmres_iters";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]gmres_last_res";
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
 }
 
 template <class Vec, class Ops>
@@ -212,6 +243,7 @@ void NewtonSolver<Vec,Ops>::ParseParameters ()
     pp_newton.query("relative_tolerance",  m_rtol);
     pp_newton.query("max_iterations",      m_maxits);
     pp_newton.query("require_convergence", m_require_convergence);
+    pp_newton.query("diagnostic_file",     this->m_diagnostic_file);
 
     const amrex::ParmParse pp_gmres("gmres");
     pp_gmres.query("verbose_int",         m_gmres_verbose_int);
@@ -228,7 +260,8 @@ template <class Vec, class Ops>
 void NewtonSolver<Vec,Ops>::Solve ( Vec&         a_U,
                               const Vec&         a_b,
                                     amrex::Real  a_time,
-                                    amrex::Real  a_dt ) const
+                                    amrex::Real  a_dt,
+                                    int          a_step) const
 {
     BL_PROFILE("NewtonSolver::Solve()");
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
@@ -249,6 +282,7 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec&         a_U,
     amrex::Real norm_rel = 0.;
 
     int iter;
+    int linear_solver_iters = 0;
     for (iter = 0; iter < m_maxits;) {
 
         // Compute residual: F(U) = U - b - R(U)
@@ -294,6 +328,7 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec&         a_U,
         // Solve linear system for Newton step [Jac]*dU = F
         m_dU.zero();
         m_linear_solver->solve( m_dU, m_F, m_gmres_rtol, m_gmres_atol );
+        linear_solver_iters += m_linear_solver->getNumIters();
 
         // Update solution
         a_U -= m_dU;
@@ -322,6 +357,26 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec&         a_U,
        }
     }
 
+    if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) {
+        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 << " ";;
+        diagnostic_file << a_time;
+        diagnostic_file << " ";;
+        diagnostic_file << iter;
+        diagnostic_file << " ";;
+        diagnostic_file << norm_abs;
+        diagnostic_file << " ";;
+        diagnostic_file << norm_rel;
+        diagnostic_file << " ";;
+        diagnostic_file << linear_solver_iters;
+        diagnostic_file << " ";;
+        diagnostic_file << m_linear_solver->getResidualNorm();
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
+
 }
 
 template <class Vec, class Ops>

Source/NonlinearSolvers/NonlinearSolver.H

@@ -59,7 +59,8 @@ public:
     virtual void Solve ( Vec&,
                    const Vec&,
                          amrex::Real,
-                         amrex::Real ) const = 0;
+                         amrex::Real,
+                         int) const = 0;
 
     /**
      * \brief Print parameters used by the nonlinear solver.
@@ -81,6 +82,7 @@ protected:
 
     bool m_is_defined = false;
     mutable bool m_verbose = true;
+    std::string m_diagnostic_file;
 
 };
 

Source/NonlinearSolvers/PicardSolver.H

@@ -13,6 +13,8 @@
 #include "Utils/TextMsg.H"
 
 #include <vector>
+#include <istream>
+#include <filesystem>
 
 /**
  * \brief Picard fixed-point iteration method to solve nonlinear
@@ -42,7 +44,8 @@ public:
     void Solve ( Vec&   a_U,
            const Vec&   a_b,
            amrex::Real  a_time,
-           amrex::Real  a_dt ) const override;
+           amrex::Real  a_dt,
+           int          a_step) const override;
 
     void GetSolverParams ( amrex::Real&  a_rtol,
                            amrex::Real&  a_atol,
@@ -114,6 +117,31 @@ void PicardSolver<Vec,Ops>::Define ( const Vec&  a_U,
 
     this->m_is_defined = true;
 
+    // Create diagnostic file and write header
+    if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) {
+
+        std::filesystem::path diagnostic_path(this->m_diagnostic_file);
+        std::filesystem::path diagnostic_dir = diagnostic_path.parent_path();
+        if (!diagnostic_dir.empty()) {
+            std::filesystem::create_directories(diagnostic_dir);
+        }
+
+        std::ofstream diagnostic_file{this->m_diagnostic_file, std::ofstream::out | std::ofstream::trunc};
+        int c = 0;
+        diagnostic_file << "#";
+        diagnostic_file << "[" << c++ << "]step()";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]time(s)";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]iters";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]norm_abs";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]norm_rel";
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
+
 }
 
 template <class Vec, class Ops>
@@ -125,14 +153,15 @@ void PicardSolver<Vec,Ops>::ParseParameters ()
     pp_picard.query("relative_tolerance",  m_rtol);
     pp_picard.query("max_iterations",      m_maxits);
     pp_picard.query("require_convergence", m_require_convergence);
-
+    pp_picard.query("diagnostic_file",     this->m_diagnostic_file);
 }
 
 template <class Vec, class Ops>
 void PicardSolver<Vec,Ops>::Solve ( Vec&         a_U,
                               const Vec&         a_b,
                                     amrex::Real  a_time,
-                                    amrex::Real  a_dt ) const
+                                    amrex::Real  a_dt,
+                                    int          a_step) const
 {
     BL_PROFILE("PicardSolver::Solve()");
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
@@ -212,6 +241,22 @@ void PicardSolver<Vec,Ops>::Solve ( Vec&         a_U,
        }
     }
 
+    if (!this->m_diagnostic_file.empty() && amrex::ParallelDescriptor::IOProcessor()) {
+        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 << " ";
+        diagnostic_file << a_time;
+        diagnostic_file << " ";
+        diagnostic_file << iter;
+        diagnostic_file << " ";
+        diagnostic_file << norm_abs;
+        diagnostic_file << " ";
+        diagnostic_file << norm_rel;
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
+
 }
 
 #endif

Tools/machines/desktop/spack-ubuntu-cuda.yaml

@@ -121,19 +121,6 @@ spack:
       modules: []
       environment: {}
       extra_rpaths: []
-  - compiler:
-      spec: [email protected]
-      paths:
-        cc: /usr/bin/gcc
-        cxx: /usr/bin/g++
-        f77: /usr/bin/gfortran
-        fc: /usr/bin/gfortran
-      flags: {}
-      operating_system: debian10
-      target: x86_64
-      modules: []
-      environment: {}
-      extra_rpaths: []
 
   # binary caches
   mirrors:

Tools/machines/desktop/spack-ubuntu-openmp.yaml

@@ -111,19 +111,6 @@ spack:
       modules: []
       environment: {}
       extra_rpaths: []
-  - compiler:
-      spec: [email protected]
-      paths:
-        cc: /usr/bin/gcc
-        cxx: /usr/bin/g++
-        f77: /usr/bin/gfortran
-        fc: /usr/bin/gfortran
-      flags: {}
-      operating_system: debian10
-      target: x86_64
-      modules: []
-      environment: {}
-      extra_rpaths: []
 
   # binary caches
   mirrors:

Tools/machines/desktop/spack-ubuntu-rocm.yaml

@@ -115,19 +115,6 @@ spack:
       modules: []
       environment: {}
       extra_rpaths: []
-  - compiler:
-      spec: [email protected]
-      paths:
-        cc: /usr/bin/gcc
-        cxx: /usr/bin/g++
-        f77: /usr/bin/gfortran
-        fc: /usr/bin/gfortran
-      flags: {}
-      operating_system: debian10
-      target: x86_64
-      modules: []
-      environment: {}
-      extra_rpaths: []
 
   # binary caches
   mirrors: