From: Laura Prieto Saavedra Date: Thu, 14 Jan 2021 23:33:10 +0000 (-0500) Subject: Convert step-20 to a simplex test X-Git-Tag: v9.3.0-rc1~624^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11545%2Fhead;p=dealii.git Convert step-20 to a simplex test --- diff --git a/tests/simplex/step-20.cc b/tests/simplex/step-20.cc new file mode 100644 index 0000000000..34e482fe7c --- /dev/null +++ b/tests/simplex/step-20.cc @@ -0,0 +1,549 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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. +// +// --------------------------------------------------------------------- + + +// Step-20 on a simplex mesh. Following modifications had to be made: +// - Create a hypercube with simplices +// - Change the FE_RT elements to Simplex::FE_P (2nd degree) and FE_DGQ to +// Simplex::FE_DGP (1st degree) +// These spaces are an alternative to solve the mixed Laplacian +// (https://link.springer.com/article/10.1007/s10092-009-0009-6) +// - Change QGauss to Simplex::QGauss and use MappingFE instead of default +// mapping. + + +#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 "../tests.h" + +namespace Step20 +{ + using namespace dealii; + + template + class MixedLaplaceProblem + { + public: + MixedLaplaceProblem(const unsigned int degree); + void + run(); + + private: + void + make_grid_and_dofs(); + void + assemble_system(); + void + solve(); + void + compute_errors() const; + void + output_results() const; + + const unsigned int degree; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + }; + + namespace PrescribedSolution + { + constexpr double alpha = 0.3; + constexpr double beta = 1; + + + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function(1) + {} + + virtual double + value(const Point & p, + const unsigned int component = 0) const override; + }; + + + + template + class PressureBoundaryValues : public Function + { + public: + PressureBoundaryValues() + : Function(1) + {} + + virtual double + value(const Point & p, + const unsigned int component = 0) const override; + }; + + + template + class ExactSolution : public Function + { + public: + ExactSolution() + : Function(dim + 1) + {} + + virtual void + vector_value(const Point &p, Vector &value) const override; + }; + + + template + double + RightHandSide::value(const Point & /*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + + template + double + PressureBoundaryValues::value(const Point &p, + const unsigned int /*component*/) const + { + return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] - + alpha * p[0] * p[0] * p[0] / 6); + } + + + + template + void + ExactSolution::vector_value(const Point &p, + Vector & values) const + { + Assert(values.size() == dim + 1, + ExcDimensionMismatch(values.size(), dim + 1)); + + values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2; + values(1) = alpha * p[0] * p[1]; + values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] - + alpha * p[0] * p[0] * p[0] / 6); + } + + template + class KInverse : public TensorFunction<2, dim> + { + public: + KInverse() + : TensorFunction<2, dim>() + {} + + virtual void + value_list(const std::vector> &points, + std::vector> & values) const override; + }; + + template + void + KInverse::value_list(const std::vector> &points, + std::vector> & values) const + { + (void)points; + AssertDimension(points.size(), values.size()); + + for (auto &value : values) + value = unit_symmetric_tensor(); + } + } // namespace PrescribedSolution + + template + MixedLaplaceProblem::MixedLaplaceProblem(const unsigned int degree) + : degree(degree) + , fe(FESystem(Simplex::FE_P(degree), dim), + 1, + Simplex::FE_DGP(degree - 1), + 1) + , dof_handler(triangulation) + {} + + template + void + MixedLaplaceProblem::make_grid_and_dofs() + { + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + 2, + -1, + 1); + triangulation.refine_global(3); + + dof_handler.distribute_dofs(fe); + + DoFRenumbering::component_wise(dof_handler); + + const std::vector dofs_per_component = + DoFTools::count_dofs_per_fe_component(dof_handler); + const unsigned int n_u = dofs_per_component[0] + dofs_per_component[1], + n_p = dofs_per_component[dim]; + + deallog << "Number of active cells: " << triangulation.n_active_cells() + << std::endl + << "Total number of cells: " << triangulation.n_cells() << std::endl + << "Number of degrees of freedom: " << dof_handler.n_dofs() << " (" + << n_u << '+' << n_p << ')' << std::endl; + + BlockDynamicSparsityPattern dsp(2, 2); + dsp.block(0, 0).reinit(n_u, n_u); + dsp.block(1, 0).reinit(n_p, n_u); + dsp.block(0, 1).reinit(n_u, n_p); + dsp.block(1, 1).reinit(n_p, n_p); + dsp.collect_sizes(); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); + + solution.reinit(2); + solution.block(0).reinit(n_u); + solution.block(1).reinit(n_p); + solution.collect_sizes(); + + system_rhs.reinit(2); + system_rhs.block(0).reinit(n_u); + system_rhs.block(1).reinit(n_p); + system_rhs.collect_sizes(); + } + + template + void + MixedLaplaceProblem::assemble_system() + { + Simplex::QGauss quadrature_formula(degree + 1); + Simplex::QGauss face_quadrature_formula(degree + 1); + MappingFE mapping(Simplex::FE_DGP(1)); + + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEFaceValues fe_face_values(mapping, + fe, + face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_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 PrescribedSolution::RightHandSide right_hand_side; + const PrescribedSolution::PressureBoundaryValues + pressure_boundary_values; + const PrescribedSolution::KInverse k_inverse; + + std::vector rhs_values(n_q_points); + std::vector boundary_values(n_face_q_points); + std::vector> k_inverse_values(n_q_points); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.value_list(fe_values.get_quadrature_points(), + rhs_values); + k_inverse.value_list(fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q); + const double div_phi_i_u = fe_values[velocities].divergence(i, q); + const double phi_i_p = fe_values[pressure].value(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j_u = + fe_values[velocities].value(j, q); + const double div_phi_j_u = + fe_values[velocities].divergence(j, q); + const double phi_j_p = fe_values[pressure].value(j, q); + + local_matrix(i, j) += + (phi_i_u * k_inverse_values[q] * phi_j_u // + - phi_i_p * div_phi_j_u // + - div_phi_i_u * phi_j_p) // + * fe_values.JxW(q); + } + + local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q); + } + + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + fe_face_values.reinit(cell, face); + + pressure_boundary_values.value_list( + fe_face_values.get_quadrature_points(), boundary_values); + + for (unsigned int q = 0; q < n_face_q_points; ++q) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + local_rhs(i) += -(fe_face_values[velocities].value(i, q) * // + fe_face_values.normal_vector(q) * // + boundary_values[q] * // + fe_face_values.JxW(q)); + } + + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + local_matrix(i, j)); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + system_rhs(local_dof_indices[i]) += local_rhs(i); + } + } + + template + void + MixedLaplaceProblem::solve() + { + const auto &M = system_matrix.block(0, 0); + const auto &B = system_matrix.block(0, 1); + + const auto &F = system_rhs.block(0); + const auto &G = system_rhs.block(1); + + auto &U = solution.block(0); + auto &P = solution.block(1); + + const auto op_M = linear_operator(M); + const auto op_B = linear_operator(B); + + ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); + SolverCG> solver_M(reduction_control_M); + PreconditionJacobi> preconditioner_M; + + preconditioner_M.initialize(M); + + const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M); + + const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; + const auto op_aS = + transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B; + + IterationNumberControl iteration_number_control_aS(30, 1.e-18); + SolverCG> solver_aS(iteration_number_control_aS); + + const auto preconditioner_S = + inverse_operator(op_aS, solver_aS, PreconditionIdentity()); + + const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G; + + SolverControl solver_control_S(2000, 1.e-12); + SolverCG> solver_S(solver_control_S); + + const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S); + + P = op_S_inv * schur_rhs; + + deallog << solver_control_S.last_step() + << " CG Schur complement iterations to obtain convergence." + << std::endl; + + U = op_M_inv * (F - op_B * P); + } + + template + void + MixedLaplaceProblem::compute_errors() const + { + const ComponentSelectFunction pressure_mask(dim, dim + 1); + const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), + dim + 1); + + PrescribedSolution::ExactSolution exact_solution; + Vector cellwise_errors(triangulation.n_active_cells()); + + QTrapezoid<1> q_trapez; + QIterated quadrature(q_trapez, degree + 2); + MappingFE mapping(Simplex::FE_DGP(1)); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + exact_solution, + cellwise_errors, + quadrature, + VectorTools::L2_norm, + &pressure_mask); + const double p_l2_error = + VectorTools::compute_global_error(triangulation, + cellwise_errors, + VectorTools::L2_norm); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + exact_solution, + cellwise_errors, + quadrature, + VectorTools::L2_norm, + &velocity_mask); + const double u_l2_error = + VectorTools::compute_global_error(triangulation, + cellwise_errors, + VectorTools::L2_norm); + + deallog << "Errors: ||e_p||_L2 = " << p_l2_error + << ", ||e_u||_L2 = " << u_l2_error << std::endl; + } + + template + void + MixedLaplaceProblem::output_results() const + { + MappingFE mapping(Simplex::FE_DGP(1)); + std::vector solution_names(dim, "u"); + solution_names.emplace_back("p"); + std::vector + interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + interpretation.push_back(DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.add_data_vector(dof_handler, + solution, + solution_names, + interpretation); + + data_out.build_patches(mapping, degree); + + std::ofstream output("solution.vtu"); + data_out.write_vtu(output); + } + + + template + void + MixedLaplaceProblem::run() + { + make_grid_and_dofs(); + assemble_system(); + solve(); + compute_errors(); + output_results(); + } +} // namespace Step20 + +int +main() +{ + initlog(); + + try + { + using namespace Step20; + + const unsigned int fe_degree = 2; + MixedLaplaceProblem<2> mixed_laplace_problem(fe_degree); + mixed_laplace_problem.run(); + } + 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/simplex/step-20.with_simplex_support=on.output b/tests/simplex/step-20.with_simplex_support=on.output new file mode 100644 index 0000000000..157e9c1463 --- /dev/null +++ b/tests/simplex/step-20.with_simplex_support=on.output @@ -0,0 +1,144 @@ + +DEAL::Number of active cells: 512 +DEAL::Total number of cells: 680 +DEAL::Number of degrees of freedom: 3714 (2178+1536) +DEAL:cg::Starting value 0.624584 +DEAL:cg::Convergence step 25 value 4.62032e-11 +DEAL:cg::Starting value 37.4651 +DEAL:cg:cg::Starting value 37.4651 +DEAL:cg:cg::Convergence step 30 value 0.0848672 +DEAL:cg:cg::Starting value 0.826345 +DEAL:cg:cg::Convergence step 25 value 5.36874e-11 +DEAL:cg:cg::Starting value 8.99626 +DEAL:cg:cg::Convergence step 30 value 0.108473 +DEAL:cg:cg::Starting value 0.257746 +DEAL:cg:cg::Convergence step 25 value 1.58401e-11 +DEAL:cg:cg::Starting value 3.18381 +DEAL:cg:cg::Convergence step 30 value 0.0283153 +DEAL:cg:cg::Starting value 0.0877011 +DEAL:cg:cg::Convergence step 25 value 6.78488e-12 +DEAL:cg:cg::Starting value 1.08833 +DEAL:cg:cg::Convergence step 30 value 0.0265170 +DEAL:cg:cg::Starting value 0.0390461 +DEAL:cg:cg::Convergence step 25 value 3.22787e-12 +DEAL:cg:cg::Starting value 0.408285 +DEAL:cg:cg::Convergence step 30 value 0.00738618 +DEAL:cg:cg::Starting value 0.0169890 +DEAL:cg:cg::Convergence step 25 value 1.31488e-12 +DEAL:cg:cg::Starting value 0.205984 +DEAL:cg:cg::Convergence step 30 value 0.00300423 +DEAL:cg:cg::Starting value 0.00694787 +DEAL:cg:cg::Convergence step 25 value 5.53097e-13 +DEAL:cg:cg::Starting value 0.0555164 +DEAL:cg:cg::Convergence step 30 value 0.00172572 +DEAL:cg:cg::Starting value 0.00208394 +DEAL:cg:cg::Convergence step 25 value 1.92428e-13 +DEAL:cg:cg::Starting value 0.0202742 +DEAL:cg:cg::Convergence step 30 value 0.000727924 +DEAL:cg:cg::Starting value 0.000839806 +DEAL:cg:cg::Convergence step 25 value 7.63218e-14 +DEAL:cg:cg::Starting value 0.00899883 +DEAL:cg:cg::Convergence step 30 value 0.000386101 +DEAL:cg:cg::Starting value 0.000366998 +DEAL:cg:cg::Convergence step 25 value 2.97977e-14 +DEAL:cg:cg::Starting value 0.00371456 +DEAL:cg:cg::Convergence step 30 value 0.000143672 +DEAL:cg:cg::Starting value 0.000149522 +DEAL:cg:cg::Convergence step 25 value 1.25747e-14 +DEAL:cg:cg::Starting value 0.00151238 +DEAL:cg:cg::Convergence step 30 value 2.76485e-05 +DEAL:cg:cg::Starting value 5.80561e-05 +DEAL:cg:cg::Convergence step 25 value 4.95193e-15 +DEAL:cg:cg::Starting value 0.000568491 +DEAL:cg:cg::Convergence step 30 value 1.51566e-05 +DEAL:cg:cg::Starting value 2.14606e-05 +DEAL:cg:cg::Convergence step 25 value 1.85827e-15 +DEAL:cg:cg::Starting value 0.000231159 +DEAL:cg:cg::Convergence step 30 value 4.26986e-06 +DEAL:cg:cg::Starting value 8.83789e-06 +DEAL:cg:cg::Convergence step 25 value 7.59785e-16 +DEAL:cg:cg::Starting value 9.93855e-05 +DEAL:cg:cg::Convergence step 30 value 1.66407e-06 +DEAL:cg:cg::Starting value 3.98372e-06 +DEAL:cg:cg::Convergence step 25 value 3.25992e-16 +DEAL:cg:cg::Starting value 3.74185e-05 +DEAL:cg:cg::Convergence step 30 value 1.12542e-06 +DEAL:cg:cg::Starting value 1.44245e-06 +DEAL:cg:cg::Convergence step 25 value 1.13963e-16 +DEAL:cg:cg::Starting value 1.41827e-05 +DEAL:cg:cg::Convergence step 30 value 3.57707e-07 +DEAL:cg:cg::Starting value 5.33178e-07 +DEAL:cg:cg::Convergence step 25 value 4.28754e-17 +DEAL:cg:cg::Starting value 5.62814e-06 +DEAL:cg:cg::Convergence step 30 value 1.27950e-07 +DEAL:cg:cg::Starting value 2.23213e-07 +DEAL:cg:cg::Convergence step 25 value 1.81003e-17 +DEAL:cg:cg::Starting value 2.12806e-06 +DEAL:cg:cg::Convergence step 30 value 4.90715e-08 +DEAL:cg:cg::Starting value 8.59275e-08 +DEAL:cg:cg::Convergence step 25 value 6.89855e-18 +DEAL:cg:cg::Starting value 8.70957e-07 +DEAL:cg:cg::Convergence step 30 value 1.25220e-08 +DEAL:cg:cg::Starting value 3.33966e-08 +DEAL:cg:cg::Convergence step 25 value 2.74359e-18 +DEAL:cg:cg::Starting value 2.94918e-07 +DEAL:cg:cg::Convergence step 30 value 7.19353e-09 +DEAL:cg:cg::Starting value 1.20206e-08 +DEAL:cg:cg::Convergence step 25 value 9.72271e-19 +DEAL:cg:cg::Starting value 1.20514e-07 +DEAL:cg:cg::Convergence step 30 value 2.48469e-09 +DEAL:cg:cg::Starting value 4.68186e-09 +DEAL:cg:cg::Convergence step 25 value 3.97912e-19 +DEAL:cg:cg::Starting value 5.18972e-08 +DEAL:cg:cg::Convergence step 30 value 1.83964e-09 +DEAL:cg:cg::Starting value 2.00490e-09 +DEAL:cg:cg::Convergence step 24 value 4.00689e-19 +DEAL:cg:cg::Starting value 1.99515e-08 +DEAL:cg:cg::Convergence step 30 value 9.89302e-10 +DEAL:cg:cg::Starting value 7.90511e-10 +DEAL:cg:cg::Convergence step 23 value 4.17184e-19 +DEAL:cg:cg::Starting value 7.90659e-09 +DEAL:cg:cg::Convergence step 30 value 3.01626e-10 +DEAL:cg:cg::Starting value 3.15023e-10 +DEAL:cg:cg::Convergence step 22 value 4.65119e-19 +DEAL:cg:cg::Starting value 3.05591e-09 +DEAL:cg:cg::Convergence step 30 value 8.78930e-11 +DEAL:cg:cg::Starting value 1.23874e-10 +DEAL:cg:cg::Convergence step 21 value 4.01642e-19 +DEAL:cg:cg::Starting value 1.44274e-09 +DEAL:cg:cg::Convergence step 30 value 3.28342e-11 +DEAL:cg:cg::Starting value 5.36043e-11 +DEAL:cg:cg::Convergence step 20 value 4.90919e-19 +DEAL:cg:cg::Starting value 5.39812e-10 +DEAL:cg:cg::Convergence step 30 value 1.21691e-11 +DEAL:cg:cg::Starting value 2.06924e-11 +DEAL:cg:cg::Convergence step 19 value 5.18474e-19 +DEAL:cg:cg::Starting value 1.95959e-10 +DEAL:cg:cg::Convergence step 30 value 3.48296e-12 +DEAL:cg:cg::Starting value 7.18708e-12 +DEAL:cg:cg::Convergence step 18 value 4.56473e-19 +DEAL:cg:cg::Starting value 7.93580e-11 +DEAL:cg:cg::Convergence step 30 value 1.28153e-12 +DEAL:cg:cg::Starting value 2.99859e-12 +DEAL:cg:cg::Convergence step 17 value 4.66330e-19 +DEAL:cg:cg::Starting value 3.15137e-11 +DEAL:cg:cg::Convergence step 30 value 5.03608e-13 +DEAL:cg:cg::Starting value 1.21338e-12 +DEAL:cg:cg::Convergence step 16 value 4.39333e-19 +DEAL:cg:cg::Starting value 1.28779e-11 +DEAL:cg:cg::Convergence step 30 value 3.26557e-13 +DEAL:cg:cg::Starting value 4.68242e-13 +DEAL:cg:cg::Convergence step 15 value 4.13608e-19 +DEAL:cg:cg::Starting value 4.57873e-12 +DEAL:cg:cg::Convergence step 30 value 1.99147e-13 +DEAL:cg:cg::Starting value 1.70285e-13 +DEAL:cg:cg::Convergence step 14 value 4.69229e-19 +DEAL:cg:cg::Starting value 1.63301e-12 +DEAL:cg:cg::Convergence step 30 value 5.33467e-14 +DEAL:cg:cg::Starting value 6.71379e-14 +DEAL:cg:cg::Convergence step 13 value 4.72485e-19 +DEAL:cg::Convergence step 33 value 7.32840e-13 +DEAL::33 CG Schur complement iterations to obtain convergence. +DEAL:cg::Starting value 0.143776 +DEAL:cg::Convergence step 19 value 1.12728e-11 +DEAL::Errors: ||e_p||_L2 = 0.00221522, ||e_u||_L2 = 2.44171e-09 \ No newline at end of file