Add information linear and non-linear solver outputs (#1076)

Description of the problem The console output of the non-linear solver in the NS solver was showing the the norms of the newton update without distinction between the velocity and pressure contributions. The final residual norm of the linear solver was not outputted, only the tolerance. Description of the solution Add the distinction between the velocity and pressure contributions in the non-linear solver output. Add the final value of the linear residual when verbosity is at verbose. If it is at extra verbose, the solver output the residual for each iteration (already implemented but not documented). How Has This Been Tested? A unit test was added. vector_problem.cc Co-authored-by: Pierre Laurentin <[email protected]> Co-authored-by: Amishga Alphonius <[email protected]>
chaos-polymtl · Apr 6, 2024 · 4632e53 · 4632e53
1 parent 95dbbc5
commit 4632e53
Show file tree

Hide file tree

Showing 26 changed files with 470 additions and 29 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,6 +3,12 @@
 All notable changes to the Lethe project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
+## [Master] - 2024-04-05
+
+### Changed
+
+- MINOR The output at every iteration in the terminal of the non-linear solver now displays the norms of the correction for both the velocity and pressure contributions for the Navier-Stokes solver when `verbosity` in the `non-linear solver` subsection is set to `verbose`. [#1076](https://github.com/lethe-cfd/lethe/pull/1076)
+
 ## [Master] - 2024-04-04
 
 ### Added

diff --git a/applications_tests/lethe-fluid/heat_transfer_heat-flux.output b/applications_tests/lethe-fluid/heat_transfer_heat-flux.output
@@ -7,11 +7,6 @@ Running on 1 MPI rank(s)...
 *****************************
 Steady iteration:        1/1
 *****************************
----------------
-Fluid Dynamics
----------------
-  -Tolerance of iterative solver is : 1e-12
-  -Iterative solver took : 0 steps 
 L2 error velocity : 0
 L2 error temperature : 1.0986e-29
 cells error_velocity error_pressure 

diff --git a/applications_tests/lethe-fluid/heat_transfer_heat-flux.prm b/applications_tests/lethe-fluid/heat_transfer_heat-flux.prm
@@ -114,4 +114,7 @@ subsection linear solver
     set relative residual = 1e-3
     set minimum residual  = 1e-8
   end
+  subsection fluid dynamics
+    set verbosity = quiet
+  end
 end
diff --git a/doc/source/parameters/cfd/linear_solver_control.rst b/doc/source/parameters/cfd/linear_solver_control.rst
@@ -76,6 +76,7 @@ In this subsection, the control options of the linear solvers are specified. The
 		-Tolerance of iterative solver is : 1e-13
 		-Iterative solver took : 0 steps 
 
+When set to ``extra verbose``, the residual at each iteration of the linear solver is printed.
 
 * The ``minimum residual`` for the linear solver.
 
@@ -231,4 +232,4 @@ Different parameters for the main components of the two geometric multigrid algo
   The default algorithms build and use ALL the multigrid levels. There are two ways to change the number of levels, either by setting the ``mg min level`` parameter OR the ``mg level min cells`` parameter. For ``lsmg`` the coarsest mesh should cover the whole domain, i.e., no hanging nodes are allowed.
 
 .. tip::
-  If ``mg verbosity`` is set to ``verbose``, the information about the levels (cells and degrees of freedom) and the number of iterations of the coarse grid solver are displayed. If this parameter is set to ``extra verbose``, apart from all the previous information, an additional table with the time it took to set up the different components of the multigrid preconditioners is also displayed. 
+  If ``mg verbosity`` is set to ``verbose``, the information about the levels (cells and degrees of freedom) and the number of iterations of the coarse grid solver are displayed. If this parameter is set to ``extra verbose``, apart from all the previous information, an additional table with the time it took to set up the different components of the multigrid preconditioners is also displayed. 
diff --git a/include/core/inexact_newton_non_linear_solver.h b/include/core/inexact_newton_non_linear_solver.h
@@ -142,12 +142,10 @@ InexactNewtonNonLinearSolver<VectorType>::solve(const bool is_initial_step)
               solver->pcout << "\talpha = " << std::setw(6) << alpha
                             << std::setw(0) << " res = "
                             << std::setprecision(this->params.display_precision)
-                            << std::setw(6) << current_res << std::setw(6)
-                            << "\tL^2(dx) = " << std::setw(6)
-                            << newton_update.l2_norm() << std::setw(6)
-                            << "\tL^infty(dx) = "
-                            << std::setprecision(this->params.display_precision)
-                            << newton_update.linfty_norm() << std::endl;
+                            << std::setw(6) << current_res;
+
+              solver->output_newton_update_norms(
+                this->params.display_precision);
             }
 
           // If it's not the first iteration of alpha check if the residual is

diff --git a/include/core/newton_non_linear_solver.h b/include/core/newton_non_linear_solver.h
@@ -133,12 +133,10 @@ NewtonNonLinearSolver<VectorType>::solve(const bool is_initial_step)
               solver->pcout << "\talpha = " << std::setw(6) << alpha
                             << std::setw(0) << " res = "
                             << std::setprecision(this->params.display_precision)
-                            << std::setw(6) << current_res << std::setw(6)
-                            << "\tL^2(dx) = " << std::setw(6)
-                            << newton_update.l2_norm() << std::setw(6)
-                            << "\tL^infty(dx) = "
-                            << std::setprecision(this->params.display_precision)
-                            << newton_update.linfty_norm() << std::endl;
+                            << std::setw(6) << current_res;
+
+              solver->output_newton_update_norms(
+                this->params.display_precision);
             }
 
           // If it's not the first iteration of alpha check if the residual is

diff --git a/include/core/physics_solver.h b/include/core/physics_solver.h
@@ -105,6 +105,14 @@ class PhysicsSolver
   virtual AffineConstraints<double> &
   get_nonzero_constraints() = 0;
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  virtual void
+  output_newton_update_norms(const unsigned int display_precision) = 0;
+
   /**
    * @brief Default way to evaluate the residual for the nonlinear solver.
    * Some application may use more complex evaluation of the residual and

diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h
@@ -309,6 +309,14 @@ class CahnHilliard : public AuxiliaryPhysics<dim, GlobalVectorType>
     return nonzero_constraints;
   }
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  void
+  output_newton_update_norms(const unsigned int display_precision) override;
+
 
 private:
   /**

diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h
@@ -382,6 +382,20 @@ class HeatTransfer : public AuxiliaryPhysics<dim, GlobalVectorType>
     return nonzero_constraints;
   }
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  void
+  output_newton_update_norms(const unsigned int display_precision) override
+  {
+    this->pcout << std::setprecision(display_precision)
+                << "\t||dT||_L2 = " << std::setw(6) << newton_update.l2_norm()
+                << std::setw(6)
+                << "\t||dT||_Linfty = " << std::setprecision(display_precision)
+                << newton_update.linfty_norm() << std::endl;
+  }
 
 private:
   /**

diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h
@@ -46,6 +46,7 @@
 
 #include <deal.II/dofs/dof_handler.h>
 
+#include <deal.II/fe/component_mask.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/mapping_fe.h>
 
@@ -163,6 +164,14 @@ class NavierStokesBase : public PhysicsSolver<VectorType>
     return nonzero_constraints;
   };
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  virtual void
+  output_newton_update_norms(const unsigned int display_precision) override;
+
   /**
    *  Generic interface routine to allow the CFD solver
    *  to cooperate with the multiphysics modules

diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h
@@ -271,6 +271,21 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   }
 
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  void
+  output_newton_update_norms(const unsigned int display_precision) override
+  {
+    this->pcout << std::setprecision(display_precision)
+                << "\t||dx||_L2 = " << std::setw(6) << newton_update.l2_norm()
+                << std::setw(6)
+                << "\t||dx||_Linfty = " << std::setprecision(display_precision)
+                << newton_update.linfty_norm() << std::endl;
+  }
+
 private:
   /**
    *  @brief Assembles the matrix associated with the solver

diff --git a/include/solvers/vof.h b/include/solvers/vof.h
@@ -415,6 +415,21 @@ class VolumeOfFluid : public AuxiliaryPhysics<dim, GlobalVectorType>
     return &present_curvature_solution;
   }
 
+  /**
+   * @brief Output the L2 and Linfty norms of the correction vector.
+   *
+   * @param[in] display_precision Number of outputted digits.
+   */
+  void
+  output_newton_update_norms(const unsigned int display_precision) override
+  {
+    this->pcout << std::setprecision(display_precision)
+                << "\t||dphi||_L2 = " << std::setw(6) << newton_update.l2_norm()
+                << std::setw(6) << "\t||dphi||_Linfty = "
+                << std::setprecision(display_precision)
+                << newton_update.linfty_norm() << std::endl;
+  }
+
 private:
   /**
    *  @brief Assembles the matrix associated with the solver

diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc
@@ -1146,7 +1146,8 @@ CahnHilliard<dim>::solve_linear_system(const bool initial_step,
         .verbosity != Parameters::Verbosity::quiet)
     {
       this->pcout << "  -Iterative solver took : " << solver_control.last_step()
-                  << " steps " << std::endl;
+                  << " steps to reach a residual norm of "
+                  << solver_control.last_value() << std::endl;
     }
 
   constraints_used.distribute(completely_distributed_solution);
@@ -1321,6 +1322,78 @@ CahnHilliard<dim>::apply_phase_filter()
     }
 }
 
+template <int dim>
+void
+CahnHilliard<dim>::output_newton_update_norms(
+  const unsigned int display_precision)
+{
+  auto mpi_communicator = triangulation->get_communicator();
+
+  FEValuesExtractors::Scalar phase_order(0);
+  FEValuesExtractors::Scalar chemical_potential(1);
+
+  ComponentMask phase_order_mask = fe->component_mask(phase_order);
+  ComponentMask chemical_potential_mask =
+    fe->component_mask(chemical_potential);
+
+  const std::vector<IndexSet> index_set_phase_order =
+    DoFTools::locally_owned_dofs_per_component(dof_handler, phase_order_mask);
+  const std::vector<IndexSet> index_set_chemical_potential =
+    DoFTools::locally_owned_dofs_per_component(dof_handler,
+                                               chemical_potential_mask);
+
+  double local_sum = 0.0;
+  double local_max = std::numeric_limits<double>::lowest();
+
+
+  for (auto j = index_set_phase_order[0].begin();
+       j != index_set_phase_order[0].end();
+       j++)
+    {
+      double dof_newton_update = newton_update[*j];
+
+      local_sum += dof_newton_update * dof_newton_update;
+
+      local_max = std::max(local_max, std::abs(dof_newton_update));
+    }
+
+
+  double global_phase_order_l2_norm =
+    std::sqrt(Utilities::MPI::sum(local_sum, mpi_communicator));
+  double global_phase_order_linfty_norm =
+    Utilities::MPI::max(local_max, mpi_communicator);
+
+  local_sum = 0.0;
+  local_max = std::numeric_limits<double>::lowest();
+
+  for (auto j = index_set_chemical_potential[1].begin();
+       j != index_set_chemical_potential[1].end();
+       j++)
+    {
+      double dof_newton_update = newton_update[*j];
+
+      local_sum += dof_newton_update * dof_newton_update;
+
+      local_max = std::max(local_max, std::abs(dof_newton_update));
+    }
+
+  double global_chemical_potential_l2_norm =
+    std::sqrt(Utilities::MPI::sum(local_sum, mpi_communicator));
+  double global_chemical_potential_linfty_norm =
+    Utilities::MPI::max(local_max, mpi_communicator);
+
+  this->pcout << std::setprecision(display_precision)
+              << "\n\t||dphi||_L2 = " << std::setw(6)
+              << global_phase_order_l2_norm << std::setw(6)
+              << "\t||dphi||_Linfty = " << std::setprecision(display_precision)
+              << global_phase_order_linfty_norm << std::endl;
+  this->pcout << std::setprecision(display_precision)
+              << "\t||deta||_L2 = " << std::setw(6)
+              << global_chemical_potential_l2_norm << std::setw(6)
+              << "\t||deta||_Linfty = " << std::setprecision(display_precision)
+              << global_chemical_potential_linfty_norm << std::endl;
+}
+
 template std::pair<Tensor<1, 2>, Tensor<1, 2>>
 CahnHilliard<2>::calculate_barycenter<GlobalVectorType>(
   const GlobalVectorType &solution,

diff --git a/source/solvers/gd_navier_stokes.cc b/source/solvers/gd_navier_stokes.cc
@@ -1154,7 +1154,9 @@ GDNavierStokesSolver<dim>::solve_system_GMRES(const bool   initial_step,
           .verbosity != Parameters::Verbosity::quiet)
       {
         this->pcout << "  -Iterative solver took : "
-                    << solver_control.last_step() << " steps " << std::endl;
+                    << solver_control.last_step()
+                    << " steps to reach a residual norm of "
+                    << solver_control.last_value() << std::endl;
       }
 
     constraints_used.distribute(this->newton_update);

diff --git a/source/solvers/gls_navier_stokes.cc b/source/solvers/gls_navier_stokes.cc
@@ -1611,7 +1611,8 @@ GLSNavierStokesSolver<dim>::solve_system_GMRES(const bool   initial_step,
               {
                 this->pcout
                   << "  -Iterative solver took : " << solver_control.last_step()
-                  << " steps " << std::endl;
+                  << " steps to reach a residual norm of "
+                  << solver_control.last_value() << std::endl;
               }
           }
           constraints_used.distribute(completely_distributed_solution);
@@ -1723,7 +1724,8 @@ GLSNavierStokesSolver<dim>::solve_system_BiCGStab(
               {
                 this->pcout
                   << "  -Iterative solver took : " << solver_control.last_step()
-                  << " steps " << std::endl;
+                  << " steps to reach a residual norm of "
+                  << solver_control.last_value() << std::endl;
               }
             constraints_used.distribute(completely_distributed_solution);
             this->newton_update = completely_distributed_solution;

diff --git a/source/solvers/heat_transfer.cc b/source/solvers/heat_transfer.cc
@@ -1390,7 +1390,8 @@ HeatTransfer<dim>::solve_linear_system(const bool initial_step,
         .verbosity != Parameters::Verbosity::quiet)
     {
       this->pcout << "  -Iterative solver took : " << solver_control.last_step()
-                  << " steps " << std::endl;
+                  << " steps to reach a residual norm of "
+                  << solver_control.last_value() << std::endl;
     }
 
   constraints_used.distribute(completely_distributed_solution);

diff --git a/source/solvers/mf_navier_stokes.cc b/source/solvers/mf_navier_stokes.cc
@@ -1782,7 +1782,8 @@ MFNavierStokesSolver<dim>::solve_system_GMRES(const bool   initial_step,
         .verbosity != Parameters::Verbosity::quiet)
     {
       this->pcout << "  -Iterative solver took : " << solver_control.last_step()
-                  << " steps " << std::endl;
+                  << " steps to reach a residual norm of "
+                  << solver_control.last_value() << std::endl;
     }
 
   constraints_used.distribute(this->newton_update);