Add information linear and non-linear solver outputs

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

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 

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

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. 

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

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

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

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:
   /**

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:
   /**

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

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

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

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,

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);

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;

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);

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);