From bff2529f6060579e6ebe3525de9336c7e672a17c Mon Sep 17 00:00:00 2001 From: hepap Date: Fri, 15 Mar 2024 15:37:15 -0400 Subject: [PATCH 01/18] Add unit test for to split a dummy solution vector into its velocity and pressure components --- tests/core/vector_problem.cc | 168 +++++++++++++++++++++++++++++++ tests/core/vector_problem.output | 0 2 files changed, 168 insertions(+) create mode 100644 tests/core/vector_problem.cc create mode 100644 tests/core/vector_problem.output diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc new file mode 100644 index 0000000000..d21aaa4886 --- /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. + * + * --------------------------------------------------------------------- + +* +* Author: Audrey Collard-Daigneault, Polytechnique Montreal, 2020- +*/ + +/** + * @brief This code tests averaging values in time with Trilinos vectors. + */ + +// Deal.II includes +#include + +#include + +#include +#include +#include + +#include + +#include + +#include + +// Lethe +#include +#include +#include + +#include + +// Tests +#include <../tests/tests.h> + +void +test() +{ + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + + 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>(2), 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; + + if (comp_i < 3) + { + 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 = DoFTools::locally_owned_dofs_per_component(dof_handler, velocity_mask); + + Vector velocity_solution(3*index_set[0].n_elements()); + + unsigned int k = 0; + for (unsigned int d = 0; d < 3; ++d) + { + + unsigned int index_set_size = index_set[d].n_elements(); + + std::cout << "IndexSet comp " << d << " size = " << index_set_size << std::endl; + + for (auto j = index_set[d].begin(); j !=index_set[d].end(); j++, k++) + { + velocity_solution[k] = dummy_solution[*j]; + + std::cout << "DoF index = " << *j << " dummy_velocity = " << velocity_solution[k] <<'\n'; + } + } + + for (unsigned int i = 0; i < 3*index_set[0].n_elements(); ++i) + { + auto comp_i = i/index_set[0].n_elements(); + + std::cout << "velocity comp_i " << comp_i << " = " << velocity_solution[i] <<'\n'; + } +} + +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..e69de29bb2 From 514565e4bd2274507cd63cc11f823168b217a5c9 Mon Sep 17 00:00:00 2001 From: hepap Date: Mon, 18 Mar 2024 10:06:52 -0400 Subject: [PATCH 02/18] Complete unit test --- tests/core/vector_problem.cc | 45 ++++++++++++++++---------------- tests/core/vector_problem.output | 5 ++++ 2 files changed, 27 insertions(+), 23 deletions(-) diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc index d21aaa4886..c3ee9250a0 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -49,8 +49,6 @@ void test() { - MPI_Comm mpi_communicator(MPI_COMM_WORLD); - Triangulation<3> tria( typename Triangulation<3>::MeshSmoothing( Triangulation<3>::smoothing_on_refinement | @@ -59,7 +57,7 @@ test() GridGenerator::hyper_cube(tria, -1, 1); DoFHandler<3> dof_handler(tria); - FESystem<3> fe(FE_Q<3>(2), 3, FE_Q<3>(1), 1); + FESystem<3> fe(FE_Q<3>(1), 3, FE_Q<3>(1), 1); dof_handler.distribute_dofs(fe); @@ -79,11 +77,9 @@ test() 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); - if (comp_i < 3) - { - cell_dummy_solution(i) = 1.0*(float(comp_i)+1); - } } cell->get_dof_indices(local_dof_indices); @@ -100,32 +96,35 @@ test() ComponentMask pressure_mask = fe.component_mask(pressure); - std::vector index_set = DoFTools::locally_owned_dofs_per_component(dof_handler, velocity_mask); + 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); - Vector velocity_solution(3*index_set[0].n_elements()); + Vector velocity_solution(3*index_set_velocity[0].n_elements()); + Vector pressure_solution(index_set_velocity[0].n_elements()); unsigned int k = 0; for (unsigned int d = 0; d < 3; ++d) - { - - unsigned int index_set_size = index_set[d].n_elements(); + { + unsigned int index_set_velocity_size = index_set_velocity[d].n_elements(); - std::cout << "IndexSet comp " << d << " size = " << index_set_size << std::endl; - - for (auto j = index_set[d].begin(); j !=index_set[d].end(); j++, k++) + for (auto j = index_set_velocity[d].begin(); j !=index_set_velocity[d].end(); j++, k++) { velocity_solution[k] = dummy_solution[*j]; - - std::cout << "DoF index = " << *j << " dummy_velocity = " << velocity_solution[k] <<'\n'; } } - for (unsigned int i = 0; i < 3*index_set[0].n_elements(); ++i) - { - auto comp_i = i/index_set[0].n_elements(); - - std::cout << "velocity comp_i " << comp_i << " = " << velocity_solution[i] <<'\n'; - } + k = 0; + for (auto j = index_set_pressure[3].begin(); j !=index_set_pressure[3].end(); j++, k++) + { + pressure_solution[k] = dummy_solution[*j]; + } + + deallog << "||u||_L2 : " << velocity_solution.l2_norm() << std::endl; + deallog << "||u||_Linfty : " << velocity_solution.linfty_norm() << std::endl; + + deallog << "||p||_L2 : " << pressure_solution.l2_norm() << std::endl; + deallog << "||p||_Linfty : " << pressure_solution.linfty_norm() << std::endl; + } int diff --git a/tests/core/vector_problem.output b/tests/core/vector_problem.output index e69de29bb2..2d4225a148 100644 --- a/tests/core/vector_problem.output +++ 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 From b043da098d01e6b62a78652493ef27ad547392e0 Mon Sep 17 00:00:00 2001 From: hepap Date: Tue, 19 Mar 2024 17:03:19 -0400 Subject: [PATCH 03/18] Add NS correction norm calculation --- include/core/newton_non_linear_solver.h | 1 + include/core/physics_solver.h | 2 + include/solvers/cahn_hilliard.h | 5 ++ include/solvers/heat_transfer.h | 6 ++ include/solvers/navier_stokes_base.h | 3 + include/solvers/tracer.h | 5 ++ include/solvers/vof.h | 5 ++ source/solvers/navier_stokes_base.cc | 77 +++++++++++++++++++++++++ tests/core/vector_problem.cc | 11 ++-- 9 files changed, 109 insertions(+), 6 deletions(-) diff --git a/include/core/newton_non_linear_solver.h b/include/core/newton_non_linear_solver.h index 2f1840b0e7..dc27d66fc8 100644 --- a/include/core/newton_non_linear_solver.h +++ b/include/core/newton_non_linear_solver.h @@ -139,6 +139,7 @@ NewtonNonLinearSolver::solve(const bool is_initial_step) << "\tL^infty(dx) = " << std::setprecision(this->params.display_precision) << newton_update.linfty_norm() << std::endl; + solver->pcout << solver->get_newton_update_norms_output(this->params.display_precision) << std::endl; } // 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..85a4893ab9 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -98,6 +98,8 @@ class PhysicsSolver get_local_evaluation_point() = 0; virtual VectorType & get_newton_update() = 0; + virtual std::string + get_newton_update_norms_output(const unsigned int display_precision) = 0; virtual VectorType & get_present_solution() = 0; virtual VectorType & diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index 8882f66100..7850d2f691 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -293,6 +293,11 @@ class CahnHilliard : public AuxiliaryPhysics { return newton_update; } + std::string + get_newton_update_norms_output(const unsigned int display_precision) override + { + return "boop"; + } GlobalVectorType & get_present_solution() override { diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index 7fa5b2ffce..a61ecfaa61 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -341,6 +341,12 @@ class HeatTransfer : public AuxiliaryPhysics { return newton_update; } + + std::string + get_newton_update_norms_output(const unsigned int display_precision) override + { + return "boop"; + } /** * @brief Getter method to access the private attribute diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index 02aa4c3cc3..a18ab47245 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -48,6 +48,7 @@ #include #include +#include #include @@ -147,6 +148,8 @@ class NavierStokesBase : public PhysicsSolver { return newton_update; }; + virtual std::string + get_newton_update_norms_output(const unsigned int display_precision) override; virtual VectorType & get_present_solution() override { diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index e15d64eceb..0baf5d5f03 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -254,6 +254,11 @@ class Tracer : public AuxiliaryPhysics { return newton_update; } + std::string + get_newton_update_norms_output(const unsigned int display_precision) override + { + return "boop"; + } GlobalVectorType & get_present_solution() override { diff --git a/include/solvers/vof.h b/include/solvers/vof.h index a455041898..0e7d6bce49 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -379,6 +379,11 @@ class VolumeOfFluid : public AuxiliaryPhysics { return newton_update; } + std::string + get_newton_update_norms_output(const unsigned int display_precision) override + { + return "boop"; + } GlobalVectorType & get_present_solution() override { diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 352743e172..168045b8aa 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2483,6 +2483,83 @@ NavierStokesBase:: } } +template +std::string +NavierStokesBase::get_newton_update_norms_output(const unsigned int display_precision) +{ + if constexpr (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); + + VectorType velocity_correction = init_temporary_vector(); + VectorType pressure_correction = init_temporary_vector(); + + double local_sum = 0.0; + double local_max = DBL_MIN; + + 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]; + + velocity_correction[*j] = dof_newton_update; + + local_sum += dof_newton_update*dof_newton_update; + + if (dof_newton_update > local_max) + { + local_max = 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 = DBL_MIN; + + for (auto j = index_set_pressure[dim].begin(); j !=index_set_pressure[dim].end(); j++) + { + double dof_newton_update = newton_update[*j]; + + pressure_correction[*j] = dof_newton_update; + + local_sum += dof_newton_update*dof_newton_update; + + if (dof_newton_update > local_max) + { + local_max = 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); + + std::string velocity_norms_output = "\t||du||_L2 = " + Utilities::to_string(velocity_correction.l2_norm()) + "\t||du||_Linfty = " + Utilities::to_string(velocity_correction.linfty_norm()) + "\n" + "\t||du||_L2 = " + Utilities::to_string(global_velocity_l2_norm) + "\t||du||_Linfty = " + Utilities::to_string(global_velocity_linfty_norm) + "\n"; + + std::string pressure_norms_output = "\t||dp||_L2 = " + Utilities::to_string(pressure_correction.l2_norm(), 6) + "\t||dp||_Linfty = " + Utilities::to_string(pressure_correction.linfty_norm(), 6) + "\n" + "\t||dp||_L2 = " + Utilities::to_string(global_pressure_l2_norm, 6) + "\t||dp||_Linfty = " + Utilities::to_string(global_pressure_linfty_norm, 6); + + return velocity_norms_output + pressure_norms_output; + + } + else + { + std::string correction_norms_output = "\t||dx||_L2 = " + Utilities::to_string(newton_update.l2_norm(), 6) + "\t||dx||_Linfty = " + Utilities::to_string(newton_update.linfty_norm(), 6); + + return correction_norms_output; + } +}; + template inline VectorType NavierStokesBase::init_temporary_vector() diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc index c3ee9250a0..8ef7d30c8f 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -98,25 +98,24 @@ test() 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); - - Vector velocity_solution(3*index_set_velocity[0].n_elements()); - Vector pressure_solution(index_set_velocity[0].n_elements()); + + Vector velocity_solution(dof_handler.n_locally_owned_dofs()); + Vector pressure_solution(dof_handler.n_locally_owned_dofs()); unsigned int k = 0; for (unsigned int d = 0; d < 3; ++d) { - unsigned int index_set_velocity_size = index_set_velocity[d].n_elements(); for (auto j = index_set_velocity[d].begin(); j !=index_set_velocity[d].end(); j++, k++) { - velocity_solution[k] = dummy_solution[*j]; + velocity_solution[*j] = dummy_solution[*j]; } } k = 0; for (auto j = index_set_pressure[3].begin(); j !=index_set_pressure[3].end(); j++, k++) { - pressure_solution[k] = dummy_solution[*j]; + pressure_solution[*j] = dummy_solution[*j]; } deallog << "||u||_L2 : " << velocity_solution.l2_norm() << std::endl; From 2153fc0d57876a1ee48dcf361463d37af550e7f2 Mon Sep 17 00:00:00 2001 From: hepap Date: Tue, 19 Mar 2024 17:24:00 -0400 Subject: [PATCH 04/18] Modify console output --- include/core/newton_non_linear_solver.h | 10 +++----- include/core/physics_solver.h | 2 +- include/solvers/cahn_hilliard.h | 4 ++-- include/solvers/heat_transfer.h | 4 ++-- include/solvers/navier_stokes_base.h | 2 +- include/solvers/tracer.h | 4 ++-- include/solvers/vof.h | 4 ++-- source/solvers/navier_stokes_base.cc | 31 +++++++++---------------- 8 files changed, 24 insertions(+), 37 deletions(-) diff --git a/include/core/newton_non_linear_solver.h b/include/core/newton_non_linear_solver.h index dc27d66fc8..0df7972155 100644 --- a/include/core/newton_non_linear_solver.h +++ b/include/core/newton_non_linear_solver.h @@ -133,13 +133,9 @@ 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; - solver->pcout << solver->get_newton_update_norms_output(this->params.display_precision) << std::endl; + << std::setw(6) << current_res; + + solver->get_newton_update_norms_output(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 85a4893ab9..6d210063b4 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -98,7 +98,7 @@ class PhysicsSolver get_local_evaluation_point() = 0; virtual VectorType & get_newton_update() = 0; - virtual std::string + virtual void get_newton_update_norms_output(const unsigned int display_precision) = 0; virtual VectorType & get_present_solution() = 0; diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index 7850d2f691..22e93d692c 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -293,10 +293,10 @@ class CahnHilliard : public AuxiliaryPhysics { return newton_update; } - std::string + void get_newton_update_norms_output(const unsigned int display_precision) override { - return "boop"; + this->pcout << "boop" << std::endl; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index a61ecfaa61..0e857bceaa 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -342,10 +342,10 @@ class HeatTransfer : public AuxiliaryPhysics return newton_update; } - std::string + void get_newton_update_norms_output(const unsigned int display_precision) override { - return "boop"; + this->pcout << "boop" << std::endl; } /** diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index a18ab47245..b88ecca6d8 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -148,7 +148,7 @@ class NavierStokesBase : public PhysicsSolver { return newton_update; }; - virtual std::string + virtual void get_newton_update_norms_output(const unsigned int display_precision) override; virtual VectorType & get_present_solution() override diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index 0baf5d5f03..4d76a3e0f7 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -254,10 +254,10 @@ class Tracer : public AuxiliaryPhysics { return newton_update; } - std::string + void get_newton_update_norms_output(const unsigned int display_precision) override { - return "boop"; + this->pcout << "boop" << std::endl; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/vof.h b/include/solvers/vof.h index 0e7d6bce49..920e7de89c 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -379,10 +379,10 @@ class VolumeOfFluid : public AuxiliaryPhysics { return newton_update; } - std::string + void get_newton_update_norms_output(const unsigned int display_precision) override { - return "boop"; + this->pcout << "boop" << std::endl; } GlobalVectorType & get_present_solution() override diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 168045b8aa..e5c5d8e480 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2484,7 +2484,7 @@ NavierStokesBase:: } template -std::string +void NavierStokesBase::get_newton_update_norms_output(const unsigned int display_precision) { if constexpr (std::is_same_v) @@ -2499,9 +2499,6 @@ NavierStokesBase::get_newton_update_norms_output(cons 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); - VectorType velocity_correction = init_temporary_vector(); - VectorType pressure_correction = init_temporary_vector(); - double local_sum = 0.0; double local_max = DBL_MIN; @@ -2511,8 +2508,6 @@ NavierStokesBase::get_newton_update_norms_output(cons { double dof_newton_update = newton_update[*j]; - velocity_correction[*j] = dof_newton_update; - local_sum += dof_newton_update*dof_newton_update; if (dof_newton_update > local_max) @@ -2532,8 +2527,6 @@ NavierStokesBase::get_newton_update_norms_output(cons { double dof_newton_update = newton_update[*j]; - pressure_correction[*j] = dof_newton_update; - local_sum += dof_newton_update*dof_newton_update; if (dof_newton_update > local_max) @@ -2545,18 +2538,16 @@ NavierStokesBase::get_newton_update_norms_output(cons 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); - std::string velocity_norms_output = "\t||du||_L2 = " + Utilities::to_string(velocity_correction.l2_norm()) + "\t||du||_Linfty = " + Utilities::to_string(velocity_correction.linfty_norm()) + "\n" + "\t||du||_L2 = " + Utilities::to_string(global_velocity_l2_norm) + "\t||du||_Linfty = " + Utilities::to_string(global_velocity_linfty_norm) + "\n"; - - std::string pressure_norms_output = "\t||dp||_L2 = " + Utilities::to_string(pressure_correction.l2_norm(), 6) + "\t||dp||_Linfty = " + Utilities::to_string(pressure_correction.linfty_norm(), 6) + "\n" + "\t||dp||_L2 = " + Utilities::to_string(global_pressure_l2_norm, 6) + "\t||dp||_Linfty = " + Utilities::to_string(global_pressure_linfty_norm, 6); - - return velocity_norms_output + pressure_norms_output; - - } - else - { - std::string correction_norms_output = "\t||dx||_L2 = " + Utilities::to_string(newton_update.l2_norm(), 6) + "\t||dx||_Linfty = " + Utilities::to_string(newton_update.linfty_norm(), 6); - - return correction_norms_output; + 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; } }; From 0aa46b98b82f39759d36681789ba69aa69e0ccc2 Mon Sep 17 00:00:00 2001 From: hepap Date: Wed, 20 Mar 2024 16:06:05 -0400 Subject: [PATCH 05/18] Finish implementation --- .../core/inexact_newton_non_linear_solver.h | 9 ++--- include/solvers/cahn_hilliard.h | 6 ++- include/solvers/heat_transfer.h | 6 ++- include/solvers/tracer.h | 6 ++- include/solvers/vof.h | 6 ++- source/solvers/cahn_hilliard.cc | 3 +- source/solvers/gd_navier_stokes.cc | 3 +- 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 | 18 +++++++-- source/solvers/tracer.cc | 3 +- source/solvers/vof.cc | 3 +- tests/core/vector_problem.cc | 40 +++++++++++++------ 14 files changed, 81 insertions(+), 34 deletions(-) diff --git a/include/core/inexact_newton_non_linear_solver.h b/include/core/inexact_newton_non_linear_solver.h index dd8ea70295..9060fbd1a3 100644 --- a/include/core/inexact_newton_non_linear_solver.h +++ b/include/core/inexact_newton_non_linear_solver.h @@ -142,12 +142,9 @@ 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->get_newton_update_norms_output(this->params.display_precision); } // If it's not the first iteration of alpha check if the residual is diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index 22e93d692c..046c8203a0 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -296,7 +296,11 @@ class CahnHilliard : public AuxiliaryPhysics void get_newton_update_norms_output(const unsigned int display_precision) override { - this->pcout << "boop" << std::endl; + 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; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index 0e857bceaa..c4e05f8d00 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -345,7 +345,11 @@ class HeatTransfer : public AuxiliaryPhysics void get_newton_update_norms_output(const unsigned int display_precision) override { - this->pcout << "boop" << std::endl; + 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; } /** diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index 4d76a3e0f7..153bf0149d 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -257,7 +257,11 @@ class Tracer : public AuxiliaryPhysics void get_newton_update_norms_output(const unsigned int display_precision) override { - this->pcout << "boop" << std::endl; + 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; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/vof.h b/include/solvers/vof.h index 920e7de89c..db4b5ab8a1 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -382,7 +382,11 @@ class VolumeOfFluid : public AuxiliaryPhysics void get_newton_update_norms_output(const unsigned int display_precision) override { - this->pcout << "boop" << std::endl; + 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; } GlobalVectorType & get_present_solution() override diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc index d7cbf7fd9a..3041f5c44c 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); diff --git a/source/solvers/gd_navier_stokes.cc b/source/solvers/gd_navier_stokes.cc index 248fb786fc..228309c7c4 100644 --- a/source/solvers/gd_navier_stokes.cc +++ b/source/solvers/gd_navier_stokes.cc @@ -1154,7 +1154,8 @@ 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..25564412d0 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..4b11f1c421 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..6797572668 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 e5c5d8e480..7845d3c23c 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2487,15 +2487,14 @@ template void NavierStokesBase::get_newton_update_norms_output(const unsigned int display_precision) { - if constexpr (std::is_same_v) + 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); @@ -2549,6 +2548,19 @@ NavierStokesBase::get_newton_update_norms_output(cons << 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 diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index a5ffc8f6ca..0b0540a2e1 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..e177dea3a7 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/vector_problem.cc b/tests/core/vector_problem.cc index 8ef7d30c8f..4b312852d3 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -99,30 +99,44 @@ test() 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); - Vector velocity_solution(dof_handler.n_locally_owned_dofs()); - Vector pressure_solution(dof_handler.n_locally_owned_dofs()); + // Vector velocity_solution(dof_handler.n_locally_owned_dofs()); + // Vector pressure_solution(dof_handler.n_locally_owned_dofs()); + + double correction_norm = 0.0; + double max_correction = DBL_MIN; - unsigned int k = 0; for (unsigned int d = 0; d < 3; ++d) { - for (auto j = index_set_velocity[d].begin(); j !=index_set_velocity[d].end(); j++, k++) + for (auto j = index_set_velocity[d].begin(); j !=index_set_velocity[d].end(); j++) { - velocity_solution[*j] = dummy_solution[*j]; + correction_norm += dummy_solution[*j]*dummy_solution[*j]; + + if (dummy_solution[*j]>max_correction) + { + max_correction = dummy_solution[*j]; + } } } - k = 0; - for (auto j = index_set_pressure[3].begin(); j !=index_set_pressure[3].end(); j++, k++) + 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++) { - pressure_solution[*j] = dummy_solution[*j]; + correction_norm += dummy_solution[*j]*dummy_solution[*j]; + + if (dummy_solution[*j]>max_correction) + { + max_correction = dummy_solution[*j]; + } } - deallog << "||u||_L2 : " << velocity_solution.l2_norm() << std::endl; - deallog << "||u||_Linfty : " << velocity_solution.linfty_norm() << std::endl; - - deallog << "||p||_L2 : " << pressure_solution.l2_norm() << std::endl; - deallog << "||p||_Linfty : " << pressure_solution.linfty_norm() << std::endl; + deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; + deallog << "||p||_Linfty : " << max_correction << std::endl; } From df0bf54734fdf62b1adc555bd72bc7af65a31fbd Mon Sep 17 00:00:00 2001 From: hepap Date: Wed, 20 Mar 2024 16:06:29 -0400 Subject: [PATCH 06/18] Indent --- .../core/inexact_newton_non_linear_solver.h | 5 +- include/core/newton_non_linear_solver.h | 5 +- include/solvers/cahn_hilliard.h | 10 +- include/solvers/heat_transfer.h | 12 +- include/solvers/navier_stokes_base.h | 2 +- include/solvers/tracer.h | 10 +- include/solvers/vof.h | 10 +- source/solvers/cahn_hilliard.cc | 4 +- source/solvers/gd_navier_stokes.cc | 5 +- source/solvers/gls_navier_stokes.cc | 8 +- source/solvers/heat_transfer.cc | 4 +- source/solvers/mf_navier_stokes.cc | 4 +- source/solvers/navier_stokes_base.cc | 156 ++++++++++-------- source/solvers/tracer.cc | 4 +- source/solvers/vof.cc | 4 +- tests/core/vector_problem.cc | 116 ++++++------- 16 files changed, 189 insertions(+), 170 deletions(-) diff --git a/include/core/inexact_newton_non_linear_solver.h b/include/core/inexact_newton_non_linear_solver.h index 9060fbd1a3..43bb395e9c 100644 --- a/include/core/inexact_newton_non_linear_solver.h +++ b/include/core/inexact_newton_non_linear_solver.h @@ -143,8 +143,9 @@ InexactNewtonNonLinearSolver::solve(const bool is_initial_step) << std::setw(0) << " res = " << std::setprecision(this->params.display_precision) << std::setw(6) << current_res; - - solver->get_newton_update_norms_output(this->params.display_precision); + + solver->get_newton_update_norms_output( + 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 0df7972155..a1e86bb6f7 100644 --- a/include/core/newton_non_linear_solver.h +++ b/include/core/newton_non_linear_solver.h @@ -134,8 +134,9 @@ NewtonNonLinearSolver::solve(const bool is_initial_step) << std::setw(0) << " res = " << std::setprecision(this->params.display_precision) << std::setw(6) << current_res; - - solver->get_newton_update_norms_output(this->params.display_precision); + + solver->get_newton_update_norms_output( + this->params.display_precision); } // If it's not the first iteration of alpha check if the residual is diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index 046c8203a0..cb2e64d745 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -296,11 +296,11 @@ class CahnHilliard : public AuxiliaryPhysics void get_newton_update_norms_output(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; + 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; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index c4e05f8d00..f842cc5267 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -341,15 +341,15 @@ class HeatTransfer : public AuxiliaryPhysics { return newton_update; } - + void get_newton_update_norms_output(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; + 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; } /** diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index b88ecca6d8..8f22d144f6 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -46,9 +46,9 @@ #include +#include #include #include -#include #include diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index 153bf0149d..ae4f3cc0b8 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -257,11 +257,11 @@ class Tracer : public AuxiliaryPhysics void get_newton_update_norms_output(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; + 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; } GlobalVectorType & get_present_solution() override diff --git a/include/solvers/vof.h b/include/solvers/vof.h index db4b5ab8a1..4c3cfb66c0 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -382,11 +382,11 @@ class VolumeOfFluid : public AuxiliaryPhysics void get_newton_update_norms_output(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; + 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; } GlobalVectorType & get_present_solution() override diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc index 3041f5c44c..34d14efaad 100644 --- a/source/solvers/cahn_hilliard.cc +++ b/source/solvers/cahn_hilliard.cc @@ -1146,8 +1146,8 @@ CahnHilliard::solve_linear_system(const bool initial_step, .verbosity != Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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/gd_navier_stokes.cc b/source/solvers/gd_navier_stokes.cc index 228309c7c4..a6f6b7ac23 100644 --- a/source/solvers/gd_navier_stokes.cc +++ b/source/solvers/gd_navier_stokes.cc @@ -1154,8 +1154,9 @@ GDNavierStokesSolver::solve_system_GMRES(const bool initial_step, .verbosity != Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " - << solver_control.last_step() << " steps to reach a residual norm of " - << solver_control.last_value() << 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 25564412d0..461e07aebb 100644 --- a/source/solvers/gls_navier_stokes.cc +++ b/source/solvers/gls_navier_stokes.cc @@ -1611,8 +1611,8 @@ GLSNavierStokesSolver::solve_system_GMRES(const bool initial_step, { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << std::endl; + << " steps to reach a residual norm of " + << solver_control.last_value() << std::endl; } } constraints_used.distribute(completely_distributed_solution); @@ -1724,8 +1724,8 @@ GLSNavierStokesSolver::solve_system_BiCGStab( { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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 4b11f1c421..1d4f6a58e1 100644 --- a/source/solvers/heat_transfer.cc +++ b/source/solvers/heat_transfer.cc @@ -1390,8 +1390,8 @@ HeatTransfer::solve_linear_system(const bool initial_step, .verbosity != Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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 6797572668..25691f79b7 100644 --- a/source/solvers/mf_navier_stokes.cc +++ b/source/solvers/mf_navier_stokes.cc @@ -1782,8 +1782,8 @@ MFNavierStokesSolver::solve_system_GMRES(const bool initial_step, .verbosity != Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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 7845d3c23c..680c6f9a3e 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2485,82 +2485,98 @@ NavierStokesBase:: template void -NavierStokesBase::get_newton_update_norms_output(const unsigned int display_precision) +NavierStokesBase::get_newton_update_norms_output( + 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 = DBL_MIN; - - 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; - - if (dof_newton_update > local_max) + 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 = DBL_MIN; + + for (unsigned int d = 0; d < dim; ++d) { - local_max = dof_newton_update; + 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; + + if (dof_newton_update > local_max) + { + local_max = 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 = DBL_MIN; - - 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; - - if (dof_newton_update > local_max) + + 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 = DBL_MIN; + + for (auto j = index_set_pressure[dim].begin(); + j != index_set_pressure[dim].end(); + j++) { - local_max = dof_newton_update; + double dof_newton_update = newton_update[*j]; + + local_sum += dof_newton_update * dof_newton_update; + + if (dof_newton_update > local_max) + { + local_max = 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; - } + + 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; - } + { + 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 diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index 0b0540a2e1..d96429fe0a 100644 --- a/source/solvers/tracer.cc +++ b/source/solvers/tracer.cc @@ -863,8 +863,8 @@ Tracer::solve_linear_system(const bool initial_step, Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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 e177dea3a7..5888fb1b8c 100644 --- a/source/solvers/vof.cc +++ b/source/solvers/vof.cc @@ -2558,8 +2558,8 @@ VolumeOfFluid::solve_linear_system(const bool initial_step, Parameters::Verbosity::quiet) { this->pcout << " -Iterative solver took : " << solver_control.last_step() - << " steps to reach a residual norm of " - << solver_control.last_value() << 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/vector_problem.cc b/tests/core/vector_problem.cc index 4b312852d3..4a177be254 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -26,9 +26,9 @@ #include +#include #include #include -#include #include @@ -49,95 +49,95 @@ void test() { - Triangulation<3> tria( - typename Triangulation<3>::MeshSmoothing( - Triangulation<3>::smoothing_on_refinement | - Triangulation<3>::smoothing_on_coarsening)); - + 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); + + 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) + 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_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); - + + 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); - // Vector velocity_solution(dof_handler.n_locally_owned_dofs()); + + 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); + + // Vector velocity_solution(dof_handler.n_locally_owned_dofs()); // Vector pressure_solution(dof_handler.n_locally_owned_dofs()); - + double correction_norm = 0.0; - double max_correction = DBL_MIN; - + double max_correction = DBL_MIN; + 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]; - - if (dummy_solution[*j]>max_correction) - { - max_correction = dummy_solution[*j]; - } + for (auto j = index_set_velocity[d].begin(); + j != index_set_velocity[d].end(); + j++) + { + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + if (dummy_solution[*j] > max_correction) + { + max_correction = 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++) + 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]; - - if (dummy_solution[*j]>max_correction) - { - max_correction = dummy_solution[*j]; - } + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + if (dummy_solution[*j] > max_correction) + { + max_correction = dummy_solution[*j]; + } } - + deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; deallog << "||p||_Linfty : " << max_correction << std::endl; - } int From e624210e86e54c9f09987f10125f0b04674c5703 Mon Sep 17 00:00:00 2001 From: hepap Date: Wed, 20 Mar 2024 16:36:45 -0400 Subject: [PATCH 07/18] Update test --- .../lethe-fluid/heat_transfer_heat-flux.output | 5 ----- .../lethe-fluid/heat_transfer_heat-flux.prm | 3 +++ tests/core/non_linear_test_system_01.h | 9 +++++++++ 3 files changed, 12 insertions(+), 5 deletions(-) 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/tests/core/non_linear_test_system_01.h b/tests/core/non_linear_test_system_01.h index 92657c34ff..973f439f6a 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 + get_newton_update_norms_output(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; + }; virtual Vector & get_present_solution() override { From 74ee0ff5700d52b9e1b58baee1312849b1db83f9 Mon Sep 17 00:00:00 2001 From: hepap Date: Wed, 20 Mar 2024 16:51:12 -0400 Subject: [PATCH 08/18] Change method name --- include/core/inexact_newton_non_linear_solver.h | 2 +- include/core/newton_non_linear_solver.h | 2 +- include/core/physics_solver.h | 2 +- include/solvers/cahn_hilliard.h | 2 +- include/solvers/heat_transfer.h | 2 +- include/solvers/navier_stokes_base.h | 2 +- include/solvers/tracer.h | 2 +- include/solvers/vof.h | 2 +- source/solvers/navier_stokes_base.cc | 2 +- tests/core/non_linear_test_system_01.h | 2 +- 10 files changed, 10 insertions(+), 10 deletions(-) diff --git a/include/core/inexact_newton_non_linear_solver.h b/include/core/inexact_newton_non_linear_solver.h index 43bb395e9c..43123d83b4 100644 --- a/include/core/inexact_newton_non_linear_solver.h +++ b/include/core/inexact_newton_non_linear_solver.h @@ -144,7 +144,7 @@ InexactNewtonNonLinearSolver::solve(const bool is_initial_step) << std::setprecision(this->params.display_precision) << std::setw(6) << current_res; - solver->get_newton_update_norms_output( + solver->output_newton_update_norms( this->params.display_precision); } diff --git a/include/core/newton_non_linear_solver.h b/include/core/newton_non_linear_solver.h index a1e86bb6f7..d77f90298c 100644 --- a/include/core/newton_non_linear_solver.h +++ b/include/core/newton_non_linear_solver.h @@ -135,7 +135,7 @@ NewtonNonLinearSolver::solve(const bool is_initial_step) << std::setprecision(this->params.display_precision) << std::setw(6) << current_res; - solver->get_newton_update_norms_output( + solver->output_newton_update_norms( this->params.display_precision); } diff --git a/include/core/physics_solver.h b/include/core/physics_solver.h index 6d210063b4..f0320cc34f 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -99,7 +99,7 @@ class PhysicsSolver virtual VectorType & get_newton_update() = 0; virtual void - get_newton_update_norms_output(const unsigned int display_precision) = 0; + output_newton_update_norms(const unsigned int display_precision) = 0; virtual VectorType & get_present_solution() = 0; virtual VectorType & diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index cb2e64d745..7093d634e5 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -294,7 +294,7 @@ class CahnHilliard : public AuxiliaryPhysics return newton_update; } void - get_newton_update_norms_output(const unsigned int display_precision) override + 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() diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index f842cc5267..fc87bc7faf 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -343,7 +343,7 @@ class HeatTransfer : public AuxiliaryPhysics } void - get_newton_update_norms_output(const unsigned int display_precision) override + 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() diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index 8f22d144f6..3f20ef7439 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -149,7 +149,7 @@ class NavierStokesBase : public PhysicsSolver return newton_update; }; virtual void - get_newton_update_norms_output(const unsigned int display_precision) override; + output_newton_update_norms(const unsigned int display_precision) override; virtual VectorType & get_present_solution() override { diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index ae4f3cc0b8..10b4ca13e9 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -255,7 +255,7 @@ class Tracer : public AuxiliaryPhysics return newton_update; } void - get_newton_update_norms_output(const unsigned int display_precision) override + 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() diff --git a/include/solvers/vof.h b/include/solvers/vof.h index 4c3cfb66c0..3be795c971 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -380,7 +380,7 @@ class VolumeOfFluid : public AuxiliaryPhysics return newton_update; } void - get_newton_update_norms_output(const unsigned int display_precision) override + 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() diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 680c6f9a3e..f21d9422eb 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2485,7 +2485,7 @@ NavierStokesBase:: template void -NavierStokesBase::get_newton_update_norms_output( +NavierStokesBase::output_newton_update_norms( const unsigned int display_precision) { if constexpr (std::is_same_v || diff --git a/tests/core/non_linear_test_system_01.h b/tests/core/non_linear_test_system_01.h index 973f439f6a..4b698b756f 100644 --- a/tests/core/non_linear_test_system_01.h +++ b/tests/core/non_linear_test_system_01.h @@ -127,7 +127,7 @@ class NonLinearProblemTestClass : public PhysicsSolver> return newton_update; }; virtual void - get_newton_update_norms_output(const unsigned int display_precision) override + 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() From 1cf37588f37ec04aec4479653f2ccbc43063cb98 Mon Sep 17 00:00:00 2001 From: hepap Date: Thu, 21 Mar 2024 09:00:36 -0400 Subject: [PATCH 09/18] Remove commented lines --- tests/core/vector_problem.cc | 125 +++++++++++++++++------------------ 1 file changed, 59 insertions(+), 66 deletions(-) diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc index 4a177be254..5df55168d1 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -11,14 +11,10 @@ * The full text of the license can be found in the file LICENSE at * the top level of the Lethe distribution. * - * --------------------------------------------------------------------- - -* -* Author: Audrey Collard-Daigneault, Polytechnique Montreal, 2020- -*/ + * ---------------------------------------------------------------------*/ /** - * @brief This code tests averaging values in time with Trilinos vectors. + * @brief This code tests dividing a solution vector in terms of the velocity and pressure components. */ // Deal.II includes @@ -26,9 +22,9 @@ #include -#include #include #include +#include #include @@ -49,95 +45,92 @@ void test() { - Triangulation<3> tria(typename Triangulation<3>::MeshSmoothing( - Triangulation<3>::smoothing_on_refinement | - Triangulation<3>::smoothing_on_coarsening)); - + 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); + + 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) + 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_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); + + 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); - - // Vector velocity_solution(dof_handler.n_locally_owned_dofs()); - // Vector pressure_solution(dof_handler.n_locally_owned_dofs()); - + + + 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 = DBL_MIN; - + double max_correction = DBL_MIN; + for (unsigned int d = 0; d < 3; ++d) + { + + for (auto j = index_set_velocity[d].begin(); j !=index_set_velocity[d].end(); j++) { - for (auto j = index_set_velocity[d].begin(); - j != index_set_velocity[d].end(); - j++) - { - correction_norm += dummy_solution[*j] * dummy_solution[*j]; - - if (dummy_solution[*j] > max_correction) - { - max_correction = dummy_solution[*j]; - } - } + correction_norm += dummy_solution[*j]*dummy_solution[*j]; + + if (dummy_solution[*j]>max_correction) + { + max_correction = 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++) + 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]; - - if (dummy_solution[*j] > max_correction) - { - max_correction = dummy_solution[*j]; - } + correction_norm += dummy_solution[*j]*dummy_solution[*j]; + + if (dummy_solution[*j]>max_correction) + { + max_correction = dummy_solution[*j]; + } } - + deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; deallog << "||p||_Linfty : " << max_correction << std::endl; + } int From 818d990256fda79fa17a66e97b15a5487c7b879a Mon Sep 17 00:00:00 2001 From: hepap Date: Thu, 21 Mar 2024 09:10:17 -0400 Subject: [PATCH 10/18] Add brief --- include/core/physics_solver.h | 10 ++- include/solvers/cahn_hilliard.h | 24 +++--- include/solvers/heat_transfer.h | 24 +++--- include/solvers/navier_stokes_base.h | 10 ++- include/solvers/tracer.h | 24 +++--- include/solvers/vof.h | 24 +++--- tests/core/vector_problem.cc | 114 +++++++++++++-------------- 7 files changed, 132 insertions(+), 98 deletions(-) diff --git a/include/core/physics_solver.h b/include/core/physics_solver.h index f0320cc34f..23d72e22ce 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -98,8 +98,6 @@ class PhysicsSolver get_local_evaluation_point() = 0; virtual VectorType & get_newton_update() = 0; - virtual void - output_newton_update_norms(const unsigned int display_precision) = 0; virtual VectorType & get_present_solution() = 0; virtual VectorType & @@ -107,6 +105,14 @@ class PhysicsSolver virtual AffineConstraints & get_nonzero_constraints() = 0; + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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 7093d634e5..012f006238 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -293,15 +293,6 @@ class CahnHilliard : public AuxiliaryPhysics { return newton_update; } - 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; - } GlobalVectorType & get_present_solution() override { @@ -318,6 +309,21 @@ class CahnHilliard : public AuxiliaryPhysics return nonzero_constraints; } + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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: /** diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index fc87bc7faf..df2d4e55ec 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -342,16 +342,6 @@ class HeatTransfer : public AuxiliaryPhysics return newton_update; } - 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; - } - /** * @brief Getter method to access the private attribute * present_solution for the physic currently solved. NB : present_solution is @@ -392,6 +382,20 @@ class HeatTransfer : public AuxiliaryPhysics return nonzero_constraints; } + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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 3f20ef7439..7dd5fe836d 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -148,8 +148,6 @@ class NavierStokesBase : public PhysicsSolver { return newton_update; }; - virtual void - output_newton_update_norms(const unsigned int display_precision) override; virtual VectorType & get_present_solution() override { @@ -166,6 +164,14 @@ class NavierStokesBase : public PhysicsSolver return nonzero_constraints; }; + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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 10b4ca13e9..dcc425ce66 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -254,15 +254,6 @@ class Tracer : public AuxiliaryPhysics { return newton_update; } - 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; - } GlobalVectorType & get_present_solution() override { @@ -280,6 +271,21 @@ class Tracer : public AuxiliaryPhysics } + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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 3be795c971..61b8396041 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -379,15 +379,6 @@ class VolumeOfFluid : public AuxiliaryPhysics { return newton_update; } - 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; - } GlobalVectorType & get_present_solution() override { @@ -424,6 +415,21 @@ class VolumeOfFluid : public AuxiliaryPhysics return &present_curvature_solution; } + /** + * @brief Output the L2 and Linfty norms of the correction vector. + * + * @param display_precision [in] Number of outputed 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/tests/core/vector_problem.cc b/tests/core/vector_problem.cc index 5df55168d1..60a3d18cf9 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -22,9 +22,9 @@ #include +#include #include #include -#include #include @@ -45,92 +45,92 @@ void test() { - Triangulation<3> tria( - typename Triangulation<3>::MeshSmoothing( - Triangulation<3>::smoothing_on_refinement | - Triangulation<3>::smoothing_on_coarsening)); - + 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); + + 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) + 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_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); - + + 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); - + + + 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 = DBL_MIN; - + double max_correction = DBL_MIN; + 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]; - - if (dummy_solution[*j]>max_correction) - { - max_correction = dummy_solution[*j]; - } + for (auto j = index_set_velocity[d].begin(); + j != index_set_velocity[d].end(); + j++) + { + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + if (dummy_solution[*j] > max_correction) + { + max_correction = 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++) + 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]; - - if (dummy_solution[*j]>max_correction) - { - max_correction = dummy_solution[*j]; - } + correction_norm += dummy_solution[*j] * dummy_solution[*j]; + + if (dummy_solution[*j] > max_correction) + { + max_correction = dummy_solution[*j]; + } } - + deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; deallog << "||p||_Linfty : " << max_correction << std::endl; - } int From db77b89fac4c7532fc2a64fe70b05a17675084f9 Mon Sep 17 00:00:00 2001 From: hepap Date: Thu, 21 Mar 2024 10:03:41 -0400 Subject: [PATCH 11/18] Add comment in the doc on the extra verbose verbosity of the linear solver --- doc/source/parameters/cfd/linear_solver_control.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/source/parameters/cfd/linear_solver_control.rst b/doc/source/parameters/cfd/linear_solver_control.rst index 524521b288..3a63b7c1b4 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 at ``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. From 20f86768577c536a876323536687b83f6c1a4ff3 Mon Sep 17 00:00:00 2001 From: hepap Date: Thu, 21 Mar 2024 10:09:53 -0400 Subject: [PATCH 12/18] Add chanhelog entry --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) 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 From 6d0b5651e5e674a4065a24a78ad0466add485db8 Mon Sep 17 00:00:00 2001 From: hepap Date: Thu, 21 Mar 2024 10:41:14 -0400 Subject: [PATCH 13/18] Remove warning because of an extra ; --- source/solvers/navier_stokes_base.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index f21d9422eb..a1c0889cef 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2577,7 +2577,7 @@ NavierStokesBase::output_newton_update_norms( << std::setprecision(display_precision) << newton_update.block(1).linfty_norm() << std::endl; } -}; +} template inline VectorType From 9dd38be682cae329fcb1e424f48f5727937b557b Mon Sep 17 00:00:00 2001 From: hepap <47506601+hepap@users.noreply.github.com> Date: Tue, 26 Mar 2024 14:03:38 -0400 Subject: [PATCH 14/18] Apply suggestions from code review Co-authored-by: Pierre Laurentin <81017102+PierreLaurentinCS@users.noreply.github.com> Co-authored-by: Amishga Alphonius <107414376+AmishgaAlphonius@users.noreply.github.com> --- doc/source/parameters/cfd/linear_solver_control.rst | 2 +- include/core/physics_solver.h | 2 +- include/solvers/cahn_hilliard.h | 2 +- include/solvers/heat_transfer.h | 2 +- include/solvers/navier_stokes_base.h | 2 +- include/solvers/tracer.h | 2 +- include/solvers/vof.h | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/source/parameters/cfd/linear_solver_control.rst b/doc/source/parameters/cfd/linear_solver_control.rst index 3a63b7c1b4..b6e2537468 100644 --- a/doc/source/parameters/cfd/linear_solver_control.rst +++ b/doc/source/parameters/cfd/linear_solver_control.rst @@ -76,7 +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 at ``extra verbose``, the residual at each iteration of the linear solver is printed. +When set to ``extra verbose``, the residual at each iteration of the linear solver is printed. * The ``minimum residual`` for the linear solver. diff --git a/include/core/physics_solver.h b/include/core/physics_solver.h index 23d72e22ce..3fa340a7ae 100644 --- a/include/core/physics_solver.h +++ b/include/core/physics_solver.h @@ -108,7 +108,7 @@ class PhysicsSolver /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ virtual void output_newton_update_norms(const unsigned int display_precision) = 0; diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index 012f006238..fd3f3b1089 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -312,7 +312,7 @@ class CahnHilliard : public AuxiliaryPhysics /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ void output_newton_update_norms(const unsigned int display_precision) override diff --git a/include/solvers/heat_transfer.h b/include/solvers/heat_transfer.h index df2d4e55ec..021a00ce39 100644 --- a/include/solvers/heat_transfer.h +++ b/include/solvers/heat_transfer.h @@ -385,7 +385,7 @@ class HeatTransfer : public AuxiliaryPhysics /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ void output_newton_update_norms(const unsigned int display_precision) override diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index 7dd5fe836d..eae196c511 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -167,7 +167,7 @@ class NavierStokesBase : public PhysicsSolver /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ virtual void output_newton_update_norms(const unsigned int display_precision) override; diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index dcc425ce66..d7818d2cc6 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -274,7 +274,7 @@ class Tracer : public AuxiliaryPhysics /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ void output_newton_update_norms(const unsigned int display_precision) override diff --git a/include/solvers/vof.h b/include/solvers/vof.h index 61b8396041..4bb0b3e646 100644 --- a/include/solvers/vof.h +++ b/include/solvers/vof.h @@ -418,7 +418,7 @@ class VolumeOfFluid : public AuxiliaryPhysics /** * @brief Output the L2 and Linfty norms of the correction vector. * - * @param display_precision [in] Number of outputed digits. + * @param[in] display_precision Number of outputted digits. */ void output_newton_update_norms(const unsigned int display_precision) override From 3c693ebc8e51bb9db3f3154096004dec51981b23 Mon Sep 17 00:00:00 2001 From: hepap Date: Tue, 26 Mar 2024 16:52:48 -0400 Subject: [PATCH 15/18] Apply comments from code review --- source/solvers/navier_stokes_base.cc | 14 ++++---------- tests/core/non_linear_test_system_01.h | 2 +- tests/core/vector_problem.cc | 7 ++----- 3 files changed, 7 insertions(+), 16 deletions(-) diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index a1c0889cef..4e4b781230 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2504,7 +2504,7 @@ NavierStokesBase::output_newton_update_norms( DoFTools::locally_owned_dofs_per_component(dof_handler, pressure_mask); double local_sum = 0.0; - double local_max = DBL_MIN; + double local_max = std::numeric_limits::lowest(); for (unsigned int d = 0; d < dim; ++d) { @@ -2516,10 +2516,7 @@ NavierStokesBase::output_newton_update_norms( local_sum += dof_newton_update * dof_newton_update; - if (dof_newton_update > local_max) - { - local_max = dof_newton_update; - } + local_max = std::max(local_max, dof_newton_update); } } @@ -2529,7 +2526,7 @@ NavierStokesBase::output_newton_update_norms( Utilities::MPI::max(local_max, this->mpi_communicator); local_sum = 0.0; - local_max = DBL_MIN; + local_max = std::numeric_limits::lowest(); for (auto j = index_set_pressure[dim].begin(); j != index_set_pressure[dim].end(); @@ -2539,10 +2536,7 @@ NavierStokesBase::output_newton_update_norms( local_sum += dof_newton_update * dof_newton_update; - if (dof_newton_update > local_max) - { - local_max = dof_newton_update; - } + local_max = std::max(local_max, dof_newton_update); } double global_pressure_l2_norm = diff --git a/tests/core/non_linear_test_system_01.h b/tests/core/non_linear_test_system_01.h index 4b698b756f..d654709f00 100644 --- a/tests/core/non_linear_test_system_01.h +++ b/tests/core/non_linear_test_system_01.h @@ -129,7 +129,7 @@ class NonLinearProblemTestClass : public PhysicsSolver> virtual void output_newton_update_norms(const unsigned int display_precision) override { - this->pcout << std::setprecision(display_precision) + 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) diff --git a/tests/core/vector_problem.cc b/tests/core/vector_problem.cc index 60a3d18cf9..356e34ceb5 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -95,7 +95,7 @@ test() DoFTools::locally_owned_dofs_per_component(dof_handler, pressure_mask); double correction_norm = 0.0; - double max_correction = DBL_MIN; + double max_correction = std::numeric_limits::lowest(); for (unsigned int d = 0; d < 3; ++d) { @@ -105,10 +105,7 @@ test() { correction_norm += dummy_solution[*j] * dummy_solution[*j]; - if (dummy_solution[*j] > max_correction) - { - max_correction = dummy_solution[*j]; - } + max_correction = std::max(max_correction, dummy_solution[*j]); } } From c186752beab25b30939d63231359ccde1256a77b Mon Sep 17 00:00:00 2001 From: hepap Date: Tue, 26 Mar 2024 17:12:48 -0400 Subject: [PATCH 16/18] Add CH output, in progress --- include/solvers/cahn_hilliard.h | 9 +--- source/solvers/cahn_hilliard.cc | 79 +++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+), 8 deletions(-) diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h index fd3f3b1089..2e1da183ad 100644 --- a/include/solvers/cahn_hilliard.h +++ b/include/solvers/cahn_hilliard.h @@ -315,14 +315,7 @@ class CahnHilliard : public AuxiliaryPhysics * @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; - } + output_newton_update_norms(const unsigned int display_precision) override; private: diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc index 34d14efaad..8c849df67d 100644 --- a/source/solvers/cahn_hilliard.cc +++ b/source/solvers/cahn_hilliard.cc @@ -1322,6 +1322,85 @@ 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, 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, 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||du||_L2 = " << std::setw(6) + << newton_update.l2_norm() << std::setw(6) + << "\t||du||_Linfty = " + << std::setprecision(display_precision) + << newton_update.linfty_norm() << std::endl; + + 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, From f533a1f4eedd8518115fddb38525c8c3bd1960ef Mon Sep 17 00:00:00 2001 From: hepap Date: Fri, 5 Apr 2024 15:14:15 -0400 Subject: [PATCH 17/18] Correct bug --- source/solvers/cahn_hilliard.cc | 121 ++++++++++++------------- source/solvers/navier_stokes_base.cc | 4 +- tests/core/non_linear_test_system_01.h | 8 +- tests/core/vector_problem.cc | 8 +- 4 files changed, 66 insertions(+), 75 deletions(-) diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc index 8c849df67d..d7e492bfdb 100644 --- a/source/solvers/cahn_hilliard.cc +++ b/source/solvers/cahn_hilliard.cc @@ -1327,78 +1327,71 @@ 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]; + auto mpi_communicator = triangulation->get_communicator(); - local_sum += dof_newton_update * dof_newton_update; + FEValuesExtractors::Scalar phase_order(0); + FEValuesExtractors::Scalar chemical_potential(1); - local_max = std::max(local_max, dof_newton_update); - } - + ComponentMask phase_order_mask = fe->component_mask(phase_order); + ComponentMask chemical_potential_mask = + fe->component_mask(chemical_potential); - 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); + 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); - local_sum = 0.0; - local_max = std::numeric_limits::lowest(); + double local_sum = 0.0; + double 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; + 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, 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||du||_L2 = " << std::setw(6) - << newton_update.l2_norm() << std::setw(6) - << "\t||du||_Linfty = " - << std::setprecision(display_precision) - << newton_update.linfty_norm() << std::endl; - - 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; - + 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>> diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 4e4b781230..2c41a40b90 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -2516,7 +2516,7 @@ NavierStokesBase::output_newton_update_norms( local_sum += dof_newton_update * dof_newton_update; - local_max = std::max(local_max, dof_newton_update); + local_max = std::max(local_max, std::abs(dof_newton_update)); } } @@ -2536,7 +2536,7 @@ NavierStokesBase::output_newton_update_norms( local_sum += dof_newton_update * dof_newton_update; - local_max = std::max(local_max, dof_newton_update); + local_max = std::max(local_max, std::abs(dof_newton_update)); } double global_pressure_l2_norm = diff --git a/tests/core/non_linear_test_system_01.h b/tests/core/non_linear_test_system_01.h index d654709f00..5d5599cb2b 100644 --- a/tests/core/non_linear_test_system_01.h +++ b/tests/core/non_linear_test_system_01.h @@ -130,10 +130,10 @@ class NonLinearProblemTestClass : public PhysicsSolver> 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; + << "\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 index 356e34ceb5..bcaedee2e3 100644 --- a/tests/core/vector_problem.cc +++ b/tests/core/vector_problem.cc @@ -105,7 +105,8 @@ test() { correction_norm += dummy_solution[*j] * dummy_solution[*j]; - max_correction = std::max(max_correction, dummy_solution[*j]); + max_correction = + std::max(max_correction, std::abs(dummy_solution[*j])); } } @@ -120,10 +121,7 @@ test() { correction_norm += dummy_solution[*j] * dummy_solution[*j]; - if (dummy_solution[*j] > max_correction) - { - max_correction = dummy_solution[*j]; - } + max_correction = std::max(max_correction, std::abs(dummy_solution[*j])); } deallog << "||p||_L2 : " << std::sqrt(correction_norm) << std::endl; From 7885817868d4659655af613b46e207c0c341e767 Mon Sep 17 00:00:00 2001 From: hepap Date: Fri, 5 Apr 2024 16:20:30 -0400 Subject: [PATCH 18/18] Fix tests --- tests/core/inexact_newton_non_linear_solver_01.output | 7 ++++++- tests/core/inexact_newton_non_linear_solver_02.output | 11 ++++++++++- tests/core/newton_non_linear_solver_01.output | 6 +++++- 3 files changed, 21 insertions(+), 3 deletions(-) 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