From 4e072de184f35c03fb7fb8d5b53c78249aa68e6c Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sat, 1 Jul 2023 07:04:37 -0400 Subject: [PATCH] test examples --- .gitignore | 1 + tests/examples/CMakeLists.txt | 7 + tests/examples/create_diffs.sh | 21 + tests/examples/step-1.diff | 0 tests/examples/step-1.output | 2 + tests/examples/step-2.diff | 0 tests/examples/step-2.output | 0 tests/examples/step-56.cc | 1171 -------------------------------- tests/examples/step-56.diff | 29 + tests/examples/step-56.output | 104 ++- tests/examples/update.sh | 21 + 11 files changed, 149 insertions(+), 1207 deletions(-) create mode 100755 tests/examples/create_diffs.sh create mode 100644 tests/examples/step-1.diff create mode 100644 tests/examples/step-1.output create mode 100644 tests/examples/step-2.diff create mode 100644 tests/examples/step-2.output delete mode 100644 tests/examples/step-56.cc create mode 100644 tests/examples/step-56.diff create mode 100755 tests/examples/update.sh diff --git a/.gitignore b/.gitignore index 4d186c07f0..d5f9ed9011 100644 --- a/.gitignore +++ b/.gitignore @@ -55,3 +55,4 @@ examples/step-*/CMakeCache.txt examples/step-*/cmake_install.cmake examples/step-*/step-? examples/step-*/step-?? +/tests/examples/*.cc diff --git a/tests/examples/CMakeLists.txt b/tests/examples/CMakeLists.txt index db9c1adcf2..1ff5768966 100644 --- a/tests/examples/CMakeLists.txt +++ b/tests/examples/CMakeLists.txt @@ -1,4 +1,11 @@ cmake_minimum_required(VERSION 3.13.4) include(../scripts/setup_testsubproject.cmake) project(testsuite CXX) + +# create tests from diffs: +execute_process( + COMMAND "./update.sh" + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} +) + deal_ii_pickup_tests() diff --git a/tests/examples/create_diffs.sh b/tests/examples/create_diffs.sh new file mode 100755 index 0000000000..515ad3f057 --- /dev/null +++ b/tests/examples/create_diffs.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +# create .diff files from current .cc files + +# Find all *.diff files in the current directory +diff_files=$(find . -maxdepth 1 -type f -name "*.diff") + +for file in $diff_files; do + # Extract the filename without the extension + filename=$(basename "$file" .diff) + src="../../examples/$filename/$filename.cc" + + # Check if the corresponding .cc file exists + if [[ -f "$src" ]]; then + echo "creating $file from $src and $filename.cc" + diff "$src" "$filename.cc" > $file + else + echo "No matching .cc file found for $diff_file" + exit -1 + fi +done diff --git a/tests/examples/step-1.diff b/tests/examples/step-1.diff new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/examples/step-1.output b/tests/examples/step-1.output new file mode 100644 index 0000000000..6aa9888604 --- /dev/null +++ b/tests/examples/step-1.output @@ -0,0 +1,2 @@ +Grid written to grid-1.svg +Grid written to grid-2.svg diff --git a/tests/examples/step-2.diff b/tests/examples/step-2.diff new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/examples/step-2.output b/tests/examples/step-2.output new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/examples/step-56.cc b/tests/examples/step-56.cc deleted file mode 100644 index 8f1ae4877f..0000000000 --- a/tests/examples/step-56.cc +++ /dev/null @@ -1,1171 +0,0 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2016 - 2022 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II 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 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE.md at - * the top level directory of deal.II. - * - * --------------------------------------------------------------------- - - * Author: Ryan Grove, Clemson University - * Timo Heister, Clemson University - */ - -// @sect3{Include files} - -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// We need to include the following file to do timings: -#include - -// This includes the files necessary for us to use geometric Multigrid -#include -#include -#include -#include -#include -#include - -#include -#include - -#include "../tests.h" - -namespace Step56 -{ - // In order to make it easy to switch between the different solvers that are - // being used, we declare an enum that can be passed as an argument to the - // constructor of the main class. - struct SolverType - { - enum type - { - FGMRES_ILU, - FGMRES_GMG, - UMFPACK - }; - }; - - // @sect3{Functions for Solution and Righthand side} - // - // The class Solution is used to define the boundary conditions and to - // compute errors of the numerical solution. Note that we need to define the - // values and gradients in order to compute L2 and H1 errors. Here we - // decided to separate the implementations for 2d and 3d using template - // specialization. - // - // Note that the first dim components are the velocity components - // and the last is the pressure. - template - class Solution : public Function - { - public: - Solution() - : Function(dim + 1) - {} - virtual double - value(const Point &p, const unsigned int component = 0) const override; - virtual Tensor<1, dim> - gradient(const Point & p, - const unsigned int component = 0) const override; - }; - - template <> - double - Solution<2>::value(const Point<2> &p, const unsigned int component) const - { - Assert(component <= 2 + 1, ExcIndexRange(component, 0, 2 + 1)); - - using numbers::PI; - const double x = p(0); - const double y = p(1); - - if (component == 0) - return sin(PI * x); - if (component == 1) - return -PI * y * cos(PI * x); - if (component == 2) - return sin(PI * x) * cos(PI * y); - - return 0; - } - - template <> - double - Solution<3>::value(const Point<3> &p, const unsigned int component) const - { - Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1)); - - using numbers::PI; - const double x = p(0); - const double y = p(1); - const double z = p(2); - - if (component == 0) - return 2.0 * sin(PI * x); - if (component == 1) - return -PI * y * cos(PI * x); - if (component == 2) - return -PI * z * cos(PI * x); - if (component == 3) - return sin(PI * x) * cos(PI * y) * sin(PI * z); - - return 0; - } - - // Note that for the gradient we need to return a Tensor<1,dim> - template <> - Tensor<1, 2> - Solution<2>::gradient(const Point<2> &p, const unsigned int component) const - { - Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); - - using numbers::PI; - const double x = p(0); - const double y = p(1); - - Tensor<1, 2> return_value; - if (component == 0) - { - return_value[0] = PI * cos(PI * x); - return_value[1] = 0.0; - } - else if (component == 1) - { - return_value[0] = y * PI * PI * sin(PI * x); - return_value[1] = -PI * cos(PI * x); - } - else if (component == 2) - { - return_value[0] = PI * cos(PI * x) * cos(PI * y); - return_value[1] = -PI * sin(PI * x) * sin(PI * y); - } - - return return_value; - } - - template <> - Tensor<1, 3> - Solution<3>::gradient(const Point<3> &p, const unsigned int component) const - { - Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); - - using numbers::PI; - const double x = p(0); - const double y = p(1); - const double z = p(2); - - Tensor<1, 3> return_value; - if (component == 0) - { - return_value[0] = 2 * PI * cos(PI * x); - return_value[1] = 0.0; - return_value[2] = 0.0; - } - else if (component == 1) - { - return_value[0] = y * PI * PI * sin(PI * x); - return_value[1] = -PI * cos(PI * x); - return_value[2] = 0.0; - } - else if (component == 2) - { - return_value[0] = z * PI * PI * sin(PI * x); - return_value[1] = 0.0; - return_value[2] = -PI * cos(PI * x); - } - else if (component == 3) - { - return_value[0] = PI * cos(PI * x) * cos(PI * y) * sin(PI * z); - return_value[1] = -PI * sin(PI * x) * sin(PI * y) * sin(PI * z); - return_value[2] = PI * sin(PI * x) * cos(PI * y) * cos(PI * z); - } - - return return_value; - } - - // Implementation of $f$. See the introduction for more information. - template - class RightHandSide : public Function - { - public: - RightHandSide() - : Function(dim + 1) - {} - - virtual double - value(const Point &p, const unsigned int component = 0) const override; - }; - - template <> - double - RightHandSide<2>::value(const Point<2> &p, const unsigned int component) const - { - Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); - - using numbers::PI; - double x = p(0); - double y = p(1); - if (component == 0) - return PI * PI * sin(PI * x) + PI * cos(PI * x) * cos(PI * y); - if (component == 1) - return -PI * PI * PI * y * cos(PI * x) - PI * sin(PI * y) * sin(PI * x); - if (component == 2) - return 0; - - return 0; - } - - template <> - double - RightHandSide<3>::value(const Point<3> &p, const unsigned int component) const - { - Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); - - using numbers::PI; - double x = p(0); - double y = p(1); - double z = p(2); - if (component == 0) - return 2 * PI * PI * sin(PI * x) + - PI * cos(PI * x) * cos(PI * y) * sin(PI * z); - if (component == 1) - return -PI * PI * PI * y * cos(PI * x) + - PI * (-1) * sin(PI * y) * sin(PI * x) * sin(PI * z); - if (component == 2) - return -PI * PI * PI * z * cos(PI * x) + - PI * cos(PI * z) * sin(PI * x) * cos(PI * y); - if (component == 3) - return 0; - - return 0; - } - - - - // @sect3{ASPECT BlockSchurPreconditioner} - - // In the following, we will implement a preconditioner that expands - // on the ideas discussed in the Results section of step-22. - // Specifically, we - // 1. use an upper block-triangular preconditioner because we want to - // use right preconditioning. - // 2. optionally allow using an inner solver for the velocity block instead - // of a single preconditioner application. - // 3. do not use InverseMatrix but explicitly call SolverCG. - // This approach is also used in the ASPECT code - // (see http://aspect.dealii.org) that solves the Stokes equations in - // the context of simulating convection in the earth mantle, and which - // has been used to solve problems on many thousands of processors. - // - // The bool flag @p do_solve_A in the constructor allows us to either - // apply the preconditioner for the velocity block once or use an inner - // iterative solver for a more accurate approximation instead. - // - // Notice how we keep track of the sum of the inner iterations - // (preconditioner applications). - template - class BlockSchurPreconditioner : public Subscriptor - { - public: - BlockSchurPreconditioner( - const BlockSparseMatrix &system_matrix, - const SparseMatrix & schur_complement_matrix, - const PreconditionerAType & preconditioner_A, - const PreconditionerSType & preconditioner_S, - const bool do_solve_A); - - void - vmult(BlockVector &dst, const BlockVector &src) const; - - mutable unsigned int n_iterations_A; - mutable unsigned int n_iterations_S; - - private: - const BlockSparseMatrix &system_matrix; - const SparseMatrix & schur_complement_matrix; - const PreconditionerAType & preconditioner_A; - const PreconditionerSType & preconditioner_S; - - const bool do_solve_A; - }; - - template - BlockSchurPreconditioner:: - BlockSchurPreconditioner( - const BlockSparseMatrix &system_matrix, - const SparseMatrix & schur_complement_matrix, - const PreconditionerAType & preconditioner_A, - const PreconditionerSType & preconditioner_S, - const bool do_solve_A) - : n_iterations_A(0) - , n_iterations_S(0) - , system_matrix(system_matrix) - , schur_complement_matrix(schur_complement_matrix) - , preconditioner_A(preconditioner_A) - , preconditioner_S(preconditioner_S) - , do_solve_A(do_solve_A) - {} - - - - template - void - BlockSchurPreconditioner::vmult( - BlockVector & dst, - const BlockVector &src) const - { - Vector utmp(src.block(0)); - - // First solve with the approximation for S - { - SolverControl solver_control(1000, 1e-6 * src.block(1).l2_norm()); - SolverCG<> cg(solver_control); - - dst.block(1) = 0.0; - cg.solve(schur_complement_matrix, - dst.block(1), - src.block(1), - preconditioner_S); - - n_iterations_S += solver_control.last_step(); - dst.block(1) *= -1.0; - } - - // Second, apply the top right block (B^T) - { - system_matrix.block(0, 1).vmult(utmp, dst.block(1)); - utmp *= -1.0; - utmp += src.block(0); - } - - // Finally, either solve with the top left block - // or just apply one preconditioner sweep - if (do_solve_A == true) - { - SolverControl solver_control(10000, utmp.l2_norm() * 1e-4); - SolverCG<> cg(solver_control); - - dst.block(0) = 0.0; - cg.solve(system_matrix.block(0, 0), - dst.block(0), - utmp, - preconditioner_A); - - n_iterations_A += solver_control.last_step(); - } - else - { - preconditioner_A.vmult(dst.block(0), utmp); - n_iterations_A += 1; - } - } - - // @sect3{The StokesProblem class} - // - // This is the main class of the problem. - template - class StokesProblem - { - public: - StokesProblem(const unsigned int pressure_degree, - SolverType::type solver_type); - void - run(); - - private: - void - setup_dofs(); - void - assemble_system(); - void - assemble_multigrid(); - void - solve(); - void - compute_errors(); - void - output_results(const unsigned int refinement_cycle) const; - - const unsigned int pressure_degree; - SolverType::type solver_type; - - Triangulation triangulation; - FESystem velocity_fe; - FESystem fe; - DoFHandler dof_handler; - DoFHandler velocity_dof_handler; - - AffineConstraints constraints; - - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - SparseMatrix pressure_mass_matrix; - - BlockVector solution; - BlockVector system_rhs; - - MGLevelObject mg_sparsity_patterns; - MGLevelObject> mg_matrices; - MGLevelObject> mg_interface_matrices; - MGConstrainedDoFs mg_constrained_dofs; - - TimerOutput computing_timer; - }; - - - - template - StokesProblem::StokesProblem(const unsigned int pressure_degree, - SolverType::type solver_type) - : pressure_degree(pressure_degree) - , solver_type(solver_type) - , triangulation(Triangulation::maximum_smoothing) - , - // Finite element for the velocity only: - velocity_fe(FE_Q(pressure_degree + 1), dim) - , - // Finite element for the whole system: - fe(velocity_fe, 1, FE_Q(pressure_degree), 1) - , dof_handler(triangulation) - , velocity_dof_handler(triangulation) - , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times) - {} - - - - // @sect4{StokesProblem::setup_dofs} - - // This function sets up the DoFHandler, matrices, vectors, and Multigrid - // structures (if needed). - template - void - StokesProblem::setup_dofs() - { - TimerOutput::Scope scope(computing_timer, "Setup"); - - system_matrix.clear(); - pressure_mass_matrix.clear(); - - // The main DoFHandler only needs active DoFs, so we are not calling - // distribute_mg_dofs(); - dof_handler.distribute_dofs(fe); - - // This block structure separates the dim velocity components from - // the pressure component (used for reordering). Note that we have - // 2 instead of dim+1 blocks like in step-22, because our FESystem - // is nested and the dim velocity components appear as one block. - std::vector block_component(2); - block_component[0] = 0; - block_component[1] = 1; - - // Velocities start at component 0: - const FEValuesExtractors::Vector velocities(0); - - // ILU behaves better if we apply a reordering to reduce fillin. There - // is no advantage in doing this for the other solvers. - if (solver_type == SolverType::FGMRES_ILU) - { - TimerOutput::Scope ilu_specific(computing_timer, "(ILU specific)"); - DoFRenumbering::Cuthill_McKee(dof_handler); - } - - // This ensures that all velocities DoFs are enumerated before the - // pressure unknowns. This allows us to use blocks for vectors and - // matrices and allows us to get the same DoF numbering for - // dof_handler and velocity_dof_handler. - DoFRenumbering::block_wise(dof_handler); - - if (solver_type == SolverType::FGMRES_GMG) - { - TimerOutput::Scope multigrid_specific(computing_timer, - "(Multigrid specific)"); - TimerOutput::Scope setup_multigrid(computing_timer, - "Setup - Multigrid"); - - // This distributes the active dofs and multigrid dofs for the - // velocity space in a separate DoFHandler as described in the - // introduction. - velocity_dof_handler.distribute_dofs(velocity_fe); - velocity_dof_handler.distribute_mg_dofs(); - - // The following block of code initializes the MGConstrainedDofs - // (using the boundary conditions for the velocity), and the - // sparsity patterns and matrices for each level. The resize() - // function of MGLevelObject will destroy all existing contained - // objects. - const std::set zero_boundary_ids = {0}; - - mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(velocity_dof_handler); - mg_constrained_dofs.make_zero_boundary_constraints(velocity_dof_handler, - zero_boundary_ids); - const unsigned int n_levels = triangulation.n_levels(); - - mg_interface_matrices.resize(0, n_levels - 1); - mg_matrices.resize(0, n_levels - 1); - mg_sparsity_patterns.resize(0, n_levels - 1); - - for (unsigned int level = 0; level < n_levels; ++level) - { - DynamicSparsityPattern csp(velocity_dof_handler.n_dofs(level), - velocity_dof_handler.n_dofs(level)); - MGTools::make_sparsity_pattern(velocity_dof_handler, csp, level); - mg_sparsity_patterns[level].copy_from(csp); - - mg_matrices[level].reinit(mg_sparsity_patterns[level]); - mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]); - } - } - - const std::vector dofs_per_block = - DoFTools::count_dofs_per_fe_block(dof_handler, block_component); - const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; - - { - constraints.clear(); - // The following makes use of a component mask for interpolation of the - // boundary values for the velocity only, which is further explained in - // the vector valued dealii step-20 tutorial. - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - VectorTools::interpolate_boundary_values(dof_handler, - 0, - Solution(), - constraints, - fe.component_mask(velocities)); - - // As discussed in the introduction, we need to fix one degree of freedom - // of the pressure variable to ensure solvability of the problem. We do - // this here by marking the first pressure dof, which has index n_u as a - // constrained dof. - if (solver_type == SolverType::UMFPACK) - constraints.add_line(n_u); - - constraints.close(); - } - - std::cout << "\tNumber of active cells: " << triangulation.n_active_cells() - << std::endl - << "\tNumber of degrees of freedom: " << dof_handler.n_dofs() - << " (" << n_u << '+' << n_p << ')' << std::endl; - - { - BlockDynamicSparsityPattern csp(dofs_per_block, dofs_per_block); - DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false); - sparsity_pattern.copy_from(csp); - } - system_matrix.reinit(sparsity_pattern); - - solution.reinit(dofs_per_block); - system_rhs.reinit(dofs_per_block); - } - - - // @sect4{StokesProblem::assemble_system} - - // In this function, the system matrix is assembled. We assemble the pressure - // mass matrix in the (1,1) block (if needed) and move it out of this location - // at the end of this function. - template - void - StokesProblem::assemble_system() - { - TimerOutput::Scope assemble(computing_timer, "Assemble"); - system_matrix = 0; - system_rhs = 0; - - // If true, we will assemble the pressure mass matrix in the (1,1) block: - const bool assemble_pressure_mass_matrix = - (solver_type == SolverType::UMFPACK) ? false : true; - - QGauss quadrature_formula(pressure_degree + 2); - - FEValues fe_values(fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values | update_gradients); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); - Vector local_rhs(dofs_per_cell); - - std::vector local_dof_indices(dofs_per_cell); - - const RightHandSide right_hand_side; - std::vector> rhs_values(n_q_points, Vector(dim + 1)); - - const FEValuesExtractors::Vector velocities(0); - const FEValuesExtractors::Scalar pressure(dim); - - std::vector> symgrad_phi_u(dofs_per_cell); - std::vector div_phi_u(dofs_per_cell); - std::vector phi_p(dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - { - fe_values.reinit(cell); - local_matrix = 0; - local_rhs = 0; - - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - rhs_values); - - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int k = 0; k < dofs_per_cell; ++k) - { - symgrad_phi_u[k] = - fe_values[velocities].symmetric_gradient(k, q); - div_phi_u[k] = fe_values[velocities].divergence(k, q); - phi_p[k] = fe_values[pressure].value(k, q); - } - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j <= i; ++j) - { - local_matrix(i, j) += - (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) - - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + - (assemble_pressure_mass_matrix ? phi_p[i] * phi_p[j] : - 0)) * - fe_values.JxW(q); - } - - const unsigned int component_i = - fe.system_to_component_index(i).first; - local_rhs(i) += fe_values.shape_value(i, q) * - rhs_values[q](component_i) * fe_values.JxW(q); - } - } - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = i + 1; j < dofs_per_cell; ++j) - local_matrix(i, j) = local_matrix(j, i); - - cell->get_dof_indices(local_dof_indices); - constraints.distribute_local_to_global(local_matrix, - local_rhs, - local_dof_indices, - system_matrix, - system_rhs); - } - - if (solver_type != SolverType::UMFPACK) - { - pressure_mass_matrix.reinit(sparsity_pattern.block(1, 1)); - pressure_mass_matrix.copy_from(system_matrix.block(1, 1)); - system_matrix.block(1, 1) = 0; - } - } - - // @sect4{StokesProblem::assemble_multigrid} - - // Here, like in step-16, we have a function that assembles the level - // and interface matrices necessary for the multigrid preconditioner. - template - void - StokesProblem::assemble_multigrid() - { - TimerOutput::Scope multigrid_specific(computing_timer, - "(Multigrid specific)"); - TimerOutput::Scope assemble_multigrid(computing_timer, - "Assemble Multigrid"); - - mg_matrices = 0.; - - QGauss quadrature_formula(pressure_degree + 2); - - FEValues fe_values(velocity_fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values | update_gradients); - - const unsigned int dofs_per_cell = velocity_fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices(dofs_per_cell); - - const FEValuesExtractors::Vector velocities(0); - - std::vector> symgrad_phi_u(dofs_per_cell); - - std::vector> boundary_constraints( - triangulation.n_levels()); - std::vector> boundary_interface_constraints( - triangulation.n_levels()); - for (unsigned int level = 0; level < triangulation.n_levels(); ++level) - { - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_refinement_edge_indices(level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_boundary_indices(level)); - boundary_constraints[level].close(); - - IndexSet idx = mg_constrained_dofs.get_refinement_edge_indices(level) & - mg_constrained_dofs.get_boundary_indices(level); - - boundary_interface_constraints[level].add_lines(idx); - boundary_interface_constraints[level].close(); - } - - // This iterator goes over all cells (not just active) - typename DoFHandler::cell_iterator cell = velocity_dof_handler.begin(), - endc = velocity_dof_handler.end(); - - for (; cell != endc; ++cell) - { - fe_values.reinit(cell); - cell_matrix = 0; - - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int k = 0; k < dofs_per_cell; ++k) - symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient(k, q); - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j <= i; ++j) - { - cell_matrix(i, j) += - (symgrad_phi_u[i] * symgrad_phi_u[j]) * fe_values.JxW(q); - } - } - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = i + 1; j < dofs_per_cell; ++j) - cell_matrix(i, j) = cell_matrix(j, i); - - cell->get_mg_dof_indices(local_dof_indices); - - boundary_constraints[cell->level()].distribute_local_to_global( - cell_matrix, local_dof_indices, mg_matrices[cell->level()]); - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - if (!mg_constrained_dofs.at_refinement_edge(cell->level(), - local_dof_indices[i]) || - mg_constrained_dofs.at_refinement_edge(cell->level(), - local_dof_indices[j])) - cell_matrix(i, j) = 0; - - boundary_interface_constraints[cell->level()] - .distribute_local_to_global(cell_matrix, - local_dof_indices, - mg_interface_matrices[cell->level()]); - } - } - - // @sect4{StokesProblem::solve} - - // This function sets up things differently based on if you want to use ILU - // or GMG as a preconditioner. Both methods share the same solver (FGMRES) - // but require a different preconditioner to be initialized. Here we time not - // only the entire solve function, but we separately time the setup of the - // preconditioner as well as the solve itself. - template - void - StokesProblem::solve() - { - TimerOutput::Scope solve(computing_timer, "Solve"); - constraints.set_zero(solution); - - if (solver_type == SolverType::UMFPACK) - { - computing_timer.enter_subsection("(UMFPACK specific)"); - computing_timer.enter_subsection("Solve - Initialize"); - - SparseDirectUMFPACK A_direct; - A_direct.initialize(system_matrix); - - computing_timer.leave_subsection(); - computing_timer.leave_subsection(); - - { - TimerOutput::Scope solve_backslash(computing_timer, - "Solve - Backslash"); - A_direct.vmult(solution, system_rhs); - } - - constraints.distribute(solution); - return; - } - - // Here we must make sure to solve for the residual with "good enough" - // accuracy - SolverControl solver_control(system_matrix.m(), - 1e-10 * system_rhs.l2_norm()); - unsigned int n_iterations_A; - unsigned int n_iterations_S; - - // This is used to pass whether or not we want to solve for A inside - // the preconditioner. One could change this to false to see if - // there is still convergence and if so does the program then run - // faster or slower - const bool use_expensive = true; - - SolverFGMRES> solver(solver_control); - - if (solver_type == SolverType::FGMRES_ILU) - { - computing_timer.enter_subsection("(ILU specific)"); - computing_timer.enter_subsection("Solve - Set-up Preconditioner"); - - std::cout << " Computing preconditioner..." << std::endl - << std::flush; - - SparseILU A_preconditioner; - A_preconditioner.initialize(system_matrix.block(0, 0)); - - SparseILU S_preconditioner; - S_preconditioner.initialize(pressure_mass_matrix); - - const BlockSchurPreconditioner, SparseILU> - preconditioner(system_matrix, - pressure_mass_matrix, - A_preconditioner, - S_preconditioner, - use_expensive); - - computing_timer.leave_subsection(); - computing_timer.leave_subsection(); - - { - TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES"); - - solver.solve(system_matrix, solution, system_rhs, preconditioner); - n_iterations_A = preconditioner.n_iterations_A; - n_iterations_S = preconditioner.n_iterations_S; - } - } - else - { - computing_timer.enter_subsection("(Multigrid specific)"); - computing_timer.enter_subsection("Solve - Set-up Preconditioner"); - - // Transfer operators between levels - MGTransferPrebuilt> mg_transfer(mg_constrained_dofs); - mg_transfer.build(velocity_dof_handler); - - // Setup coarse grid solver - FullMatrix coarse_matrix; - coarse_matrix.copy_from(mg_matrices[0]); - MGCoarseGridHouseholder<> coarse_grid_solver; - coarse_grid_solver.initialize(coarse_matrix); - - using Smoother = PreconditionSOR>; - mg::SmootherRelaxation> mg_smoother; - mg_smoother.initialize(mg_matrices); - mg_smoother.set_steps(2); - - // Multigrid, when used as a preconditioner for CG, needs to be a - // symmetric operator, so the smoother must be symmetric - mg_smoother.set_symmetric(true); - - mg::Matrix> mg_matrix(mg_matrices); - mg::Matrix> mg_interface_up(mg_interface_matrices); - mg::Matrix> mg_interface_down(mg_interface_matrices); - - // Now, we are ready to set up the V-cycle operator and the multilevel - // preconditioner. - Multigrid> mg( - mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother); - mg.set_edge_matrices(mg_interface_down, mg_interface_up); - - PreconditionMG, MGTransferPrebuilt>> - A_Multigrid(velocity_dof_handler, mg, mg_transfer); - - SparseILU S_preconditioner; - S_preconditioner.initialize(pressure_mass_matrix, - SparseILU::AdditionalData()); - - const BlockSchurPreconditioner< - PreconditionMG, - MGTransferPrebuilt>>, - SparseILU> - preconditioner(system_matrix, - pressure_mass_matrix, - A_Multigrid, - S_preconditioner, - use_expensive); - - computing_timer.leave_subsection(); - computing_timer.leave_subsection(); - - { - TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES"); - solver.solve(system_matrix, solution, system_rhs, preconditioner); - n_iterations_A = preconditioner.n_iterations_A; - n_iterations_S = preconditioner.n_iterations_S; - } - } - - constraints.distribute(solution); - - std::cout - << std::endl - << "\tNumber of FGMRES iterations: " << solver_control.last_step() - << std::endl - << "\tTotal number of iterations used for approximation of A inverse: " - << n_iterations_A << std::endl - << "\tTotal number of iterations used for approximation of S inverse: " - << n_iterations_S << std::endl - << std::endl; - } - - - // @sect4{StokesProblem::process_solution} - - // This function computes the L2 and H1 errors of the solution. For this, - // we need to make sure the pressure has mean zero. - template - void - StokesProblem::compute_errors() - { - // Compute the mean pressure $\frac{1}{\Omega} \int_{\Omega} p(x) dx $ - // and then subtract it from each pressure coefficient. This will result - // in a pressure with mean value zero. Here we make use of the fact that - // the pressure is component $dim$ and that the finite element space - // is nodal. - const double mean_pressure = VectorTools::compute_mean_value( - dof_handler, QGauss(pressure_degree + 2), solution, dim); - solution.block(1).add(-mean_pressure); - std::cout << " Note: The mean value was adjusted by " << -mean_pressure - << std::endl; - - const ComponentSelectFunction pressure_mask(dim, dim + 1); - const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), - dim + 1); - - Vector difference_per_cell(triangulation.n_active_cells()); - VectorTools::integrate_difference(dof_handler, - solution, - Solution(), - difference_per_cell, - QGauss(pressure_degree + 2), - VectorTools::L2_norm, - &velocity_mask); - - const double Velocity_L2_error = - VectorTools::compute_global_error(triangulation, - difference_per_cell, - VectorTools::L2_norm); - - VectorTools::integrate_difference(dof_handler, - solution, - Solution(), - difference_per_cell, - QGauss(pressure_degree + 2), - VectorTools::L2_norm, - &pressure_mask); - - const double Pressure_L2_error = - VectorTools::compute_global_error(triangulation, - difference_per_cell, - VectorTools::L2_norm); - - VectorTools::integrate_difference(dof_handler, - solution, - Solution(), - difference_per_cell, - QGauss(pressure_degree + 2), - VectorTools::H1_norm, - &velocity_mask); - - const double Velocity_H1_error = - VectorTools::compute_global_error(triangulation, - difference_per_cell, - VectorTools::H1_norm); - - std::cout << std::endl - << " Velocity L2 Error: " << Velocity_L2_error << std::endl - << " Pressure L2 Error: " << Pressure_L2_error << std::endl - << " Velocity H1 Error: " << Velocity_H1_error << std::endl; - } - - - // @sect4{StokesProblem::output_results} - - // This function generates graphical output like it is done in step-22. - template - void - StokesProblem::output_results(const unsigned int refinement_cycle) const - { - std::vector solution_names(dim, "velocity"); - solution_names.emplace_back("pressure"); - - std::vector - data_component_interpretation( - dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation.push_back( - DataComponentInterpretation::component_is_scalar); - - DataOut data_out; - data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, - solution_names, - DataOut::type_dof_data, - data_component_interpretation); - data_out.build_patches(); - - std::ofstream output( - "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk"); - data_out.write_vtk(output); - } - - - - // @sect4{StokesProblem::run} - - // The last step in the Stokes class is, as usual, the function that - // generates the initial grid and calls the other functions in the - // respective order. - template - void - StokesProblem::run() - { - GridGenerator::hyper_cube(triangulation); - triangulation.refine_global(2); - - if (solver_type == SolverType::FGMRES_ILU) - std::cout << "Now running with ILU" << std::endl; - else if (solver_type == SolverType::FGMRES_GMG) - std::cout << "Now running with Multigrid" << std::endl; - else - std::cout << "Now running with UMFPACK" << std::endl; - - - for (unsigned int refinement_cycle = 0; refinement_cycle < 2; - ++refinement_cycle) - { - std::cout << "Refinement cycle " << refinement_cycle << std::endl; - - if (refinement_cycle > 0) - triangulation.refine_global(1); - - std::cout << " Set-up..." << std::endl; - setup_dofs(); - - std::cout << " Assembling..." << std::endl; - assemble_system(); - - if (solver_type == SolverType::FGMRES_GMG) - { - std::cout << " Assembling Multigrid..." << std::endl; - - assemble_multigrid(); - } - - std::cout << " Solving..." << std::flush; - solve(); - - compute_errors(); - - output_results(refinement_cycle); - - computing_timer.reset(); - } - } -} // namespace Step56 - -// @sect3{The main function} -int -main() -{ - try - { - using namespace Step56; - - const int degree = 1; - const int dim = 3; - // options for SolverType: UMFPACK FGMRES_ILU FGMRES_GMG - { - StokesProblem flow_problem(degree, SolverType::FGMRES_GMG); - flow_problem.run(); - } - { - StokesProblem flow_problem(degree, SolverType::FGMRES_ILU); - flow_problem.run(); - } - } - catch (const 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/examples/step-56.diff b/tests/examples/step-56.diff new file mode 100644 index 0000000000..cb1bf1dbc7 --- /dev/null +++ b/tests/examples/step-56.diff @@ -0,0 +1,29 @@ +1057c1057 +< triangulation.refine_global(6 - dim); +--- +> triangulation.refine_global(2); +1097c1097 +< std::cout << " VM Peak: " << mem.VmPeak << std::endl; +--- +> // std::cout << " VM Peak: " << mem.VmPeak << std::endl; +1099c1099 +< computing_timer.print_summary(); +--- +> // computing_timer.print_summary(); +1113c1113 +< const int dim = 3; +--- +> const int dim = 2; +1115,1117c1115,1122 +< StokesProblem flow_problem(degree, SolverType::FGMRES_GMG); +< +< flow_problem.run(); +--- +> { +> StokesProblem flow_problem(degree, SolverType::FGMRES_GMG); +> flow_problem.run(); +> } +> { +> StokesProblem flow_problem(degree, SolverType::FGMRES_ILU); +> flow_problem.run(); +> } diff --git a/tests/examples/step-56.output b/tests/examples/step-56.output index b4f4b729b1..5edf0cac08 100644 --- a/tests/examples/step-56.output +++ b/tests/examples/step-56.output @@ -1,66 +1,98 @@ Now running with Multigrid Refinement cycle 0 Set-up... - Number of active cells: 64 - Number of degrees of freedom: 2312 (2187+125) + Number of active cells: 16 + Number of degrees of freedom: 187 (162+25) Assembling... Assembling Multigrid... Solving... - Number of FGMRES iterations: 22 - Total number of iterations used for approximation of A inverse: 69 - Total number of iterations used for approximation of S inverse: 23 + Number of FGMRES iterations: 14 + Total number of iterations used for approximation of A inverse: 45 + Total number of iterations used for approximation of S inverse: 15 - Note: The mean value was adjusted by -3.36485e-15 + Note: The mean value was adjusted by 8.38229e-15 - Velocity L2 Error: 0.00539774 - Pressure L2 Error: 0.019341 - Velocity H1 Error: 0.165265 + Velocity L2 Error: 0.00340778 + Pressure L2 Error: 0.017624 + Velocity H1 Error: 0.105243 Refinement cycle 1 Set-up... - Number of active cells: 512 - Number of degrees of freedom: 15468 (14739+729) + Number of active cells: 64 + Number of degrees of freedom: 659 (578+81) + Assembling... + Assembling Multigrid... + Solving... + Number of FGMRES iterations: 16 + Total number of iterations used for approximation of A inverse: 52 + Total number of iterations used for approximation of S inverse: 17 + + Note: The mean value was adjusted by -1.75663e-14 + + Velocity L2 Error: 0.000426317 + Pressure L2 Error: 0.00414007 + Velocity H1 Error: 0.0264094 +Refinement cycle 2 + Set-up... + Number of active cells: 256 + Number of degrees of freedom: 2467 (2178+289) Assembling... Assembling Multigrid... Solving... - Number of FGMRES iterations: 21 - Total number of iterations used for approximation of A inverse: 67 - Total number of iterations used for approximation of S inverse: 22 + Number of FGMRES iterations: 16 + Total number of iterations used for approximation of A inverse: 52 + Total number of iterations used for approximation of S inverse: 17 - Note: The mean value was adjusted by -1.51964e-13 + Note: The mean value was adjusted by -2.40243e-14 - Velocity L2 Error: 0.000670888 - Pressure L2 Error: 0.0036533 - Velocity H1 Error: 0.0414704 + Velocity L2 Error: 5.33246e-05 + Pressure L2 Error: 0.00102063 + Velocity H1 Error: 0.00660864 Now running with ILU Refinement cycle 0 Set-up... - Number of active cells: 64 - Number of degrees of freedom: 2312 (2187+125) + Number of active cells: 16 + Number of degrees of freedom: 187 (162+25) Assembling... Solving... Computing preconditioner... - Number of FGMRES iterations: 22 - Total number of iterations used for approximation of A inverse: 113 - Total number of iterations used for approximation of S inverse: 23 + Number of FGMRES iterations: 14 + Total number of iterations used for approximation of A inverse: 75 + Total number of iterations used for approximation of S inverse: 15 - Note: The mean value was adjusted by 1.55695e-14 + Note: The mean value was adjusted by 3.10407e-15 - Velocity L2 Error: 0.00539774 - Pressure L2 Error: 0.019341 - Velocity H1 Error: 0.165265 + Velocity L2 Error: 0.00340778 + Pressure L2 Error: 0.017624 + Velocity H1 Error: 0.105243 Refinement cycle 1 Set-up... - Number of active cells: 512 - Number of degrees of freedom: 15468 (14739+729) + Number of active cells: 64 + Number of degrees of freedom: 659 (578+81) + Assembling... + Solving... Computing preconditioner... + + Number of FGMRES iterations: 16 + Total number of iterations used for approximation of A inverse: 134 + Total number of iterations used for approximation of S inverse: 17 + + Note: The mean value was adjusted by -3.36177e-15 + + Velocity L2 Error: 0.000426317 + Pressure L2 Error: 0.00414007 + Velocity H1 Error: 0.0264094 +Refinement cycle 2 + Set-up... + Number of active cells: 256 + Number of degrees of freedom: 2467 (2178+289) Assembling... Solving... Computing preconditioner... - Number of FGMRES iterations: 21 - Total number of iterations used for approximation of A inverse: 180 - Total number of iterations used for approximation of S inverse: 22 + Number of FGMRES iterations: 16 + Total number of iterations used for approximation of A inverse: 229 + Total number of iterations used for approximation of S inverse: 17 - Note: The mean value was adjusted by -1.52194e-13 + Note: The mean value was adjusted by -1.16485e-15 - Velocity L2 Error: 0.000670888 - Pressure L2 Error: 0.0036533 - Velocity H1 Error: 0.0414704 + Velocity L2 Error: 5.33247e-05 + Pressure L2 Error: 0.00102063 + Velocity H1 Error: 0.00660864 diff --git a/tests/examples/update.sh b/tests/examples/update.sh new file mode 100755 index 0000000000..cede5e4bd4 --- /dev/null +++ b/tests/examples/update.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +# use .diff files to create .cc files + +# Find all *.diff files in the current directory +diff_files=$(find . -maxdepth 1 -type f -name "*.diff") + +for file in $diff_files; do + # Extract the filename without the extension + filename=$(basename "$file" .diff) + src="../../examples/$filename/$filename.cc" + + # Check if the corresponding .cc file exists + if [[ -f "$src" ]]; then + # Apply the diff and save it to a new file + patch "$src" < "$file" -o "$filename.cc" + else + echo "No matching .cc file found for $diff_file" + exit -1 + fi +done -- 2.39.5