From 4632e531764aa19ddd111f8c3d7d859fdd90f3f6 Mon Sep 17 00:00:00 2001 From: hepap <47506601+hepap@users.noreply.github.com> Date: Sat, 6 Apr 2024 03:17:34 -0400 Subject: [PATCH] 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 <81017102+PierreLaurentinCS@users.noreply.github.com> Co-authored-by: Amishga Alphonius <107414376+AmishgaAlphonius@users.noreply.github.com> --- CHANGELOG.md | 6 + .../heat_transfer_heat-flux.output | 5 - .../lethe-fluid/heat_transfer_heat-flux.prm | 3 + .../parameters/cfd/linear_solver_control.rst | 3 +- .../core/inexact_newton_non_linear_solver.h | 10 +- include/core/newton_non_linear_solver.h | 10 +- include/core/physics_solver.h | 8 + include/solvers/cahn_hilliard.h | 8 + include/solvers/heat_transfer.h | 14 ++ include/solvers/navier_stokes_base.h | 9 + include/solvers/tracer.h | 15 ++ include/solvers/vof.h | 15 ++ source/solvers/cahn_hilliard.cc | 75 +++++++- source/solvers/gd_navier_stokes.cc | 4 +- source/solvers/gls_navier_stokes.cc | 6 +- source/solvers/heat_transfer.cc | 3 +- source/solvers/mf_navier_stokes.cc | 3 +- source/solvers/navier_stokes_base.cc | 90 ++++++++++ source/solvers/tracer.cc | 3 +- source/solvers/vof.cc | 3 +- ...inexact_newton_non_linear_solver_01.output | 7 +- ...inexact_newton_non_linear_solver_02.output | 11 +- tests/core/newton_non_linear_solver_01.output | 6 +- tests/core/non_linear_test_system_01.h | 9 + tests/core/vector_problem.cc | 168 ++++++++++++++++++ tests/core/vector_problem.output | 5 + 26 files changed, 470 insertions(+), 29 deletions(-) create mode 100644 tests/core/vector_problem.cc create mode 100644 tests/core/vector_problem.output diff --git a/CHANGELOG.md b/CHANGELOG.md index c5d2b7d9f2..b687616228 100644 --- 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 index 7ea129a173..74d5a8c4af 100755 --- 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 index 2b24f0429f..bdb5a0a0e3 100755 --- 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 index 524521b288..b6e2537468 100644 --- 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. \ No newline at end of file + 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 index dd8ea70295..43123d83b4 100644 --- a/include/core/inexact_newton_non_linear_solver.h +++ b/include/core/inexact_newton_non_linear_solver.h @@ -142,12 +142,10 @@ InexactNewtonNonLinearSolver::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 index 2f1840b0e7..d77f90298c 100644 --- a/include/core/newton_non_linear_solver.h +++ b/include/core/newton_non_linear_solver.h @@ -133,12 +133,10 @@ NewtonNonLinearSolver::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 index 3d40c2eb42..3fa340a7ae 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -105,6 +105,14 @@ class PhysicsSolver virtual AffineConstraints & 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 index 8882f66100..2e1da183ad 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -309,6 +309,14 @@ class CahnHilliard : public AuxiliaryPhysics 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 index 7fa5b2ffce..021a00ce39 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -382,6 +382,20 @@ class HeatTransfer : public AuxiliaryPhysics 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 index 02aa4c3cc3..eae196c511 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -46,6 +46,7 @@ #include +#include #include #include @@ -163,6 +164,14 @@ class NavierStokesBase : public PhysicsSolver 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 index e15d64eceb..d7818d2cc6 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -271,6 +271,21 @@ class Tracer : public AuxiliaryPhysics } + /** + * @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 index a455041898..4bb0b3e646 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -415,6 +415,21 @@ class VolumeOfFluid : public AuxiliaryPhysics 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 index d7cbf7fd9a..d7e492bfdb 100644 --- a/source/solvers/cahn_hilliard.cc +++ b/source/solvers/cahn_hilliard.cc @@ -1146,7 +1146,8 @@ CahnHilliard::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::apply_phase_filter() } } +template +void +CahnHilliard::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 index_set_phase_order = + DoFTools::locally_owned_dofs_per_component(dof_handler, phase_order_mask); + const std::vector 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::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::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>> CahnHilliard<2>::calculate_barycenter( const GlobalVectorType &solution, diff --git a/source/solvers/gd_navier_stokes.cc b/source/solvers/gd_navier_stokes.cc index 248fb786fc..a6f6b7ac23 100644 --- a/source/solvers/gd_navier_stokes.cc +++ b/source/solvers/gd_navier_stokes.cc @@ -1154,7 +1154,9 @@ GDNavierStokesSolver::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 index 9602e73fe8..461e07aebb 100644 --- a/source/solvers/gls_navier_stokes.cc +++ b/source/solvers/gls_navier_stokes.cc @@ -1611,7 +1611,8 @@ GLSNavierStokesSolver::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::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 index f9aa19347c..1d4f6a58e1 100644 --- a/source/solvers/heat_transfer.cc +++ b/source/solvers/heat_transfer.cc @@ -1390,7 +1390,8 @@ HeatTransfer::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 index eacd1a8573..25691f79b7 100644 --- a/source/solvers/mf_navier_stokes.cc +++ b/source/solvers/mf_navier_stokes.cc @@ -1782,7 +1782,8 @@ MFNavierStokesSolver::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); diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 352743e172..2c41a40b90 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2483,6 +2483,96 @@ NavierStokesBase:: } } +template +void +NavierStokesBase::output_newton_update_norms( + const unsigned int display_precision) +{ + if constexpr (std::is_same_v || + std::is_same_v>) + { + FEValuesExtractors::Vector velocities(0); + FEValuesExtractors::Scalar pressure(dim); + + ComponentMask velocity_mask = fe->component_mask(velocities); + ComponentMask pressure_mask = fe->component_mask(pressure); + + const std::vector index_set_velocity = + DoFTools::locally_owned_dofs_per_component(dof_handler, velocity_mask); + const std::vector index_set_pressure = + DoFTools::locally_owned_dofs_per_component(dof_handler, pressure_mask); + + double local_sum = 0.0; + double local_max = std::numeric_limits::lowest(); + + for (unsigned int d = 0; d < dim; ++d) + { + for (auto j = index_set_velocity[d].begin(); + j != index_set_velocity[d].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_velocity_l2_norm = + std::sqrt(Utilities::MPI::sum(local_sum, this->mpi_communicator)); + double global_velocity_linfty_norm = + Utilities::MPI::max(local_max, this->mpi_communicator); + + local_sum = 0.0; + local_max = std::numeric_limits::lowest(); + + for (auto j = index_set_pressure[dim].begin(); + j != index_set_pressure[dim].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_pressure_l2_norm = + std::sqrt(Utilities::MPI::sum(local_sum, this->mpi_communicator)); + double global_pressure_linfty_norm = + Utilities::MPI::max(local_max, this->mpi_communicator); + + this->pcout << std::setprecision(display_precision) + << "\n\t||du||_L2 = " << std::setw(6) + << global_velocity_l2_norm << std::setw(6) + << "\t||du||_Linfty = " + << std::setprecision(display_precision) + << global_velocity_linfty_norm << std::endl; + this->pcout << std::setprecision(display_precision) + << "\t||dp||_L2 = " << std::setw(6) << global_pressure_l2_norm + << std::setw(6) << "\t||dp||_Linfty = " + << std::setprecision(display_precision) + << global_pressure_linfty_norm << std::endl; + } + if constexpr (std::is_same_v) + { + this->pcout << std::setprecision(display_precision) + << "\t||du||_L2 = " << std::setw(6) + << newton_update.block(0).l2_norm() << std::setw(6) + << "\t||du||_Linfty = " + << std::setprecision(display_precision) + << newton_update.block(0).linfty_norm() << std::endl; + this->pcout << std::setprecision(display_precision) + << "\t||dp||_L2 = " << std::setw(6) + << newton_update.block(1).l2_norm() << std::setw(6) + << "\t||dp||_Linfty = " + << std::setprecision(display_precision) + << newton_update.block(1).linfty_norm() << std::endl; + } +} + template inline VectorType NavierStokesBase::init_temporary_vector() diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index a5ffc8f6ca..d96429fe0a 100644 --- a/source/solvers/tracer.cc +++ b/source/solvers/tracer.cc @@ -863,7 +863,8 @@ Tracer::solve_linear_system(const bool initial_step, 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/vof.cc b/source/solvers/vof.cc index cae784a302..5888fb1b8c 100644 --- a/source/solvers/vof.cc +++ b/source/solvers/vof.cc @@ -2558,7 +2558,8 @@ VolumeOfFluid::solve_linear_system(const bool initial_step, 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; } // Update constraints and newton vectors diff --git a/tests/core/inexact_newton_non_linear_solver_01.output b/tests/core/inexact_newton_non_linear_solver_01.output index fcbfba0964..5e0a806fa8 100644 --- a/tests/core/inexact_newton_non_linear_solver_01.output +++ b/tests/core/inexact_newton_non_linear_solver_01.output @@ -1,4 +1,9 @@ DEAL::Creating solver DEAL::Solving non-linear system -DEAL::The final solution is : 1.22474 -1.50000 +DEAL:: ||dx||_L2 = 1.521 ||dx||_Linfty = 1.500 +DEAL:: ||dx||_L2 = 0.03125 ||dx||_Linfty = 0.03125 +DEAL:: ||dx||_L2 = 0.006010 ||dx||_Linfty = 0.006010 +DEAL:: ||dx||_L2 = 1.482e-05 ||dx||_Linfty = 1.482e-05 +DEAL:: ||dx||_L2 = 7.297e-08 ||dx||_Linfty = 7.297e-08 +DEAL::The final solution is : 1.225 -1.500 diff --git a/tests/core/inexact_newton_non_linear_solver_02.output b/tests/core/inexact_newton_non_linear_solver_02.output index fcbfba0964..3ce63a7073 100644 --- a/tests/core/inexact_newton_non_linear_solver_02.output +++ b/tests/core/inexact_newton_non_linear_solver_02.output @@ -1,4 +1,13 @@ DEAL::Creating solver DEAL::Solving non-linear system -DEAL::The final solution is : 1.22474 -1.50000 +DEAL:: ||dx||_L2 = 1.521 ||dx||_Linfty = 1.500 +DEAL:: ||dx||_L2 = 0.03125 ||dx||_Linfty = 0.03125 +DEAL:: ||dx||_L2 = 0.006010 ||dx||_Linfty = 0.006010 +DEAL:: ||dx||_L2 = 1.482e-05 ||dx||_Linfty = 1.482e-05 +DEAL:: ||dx||_L2 = 7.297e-08 ||dx||_Linfty = 7.297e-08 +DEAL:: ||dx||_L2 = 1.514 ||dx||_Linfty = 1.500 +DEAL:: ||dx||_L2 = 0.01956 ||dx||_Linfty = 0.01956 +DEAL:: ||dx||_L2 = 6.168e-05 ||dx||_Linfty = 6.168e-05 +DEAL:: ||dx||_L2 = 3.019e-07 ||dx||_Linfty = 3.019e-07 +DEAL::The final solution is : 1.225 -1.500 diff --git a/tests/core/newton_non_linear_solver_01.output b/tests/core/newton_non_linear_solver_01.output index fcbfba0964..a6ea433303 100644 --- a/tests/core/newton_non_linear_solver_01.output +++ b/tests/core/newton_non_linear_solver_01.output @@ -1,4 +1,8 @@ DEAL::Creating solver DEAL::Solving non-linear system -DEAL::The final solution is : 1.22474 -1.50000 +DEAL:: ||dx||_L2 = 1.521 ||dx||_Linfty = 1.500 +DEAL:: ||dx||_L2 = 0.02500 ||dx||_Linfty = 0.02500 +DEAL:: ||dx||_L2 = 0.0002551 ||dx||_Linfty = 0.0002551 +DEAL:: ||dx||_L2 = 2.657e-08 ||dx||_Linfty = 2.657e-08 +DEAL::The final solution is : 1.225 -1.500 diff --git a/tests/core/non_linear_test_system_01.h b/tests/core/non_linear_test_system_01.h index 92657c34ff..5d5599cb2b 100644 --- a/tests/core/non_linear_test_system_01.h +++ b/tests/core/non_linear_test_system_01.h @@ -126,6 +126,15 @@ class NonLinearProblemTestClass : public PhysicsSolver> { return newton_update; }; + virtual void + output_newton_update_norms(const unsigned int display_precision) override + { + deallog << 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; + }; virtual Vector & get_present_solution() override { diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc new file mode 100644 index 0000000000..bcaedee2e3 --- /dev/null +++ b/tests/core/vector_problem.cc @@ -0,0 +1,168 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2019 - by the Lethe authors + * + * This file is part of the Lethe library + * + * The Lethe library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 3.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the Lethe distribution. + * + * ---------------------------------------------------------------------*/ + +/** + * @brief This code tests dividing a solution vector in terms of the velocity and pressure components. + */ + +// Deal.II includes +#include + +#include + +#include +#include +#include + +#include + +#include + +#include + +// Lethe +#include +#include +#include + +#include + +// Tests +#include <../tests/tests.h> + +void +test() +{ + Triangulation<3> tria(typename Triangulation<3>::MeshSmoothing( + Triangulation<3>::smoothing_on_refinement | + Triangulation<3>::smoothing_on_coarsening)); + + GridGenerator::hyper_cube(tria, -1, 1); + DoFHandler<3> dof_handler(tria); + + FESystem<3> fe(FE_Q<3>(1), 3, FE_Q<3>(1), 1); + + dof_handler.distribute_dofs(fe); + + Vector dummy_solution; + + dummy_solution.reinit(dof_handler.n_dofs()); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + Vector cell_dummy_solution(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_dummy_solution = 0; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto comp_i = fe.system_to_component_index(i).first; + + cell_dummy_solution(i) = 1.0 * (float(comp_i) + 1); + } + + cell->get_dof_indices(local_dof_indices); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + dummy_solution(local_dof_indices[i]) += cell_dummy_solution(i); + } + + FEValuesExtractors::Vector velocities(0); + FEValuesExtractors::Scalar pressure(3); + + ComponentMask velocity_mask = fe.component_mask(velocities); + ComponentMask pressure_mask = fe.component_mask(pressure); + + + std::vector index_set_velocity = + DoFTools::locally_owned_dofs_per_component(dof_handler, velocity_mask); + std::vector index_set_pressure = + DoFTools::locally_owned_dofs_per_component(dof_handler, pressure_mask); + + double correction_norm = 0.0; + double max_correction = std::numeric_limits::lowest(); + + for (unsigned int d = 0; d < 3; ++d) + { + for (auto j = index_set_velocity[d].begin(); + j != index_set_velocity[d].end(); + j++) + { + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + max_correction = + std::max(max_correction, std::abs(dummy_solution[*j])); + } + } + + deallog << "||u||_L2 : " << std::sqrt(correction_norm) << std::endl; + deallog << "||u||_Linfty : " << max_correction << std::endl; + + correction_norm = 0.0; + max_correction = DBL_MIN; + + for (auto j = index_set_pressure[3].begin(); j != index_set_pressure[3].end(); + j++) + { + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + max_correction = std::max(max_correction, std::abs(dummy_solution[*j])); + } + + deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; + deallog << "||p||_Linfty : " << max_correction << std::endl; +} + +int +main(int argc, char **argv) +{ + try + { + initlog(); + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + test(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/core/vector_problem.output b/tests/core/vector_problem.output new file mode 100644 index 0000000000..2d4225a148 --- /dev/null +++ b/tests/core/vector_problem.output @@ -0,0 +1,5 @@ + +DEAL::||u||_L2 : 10.5830 +DEAL::||u||_Linfty : 3.00000 +DEAL::||p||_L2 : 11.3137 +DEAL::||p||_Linfty : 4.00000