From 667c59f573731fc8b55ae935d56641e57f3012c6 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 9 Sep 2024 16:12:45 -0400 Subject: [PATCH] add Scott-Vogelius test --- tests/simplex/stokes-sv.cc | 726 +++++++++++++++++++++++++++++++++ tests/simplex/stokes-sv.output | 31 ++ 2 files changed, 757 insertions(+) create mode 100644 tests/simplex/stokes-sv.cc create mode 100644 tests/simplex/stokes-sv.output diff --git a/tests/simplex/stokes-sv.cc b/tests/simplex/stokes-sv.cc new file mode 100644 index 0000000000..7e2c8e8252 --- /dev/null +++ b/tests/simplex/stokes-sv.cc @@ -0,0 +1,726 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2021 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +// Stokes on a simplex mesh using barycenter refinement to show +// that Scott-Vogelius elements are pointwise divergence free. + +#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 "deal.II/numerics/vector_tools_mean_value.h" +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + +// #define HEX + + +namespace Step56 +{ + using namespace dealii; + +#ifdef HEX + template + using QuadratureT = QGauss; +#else + template + using QuadratureT = QGaussSimplex; +#endif + + 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; + } + + 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; + } + + 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; + } + + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_matrix, + const SparseMatrix &schur_complement_matrix, + const PreconditionerAType &preconditioner_A, + const PreconditionerSType &preconditioner_S); + + 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; + }; + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_matrix, + const SparseMatrix &schur_complement_matrix, + const PreconditionerAType &preconditioner_A, + const PreconditionerSType &preconditioner_S) + : 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) + {} + + + + template + void + BlockSchurPreconditioner::vmult( + BlockVector &dst, + const BlockVector &src) const + { + Vector utmp(src.block(0)); + + { + n_iterations_S += 1; + preconditioner_S.vmult(dst.block(1), src.block(1)); + dst.block(1) *= -1.0; + } + + { + system_matrix.block(0, 1).vmult(utmp, dst.block(1)); + utmp *= -1.0; + utmp += src.block(0); + } + + { + preconditioner_A.vmult(dst.block(0), utmp); + n_iterations_A += 1; + } + } + + template + class StokesProblem + { + public: + StokesProblem(const unsigned int pressure_degree); + void + run(); + + private: + void + setup_dofs(); + void + assemble_system(); + void + solve(); + void + compute_errors(); + void + output_results(const unsigned int refinement_cycle) const; + + const unsigned int pressure_degree; + + Triangulation triangulation; + +#ifdef HEX + MappingQ1 mapping; +#else + MappingFE mapping; +#endif + + 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; + }; + + + + template + StokesProblem::StokesProblem(const unsigned int pressure_degree) + + : pressure_degree(pressure_degree) +#ifdef HEX + , velocity_fe(FE_Q(pressure_degree + 1), dim) + , fe(velocity_fe, 1, FE_Q(pressure_degree), 1) +#else + , mapping(FE_SimplexP(1)) + , velocity_fe(FE_SimplexP(pressure_degree + 1), dim) + , fe(velocity_fe, 1, FE_SimplexDGP(pressure_degree), 1) +#endif + , dof_handler(triangulation) + , velocity_dof_handler(triangulation) + {} + + + template + void + StokesProblem::setup_dofs() + { + system_matrix.clear(); + pressure_mass_matrix.clear(); + + dof_handler.distribute_dofs(fe); + + std::vector block_component(2); + block_component[0] = 0; + block_component[1] = 1; + + const FEValuesExtractors::Vector velocities(0); + + DoFRenumbering::block_wise(dof_handler); + + 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]; + const unsigned int n_p = dofs_per_block[1]; + + { + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(mapping, + dof_handler, + 0, + Solution(), + constraints, + fe.component_mask(velocities)); + + constraints.close(); + } + + deallog << "\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); + } + + template + void + StokesProblem::assemble_system() + { + system_matrix = 0; + system_rhs = 0; + + const bool assemble_pressure_mass_matrix = true; + + const QuadratureT quadrature_formula(pressure_degree + 2); + + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + const unsigned int dofs_per_cell = fe.n_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); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + 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); + } + + { + 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; + } + } + + template + void + StokesProblem::solve() + { + constraints.set_zero(solution); + + SolverControl solver_control(10000, + 1e-10 * system_rhs.l2_norm(), + false, + false); + unsigned int n_iterations_A; + unsigned int n_iterations_S; + + SolverGMRES> solver( + solver_control, + SolverGMRES>::AdditionalData( + 50, true)); // right preconditioning + + { + 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); + + { + 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); + } + + template + void + StokesProblem::compute_errors() + { + const double mean_pressure = + VectorTools::compute_mean_value(mapping, + dof_handler, + QuadratureT(pressure_degree + 2), + solution, + dim); + VectorTools::add_constant(solution, dof_handler, dim, -mean_pressure); + + 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(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + QuadratureT(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(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + QuadratureT(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(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + QuadratureT(pressure_degree + 2), + VectorTools::H1_norm, + &velocity_mask); + + const double Velocity_H1_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_norm); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, + QuadratureT(pressure_degree + 2), + VectorTools::Hdiv_seminorm, + &velocity_mask); + + const double Velocity_Hdiv_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::Hdiv_seminorm); + deallog << 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 + << " Velocity Hdiv Err: " << Velocity_Hdiv_error << std::endl; + } + + 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); + } + + template + void + StokesProblem::run() + { + for (unsigned int refinement_cycle = 0; refinement_cycle < 3; + ++refinement_cycle) + { + deallog << "Refinement cycle " << refinement_cycle << std::endl; + + Triangulation s_tria; + GridGenerator::subdivided_hyper_cube_with_simplices( + s_tria, std::pow(1 + refinement_cycle, 2)); + triangulation.clear(); + GridGenerator::alfeld_split_of_simplex_mesh(s_tria, triangulation); + + deallog << " Set-up..." << std::endl; + setup_dofs(); + + deallog << " Assembling..." << std::endl; + assemble_system(); + + deallog << " Solving..." << std::flush; + solve(); + + compute_errors(); + + output_results(refinement_cycle); + } + } +} // namespace Step56 + +int +main() +{ + initlog(); + { + using namespace Step56; + + const int degree = 1; + const int dim = 2; + StokesProblem flow_problem(degree); + + flow_problem.run(); + } + return 0; +} diff --git a/tests/simplex/stokes-sv.output b/tests/simplex/stokes-sv.output new file mode 100644 index 0000000000..72f547a5ae --- /dev/null +++ b/tests/simplex/stokes-sv.output @@ -0,0 +1,31 @@ + +DEAL::Refinement cycle 0 +DEAL:: Set-up... +DEAL:: Number of active cells: 6 +DEAL:: Number of degrees of freedom: 52 (34+18) +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: Velocity L2 Error: 0.314683 +DEAL:: Pressure L2 Error: 3.40767 +DEAL:: Velocity H1 Error: 2.13131 +DEAL:: Velocity Hdiv Err: 3.33043e-09 +DEAL::Refinement cycle 1 +DEAL:: Set-up... +DEAL:: Number of active cells: 96 +DEAL:: Number of degrees of freedom: 706 (418+288) +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: Velocity L2 Error: 0.00815750 +DEAL:: Pressure L2 Error: 0.782868 +DEAL:: Velocity H1 Error: 0.275850 +DEAL:: Velocity Hdiv Err: 1.64529e-08 +DEAL::Refinement cycle 2 +DEAL:: Set-up... +DEAL:: Number of active cells: 486 +DEAL:: Number of degrees of freedom: 3476 (2018+1458) +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: Velocity L2 Error: 0.000739254 +DEAL:: Pressure L2 Error: 0.190631 +DEAL:: Velocity H1 Error: 0.0599330 +DEAL:: Velocity Hdiv Err: 6.61221e-09 -- 2.39.5