From 88d865167aaceba18ffeee838085c780eee40502 Mon Sep 17 00:00:00 2001 From: tcclevenger Date: Tue, 19 Feb 2019 16:15:20 -0700 Subject: [PATCH] add step-63.cc (no doc) --- examples/step-63/block_jacobi.prm | 14 + examples/step-63/block_sor.prm | 14 + examples/step-63/jacobi.prm | 14 + examples/step-63/sor.prm | 15 + examples/step-63/step-63.cc | 1082 ++++++++++++++++++++++++++++- 5 files changed, 1135 insertions(+), 4 deletions(-) create mode 100644 examples/step-63/block_jacobi.prm create mode 100644 examples/step-63/block_sor.prm create mode 100644 examples/step-63/jacobi.prm create mode 100644 examples/step-63/sor.prm diff --git a/examples/step-63/block_jacobi.prm b/examples/step-63/block_jacobi.prm new file mode 100644 index 0000000000..a8dc7f5566 --- /dev/null +++ b/examples/step-63/block_jacobi.prm @@ -0,0 +1,14 @@ +# Finite Element degree +set fe degree = 1 + +# Smoother Type: sor|jacobi|block sor|block jacobi +set smoother type = block jacobi + +# Dof renumbering: no|random|downstream|upstream +set dof renumbering = downstream + +# With streamline diffusion: true|false +set with sd = true + +# Output: true|false +set output = false diff --git a/examples/step-63/block_sor.prm b/examples/step-63/block_sor.prm new file mode 100644 index 0000000000..f5991834c4 --- /dev/null +++ b/examples/step-63/block_sor.prm @@ -0,0 +1,14 @@ +# Finite Element degree +set fe degree = 1 + +# Smoother Type: sor|jacobi|block sor|block jacobi +set smoother type = block sor + +# Dof renumbering: none|random|downstream|upstream +set dof renumbering = downstream + +# With streamline diffusion: true|false +set with sd = true + +# Output: true|false +set output = false diff --git a/examples/step-63/jacobi.prm b/examples/step-63/jacobi.prm new file mode 100644 index 0000000000..d6a55ff53b --- /dev/null +++ b/examples/step-63/jacobi.prm @@ -0,0 +1,14 @@ +# Finite Element degree +set fe degree = 1 + +# Smoother Type: sor|jacobi|block sor|block jacobi +set smoother type = jacobi + +# Dof renumbering: no|random|downstream|upstream +set dof renumbering = downstream + +# With streamline diffusion: true|false +set with sd = true + +# Output: true|false +set output = false diff --git a/examples/step-63/sor.prm b/examples/step-63/sor.prm new file mode 100644 index 0000000000..47ab13c3d5 --- /dev/null +++ b/examples/step-63/sor.prm @@ -0,0 +1,15 @@ +# Finite Element degree +set fe degree = 1 + +# Smoother Type: sor|jacobi|block sor|block jacobi +set smoother type = sor + +# Dof renumbering: no|random|downstream|upstream +set dof renumbering = downstream + + +# With streamline diffusion: true|false +set with sd = true + +# Output: true|false +set output = false diff --git a/examples/step-63/step-63.cc b/examples/step-63/step-63.cc index 2534596dc6..dcfcd1fa54 100644 --- a/examples/step-63/step-63.cc +++ b/examples/step-63/step-63.cc @@ -19,16 +19,1091 @@ // @note: This is work in progress and will be an example for block smoothers // in geometric multigrid. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + + + +#include #include -#include -int main() + +#include +#include + + + +namespace Step100 +{ + using namespace dealii; + + + template + struct ScratchData + { + ScratchData(const FiniteElement &fe, + const unsigned int quadrature_degree) + : fe_values(fe, + QGauss(quadrature_degree), + update_values | update_gradients | update_hessians | + update_quadrature_points | update_JxW_values) + {} + + ScratchData(const ScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + update_values | update_gradients | update_hessians | + update_quadrature_points | update_JxW_values) + {} + + FEValues fe_values; + }; + struct CopyData + { + unsigned int level; + unsigned int dofs_per_cell; + + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + }; + + + struct Settings + { + bool try_parse(const std::string &prm_filename); + + unsigned int fe_degree; + std::string smoother_type; + std::string dof_renum; + bool with_sd; + bool output; + }; + + + bool Settings::try_parse(const std::string &prm_filename) + { + ParameterHandler prm; + + prm.declare_entry("fe degree", + "1", + Patterns::Integer(0), + "Finite Element degree"); + prm.declare_entry("smoother type", + "block sor", + Patterns::Selection("sor|jacobi|block sor|block jacobi"), + "Smoother Type: sor|jacobi|block sor|block jacobi"); + prm.declare_entry("dof renumbering", + "downstream", + Patterns::Selection("none|random|downstream|upstream"), + "Dof renumbering: none|random|downstream|upstream"); + prm.declare_entry("with sd", + "true", + Patterns::Bool(), + "With streamline diffusion: true|false"); + prm.declare_entry("output", + "true", + Patterns::Bool(), + "Generate graphical output: true|false"); + + try + { + prm.parse_input(prm_filename); + } + catch (const dealii::PathSearch::ExcFileNotFound &) + { + if (prm_filename.size() > 0) + std::cerr << "ERRROR: could not open the .prm file '" << prm_filename + << "'" << std::endl; + else + std::cerr << "Usage: please pass a .prm file as the first argument" + << std::endl; + + prm.print_parameters(std::cout, ParameterHandler::Text); + return false; + } + this->fe_degree = prm.get_integer("fe degree"); + this->smoother_type = prm.get("smoother type"); + this->dof_renum = prm.get("dof renumbering"); + this->with_sd = prm.get_bool("with sd"); + this->output = prm.get_bool("output"); + + return true; + } + + + + namespace + { + template + struct CompareDownstream + { + /** + * Constructor. + */ + CompareDownstream(const Tensor<1, dim> &dir) + : dir(dir) + {} + /** + * Return true if c1 less c2. + */ + bool operator()(const Iterator &c1, const Iterator &c2) const + { + const Tensor<1, dim> diff = c2->center() - c1->center(); + return (diff * dir > 0); + } + + private: + /** + * Flow direction. + */ + const Tensor<1, dim> dir; + }; + + // Functions for creating permutation of cells for output and Block + // smoothers + template + std::vector + create_downstream_order(const DoFHandler &dof, + const Tensor<1, dim> direction, + const unsigned int level) + { + std::vector::level_cell_iterator> ordered_cells; + ordered_cells.reserve(dof.get_triangulation().n_cells(level)); + const CompareDownstream::level_cell_iterator, + dim> + comparator(direction); + + typename DoFHandler::level_cell_iterator cell = dof.begin(level); + typename DoFHandler::level_cell_iterator endc = dof.end(level); + for (; cell != endc; ++cell) + ordered_cells.push_back(cell); + + std::sort(ordered_cells.begin(), ordered_cells.end(), comparator); + + std::vector ordered_indices; + ordered_indices.reserve(dof.get_triangulation().n_cells(level)); + + for (unsigned int i = 0; i < ordered_cells.size(); ++i) + ordered_indices.push_back(ordered_cells[i]->index()); + + return ordered_indices; + } + + template + std::vector + create_downstream_order(const DoFHandler &dof, + const Tensor<1, dim> direction) + { + std::vector::active_cell_iterator> ordered_cells; + ordered_cells.reserve(dof.get_triangulation().n_active_cells()); + const CompareDownstream::active_cell_iterator, + dim> + comparator(direction); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(); + typename DoFHandler::active_cell_iterator endc = dof.end(); + for (; cell != endc; ++cell) + ordered_cells.push_back(cell); + + std::sort(ordered_cells.begin(), ordered_cells.end(), comparator); + + std::vector ordered_indices; + ordered_indices.reserve(dof.get_triangulation().n_active_cells()); + + for (unsigned int i = 0; i < ordered_cells.size(); ++i) + ordered_indices.push_back(ordered_cells[i]->index()); + + return ordered_indices; + } + + + + template + std::vector create_random_order(const DoFHandler &dof, + const unsigned int level) + { + const unsigned int n_cells = dof.get_triangulation().n_cells(level); + + std::vector ordered_cells; + ordered_cells.reserve(n_cells); + + typename DoFHandler::cell_iterator cell = dof.begin(level); + typename DoFHandler::cell_iterator endc = dof.end(level); + for (; cell != endc; ++cell) + ordered_cells.push_back(cell->index()); + + // shuffle the elements; the following is essentially std::shuffle (which + // is new in C++11) but with a boost URNG + ::boost::mt19937 random_number_generator; + for (unsigned int i = 1; i < n_cells; ++i) + { + // get a random number between 0 and i (inclusive) + const unsigned int j = + ::boost::random::uniform_int_distribution<>(0, i)( + random_number_generator); + + // if possible, swap the elements + if (i != j) + std::swap(ordered_cells[i], ordered_cells[j]); + } + + return ordered_cells; + } + + template + std::vector create_random_order(const DoFHandler &dof) + { + const unsigned int n_cells = dof.get_triangulation().n_active_cells(); + + std::vector ordered_cells; + ordered_cells.reserve(n_cells); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(); + typename DoFHandler::active_cell_iterator endc = dof.end(); + for (; cell != endc; ++cell) + ordered_cells.push_back(cell->index()); + + // shuffle the elements; the following is essentially std::shuffle (which + // is new in C++11) but with a boost URNG + ::boost::mt19937 random_number_generator; + for (unsigned int i = 1; i < n_cells; ++i) + { + // get a random number between 0 and i (inclusive) + const unsigned int j = + ::boost::random::uniform_int_distribution<>(0, i)( + random_number_generator); + + // if possible, swap the elements + if (i != j) + std::swap(ordered_cells[i], ordered_cells[j]); + } + + return ordered_cells; + } + + } // namespace + + + + // RHS and boundary is an adaptation from one in + // Finite Elements and Fast Iterative Solvers: with Applications + // in Incompressible Fluid Dynamics + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function() + {} + + virtual double value(const Point & p, + const unsigned int component = 0) const; + + virtual void value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; + + private: + static const Point center_point; + }; + + template + double RightHandSide::value(const Point & p, + const unsigned int component) const + { + Assert(component == 0, ExcIndexRange(component, 0, 1)); + (void)component; + (void)p; + + return 0.0; + } + + + + template + void RightHandSide::value_list(const std::vector> &points, + std::vector & values, + const unsigned int component) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < points.size(); ++i) + values[i] = RightHandSide::value(points[i], component); + } + + + + template + class BoundaryValues : public Function + { + public: + BoundaryValues() + : Function() + {} + + virtual double value(const Point & p, + const unsigned int component = 0) const; + + virtual void value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; + }; + + + template + double BoundaryValues::value(const Point & p, + const unsigned int component) const + { + Assert(component == 0, ExcIndexRange(component, 0, 1)); + (void)component; + + if (std::fabs(p[0] * p[0] + p[1] * p[1] - 0.3 * 0.3) < + 1e-8 // around cylinder + || std::fabs(p[0] + 1) < 1e-8 // x == -1 + || std::fabs(p[1] - 1) < 1e-8 // y == 1 + || (std::fabs(p[1] + 1) < 1e-8 && p[0] < 0.5) // y == -1, x <= 0.5 + ) + { + return 0.0; + } + else if (std::fabs(p[0] - 1) < 1e-8 // x = 1 + || (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5) // y == -1, x > 0.5 + ) + { + return 1.0; + } + else + { + return 0.0; + } + } + + + + template + void BoundaryValues::value_list(const std::vector> &points, + std::vector & values, + const unsigned int component) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < points.size(); ++i) + values[i] = BoundaryValues::value(points[i], component); + } + + template + double delta_value(const double hk, + const double eps, + const Tensor<1, dim> dir, + const double pk) + { + // Value defined in 'On discontinuity–capturing methods for + // convection–diffusion equations' + double Peclet = dir.norm() * hk / (2.0 * eps * pk); + double coth = + (1.0 + std::exp(-2.0 * Peclet)) / (1.0 - std::exp(-2.0 * Peclet)); + + return hk / (2.0 * dir.norm() * pk) * (coth - 1.0 / Peclet); + } + + + template + class AdvectionProblem + { + public: + AdvectionProblem(Settings settings); + void run(); + + private: + void setup_system(); + + template + void assemble_cell(const IteratorType &cell, + ScratchData & scratch_data, + CopyData & copy_data); + void assemble_system_and_multigrid(); + + std::unique_ptr>> create_smoother(); + + void solve(); + void refine_grid(); + void output_results(const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + + FE_Q fe; + MappingQ mapping; + unsigned int quad_degree; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject mg_interface_sparsity_patterns; + + MGLevelObject> mg_matrices; + MGLevelObject> mg_interface_in; + MGLevelObject> mg_interface_out; + + MGConstrainedDoFs mg_constrained_dofs; + + Settings settings; + const double epsilon; + Tensor<1, dim> advection_direction; + }; + + + + template + AdvectionProblem::AdvectionProblem(Settings settings) + : triangulation(Triangulation::limit_level_difference_at_vertices) + , dof_handler(triangulation) + , fe(settings.fe_degree) + , mapping(settings.fe_degree) + , quad_degree(2 * fe.degree + 2) + , settings(settings) + , epsilon(0.005) + { + // Set Advection direction (problem is an adaptation from one in + // Finite Elements and Fast Iterative Solvers: with Applications + // in Incompressible Fluid Dynamics) + advection_direction[0] = -std::sin(numbers::PI / 6.0); + if (dim > 1) + advection_direction[1] = std::cos(numbers::PI / 6.0); + if (dim > 2) + advection_direction[2] = std::sin(numbers::PI / 6.0); + } + + + + template + void AdvectionProblem::setup_system() + { + const unsigned int n_levels = triangulation.n_levels(); + + dof_handler.distribute_dofs(fe); + + // We could renumber the active dofs with DoFRenumbering::downstream() + // here, but the smoothers only act on multigrid levels and as such, this + // wouldn't matter. Instead, we will renumber the DoFs on each multigrid + // level below. + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, BoundaryValues(), constraints); + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 1, BoundaryValues(), constraints); + constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /*keep_constrained_dofs = */ false); + + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + + // Setup GMG DoFs + dof_handler.distribute_mg_dofs(); + + // Renumber DoFs on each level in downstream or upstream direction if + // needed. This is only necessary for point smoothers (SOR and Jacobi) as + // the block smoothers operate on cells (see create_smoother()): + if (settings.smoother_type == "sor" || settings.smoother_type == "jacobi") + { + if (settings.dof_renum == "downstream" || + settings.dof_renum == "upstream") + { + const Tensor<1, dim> direction = + (settings.dof_renum == "upstream" ? -1.0 : 1.0) * + advection_direction; + + for (unsigned int level = 0; level < n_levels; ++level) + DoFRenumbering::downstream(dof_handler, + level, + direction, + /*dof_wise_renumbering = */ true); + } + else if (settings.dof_renum == "random") + { + for (unsigned int level = 0; level < n_levels; ++level) + DoFRenumbering::random(dof_handler, level); + } + else + Assert(false, ExcNotImplemented()); + } + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(dof_handler); + + std::set dirichlet_boundary_ids = {0, 1}; + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary_ids); + + mg_matrices.resize(0, n_levels - 1); + mg_matrices.clear_elements(); + mg_interface_in.resize(0, n_levels - 1); + mg_interface_in.clear_elements(); + mg_interface_out.resize(0, n_levels - 1); + mg_interface_out.clear_elements(); + mg_sparsity_patterns.resize(0, n_levels - 1); + mg_interface_sparsity_patterns.resize(0, n_levels - 1); + + for (unsigned int level = 0; level < n_levels; ++level) + { + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern(dof_handler, dsp, level); + mg_sparsity_patterns[level].copy_from(dsp); + mg_matrices[level].reinit(mg_sparsity_patterns[level]); + } + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_interface_sparsity_pattern(dof_handler, + mg_constrained_dofs, + dsp, + level); + mg_interface_sparsity_patterns[level].copy_from(dsp); + + // We need both interface in and out matrices since our problem is not + // symmetric + mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]); + mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]); + } + } + } + + + template + template + void AdvectionProblem::assemble_cell(const IteratorType &cell, + ScratchData & scratch_data, + CopyData & copy_data) + { + const unsigned int level = cell->level(); + copy_data.level = level; + + const unsigned int dofs_per_cell = + scratch_data.fe_values.get_fe().dofs_per_cell; + copy_data.dofs_per_cell = dofs_per_cell; + + const unsigned int n_q_points = + scratch_data.fe_values.get_quadrature().size(); + copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + + if (!cell->is_level_cell()) + copy_data.cell_rhs.reinit(dofs_per_cell); + + copy_data.local_dof_indices.resize(dofs_per_cell); + cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices); + + scratch_data.fe_values.reinit(cell); + + const RightHandSide right_hand_side; + std::vector rhs_values(n_q_points); + + right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(), + rhs_values); + + const double delta = settings.with_sd ? delta_value(cell->diameter(), + epsilon, + advection_direction, + settings.fe_degree) : + 0.0; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + copy_data.cell_matrix(i, j) += + (epsilon * scratch_data.fe_values.shape_grad(j, q_point) * + scratch_data.fe_values.shape_grad(i, q_point) * + scratch_data.fe_values.JxW(q_point)) + + ((advection_direction * + scratch_data.fe_values.shape_grad(j, q_point)) * + scratch_data.fe_values.shape_value(i, q_point)) * + scratch_data.fe_values.JxW(q_point); + + if (settings.with_sd) + copy_data.cell_matrix(i, j) += + delta * + (advection_direction * + scratch_data.fe_values.shape_grad(j, q_point)) * + (advection_direction * + scratch_data.fe_values.shape_grad(i, q_point)) * + scratch_data.fe_values.JxW(q_point) - + delta * epsilon * + trace(scratch_data.fe_values.shape_hessian(j, q_point)) * + (advection_direction * + scratch_data.fe_values.shape_grad(i, q_point)) * + scratch_data.fe_values.JxW(q_point); + } + if (!cell->is_level_cell()) + { + copy_data.cell_rhs(i) += + scratch_data.fe_values.shape_value(i, q_point) * + rhs_values[q_point] * scratch_data.fe_values.JxW(q_point); + if (settings.with_sd) + copy_data.cell_rhs(i) += + delta * rhs_values[q_point] * advection_direction * + scratch_data.fe_values.shape_grad(i, q_point) * + scratch_data.fe_values.JxW(q_point); + } + } + } + + + template + void AdvectionProblem::assemble_system_and_multigrid() + { + auto cell_worker_active = + [&](const decltype(dof_handler.begin_active()) &cell, + ScratchData & scratch_data, + CopyData & copy_data) { + this->assemble_cell(cell, scratch_data, copy_data); + }; + + + + auto copier_active = [&](const CopyData ©_data) { + constraints.distribute_local_to_global(copy_data.cell_matrix, + copy_data.cell_rhs, + copy_data.local_dof_indices, + system_matrix, + system_rhs); + }; + + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker_active, + copier_active, + ScratchData(fe, quad_degree), + CopyData(), + MeshWorker::assemble_own_cells); + + + // Assemble GMG + std::vector> boundary_constraints( + triangulation.n_global_levels()); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + IndexSet dofset; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + level, + dofset); + boundary_constraints[level].reinit(dofset); + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_refinement_edge_indices(level)); + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_boundary_indices(level)); + boundary_constraints[level].close(); + } + + auto cell_worker_mg = [&](const decltype(dof_handler.begin_mg()) &cell, + ScratchData &scratch_data, + CopyData & copy_data) { + this->assemble_cell(cell, scratch_data, copy_data); + }; + + auto copier_mg = [&](const CopyData ©_data) { + boundary_constraints[copy_data.level].distribute_local_to_global( + copy_data.cell_matrix, + copy_data.local_dof_indices, + mg_matrices[copy_data.level]); + + // If (i,j) is an interface_out dof pair, then (j,i) is an interface_in + // dof pair. Note: for interface_in, we load the transpose of the + // interface entries, i.e., the entry for dof pair (j,i) is stored in + // interface_in(i,j). + for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i) + for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j) + if (mg_constrained_dofs.is_interface_matrix_entry( + copy_data.level, + copy_data.local_dof_indices[i], + copy_data.local_dof_indices[j])) + { + mg_interface_out[copy_data.level].add( + copy_data.local_dof_indices[i], + copy_data.local_dof_indices[j], + copy_data.cell_matrix(i, j)); + mg_interface_in[copy_data.level].add( + copy_data.local_dof_indices[i], + copy_data.local_dof_indices[j], + copy_data.cell_matrix(j, i)); + } + }; + + MeshWorker::mesh_loop(dof_handler.begin_mg(), + dof_handler.end_mg(), + cell_worker_mg, + copier_mg, + ScratchData(fe, quad_degree), + CopyData(), + MeshWorker::assemble_own_cells); + } + + + template + std::unique_ptr>> + AdvectionProblem::create_smoother() + { + if (settings.smoother_type == "sor") + { + typedef PreconditionSOR> Smoother; + + auto smoother = + std_cxx14::make_unique, + Smoother, + Vector>>(); + smoother->initialize(mg_matrices, + Smoother::AdditionalData(fe.degree == 1 ? 1.0 : + 0.62)); + smoother->set_steps(2); + return smoother; + } + else if (settings.smoother_type == "jacobi") + { + typedef PreconditionJacobi> Smoother; + auto smoother = + std_cxx14::make_unique, + Smoother, + Vector>>(); + smoother->initialize(mg_matrices, + Smoother::AdditionalData(fe.degree == 1 ? 0.6667 : + 0.47)); + smoother->set_steps(4); + return smoother; + } + else if (settings.smoother_type == "block sor") + { + typedef RelaxationBlockSOR, double, Vector> + Smoother; + + static MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_levels() - 1); + + for (unsigned int level = 0; level < triangulation.n_levels(); ++level) + { + DoFTools::make_cell_patches(smoother_data[level].block_list, + dof_handler, + level); + + smoother_data[level].relaxation = 1.0; + smoother_data[level].inversion = PreconditionBlockBase::svd; + + std::vector ordered_indices; + if (settings.dof_renum == "downstream") + ordered_indices = create_downstream_order(dof_handler, + advection_direction, + level); + else if (settings.dof_renum == "upstream") + ordered_indices = + create_downstream_order(dof_handler, + -1.0 * advection_direction, + level); + else if (settings.dof_renum == "random") + ordered_indices = create_random_order(dof_handler, level); + else if (settings.dof_renum == "none") + { + // Do nothing + } + else + AssertThrow(false, ExcNotImplemented()); + + smoother_data[level].order = + std::vector>(1, ordered_indices); + } + + auto smoother = + std_cxx14::make_unique, + Smoother, + Vector>>(); + smoother->initialize(mg_matrices, smoother_data); + smoother->set_steps(1); + return smoother; + } + else if (settings.smoother_type == "block jacobi") + { + typedef RelaxationBlockJacobi, + double, + Vector> + Smoother; + + static MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_levels() - 1); + + for (unsigned int level = 0; level < triangulation.n_levels(); ++level) + { + DoFTools::make_cell_patches(smoother_data[level].block_list, + dof_handler, + level); + + smoother_data[level].relaxation = 0.25; + smoother_data[level].inversion = PreconditionBlockBase::svd; + + std::vector ordered_indices; + if (settings.dof_renum == "downstream") + ordered_indices = create_downstream_order(dof_handler, + advection_direction, + level); + else if (settings.dof_renum == "upstream") + ordered_indices = + create_downstream_order(dof_handler, + -1.0 * advection_direction, + level); + else if (settings.dof_renum == "random") + ordered_indices = create_random_order(dof_handler, level); + else if (settings.dof_renum == "none") + { + // Do nothing + } + else + AssertThrow(false, ExcNotImplemented()); + + smoother_data[level].order = + std::vector>(1, ordered_indices); + } + + auto smoother = + std_cxx14::make_unique, + Smoother, + Vector>>(); + smoother->initialize(mg_matrices, smoother_data); + smoother->set_steps(2); + return smoother; + } + else + AssertThrow(false, ExcNotImplemented()); + } + + + template + void AdvectionProblem::solve() + { + Timer time; + + const double solve_tol = 1e-8 * system_rhs.l2_norm(); + const unsigned int max_iters = 200; + SolverControl solver_control(max_iters, solve_tol, true, true); + solver_control.enable_history_data(); + + typedef MGTransferPrebuilt> Transfer; + Transfer mg_transfer(mg_constrained_dofs); + mg_transfer.build_matrices(dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from(mg_matrices[0]); + MGCoarseGridHouseholder> coarse_grid_solver; + coarse_grid_solver.initialize(coarse_matrix); + + std::unique_ptr>> mg_smoother = create_smoother(); + + mg::Matrix> mg_matrix(mg_matrices); + mg::Matrix> mg_interface_matrix_in(mg_interface_in); + mg::Matrix> mg_interface_matrix_out(mg_interface_out); + + Multigrid> mg( + mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother); + mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in); + + PreconditionMG, Transfer> preconditioner(dof_handler, + mg, + mg_transfer); + + + std::cout << " Solving with GMRES to tol " << solve_tol << "..." + << std::endl; + SolverGMRES<> solver(solver_control); + + time.restart(); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + time.stop(); + + std::cout << " converged in " << solver_control.last_step() + << " iterations" + << " in " << time.last_wall_time() << " seconds " << std::endl; + + constraints.distribute(solution); + } + + + + template + void AdvectionProblem::output_results(const unsigned int cycle) const + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + + // Here we generate an index for each cell to visualize the ordering used + // by the smoothers. Note that we do this only for the active cells + // instead of the levels, where the smoothers are actually used. For the + // point smoothers we renumber DoFs instead of cells, so this is only an + // approximation of what happens in reality. Finally, the random ordering + // is not the random ordering we actually use (see create_smoother() for + // that). + const unsigned int n_active_cells = triangulation.n_active_cells(); + Vector cell_indices(n_active_cells); + + { + // First generate a permutation vector for the cell indices: + std::vector ordered_indices; + if (settings.dof_renum == "downstream") + { + ordered_indices = + create_downstream_order(dof_handler, advection_direction); + } + else if (settings.dof_renum == "upstream") + { + ordered_indices = + create_downstream_order(dof_handler, -1.0 * advection_direction); + } + else if (settings.dof_renum == "random") + { + ordered_indices = create_random_order(dof_handler); + } + else if (settings.dof_renum == "none") + { + ordered_indices.resize(n_active_cells); + for (unsigned int i = 0; i < n_active_cells; ++i) + ordered_indices[i] = i; + } + else + AssertThrow(false, ExcNotImplemented()); + + // Then copy the permutation in ordered_indices into an output vector: + for (unsigned int i = 0; i < n_active_cells; ++i) + cell_indices(ordered_indices[i]) = static_cast(i); + } + + data_out.add_data_vector(cell_indices, "cell_index"); + + data_out.build_patches(); + + std::string filename = + "solution-" + Utilities::int_to_string(cycle) + ".vtu"; + std::ofstream output(filename.c_str()); + data_out.write_vtu(output); + } + + + template + void AdvectionProblem::run() + { + for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7 : 5); + ++cycle) + { + std::cout << " Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube_with_cylindrical_hole( + triangulation, 0.3, 1.0, 0.5, 1, false); + static const SphericalManifold manifold_description( + Point(0, 0)); + triangulation.set_manifold(1, manifold_description); + } + + triangulation.refine_global(); + + setup_system(); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << " (" + << triangulation.n_levels() << " levels)" << std::endl; + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() << std::endl; + + assemble_system_and_multigrid(); + + solve(); + + if (settings.output) + output_results(cycle); + + std::cout << std::endl; + } + } +} // namespace Step100 + + +int main(int argc, char *argv[]) { try { - // do nothing. + Step100::Settings settings; + if (!settings.try_parse((argc > 1) ? (argv[1]) : "")) + return 0; + + Step100::AdvectionProblem<2> advection_problem_2d(settings); + advection_problem_2d.run(); } catch (std::exception &exc) { @@ -41,7 +1116,6 @@ int main() << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - return 1; } catch (...) -- 2.39.5