From 8ab7de6561ea2957986fc562c6fa98e4208bc123 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 9 May 2019 18:32:55 +0200 Subject: [PATCH] Add step-65 --- examples/step-65/CMakeLists.txt | 39 +++ examples/step-65/doc/intro.dox | 14 ++ examples/step-65/doc/results.dox | 61 +++++ examples/step-65/step-65.cc | 418 +++++++++++++++++++++++++++++++ 4 files changed, 532 insertions(+) create mode 100644 examples/step-65/CMakeLists.txt create mode 100644 examples/step-65/doc/intro.dox create mode 100644 examples/step-65/doc/results.dox create mode 100644 examples/step-65/step-65.cc diff --git a/examples/step-65/CMakeLists.txt b/examples/step-65/CMakeLists.txt new file mode 100644 index 0000000000..bcedad9028 --- /dev/null +++ b/examples/step-65/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-65 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-65") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +FIND_PACKAGE(deal.II 9.2.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-65/doc/intro.dox b/examples/step-65/doc/intro.dox new file mode 100644 index 0000000000..26b84f25d5 --- /dev/null +++ b/examples/step-65/doc/intro.dox @@ -0,0 +1,14 @@ + +
+ + +This program was contributed by Martin Kronbichler. + + + +

Introduction

+ +This tutorial program presents a way to handle expensive manifold classes, +exemplified for the TransfiniteInterpolationManifold, a special manifold class +that can blend a curved approximation on a boundary with straight-sided +descriptions elsewhere. diff --git a/examples/step-65/doc/results.dox b/examples/step-65/doc/results.dox new file mode 100644 index 0000000000..9738b1c95c --- /dev/null +++ b/examples/step-65/doc/results.dox @@ -0,0 +1,61 @@ + +

Results

+ +

Program output

+ +@code +$ make run +Scanning dependencies of target step-65 +[ 33%] Building CXX object CMakeFiles/step-65.dir/step-65.cc.o +[ 66%] Linking CXX executable step-65 +[ 66%] Built target step-65 +[100%] Run step-65 with Release configuration + +====== Running with the basic MappingQGeneric class ====== + + Number of active cells: 6656 + Number of degrees of freedom: 181609 + Number of solver iterations: 285 + L2 error vs exact solution: 8.99328e-08 + H1 error vs exact solution: 6.45341e-06 + + ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 56.3s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assemble linear system | 1 | 5.39s | 9.6% | +| Compute constraints | 1 | 0.109s | 0.19% | +| Compute error estimator | 1 | 17.3s | 31% | +| Compute error norms | 1 | 9.53s | 17% | +| Solve linear system | 1 | 10.6s | 19% | +| Write output | 2 | 10.3s | 18% | ++---------------------------------+-----------+------------+------------+ + +====== Running with the optimized MappingQCache class ====== + + Memory consumption cache: 22.9976 MB + Number of active cells: 6656 + Number of degrees of freedom: 181609 + Number of solver iterations: 285 + L2 error vs exact solution: 8.99328e-08 + H1 error vs exact solution: 6.45341e-06 + + ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 19.7s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assemble linear system | 1 | 0.876s | 4.4% | +| Compute constraints | 1 | 0.00368s | 0% | +| Compute error estimator | 1 | 0.488s | 2.5% | +| Compute error norms | 1 | 0.513s | 2.6% | +| Initialize mapping cache | 1 | 5.2s | 26% | +| Solve linear system | 1 | 10.6s | 54% | +| Write output | 2 | 1.86s | 9.5% | ++---------------------------------+-----------+------------+------------+ + +[100%] Built target run +@endcode diff --git a/examples/step-65/step-65.cc b/examples/step-65/step-65.cc new file mode 100644 index 0000000000..44bdb5abdb --- /dev/null +++ b/examples/step-65/step-65.cc @@ -0,0 +1,418 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + +// This tutorial program was contributed by Martin Kronbichler + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + + + +namespace step65 +{ + using namespace dealii; + + template + class ExactSolution : public Function + { + public: + virtual double value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + if (p.norm_square() < 0.25) + return p.norm_square(); + else + return 0.1 * p.norm_square() + (0.25 - 0.025); + } + + virtual Tensor<1, dim> + gradient(const Point &p, + const unsigned int /*component*/ = 0) const override + { + if (p.norm_square() < 0.25) + return 2. * p; + else + return 0.2 * p; + } + }; + + + template + class Coefficient : public Function + { + public: + virtual double value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + if (p.norm_square() < 0.25) + return 0.5; + else + return 5.; + } + }; + + + + template + class PoissonProblem + { + public: + PoissonProblem(); + void run(); + + private: + void create_grid(); + void setup_system(const Mapping &mapping); + void assemble_system(const Mapping &mapping); + void solve(); + void postprocess(const Mapping &mapping); + void output_results(const Mapping &mapping); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + TimerOutput timer; + }; + + + + template + PoissonProblem::PoissonProblem() + : fe(3) + , dof_handler(triangulation) + , timer(std::cout, TimerOutput::never, TimerOutput::wall_times) + {} + + + + template + void PoissonProblem::create_grid() + { + Triangulation tria_outer, tria_inner; + GridGenerator::hyper_shell( + tria_outer, Point(), 0.5, std::sqrt(dim), 2 * dim); + GridGenerator::hyper_ball(tria_inner, Point(), 0.5); + GridGenerator::merge_triangulations(tria_inner, tria_outer, triangulation); + triangulation.reset_all_manifolds(); + triangulation.set_all_manifold_ids(0); + for (const auto &cell : triangulation.cell_iterators()) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + bool face_at_sphere_boundary = true; + for (unsigned int v = 0; + v < GeometryInfo::vertices_per_cell; + ++v) + if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) > + 1e-12) + face_at_sphere_boundary = false; + if (face_at_sphere_boundary) + cell->face(f)->set_all_manifold_ids(1); + } + if (cell->center().norm_square() < 0.25) + cell->set_material_id(1); + else + cell->set_material_id(0); + } + triangulation.set_manifold(1, SphericalManifold()); + TransfiniteInterpolationManifold transfinite_manifold; + transfinite_manifold.initialize(triangulation); + triangulation.set_manifold(0, transfinite_manifold); + + triangulation.refine_global(9 - 2 * dim); + } + + + + template + void PoissonProblem::setup_system(const Mapping &mapping) + { + dof_handler.distribute_dofs(fe); + std::cout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl; + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + { + TimerOutput::Scope scope(timer, "Compute constraints"); + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, ExactSolution(), constraints); + constraints.close(); + } + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + + + template + void PoissonProblem::assemble_system(const Mapping &mapping) + { + TimerOutput::Scope scope(timer, "Assemble linear system"); + QGauss quadrature_formula(fe.degree + 1); + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + FullMatrix partial_matrix(dofs_per_cell, dim * n_q_points); + std::vector local_dof_indices(dofs_per_cell); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_rhs = 0.; + fe_values.reinit(cell); + for (unsigned int q_index = 0; q_index < n_q_points; ++q_index) + { + const double current_coefficient = + Coefficient().value(fe_values.quadrature_point(q_index)); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int d = 0; d < dim; ++d) + partial_matrix(i, q_index * dim + d) = + std::sqrt(fe_values.JxW(q_index) * current_coefficient) * + fe_values.shape_grad(i, q_index)[d]; + cell_rhs(i) += + (fe_values.shape_value(i, q_index) * // phi_i(x_q) + (-dim) * // f(x_q) + fe_values.JxW(q_index)); // dx + } + } + partial_matrix.mTmult(cell_matrix, partial_matrix); + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } + } + + + + template + void PoissonProblem::solve() + { + TimerOutput::Scope scope(timer, "Solve linear system"); + SolverControl solver_control(1000, 1e-12); + SolverCG<> solver(solver_control); + PreconditionJacobi<> preconditioner; + preconditioner.initialize(system_matrix); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + constraints.distribute(solution); + std::cout << " Number of solver iterations: " + << solver_control.last_step() << std::endl; + } + + + + template + void PoissonProblem::postprocess(const Mapping &mapping) + { + { + TimerOutput::Scope scope(timer, "Write output"); + Timer time; + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags(flags); + + Vector material_ids(triangulation.n_active_cells()); + for (const auto &cell : triangulation.active_cell_iterators()) + material_ids[cell->active_cell_index()] = cell->material_id(); + data_out.add_data_vector(material_ids, "material_ids"); + + data_out.build_patches(mapping, + fe.degree, + DataOut::curved_inner_cells); + + std::ofstream file( + ("solution-" + + std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) + + ".vtu") + .c_str()); + data_out.write_vtu(file); + } + + { + TimerOutput::Scope scope(timer, "Compute error norms"); + Vector norm_per_cell_p(triangulation.n_active_cells()); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + ExactSolution(), + norm_per_cell_p, + QGauss(fe.degree + 2), + VectorTools::L2_norm); + std::cout << " L2 error vs exact solution: " + << norm_per_cell_p.l2_norm() << std::endl; + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + ExactSolution(), + norm_per_cell_p, + QGauss(fe.degree + 2), + VectorTools::H1_norm); + std::cout << " H1 error vs exact solution: " + << norm_per_cell_p.l2_norm() << std::endl; + } + + { + TimerOutput::Scope scope(timer, "Compute error estimator"); + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + mapping, + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + solution, + estimated_error_per_cell); + } + } + + + + template + void PoissonProblem::output_results(const Mapping &mapping) + { + TimerOutput::Scope scope(timer, "Write output"); + Timer time; + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags(flags); + + Vector material_ids(triangulation.n_active_cells()); + for (const auto &cell : triangulation.active_cell_iterators()) + material_ids[cell->active_cell_index()] = cell->material_id(); + data_out.add_data_vector(material_ids, "material_ids"); + + data_out.build_patches(mapping, + fe.degree, + DataOut::curved_inner_cells); + + std::ofstream file( + ("solution-" + + std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) + ".vtu") + .c_str()); + data_out.write_vtu(file); + } + + + + template + void PoissonProblem::run() + { + create_grid(); + { + std::cout << std::endl + << "====== Running with the basic MappingQGeneric class ====== " + << std::endl + << std::endl; + + MappingQGeneric mapping(fe.degree + 1); + setup_system(mapping); + assemble_system(mapping); + solve(); + postprocess(mapping); + output_results(mapping); + + timer.print_summary(); + } + + timer.reset(); + + { + std::cout + << "====== Running with the optimized MappingQCache class ====== " + << std::endl + << std::endl; + + MappingQCache mapping(fe.degree + 1); + { + TimerOutput::Scope scope(timer, "Initialize mapping cache"); + mapping.initialize(triangulation, MappingQGeneric(fe.degree + 1)); + } + std::cout << " Memory consumption cache: " + << 1e-6 * mapping.memory_consumption() << " MB" << std::endl; + + setup_system(mapping); + assemble_system(mapping); + solve(); + postprocess(mapping); + output_results(mapping); + + timer.print_summary(); + } + } + +} // namespace step65 + + +int main() +{ + step65::PoissonProblem<3> test_program; + test_program.run(); + return 0; +} -- 2.39.5