From: Marco Feder Date: Wed, 20 Jul 2022 17:38:15 +0000 (+0200) Subject: Clean a little and fix indentation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F114%2Fhead;p=code-gallery.git Clean a little and fix indentation --- diff --git a/advection_reaction_estimator/README.md b/advection_reaction_estimator/README.md index b33d794..f243a2a 100644 --- a/advection_reaction_estimator/README.md +++ b/advection_reaction_estimator/README.md @@ -133,15 +133,15 @@ This solution has an internal layer along the line $y=\frac{1}{2} -x$, hence we The next image is the 3D view of the numerical solution: -![Screenshot](doc/images/warp_by_scalar_solution_layer.png) +![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) +![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) +![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 diff --git a/advection_reaction_estimator/include/DG_advection_reaction.h b/advection_reaction_estimator/include/DG_advection_reaction.h index 0146eab..24c5235 100644 --- a/advection_reaction_estimator/include/DG_advection_reaction.h +++ b/advection_reaction_estimator/include/DG_advection_reaction.h @@ -1,124 +1,146 @@ #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 #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 +// 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 + +// To enable parameter handling #include #include -#include #include +#include #include #include #include #include -#include #include +#include +using namespace dealii; -//This is a struct used only for throwing an exception when theta parameter is not okay. +// 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(); } + 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(); + 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 + 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_advection_reaction_estimator"; + std::map constants; + ParsedConvergenceTable error_table; + + bool use_direct_solver = true; + unsigned int n_refinement_cycles = 8; + 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 index 91eac65..97407d3 100644 --- a/advection_reaction_estimator/main.cc +++ b/advection_reaction_estimator/main.cc @@ -1,62 +1,62 @@ #include "include/DG_advection_reaction.h" -int main(int argc, char **argv) +int +main(int argc, char **argv) { - try - { - std::string par_name = ""; - if (argc > 1) { - par_name = argv[1]; + std::string par_name = ""; + if (argc > 1) + { + par_name = argv[1]; + } + else + { + par_name = "parameters.prm"; + } + deallog.depth_console(2); + AdvectionReaction<2> problem; + problem.initialize_params(par_name); + problem.run(); } - deallog.depth_console(2); //solver infos - AdvectionReaction<2> dgmethod; - if (par_name != "") + catch (std::exception &exc) { - dgmethod.initialize_params(par_name); + 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; } - - 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; - } + { + 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; - } + { + 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 index 7a69bd7..33cc939 100644 --- a/advection_reaction_estimator/source/DG_advection_reaction.cc +++ b/advection_reaction_estimator/source/DG_advection_reaction.cc @@ -12,27 +12,25 @@ * 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 +// Compute and returns the wind field b template -Tensor<1, dim> beta(const Point &p) +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; + Assert(dim > 1, 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; + return wind_field; } // @sect3{The ScratchData and CopyData classes} @@ -45,781 +43,882 @@ Tensor<1, dim> beta(const Point &p) 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; + 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; + 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); - } + 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) +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]; + 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) +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(); + Assert(dim > 1, ExcMessage("Not implemented in 1D.")); + 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) +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")); - } + 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) +void +AdvectionReaction::parse_string(const std::string ¶meters) { - ParameterAcceptor::prm.parse_input_from_string(parameters); - ParameterAcceptor::parse_all_parameters(); + ParameterAcceptor::prm.parse_input_from_string(parameters); + ParameterAcceptor::parse_all_parameters(); } + + template -void AdvectionReaction::setup_system() +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()); + // 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. + + +// 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() +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_of_shape_values(j, qpoint) * fe_iv.jump_in_shape_values(i, qpoint) + - theta * std::abs(beta_dot_n) * fe_iv.jump_in_shape_values(j, qpoint) * fe_iv.jump_in_shape_values(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); + 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_of_shape_values(j, qpoint) * + fe_iv.jump_in_shape_values(i, qpoint) + + theta * std::abs(beta_dot_n) * + fe_iv.jump_in_shape_values(j, qpoint) * + fe_iv.jump_in_shape_values(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() +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; - } + 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) +// 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() +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()); - } + 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 +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); + 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() +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 + 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 Solution doesn't + // 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. +// 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() +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; + 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. +// 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() +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); + // 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 + + +// Usual run function, which runs over several refinement cycles template -void AdvectionReaction::run() +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'; - { - } + 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; + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + setup_system(); + 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] << std::endl; } // Explicit instantiation +template class AdvectionReaction<1>; template class AdvectionReaction<2>; +template class AdvectionReaction<3>;