From: Marco Feder Date: Tue, 28 Sep 2021 17:52:20 +0000 (+0200) Subject: first version of adv reaction estimator X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F91%2Fhead;p=code-gallery.git first version of adv reaction estimator --- diff --git a/advection_reaction_estimator/CMakeLists.txt b/advection_reaction_estimator/CMakeLists.txt new file mode 100644 index 0000000..e0bef59 --- /dev/null +++ b/advection_reaction_estimator/CMakeLists.txt @@ -0,0 +1,28 @@ +# Set the name of the project and target: +SET(TARGET "DG_advection_reaction_estimator") + +# Declare all source files the target consists of: +SET(TARGET_SRC + source/DG_advection_reaction.cc + main.cc +) + +PROJECT(${TARGET} CXX) + + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.3.0 + 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/advection_reaction_estimator/README.md b/advection_reaction_estimator/README.md new file mode 100644 index 0000000..65558e5 --- /dev/null +++ b/advection_reaction_estimator/README.md @@ -0,0 +1,122 @@ +# A posteriori error estimator for first order hyperbolic problems + + +## Running the code: + +As in the tutorial programs, type + +`cmake -DDEAL_II_DIR=/path/to/deal.II .` + +on the command line to configure the program. After that you can compile with `make` and run with either `make run` or using + +`./DG_advection_reaction` + +on the command line. + +### Parameter file: + +If you run `./DG_advection_reaction parameters.prm`, an error message will tell you that a parameter file has been created for you. You can open it and change some useful parameters like the number of refinement cycles, the advection coefficient, and others. If you don't specify anything, then the default values used for the test case (see paragraph below) will be used. + + + + +## The problem: +This program solves the problem, for $\Omega \in \mathbb{R^2}$ + +$$\begin{cases} b \cdot \nabla u + c u = f \qquad \text{in } \Omega \\ +\qquad \qquad u=g \qquad \text{on } \partial_{-}\Omega \end{cases}$$ + +where $g \in L^2(\partial_{-}\Omega)$ and $\partial_{-}\Omega=\{ x \in \partial \Omega: b(x)\cdot n(x) <0\}$ is the inflow part of the boundary, with $b=(b_1,b_2) \in \mathbb{R^2}$. As we know from classical DG theory, we need to ensure that $$c(x) - \frac{1}{2}\nabla \cdot b \geq \gamma_0 >0$$for some positive $\gamma_0$ so that we have coercivity in $L^2$ at the continuous level. Discrete coercivity is achieved by using a stronger norm which takes care of jumps, see Di-Pietro Ern [1] for details. + + +## The weak formulation: + + + +As trial space we choose $V_h = \{ v_h \in L^2(\Omega): v_h \in P^1(\mathbb{T_h})\} \notin H^1(\Omega)$. If we integrate by parts and sum over all cells + +$$\sum_{T \in \mathbb{T}_h} \Bigl( (-u,\beta \cdot \nabla v_h) _T + (c u,v_h)_T + \bigl<(b \cdot n) u ,v_h \bigr>_{\partial T} \Bigr) = (f,v_h)_{\Omega}$$ + +and use the so-called DG magic formula and exploit the property $[bu]_{\mathbb{F}^i} = 0$ where $\mathbb{F}^i$ are set of internal faces we obtain the (unstable!) formulation: + +Find $u_h \in V_h$: + +$$ + a_h(u_h,v_h) + b_h(u_h,v_h)=l(v_h) \qquad \forall v_h \in V_h +$$ +where +$$ +a_h(u,v_h)=\sum_{T \in \mathbb{T}_h} \Bigl( (-u,b \cdot \nabla v_h) _T + (c u,v_h)_T \Bigr) +$$ + +$$ b_h(u,v_h)= \sum_{F \not \in \partial_{-}\Omega} \bigl< \{ b u\}, [v_h]\bigr>_F $$ + +$$ + l(v_h)= (f,v_h)_{\Omega} - \sum_{F \in \partial_{-}\Omega} \bigl< (b \cdot n) g,v_h \bigr>_F +$$ + +It's well known this formulation is coercive only in $L^2$, hence the formulation is unstable as we don't "see" the derivatives. To stabilize this, we can use a jump-penalty term, i.e. our $b_h$ is replaced by: + +$$b_h^s(u_h,v_h)=b_h(u_h,v_h)+ \sum_{F \in \mathbb{F}^i} \bigl< c_F [u_h],[v_h] \bigr> $$ + +where $c_F>0$ is a function on each edge such that $c_F \geq \theta |b \cdot n|$ for some positive $\theta$. In this program, $\theta=\frac{1}{2}$ and $c_F = \frac{1}{2} |b \cdot n|$, which corresponds to an upwind formulation. Notice that consistency is trivially achieved, as $[u]_{\mathbb{F}^i} =0$. This formulation is stable in the energy norm + +$$ |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 + \sum_{F \in \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2 \Bigr)^{\frac{1}{2}}$$ + +(well defined on $H^1(\Omega) + V_h$) and moreover we have the a-priori bound: + +$$|||u-u_h||| \leq C h^{k+\frac{1}{2}}||u||_{k+1,\Omega} $$ + +valid for $u \in H^{k+1}(\Omega)$. + +See Brezzi-Marini-Süli [3] for more details. + + + +## A-posteriori error estimator: + +The estimator is the one proposed by Georgoulis, Edward Hall and Charalambos Makridakis in [3]. This approach is quite different with respect to other works in the field, as the authors are trying to develop an estimator for the original hyperbolic problem, rather than taking the hyperbolic regime as the vanishing diffusivity limit. + +The reliability is: + +$$|||u-u_h|||^2 \leq C || \sqrt{b \cdot n}[u_h]||_{\Gamma^{-}}^2 + C \sum_{T \in \mathbb{T}_h}\Bigl( ||\beta (g-u_h^+)||_{\partial_{-}T \cap \partial_{-} \Omega}^2 +||f-c u_h - \Pi(f- cu_h)||_T^2 \Bigr)$$ + +where: + +- $\Pi$ is the (local) $L^2$ orthogonal projection onto $V_h$ + +- $\Gamma$ is the skeleton of the mesh + +- $c$ is constant + +- $\beta = |b \cdot n|$ + +- $u_h^+$ is the interior trace from the current cell $T$ of a the finite element function $u_h$. + + + +## Test case: + +The following test case has been taken from [3]. Consider: +- $c=1$ +- $b=(1,1)$ +- $f$ to be such that the exact solution is $u(x,y)=\tanh(100(x+y-\frac{1}{2}))$ +This solution has an internal layer along the line $y=\frac{1}{2} -x$, hence we would like to see that part of the domain to be much more refined than the rest. + +The next image is the 3D view of the numerical solution: + +![Screenshot](doc/images/warp_by_scalar_solution_layer.png) + +More interestingly, we see that the estimator has been able to capture the layer. Here a bulk-chasing criterion is used, with bottom fraction ´0.5´ and no coarsening. This mesh is obtained after 12 refinement cycles. +![Screenshot](doc/images/refined_mesh_internal_layer.png) + + +If we look at the decrease of the energy norm of the error in the globally refined case and in the adaptively case, with respect to the DoFs, we obtain: + +![Screenshot](doc/images/adaptive_vs_global_refinement.png) + +## References +* [1] Emmanuil H. Georgoulis, Edward Hall and Charalambos Makridakis (2013), Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems. DOI: [10.1007/978-3-319-01818-8_8 +](https://link.springer.com/chapter/10.1007%2F978-3-319-01818-8_8) +* [2] Di Pietro, Daniele Antonio and Ern, Alexandre (2012), Mathematical Aspects of Discontinuous Galerkin Methods. ISBN: [978-3-642-22980-0](https://www.springer.com/gp/book/9783642229794) +* [3] Franco Brezzi, Luisa Donatella Marini and Endre Süli (2004) Discontinuous Galerkin Methods for First-Order Hyperbolic Problems. DOI: [10.1142/S0218202504003866](https://doi.org/10.1142/S0218202504003866) \ No newline at end of file diff --git a/advection_reaction_estimator/doc/author b/advection_reaction_estimator/doc/author new file mode 100644 index 0000000..1fe9053 --- /dev/null +++ b/advection_reaction_estimator/doc/author @@ -0,0 +1 @@ +Marco Feder diff --git a/advection_reaction_estimator/doc/builds-on b/advection_reaction_estimator/doc/builds-on new file mode 100644 index 0000000..154f8bf --- /dev/null +++ b/advection_reaction_estimator/doc/builds-on @@ -0,0 +1 @@ +step-12 step-74 diff --git a/advection_reaction_estimator/doc/entry-name b/advection_reaction_estimator/doc/entry-name new file mode 100644 index 0000000..8293f59 --- /dev/null +++ b/advection_reaction_estimator/doc/entry-name @@ -0,0 +1 @@ +Adaptive advection-reaction diff --git a/advection_reaction_estimator/doc/images/adaptive_vs_global_refinement.png b/advection_reaction_estimator/doc/images/adaptive_vs_global_refinement.png new file mode 100644 index 0000000..70d3039 Binary files /dev/null and b/advection_reaction_estimator/doc/images/adaptive_vs_global_refinement.png differ diff --git a/advection_reaction_estimator/doc/images/refined_mesh_internal_layer.png b/advection_reaction_estimator/doc/images/refined_mesh_internal_layer.png new file mode 100644 index 0000000..bd38f0e Binary files /dev/null and b/advection_reaction_estimator/doc/images/refined_mesh_internal_layer.png differ diff --git a/advection_reaction_estimator/doc/images/warp_by_scalar_solution_layer.png b/advection_reaction_estimator/doc/images/warp_by_scalar_solution_layer.png new file mode 100644 index 0000000..f0c30f3 Binary files /dev/null and b/advection_reaction_estimator/doc/images/warp_by_scalar_solution_layer.png differ diff --git a/advection_reaction_estimator/doc/tooltip b/advection_reaction_estimator/doc/tooltip new file mode 100644 index 0000000..bed031f --- /dev/null +++ b/advection_reaction_estimator/doc/tooltip @@ -0,0 +1 @@ +Implementation of an a-posteriori error estimator for first order hyperbolic problems diff --git a/advection_reaction_estimator/include/DG_advection_reaction.h b/advection_reaction_estimator/include/DG_advection_reaction.h new file mode 100644 index 0000000..0146eab --- /dev/null +++ b/advection_reaction_estimator/include/DG_advection_reaction.h @@ -0,0 +1,124 @@ +#ifndef INCLUDE_DG_UPWIND_H_ +#define INCLUDE_DG_UPWIND_H_ + +// The first few files have already been covered in tutorials and will +// thus not be further commented on: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +// This header is needed for FEInterfaceValues to compute integrals on +// interfaces: +#include +//Solver +#include +#include +#include +// We are going to use gradients as refinement indicator. +#include +// Using using the mesh_loop from the MeshWorker framework +#include + +#include + +//To enable parameter handling +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +//This is a struct used only for throwing an exception when theta parameter is not okay. +struct theta_exc +{ + std::string message; + theta_exc(std::string &&s) : message{std::move(s)} {}; + const char *what() const { return message.c_str(); } +}; + +using namespace dealii; +// @sect3{Class declaration} +// In the following we have the declaration of the functions used in the program. As we want to use +// parameter files, we need to derive our class from `ParameterAcceptor`. +template +class AdvectionReaction : ParameterAcceptor +{ +public: + AdvectionReaction(); + void initialize_params(const std::string &filename); + void run(); + +private: + using Iterator = typename DoFHandler::active_cell_iterator; + void parse_string(const std::string ¶meters); + void setup_system(); + void assemble_system(); + void solve(); + void refine_grid(); + void output_results(const unsigned int cycle) const; + void compute_error(); + double compute_energy_norm(); + void compute_local_projection_and_estimate(); + + Triangulation triangulation; + const MappingQ1 mapping; + + // Furthermore we want to use DG elements. + std::unique_ptr> fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector right_hand_side; + Vector energy_norm_square_per_cell; + Vector error_indicator_per_cell; + + // So far we declared the usual objects. Hereafter we declare `FunctionParser` objects + FunctionParser exact_solution; + FunctionParser boundary_conditions; + FunctionParser rhs; + FunctionParser advection_coeff; + + unsigned int fe_degree = 1; + + // and then we define default values that will be parsed from the following strings + std::string exact_solution_expression = "tanh(100*(x+y-0.5))"; //internal layer solution + std::string rhs_expression = "-200*tanh(100*x + 100*y - 50.0)^2 + tanh(100*x + 100*y - 50.0) + 200"; + std::string advection_coefficient_expression = "1.0"; + std::string boundary_conditions_expression = "tanh(100*x + 100*y - 50.0)"; + std::string refinement = "residual"; + std::string output_filename = "DG_estimator"; + std::map constants; + ParsedConvergenceTable error_table; + + bool use_direct_solver = true; + unsigned int n_refinement_cycles = 14; + unsigned int n_global_refinements = 3; + double theta = 0.5; //default is 0.5 so that I have classical upwind flux +}; + +#endif /* INCLUDE_DG_UPWIND_H_ */ diff --git a/advection_reaction_estimator/main.cc b/advection_reaction_estimator/main.cc new file mode 100644 index 0000000..91eac65 --- /dev/null +++ b/advection_reaction_estimator/main.cc @@ -0,0 +1,62 @@ +#include "include/DG_advection_reaction.h" + +int main(int argc, char **argv) +{ + + try + { + std::string par_name = ""; + if (argc > 1) + { + par_name = argv[1]; + } + deallog.depth_console(2); //solver infos + AdvectionReaction<2> dgmethod; + if (par_name != "") + { + dgmethod.initialize_params(par_name); + } + + dgmethod.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 (const theta_exc &theta_range) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << theta_range.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/advection_reaction_estimator/source/DG_advection_reaction.cc b/advection_reaction_estimator/source/DG_advection_reaction.cc new file mode 100644 index 0000000..e3f77c2 --- /dev/null +++ b/advection_reaction_estimator/source/DG_advection_reaction.cc @@ -0,0 +1,825 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 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. + * + * --------------------------------------------------------------------- + + + + * + * Author: Marco Feder, SISSA, 2021 + * + */ + +#include "../include/DG_advection_reaction.h" + +//Compute and returns the wind field b +template +Tensor<1, dim> beta(const Point &p) +{ + Assert(dim >= 2, ExcNotImplemented()); + (void)p; //suppress warnings from p + Tensor<1, dim> wind_field; + wind_field[0] = 1.0; + wind_field[1] = 1.0; + + return wind_field; +} + +// @sect3{The ScratchData and CopyData classes} +// +// The following objects are the scratch and copy objects we use in the call +// to MeshWorker::mesh_loop(). The new object is the FEInterfaceValues object, +// that works similar to FEValues or FEFacesValues, except that it acts on +// an interface between two cells and allows us to assemble the interface +// terms in our weak form. +template +struct ScratchData +{ + ScratchData(const Mapping &mapping, const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &quadrature_face, + const UpdateFlags update_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values, + const UpdateFlags interface_update_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values | update_normal_vectors) : fe_values(mapping, fe, quadrature, update_flags), fe_interface_values(mapping, fe, quadrature_face, interface_update_flags) + { + } + + ScratchData(const ScratchData &scratch_data) : fe_values(scratch_data.fe_values.get_mapping(), + scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + scratch_data.fe_values.get_update_flags()), + fe_interface_values( + scratch_data.fe_interface_values.get_mapping(), + scratch_data.fe_interface_values.get_fe(), + scratch_data.fe_interface_values.get_quadrature(), + scratch_data.fe_interface_values.get_update_flags()) + { + } + + FEValues fe_values; + FEInterfaceValues fe_interface_values; +}; + +struct CopyDataFace +{ + FullMatrix cell_matrix; + std::vector joint_dof_indices; + std::array values; + std::array cell_indices; +}; + +struct CopyData +{ + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + + double value; + double value_estimator; + unsigned int cell_index; + + FullMatrix cell_mass_matrix; + Vector cell_mass_rhs; + + template + void reinit(const Iterator &cell, unsigned int dofs_per_cell) + { + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell); + + cell_rhs.reinit(dofs_per_cell); + cell_mass_rhs.reinit(dofs_per_cell); + + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } +}; + +// @sect3{Auxiliary function} +// This auxiliary function is taken from step-74 and it's used to +// compute the jump of the finite element function $u_h$ on a face. +template +void get_function_jump(const FEInterfaceValues &fe_iv, + const Vector &solution, + std::vector &jump) +{ + const unsigned int n_q = fe_iv.n_quadrature_points; + std::array, 2> face_values; + jump.resize(n_q); + for (unsigned int i = 0; i < 2; ++i) + { + face_values[i].resize(n_q); + fe_iv.get_fe_face_values(i).get_function_values(solution, + face_values[i]); + } + for (unsigned int q = 0; q < n_q; ++q) + jump[q] = face_values[0][q] - face_values[1][q]; +} + +// We start with the constructor. The 1 in the constructor call of +// fe is the polynomial degree. +template +AdvectionReaction::AdvectionReaction() : mapping(), + dof_handler(triangulation) +{ + + add_parameter("Finite element degree", fe_degree); + add_parameter("Problem constants", constants); + add_parameter("Output filename", output_filename); + add_parameter("Use direct solver", use_direct_solver); + add_parameter("Number of refinement cycles", n_refinement_cycles); + add_parameter("Number of global refinement", n_global_refinements); + add_parameter("Refinement", refinement); + add_parameter("Exact solution expression", exact_solution_expression); + add_parameter("Boundary conditions expression", boundary_conditions_expression); + add_parameter("Theta", theta); + add_parameter("Advection coefficient expression", advection_coefficient_expression); + add_parameter("Right hand side expression", rhs_expression); + + // + this->prm.enter_subsection("Error table"); + error_table.add_parameters(this->prm); + this->prm.leave_subsection(); +} + +template +void AdvectionReaction::initialize_params(const std::string &filename) +{ + + ParameterAcceptor::initialize(filename, "last_used_parameters.prm", ParameterHandler::Short); + if (theta < 0.0 || theta > 10.0 || std::abs(theta) < 1e-12) + { + throw(theta_exc("Theta parameter is not in a suitable range: see paper by Brezzi, Marini, Suli for an extended discussion")); + } +} + +template +void AdvectionReaction::parse_string(const std::string ¶meters) +{ + ParameterAcceptor::prm.parse_input_from_string(parameters); + ParameterAcceptor::parse_all_parameters(); +} + +template +void AdvectionReaction::setup_system() +{ + + // first need to distribute the DoFs. + if (!fe) + { + fe = std::make_unique>(fe_degree); + const auto vars = dim == 2 ? "x,y" : "x,y,z"; + exact_solution.initialize(vars, exact_solution_expression, constants); + rhs.initialize(vars, rhs_expression, constants); + advection_coeff.initialize(vars, advection_coefficient_expression, constants); + boundary_conditions.initialize(vars, boundary_conditions_expression, constants); + } + dof_handler.distribute_dofs(*fe); + + // To build the sparsity pattern for DG discretizations, we can call the + // function analogue to DoFTools::make_sparsity_pattern, which is called + // DoFTools::make_flux_sparsity_pattern: + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); //DG sparsity pattern generator + sparsity_pattern.copy_from(dsp); + + // Finally, we set up the structure of all components of the linear system. + system_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + right_hand_side.reinit(dof_handler.n_dofs()); +} + +//in the call to MeshWorker::mesh_loop() we only need to specify what should happen on +// each cell, each boundary face, and each interior face. These three tasks +// are handled by the lambda functions inside the function below. + +template +void AdvectionReaction::assemble_system() +{ + + using Iterator = typename DoFHandler::active_cell_iterator; + + const QGauss quadrature = fe->tensor_degree() + 1; + const QGauss quadrature_face = fe->tensor_degree() + 1; + + // This is the function that will be executed for each cell. + const auto cell_worker = [&](const Iterator &cell, + ScratchData &scratch_data, CopyData ©_data) + { + FEValues fe_values_continuous(*fe, + quadrature, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int n_dofs = scratch_data.fe_values.get_fe().n_dofs_per_cell(); + copy_data.reinit(cell, n_dofs); + scratch_data.fe_values.reinit(cell); + + const auto &q_points = scratch_data.fe_values.get_quadrature_points(); + + const FEValues &fe_v = scratch_data.fe_values; + const std::vector &JxW = fe_v.get_JxW_values(); + + for (unsigned int point = 0; point < fe_v.n_quadrature_points; + ++point) + { + auto beta_q = beta(q_points[point]); + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + copy_data.cell_matrix(i, j) += (-beta_q // -\beta + * fe_v.shape_grad(i, point) // \nabla \phi_i + * fe_v.shape_value(j, point) // \phi_j + + + advection_coeff.value(q_points[point]) * //gamma + fe_v.shape_value(i, point) //phi_i + * fe_v.shape_value(j, point) //phi_j + ) * + JxW[point]; // dx + } + copy_data.cell_rhs(i) += + rhs.value(q_points[point]) // f(x_q) + * fe_v.shape_value(i, point) //phi_i(x_q) + * JxW[point]; //dx + } + } + }; + + // This is the function called for boundary faces and consists of a normal + // integration using FEFaceValues. New is the logic to decide if the term + // goes into the system matrix (outflow) or the right-hand side (inflow). + const auto boundary_worker = [&](const Iterator &cell, + const unsigned int &face_no, ScratchData &scratch_data, + CopyData ©_data) + { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &fe_face = + scratch_data.fe_interface_values.get_fe_face_values(0); + + const auto &q_points = fe_face.get_quadrature_points(); + + const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell(); + const std::vector &JxW = fe_face.get_JxW_values(); + const std::vector> &normals = + fe_face.get_normal_vectors(); + + std::vector g(q_points.size()); + exact_solution.value_list(q_points, g); + + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double beta_dot_n = beta(q_points[point]) * normals[point]; + + if (beta_dot_n > 0) + { + for (unsigned int i = 0; i < n_facet_dofs; ++i) + for (unsigned int j = 0; j < n_facet_dofs; ++j) + copy_data.cell_matrix(i, j) += fe_face.shape_value(i, + point) // \phi_i + * fe_face.shape_value(j, point) // \phi_j + * beta_dot_n // \beta . n + * JxW[point]; // dx + } + else + for (unsigned int i = 0; i < n_facet_dofs; ++i) + copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i + * g[point] // g*/ + * beta_dot_n // \beta . n + * JxW[point]; // dx + } + }; + + // This is the function called on interior faces. The arguments specify + // cells, face and subface indices (for adaptive refinement). We just pass + // them along to the reinit() function of FEInterfaceValues. + const auto face_worker = [&](const Iterator &cell, const unsigned int &f, + const unsigned int &sf, const Iterator &ncell, + const unsigned int &nf, const unsigned int &nsf, + ScratchData &scratch_data, CopyData ©_data) + { + FEInterfaceValues &fe_iv = scratch_data.fe_interface_values; + fe_iv.reinit(cell, f, sf, ncell, nf, nsf); + const auto &q_points = fe_iv.get_quadrature_points(); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + + const unsigned int n_dofs = fe_iv.n_current_interface_dofs(); + copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices(); + + copy_data_face.cell_matrix.reinit(n_dofs, n_dofs); + + const std::vector &JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + { + const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint]; + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + copy_data_face.cell_matrix(i, j) += (beta(q_points[qpoint]) * normals[qpoint] * fe_iv.average(j, qpoint) * fe_iv.jump(i, qpoint) + + theta * std::abs(beta_dot_n) * fe_iv.jump(j, qpoint) * fe_iv.jump(i, qpoint)) * + JxW[qpoint]; + } + } + } + }; + + // The following lambda function will handle copying the data from the + // cell and face assembly into the global matrix and right-hand side. + // + // While we would not need an AffineConstraints object, because there are + // no hanging node constraints in DG discretizations, we use an empty + // object here as this allows us to use its `copy_local_to_global` + // functionality. + const AffineConstraints constraints; + + const auto copier = [&](const CopyData &c) + { + constraints.distribute_local_to_global(c.cell_matrix, c.cell_rhs, + c.local_dof_indices, system_matrix, right_hand_side); + + for (auto &cdf : c.face_data) + { + constraints.distribute_local_to_global(cdf.cell_matrix, + cdf.joint_dof_indices, system_matrix); + } + }; + + ScratchData scratch_data(mapping, *fe, quadrature, quadrature_face); + CopyData copy_data; + + // Here, we finally handle the assembly. We pass in ScratchData and + // CopyData objects, the lambda functions from above, an specify that we + // want to assemble interior faces once. + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), + cell_worker, copier, scratch_data, copy_data, + MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::assemble_own_interior_faces_once, + boundary_worker, face_worker); +} + +template +void AdvectionReaction::solve() +{ + + if (use_direct_solver) + { + + SparseDirectUMFPACK system_matrix_inverse; + system_matrix_inverse.initialize(system_matrix); + system_matrix_inverse.vmult(solution, right_hand_side); + } + else + { + // Here we have a classic iterative solver, as done in many tutorials: + SolverControl solver_control(1000, 1e-15); + SolverRichardson> solver(solver_control); + PreconditionBlockSSOR> preconditioner; + preconditioner.initialize(system_matrix, fe->n_dofs_per_cell()); + solver.solve(system_matrix, solution, right_hand_side, preconditioner); + std::cout << " Solver converged in " << solver_control.last_step() + << " iterations." << std::endl; + } +} + +// @sect3{Mesh refinement} +// We refine the grid according the proposed estimator or with an approximation to the gradient of the solution. +// The first option is the default one (you can see it in the header file) +template +void AdvectionReaction::refine_grid() +{ + + if (refinement == "residual") + { + + //If the `refinement` string is `"residual"`, then we first compute the local projection + compute_local_projection_and_estimate(); + //We then set the refinement fraction and as usual we execute the refinement. + const double refinement_fraction = 0.6; + GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, error_indicator_per_cell, refinement_fraction, 0.0); + triangulation.execute_coarsening_and_refinement(); + } + else if (refinement == "gradient") + { + + Vector gradient_indicator(triangulation.n_active_cells()); + + // Now the approximate gradients are computed + DerivativeApproximation::approximate_gradient(mapping, dof_handler, + solution, gradient_indicator); + + // and they are cell-wise scaled by the factor $h^{1+d/2}$ + unsigned int cell_no = 0; + for (const auto &cell : dof_handler.active_cell_iterators()) + gradient_indicator(cell_no++) *= std::pow(cell->diameter(), + 1 + 1.0 * dim / 2); + + // Finally they serve as refinement indicator. + GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, + gradient_indicator, 0.25, 0.0); + + triangulation.execute_coarsening_and_refinement(); + std::cout << gradient_indicator.l2_norm() << '\n'; + } + else if (refinement == "global") + { + triangulation.refine_global(1); //just for testing on uniformly refined meshes + } + else + { + Assert(false, ExcInternalError()); + } +} +// The output of this program consists of a vtk file of the adaptively +// refined grids and the numerical solutions. +template +void AdvectionReaction::output_results(const unsigned int cycle) const +{ + const std::string filename = "solution-" + std::to_string(cycle) + ".vtk"; + std::cout << " Writing solution to <" << filename << ">" << std::endl; + std::ofstream output(filename); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u", DataOut::type_dof_data); + + data_out.build_patches(mapping); + + data_out.write_vtk(output); +} + +template +void AdvectionReaction::compute_error() +{ + error_table.error_from_exact(mapping, dof_handler, solution, exact_solution); //be careful: a FD approximation of the gradient is used to compute the H^1 norm if you're not relying on SymbolicFunction class + // error_table.error_from_exact(mapping, dof_handler, solution, Solution()); //provided that Solution implements the Gradient function +} + +// @sect3{Compute the energy norm} +// The energy norm is defined as $ |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 + \sum_{F \in \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2 \Bigr)^{\frac{1}{2}}$ +// Notice that in the current case we have $c_f = \frac{|b \cdot n|}{2}$ +// Like in the assembly, all the contributions are handled separately by using ScratchData and CopyData objects. +template +double AdvectionReaction::compute_energy_norm() +{ + + energy_norm_square_per_cell.reinit(triangulation.n_active_cells()); + + using Iterator = typename DoFHandler::active_cell_iterator; + + // We start off by adding cell contributions + const auto cell_worker = [&](const Iterator &cell, + ScratchData &scratch_data, CopyData ©_data) + { + const unsigned int n_dofs = + scratch_data.fe_values.get_fe().n_dofs_per_cell(); + copy_data.reinit(cell, n_dofs); + scratch_data.fe_values.reinit(cell); + + copy_data.cell_index = cell->active_cell_index(); + + const auto &q_points = scratch_data.fe_values.get_quadrature_points(); + const FEValues &fe_v = scratch_data.fe_values; + const std::vector &JxW = fe_v.get_JxW_values(); + + double error_square_norm{0.0}; + std::vector sol_u(fe_v.n_quadrature_points); + fe_v.get_function_values(solution, sol_u); + + for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point) + { + const double diff = (sol_u[point] - exact_solution.value(q_points[point])); + error_square_norm += diff * diff * JxW[point]; + } + copy_data.value = error_square_norm; + }; + + // Here we add contributions coming from the internal faces + const auto face_worker = [&](const Iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData &scratch_data, + CopyData ©_data) + { + FEInterfaceValues &fe_iv = scratch_data.fe_interface_values; + fe_iv.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); + + const auto &q_points = fe_iv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + const std::vector &JxW = fe_iv.get_JxW_values(); + std::vector g(n_q_points); + + std::vector jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + const std::vector> &normals = fe_iv.get_normal_vectors(); + + double error_jump_square{0.0}; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const double beta_dot_n = theta * std::abs(beta(q_points[point]) * normals[point]); + error_jump_square += beta_dot_n * jump[point] * jump[point] * JxW[point]; + } + + copy_data.value = error_jump_square; + }; + + // Finally, we add the boundary contributions + const auto boundary_worker = [&](const Iterator &cell, + const unsigned int &face_no, + ScratchData &scratch_data, + CopyData ©_data) + { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &fe_fv = scratch_data.fe_interface_values.get_fe_face_values(0); + const auto &q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + const std::vector &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const std::vector> &normals = fe_fv.get_normal_vectors(); + + double difference_norm_square = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double beta_dot_n = theta * std::abs(beta(q_points[point]) * normals[point]); + const double diff = (boundary_conditions.value(q_points[point]) - sol_u[point]); + difference_norm_square += beta_dot_n * diff * diff * JxW[point]; + } + copy_data.value = difference_norm_square; + }; + + const auto copier = [&](const auto ©_data) + { + if (copy_data.cell_index != numbers::invalid_unsigned_int) + { + energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value; + } + for (auto &cdf : copy_data.face_data) + for (unsigned int j = 0; j < 2; ++j) + energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; + }; + + ScratchData scratch_data(mapping, *fe, QGauss{fe->tensor_degree() + 1}, + QGauss{fe->tensor_degree() + 1}); + + CopyData copy_data; + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + copy_data, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_own_interior_faces_once | + MeshWorker::assemble_boundary_faces, + boundary_worker, + face_worker); + + const double energy_error = std::sqrt(energy_norm_square_per_cell.l1_norm()); + return energy_error; +} + +// @sect3{Computing the estimator} +// In the estimator, we have to compute the term $||f- c u_h - \Pi(f- c u_h)||_{T}^{2}$ over a generic cell $T$. To achieve this, we first need to +// compute the projection involving the finite element function $u_h$. Using the definition of orthogonal projection, we're required to solve cellwise +// $(v_h,f-c u_h)_T = (v_h,\Pi)_T \qquad \forall v_h \in V_h$ for $\Pi$, which means that we have to build a mass-matrix on each cell. +// Once we have the projection, which is a finite element function, we can add its contribution in the cell_worker lambda. +// As done in step-74, the square of the error indicator is computed. +// +template +void AdvectionReaction::compute_local_projection_and_estimate() +{ + + // Compute the term $||f-c u_h - \Pi(f- cu_h)||_T^2$ + using Iterator = typename DoFHandler::active_cell_iterator; + error_indicator_per_cell.reinit(triangulation.n_active_cells()); + + const auto cell_worker = [&](const Iterator &cell, + ScratchData &scratch_data, CopyData ©_data) + { + const unsigned int n_dofs = scratch_data.fe_values.get_fe().n_dofs_per_cell(); + + copy_data.reinit(cell, n_dofs); + scratch_data.fe_values.reinit(cell); + copy_data.cell_index = cell->active_cell_index(); + + const auto &q_points = scratch_data.fe_values.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + + const FEValues &fe_v = scratch_data.fe_values; + const std::vector &JxW = fe_v.get_JxW_values(); + + std::vector sol_u_at_quadrature_points(fe_v.n_quadrature_points); + fe_v.get_function_values(solution, sol_u_at_quadrature_points); + + //Compute local L^2 projection of $f- c u_h$ over the local finite element space + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < n_dofs; ++i) + { + for (unsigned int j = 0; j < n_dofs; ++j) + { + + copy_data.cell_mass_matrix(i, j) += fe_v.shape_value(i, point) * //phi_i(x_q) + fe_v.shape_value(j, point) * //phi_j(x_q) + JxW[point]; // dx(x_q) + } + copy_data.cell_mass_rhs(i) += + (rhs.value(q_points[point]) * // f(x_q) + fe_v.shape_value(i, point) //phi_i(x_q) + - + advection_coeff.value(q_points[point]) * + fe_v.shape_value(i, point) * //c*phi_i(x_q) + sol_u_at_quadrature_points[point]) * //u_h(x_q) + JxW[point]; //dx + } + } + + FullMatrix inverse(fe_v.n_quadrature_points, fe_v.n_quadrature_points); + inverse.invert(copy_data.cell_mass_matrix); + Vector proj(fe_v.n_quadrature_points); //projection of (f-c*U_h) on the local fe_space + inverse.vmult(proj, copy_data.cell_mass_rhs); //M^{-1}*rhs = proj + + double square_norm_over_cell = 0.0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const double diff = rhs.value(q_points[point]) - sol_u_at_quadrature_points[point] - proj[point]; + square_norm_over_cell += diff * diff * JxW[point]; + } + copy_data.value_estimator = square_norm_over_cell; + }; + + // Finally we have the boundary term with $||\beta (g-u_h^+)||^2$ + const auto boundary_worker = [&](const Iterator &cell, + const unsigned int &face_no, + ScratchData &scratch_data, + CopyData ©_data) + { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &fe_fv = scratch_data.fe_interface_values.get_fe_face_values(0); + const auto &q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + const std::vector &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + exact_solution.value_list(q_points, g); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const std::vector> &normals = fe_fv.get_normal_vectors(); + + double square_norm_over_bdary_face = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double beta_dot_n = beta(q_points[point]) * normals[point]; + + if (beta_dot_n < 0) //\partial_{-T} \cap \partial_{- \Omega} + { + const double diff = std::abs(beta_dot_n) * (g[point] - sol_u[point]); + square_norm_over_bdary_face += diff * diff * JxW[point]; + } + } + copy_data.value_estimator += square_norm_over_bdary_face; + }; + + // Then compute the interior face terms with $|| \sqrt{b \cdot n}[u_h]||^2$ + const auto face_worker = [&](const Iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData &scratch_data, + CopyData ©_data) + { + FEInterfaceValues &fe_iv = scratch_data.fe_interface_values; + fe_iv.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); + + const auto &q_points = fe_iv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + + const std::vector &JxW = fe_iv.get_JxW_values(); + std::vector g(n_q_points); + + std::vector jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + const std::vector> &normals = fe_iv.get_normal_vectors(); + + double error_jump_square{0.0}; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const double beta_dot_n = beta(q_points[point]) * normals[point]; + if (beta_dot_n < 0) + { + error_jump_square += std::abs(beta_dot_n) * jump[point] * jump[point] * JxW[point]; + } + } + + copy_data_face.values[0] = error_jump_square; + copy_data_face.values[1] = copy_data_face.values[0]; + }; + + ScratchData scratch_data(mapping, *fe, QGauss{fe->tensor_degree() + 1}, + QGauss{fe->tensor_degree() + 1}); + + const auto copier = [&](const auto ©_data) + { + if (copy_data.cell_index != numbers::invalid_unsigned_int) + { + error_indicator_per_cell[copy_data.cell_index] += copy_data.value_estimator; + } + for (auto &cdf : copy_data.face_data) + { + for (unsigned int j = 0; j < 2; ++j) + { + error_indicator_per_cell[cdf.cell_indices[j]] += cdf.values[j]; + } + } + }; + + // Here, we finally handle the assembly of the Mass matrix (M)_{ij} = (\phi_j, \phi_i)_T. We pass in ScratchData and + // CopyData objects + CopyData copy_data; + MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), + cell_worker, copier, scratch_data, copy_data, + MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::assemble_own_interior_faces_once, + boundary_worker, face_worker); +} + +//Usual run function, which runs over several refinement cycles +template +void AdvectionReaction::run() +{ + std::vector energy_errors; + std::vector dofs_hist; + for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle) + { + std::cout << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(n_global_refinements); + } + else + { + refine_grid(); + } + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + assemble_system(); + solve(); + compute_error(); + output_results(cycle); + energy_errors.emplace_back(compute_energy_norm()); + dofs_hist.emplace_back(triangulation.n_active_cells()); + } + error_table.output_table(std::cout); + + for (unsigned int i = 0; i < n_refinement_cycles; ++i) + std::cout << "Cycle " << i << "\t" << energy_errors[i] << '\n'; + { + } +} +// Explicit instantiation +template class AdvectionReaction<2>;