From: tcclevenger Date: Fri, 10 May 2019 19:19:57 +0000 (-0600) Subject: Changelog, prm fixes, initial docs, add new parameters X-Git-Tag: v9.1.0-rc2~14^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1d8e80cb74c0350a421459e60b4525d060e14c9f;p=dealii.git Changelog, prm fixes, initial docs, add new parameters --- diff --git a/doc/news/changes/major/20190510Clevenger b/doc/news/changes/major/20190510Clevenger new file mode 100644 index 0000000000..a4f820fd6e --- /dev/null +++ b/doc/news/changes/major/20190510Clevenger @@ -0,0 +1,3 @@ +New: The Step-63 tutorial has been added which discusses block smoothers vs. point smoothers inside a geometric multigrid v-cycle. This tutorial also implements a GMG solver for a non-symmetric PDE, and demonstrates the effect of DoF/cell renumbering for multiplicative smoothers with advection dominated problems. +
+(Conrad Clevenger, 2019/05/10) diff --git a/examples/step-63/block_jacobi.prm b/examples/step-63/block_jacobi.prm index 9e6f98326e..941a0df9ed 100644 --- a/examples/step-63/block_jacobi.prm +++ b/examples/step-63/block_jacobi.prm @@ -1,9 +1,15 @@ +#Epsilon +set Epsilon = 0.005 + # Finite Element degree set Fe degree = 1 # Smoother Type: SOR|Jacobi|block SOR|block Jacobi set Smoother type = block Jacobi +# Number of smoothing steps +set Smoothing steps = 3 + # DoF renumbering: no|downstream|upstream|random set DoF renumbering = downstream @@ -12,3 +18,4 @@ set With streamline diffusion = true # Output: true|false set Output = false + diff --git a/examples/step-63/block_sor.prm b/examples/step-63/block_sor.prm index 375e7e63cb..3b20e403a2 100644 --- a/examples/step-63/block_sor.prm +++ b/examples/step-63/block_sor.prm @@ -1,9 +1,15 @@ +#Epsilon +set Epsilon = 0.005 + # Finite Element degree set Fe degree = 1 # Smoother Type: SOR|Jacobi|block SOR|block Jacobi set Smoother type = block SOR +# Number of smoothing steps +set Smoothing steps = 1 + # DoF renumbering: no|downstream|upstream|random set DoF renumbering = downstream @@ -11,4 +17,4 @@ set DoF renumbering = downstream set With streamline diffusion = true # Output: true|false -set Output = false \ No newline at end of file +set Output = false diff --git a/examples/step-63/doc/kind b/examples/step-63/doc/kind index 6816e9090f..387b9f817c 100644 --- a/examples/step-63/doc/kind +++ b/examples/step-63/doc/kind @@ -1 +1,2 @@ unfinished + diff --git a/examples/step-63/doc/results.dox b/examples/step-63/doc/results.dox index b5eaba9377..d2fe3c9640 100644 --- a/examples/step-63/doc/results.dox +++ b/examples/step-63/doc/results.dox @@ -1,2 +1,472 @@

Results

+

GMRES Iteration Counts

+ +The major advantage for GMG is that it is an $\mathcal{O}(n)$ method, that is, the complexity of the problem increases linearly with the problem size. To show then that the linear solver presented in this tutorial is also $\mathcal{O}(n)$, all one needs to do is show that the iteration counts for the GMRES solve stay roughly constant as we refine the mesh. + +Each of the following tables gives the GMRES iteration counts to reduce the initial residual by 1e-8. + +
DoF/Cell Renumbering:
+ +Staring with the additive smoothers, we see that renumbering the DoFs/cells has no effect on convergence speed. This is because these smoothers compute operations on each DoF (point-smoother) or cell (block-smoother) independently and add up the results. Since we can define these smoothers as an application of a sum of matrices, and matrix addition is communicative, the order at which we sum the different components will not affect the end result. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
$Q_1$Smoother (smoothing steps)
Jacobi (6)Block Jacobi (3)
Renumbering StrategyRenumbering Strategy
CellsDoFsDownstreamRandomUpstreamDownstreamRandomUpstream
32483 + 3 + 3 + 3 + 3 + 3 +
1281606 + 6 + 6 + 6 + 6 + 6 +
51257611 + 11 + 11 + 9 + 9 + 9 +
2048217615 + 15 + 15 + 13 + 13 + 13 +
8192844818 + 18 + 18 + 15 + 15 + 15 +
327683328020 + 20 + 20 + 16 + 16 + 16 +
13107213209620 + 20 + 20 + 16 + 16 + 16 +
+ +On the other hand, for multiplicative smoothers, we can speed up convergence by renumbering the DoFs/cells in the advection direction, and similarly, we can slow down convergence if we do the renumbering in the opposite direction. This is because the multiplicative smoothers transfer information from DoF to DoF, and the solutions at any one DoF of an advection-dominate problem is highly dependent on the solution at DoFs upstream of the advection direction. + +As we increase $\varepsilon$, however, the problem becomes more diffusion-dominated and the effectiveness of DoF/cell renumbering drops. On the other hand, in the extreme case of $\varepsilon=0$ (advection-only), we have a 1st-order PDE and multiplicative methods become effective solvers (Note: special care must be taken for the boundary conditions in this case). + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
$Q_1$Smoother (smoothing steps)
SOR (3)Block SOR (1)
Renumbering StrategyRenumbering Strategy
CellsDoFsDownstreamRandomUpstreamDownstreamRandomUpstream
32482 + 2 + 3 + 2 + 2 + 3 +
1281605 + 5 + 7 + 5 + 5 + 7 +
5125767 + 9 + 11 + 7 + 7 + 12 +
2048217610 + 12 + 15 + 8 + 10 + 17 +
8192844811 + 15 + 19 + 10 + 11 + 20 +
327683328012 + 16 + 20 + 10 + 12 + 21 +
13107213209612 + 16 + 19 + 11 + 12 + 21 +
+ + +
Point vs. Block Smoothers
+ +We will limit the results to runs using the downstream renumbering. Here is a cross comparison of all four smoothers for both $Q_1$ and $Q_3$ elements: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ $Q_1$Smoother (smoothing steps)$Q_3$Smoother (smoothing steps)
Cells + DoFsJacobi (6)Block Jacobi (3)SOR (3)Block SOR (1)DoFsJacobi (6)Block Jacobi (3)SOR (3)Block SOR (1)
32 + 483 + 3 + 2 + 2 + + 33615 + 14 + 15 + 6 +
128 + 1606 + 6 + 5 + 5 + + 124823 + 18 + 21 + 9 +
512 + 57611 + 9 + 7 + 7 + + 480029 + 21 + 28 + 9 +
2048 + 217615 + 13 + 10 + 8 + + 1881633 + 22 + 32 + 9 +
8192 + 844818 + 15 + 11 + 10 + + 7449635 + 22 + 34 + 10 +
32768 + 3328020 + 16 + 12 + 10 + +
131072 + 13209620 + 16 + 12 + 11 + +
+ +We see that for $Q_1$, both multiplicative smoothers require a smaller combination of smoothing steps and iteration counts than either additive smoother. However, when we in increase the degree to a $Q_3$ element, there is a clear advantage for the block smoothers in terms of the increase in smoothing steps and iterations required to solve. Specifically, the block SOR smoother gives constant iterations over degree, and the block Jacobi only sees about a 38% increase in iterations compared to 75% and 183% for Jacobi and SOR respectively. + +

Cost

+ +Iteration counts do not tell the full story in the optimality of a one smoother over another. Obviously we must examine the cost of an iteration. Block smoothers here are at a disadvantage as they are having to construct and invert a cell matrix for each cell. Here is a comparison of solve times for a $Q_3$ element with 74,496 DoFs: + + + + + + + + + + + + + + + + + + + +
$Q_3$Smoother (smoothing steps)
DoFsJacobi (6)Block Jacobi (3)SOR (3)Block SOR (1)
744960.689 s + 5.826 s + 1.189 s + 1.021 s +
+ +The smoother that requires the most iterations (Jacobi) actually take the shortest time (roughly 2/3 the time of the next fastest method). This is because all that is requires to apply a Jacobi smoothing step is multiplication by a diagonal matrix which is very cheap. On the other hand, while SOR requires over 3x more iterations (each with 3x more smoothing steps) than block SOR, the times are roughly equivalent, implying that a smoothing step of block SOR is roughly 9x slower than a smoothing step of SOR. Lastly, block Jacobi is almost 6x more expensive than block SOR, which intuitively makes sense from the fact that 1 step of each method has the same cost (inverting the cell matrices and either adding or multiply them together), and block Jacobi has 3 times the smoothing step per iteration with 2 times the iterations. + + +

Additional Points

+ +There are a few more important points to mention: + +1. For a mesh distributed in parallel, multiplicative methods cannot be executed over the entire domain. One can use a hybrid method where a multiplicative smoother is applied on each subdomain, but as you increase the number of subdomains, the method converges to that of an additive method. This is a major disadvantage to these methods. + +2. Current research into block smoothers suggest that soon we will be able to compute the inverse of the cell matrices much cheaper than what is currently being done inside deal.II. It seems that one should be able to take advantage of matrix free implementations and also the fact that, in the interior of the domain, cell matrices tend to look very similar. + +Combining 1. and 2. gives a good reason for expecting that a method like block Jacobi could become very powerful in the future, even though currently for these examples it is quite slow. + + + +

Possible Extensions

+ +
Constant iterations for $Q_5$
+ +Change the number of smoothing steps and the smoother relaxation parameter (set in Smoother::AdditionalData() inside create_smoother(), only necessary for point smoothers) so that we maintain a constant number of iterations for a $Q_5$ element. + +
Effectiveness of renumbering for changing $\varepsilon$
+ +Increase/decrease the parameter ``Epsilon" in the .prm files of the multiplicative methods and observe the relationship of downstream vs. upstream reordering. + + + + diff --git a/examples/step-63/jacobi.prm b/examples/step-63/jacobi.prm index 4c7365e06a..75314cd087 100644 --- a/examples/step-63/jacobi.prm +++ b/examples/step-63/jacobi.prm @@ -1,14 +1,20 @@ +#Epsilon +set Epsilon = 0.005 + # Finite Element degree set Fe degree = 1 # Smoother Type: SOR|Jacobi|block SOR|block Jacobi set Smoother type = Jacobi +# Number of smoothing steps +set Smoothing steps = 6 + # DoF renumbering: no|downstream|upstream|random -set DoF renumbering = downstream +set DoF renumbering = downstream # With streamline diffusion: true|false set With streamline diffusion = true # Output: true|false -set Output = false \ No newline at end of file +set Output = false diff --git a/examples/step-63/sor.prm b/examples/step-63/sor.prm index dbdc74c198..51d83bf3bb 100644 --- a/examples/step-63/sor.prm +++ b/examples/step-63/sor.prm @@ -1,14 +1,20 @@ +#Epsilon +set Epsilon = 0.005 + # Finite Element degree -set Fe degree = 1 +set Fe degree = 1 # Smoother Type: SOR|Jacobi|block SOR|block Jacobi set Smoother type = SOR +# Number of smoothing steps +set Smoothing steps = 3 + # DoF renumbering: no|downstream|upstream|random -set DoF renumbering = downstream +set DoF renumbering = downstream # With streamline diffusion: true|false set With streamline diffusion = true # Output: true|false -set Output = false \ No newline at end of file +set Output = false diff --git a/examples/step-63/step-63.cc b/examples/step-63/step-63.cc index ba6230ca31..65d2fae72e 100644 --- a/examples/step-63/step-63.cc +++ b/examples/step-63/step-63.cc @@ -17,24 +17,22 @@ * Timo Heister, University of Utah */ -// @note: This is work in progress and will be an example for block smoothers -// in geometric multigrid. +// @sect3{Include files} + +// Typical files needed for standard deal.II: #include -#include -#include #include #include -#include #include #include -#include + +#include #include #include #include #include #include #include -#include #include #include @@ -43,38 +41,51 @@ #include #include #include +#include + #include #include #include #include + +#include +#include #include + #include #include #include -#include -#include -#include -#include -#include +// Include all relevant multilevel files: #include #include #include -#include #include #include #include #include +// C++: #include #include #include +// We will be using Meshworker::mesh_loop functionality for assembling matrices: +#include + namespace Step63 { using namespace dealii; + // @sect3{MeshWorker Data} + + // The following are structure needed by assemble_cell() + // function used by Meshworker::mesh_loop(). ScratchData + // contains a FeValues object with is needed for assembling + // a cells local contribution, while CopyData contains the + // output from a cells local contribution and necessary information + // to copy that to the global system. template struct ScratchData @@ -112,6 +123,13 @@ namespace Step63 }; + + // @sect3{Problem parameters} + + // We will use ParameterHandler to pass in parameters at runtime. The + // structure Settings parses and stores these parameters to be queried + // throughout the program. + struct Settings { enum DoFRenumberingStrategy @@ -124,8 +142,10 @@ namespace Step63 void get_parameters(const std::string &prm_filename); + double epsilon; unsigned int fe_degree; std::string smoother_type; + unsigned int smoothing_steps; DoFRenumberingStrategy dof_renumbering; bool with_streamline_diffusion; bool output; @@ -135,14 +155,20 @@ namespace Step63 { ParameterHandler prm; + prm.declare_entry("Epsilon", "0.005", Patterns::Double(0), "Epsilon"); + 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"), + Patterns::Selection("SOR|Jacobi|block SOR|block Jacobi"), "Smoother Type: SOR|Jacobi|block SOR|block Jacobi"); + prm.declare_entry("Smoothing steps", + "2", + Patterns::Integer(1), + "Number of smoothing steps"); prm.declare_entry("DoF renumbering", "downstream", Patterns::Selection("none|downstream|upstream|random"), @@ -165,8 +191,10 @@ namespace Step63 prm.parse_input(prm_filename); - fe_degree = prm.get_integer("Fe degree"); - smoother_type = prm.get("Smoother type"); + epsilon = prm.get_double("Epsilon"); + fe_degree = prm.get_integer("Fe degree"); + smoother_type = prm.get("Smoother type"); + smoothing_steps = prm.get_integer("Smoothing steps"); std::string renumbering = prm.get("DoF renumbering"); if (renumbering == "none") @@ -183,8 +211,18 @@ namespace Step63 } - // Functions for creating permutation of cells for output and Block - // smoothers + // @sect1{Cell permutations} + // + // The ordering in which cells and degrees of freedom are traversed + // will play a roll in the speed of convergence for multiplicative + // methods. Here we define functions which return a specific ordering + // of cells to be used by the block smoothers. + + // For each type of cell ordering, we define a function for the active + // mesh and one for a level mesh. While the only reordering necessary + // for solving the system will be on the level meshes, we include the + // active reordering for visualization purposes in output_results(). + template std::vector create_downstream_cell_ordering(const DoFHandler &dof_handler, @@ -236,8 +274,6 @@ namespace Step63 return ordered_indices; } - - template std::vector create_random_cell_ordering(const DoFHandler &dof_handler, @@ -251,15 +287,15 @@ namespace Step63 for (const auto &cell : dof_handler.cell_iterators_on_level(level)) ordered_cells.push_back(cell->index()); - // shuffle the elements + // Shuffle the elements: std::mt19937 random_number_generator; for (unsigned int i = 1; i < n_cells; ++i) { - // get a random number between 0 and i (inclusive) + // Get a random number between 0 and i (inclusive): const unsigned int j = std::uniform_int_distribution<>(0, i)(random_number_generator); - // if possible, swap the elements + // If possible, swap the elements: if (i != j) std::swap(ordered_cells[i], ordered_cells[j]); } @@ -280,15 +316,15 @@ namespace Step63 for (const auto &cell : dof_handler.active_cell_iterators()) ordered_cells.push_back(cell->index()); - // shuffle the elements + // shuffle the elements: std::mt19937 random_number_generator; for (unsigned int i = 1; i < n_cells; ++i) { - // get a random number between 0 and i (inclusive) + // get a random number between 0 and i (inclusive): const unsigned int j = std::uniform_int_distribution<>(0, i)(random_number_generator); - // if possible, swap the elements + // if possible, swap the elements: if (i != j) std::swap(ordered_cells[i], ordered_cells[j]); } @@ -297,7 +333,13 @@ namespace Step63 } + // @sect3{Right-hand Side and Boundary Values} + + // The problem solved in this tutorial is an adaptation of Ex. ___ + // found in ______ (how to cite?), namely, we add a hole in the middle + // of our domain. + // We have a zero right-hand side. template class RightHandSide : public Function { @@ -339,7 +381,9 @@ namespace Step63 } - + // We have Dirichlet boundary conditions. On a connected portion of the + // outer, square boundary we set the value to 1, and we set the value to 0 + // everywhere else (including the inner, circular boundary). template class BoundaryValues : public Function { @@ -364,9 +408,9 @@ namespace Step63 Assert(component == 0, ExcIndexRange(component, 0, 1)); (void)component; - 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 - ) + // Set boundary to 1 if $x=1$, or if $x>0.5$ and $y=-1$. + if (std::fabs(p[0] - 1) < 1e-8 || + (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5)) { return 1.0; } @@ -390,14 +434,19 @@ namespace Step63 values[i] = BoundaryValues::value(points[i], component); } + + + // @sect3{Streamline Diffusion} + + // Streamline diffusion stabilization term. Value defined in + // 'On discontinuity–capturing methods for convection–diffusion + // equations' (cite?) template double compute_stabilization_delta(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' const double Peclet = dir.norm() * hk / (2.0 * eps * pk); const double coth = (1.0 + std::exp(-2.0 * Peclet)) / (1.0 - std::exp(-2.0 * Peclet)); @@ -406,6 +455,17 @@ namespace Step63 } + // @sect3{AdvectionProlem class} + + // This main class of the program, and should look very similar to step-16. + // The major difference is that, since we are defining our multigrid smoother + // at runtime, we choose to define a function create_smoother() and a class + // object mg_smoother which is a std::unique_ptr to a smoother that is derived + // from MGSmoother. Note that for smoother derived from RelaxationBlock, we + // must include a smoother_data object for each level. This will contain + // information about the cell ordering and the method of inverting cell + // matrices. + template class AdvectionProblem { @@ -453,16 +513,15 @@ namespace Step63 mg::Matrix> mg_interface_matrix_in; mg::Matrix> mg_interface_matrix_out; + std::unique_ptr>> mg_smoother; + using SmootherType = RelaxationBlock, double, Vector>; using SmootherAdditionalDataType = SmootherType::AdditionalData; - std::unique_ptr>> mg_smoother; - MGLevelObject smoother_data; MGConstrainedDoFs mg_constrained_dofs; - const double epsilon; Tensor<1, dim> advection_direction; const Settings settings; @@ -476,29 +535,33 @@ namespace Step63 , dof_handler(triangulation) , fe(settings.fe_degree) , mapping(settings.fe_degree) - , epsilon(0.005) , settings(settings) { 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); + AssertThrow(false, ExcNotImplemented()); } + // @sect4{AdvectionProblem::setup_system} + + // Here we setup the DoFHandler, ConstraintMatrix, and sparsity patterns for + // both active and multigrid level meshes. template void AdvectionProblem::setup_system() { const unsigned int n_levels = triangulation.n_levels(); + // Setup active DoFs: dof_handler.distribute_dofs(fe); - // We could renumber the active dofs with DoFRenumbering::downstream() + // We could renumber the active DoFs with the DoFRenumbering class // 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. + // wouldn't matter for the computations. Instead, we will renumber the + // DoFs on each multigrid level below. solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); @@ -523,7 +586,7 @@ namespace Step63 system_matrix.reinit(sparsity_pattern); - // Setup GMG DoFs + // Setup GMG DoFs: dof_handler.distribute_mg_dofs(); // Renumber DoFs on each level in downstream or upstream direction if @@ -591,8 +654,8 @@ namespace Step63 level); mg_interface_sparsity_patterns[level].copy_from(dsp); - // We need both interface in and out matrices since our problem is not - // symmetric + // Unlike the other GMG tutorials, we need both interface in and out + // matrices since our problem is non-symmetric. mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]); mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]); } @@ -600,6 +663,13 @@ namespace Step63 } + // @sect4{AdvectionProblem::assemble_cell} + + // Here we define the assembly of the linear system on each cell to be used by + // the mesh_loop() function below. This one function assembles the cell matrix + // for both and active and a level cell, and only assembles a right-hand side + // if called for an active cell. + template template void AdvectionProblem::assemble_cell(const IteratorType &cell, @@ -630,9 +700,14 @@ namespace Step63 right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(), rhs_values); + // If we are using streamline diffusion we must add its contribution + // to both the cell matrix and the cell right-handside. If we are not + // using streamline diffusion, setting $\delta=0$ negates this contribution + // below and we are left with the standard, Galerkin finite element + // assembly. const double delta = settings.with_streamline_diffusion ? compute_stabilization_delta(cell->diameter(), - epsilon, + settings.epsilon, advection_direction, settings.fe_degree) : 0.0; @@ -643,36 +718,36 @@ namespace Step63 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) * + // Galerkin contribution: + (settings.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) + + // Streamline diffusion contribution: + 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 * settings.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 (settings.with_streamline_diffusion) - 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() == false) { copy_data.cell_rhs(i) += + // Galerkin contribution: scratch_data.fe_values.shape_value(i, q_point) * - rhs_values[q_point] * scratch_data.fe_values.JxW(q_point); - if (settings.with_streamline_diffusion) - copy_data.cell_rhs(i) += - delta * rhs_values[q_point] * advection_direction * + rhs_values[q_point] * scratch_data.fe_values.JxW(q_point) + + // Streamline diffusion contribution: + delta * rhs_values[q_point] * advection_direction * scratch_data.fe_values.shape_grad(i, q_point) * scratch_data.fe_values.JxW(q_point); } @@ -680,6 +755,11 @@ namespace Step63 } + // @sect4{AdvectionProblem::assemble_system_and_multigrid} + + // Here we employ Meshworker::mesh_loop() to go over cells and assemble the + // system_matrix, system_rhs, and all mg_matrices for us. + template void AdvectionProblem::assemble_system_and_multigrid() { @@ -709,6 +789,9 @@ namespace Step63 CopyData(), MeshWorker::assemble_own_cells); + // Unlike the constraints for the active level, we choose to create + // local constraint matrices for each multigrid level since they are + // never needed elsewhere in the program. std::vector> boundary_constraints( triangulation.n_global_levels()); for (unsigned int level = 0; level < triangulation.n_global_levels(); @@ -740,7 +823,10 @@ namespace Step63 // 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). + // interface_in(i,j). This is an optimization for the symmetric case + // which allows only one matrix to be used when setting the edge_matrices + // in solve(). Here, however, since our problem is non-symmetric, we must + // store both interface_in and interface_out matrices. 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( @@ -769,6 +855,35 @@ namespace Step63 } + // @sect4{AdvectionProblem::setup_smoother} + + // Here we setup the smoother based on the settings in the .prm. The two + // options that are of significance is the number of pre- and post-smoothing + // steps on each level of the multigrid v-cycle and the relaxation parameter. + + // Since multiplicative methods tend to be more powerful than additive method, + // fewer smoothing steps are required to see convergence indepedent of mesh + // size. The same hold for block smoothers over point smoothers. This is + // reflected in the choice for the number of smoothing steps for each type of + // smoother below. + + // The relaxation parameter for point smoothers is chosen based on trial and + // error, and they reflect values necessary to keep the iteration counts in + // the GMRES solve constant (or as close as possible) as we refine the mesh. + // The two values given for both ``Jacobi" and ``SOR" are for degree 1 and + // degree 3 finite elements. If the user wants to change to another degree, + // they may need to adjust these numbers. For block smoothers, this parameter + // has a more straightforward interpretation, namely that for additive methods + // in 2D, a DoF can have a repeated contribution from up to 4 cells, + // therefore we must relax these methods by 0.25 to compensate. This is not an + // issue for multiplicative methods as each cell inverse application carries + // new information to all its DoFs. + + // Finally, as mention above, the point smoothers only operate on DoFs, and + // the block smoothers on cells, so only the block smoothers need to be given + // information regarding cell orderings. DoF ordering for point smoothers has + // already been taken care of in setup_system(). + template void AdvectionProblem::setup_smoother() { @@ -783,7 +898,7 @@ namespace Step63 smoother->initialize(mg_matrices, Smoother::AdditionalData(fe.degree == 1 ? 1.0 : 0.62)); - smoother->set_steps(2); + smoother->set_steps(settings.smoothing_steps); mg_smoother = std::move(smoother); } else if (settings.smoother_type == "Jacobi") @@ -796,7 +911,7 @@ namespace Step63 smoother->initialize(mg_matrices, Smoother::AdditionalData(fe.degree == 1 ? 0.6667 : 0.47)); - smoother->set_steps(4); + smoother->set_steps(settings.smoothing_steps); mg_smoother = std::move(smoother); } else if (settings.smoother_type == "block SOR" || @@ -817,6 +932,8 @@ namespace Step63 std::vector ordered_indices; switch (settings.dof_renumbering) { + // Order the cells downstream with respect + // to the advection direction. case Settings::DoFRenumberingStrategy::downstream: ordered_indices = create_downstream_cell_ordering(dof_handler, @@ -824,6 +941,9 @@ namespace Step63 level); break; + // Order the cells upstream with respect to the advection + // direction, i.e., downstream with respect to the negative + // of the advection direction. case Settings::DoFRenumberingStrategy::upstream: ordered_indices = create_downstream_cell_ordering(dof_handler, @@ -831,12 +951,14 @@ namespace Step63 level); break; + // Order the cells randomly. case Settings::DoFRenumberingStrategy::random: ordered_indices = create_random_cell_ordering(dof_handler, level); break; - case Settings::DoFRenumberingStrategy::none: // Do nothing + // Keep the default cell ordering (z-order, see Glossary). + case Settings::DoFRenumberingStrategy::none: break; default: @@ -855,7 +977,7 @@ namespace Step63 RelaxationBlockSOR, double, Vector>, Vector>>(); smoother->initialize(mg_matrices, smoother_data); - smoother->set_steps(1); + smoother->set_steps(settings.smoothing_steps); mg_smoother = std::move(smoother); } else if (settings.smoother_type == "block Jacobi") @@ -867,7 +989,7 @@ namespace Step63 Vector>, Vector>>(); smoother->initialize(mg_matrices, smoother_data); - smoother->set_steps(2); + smoother->set_steps(settings.smoothing_steps); mg_smoother = std::move(smoother); } } @@ -876,6 +998,29 @@ namespace Step63 } + // @sect4{AdvectionProblem::solve} + + // Before we can solve the system, we must first set up the multigrid + // preconditioner. This is requires the setup of the transfer between levels, + // the coarse matrix solver, and the smoother. This setup follows almost + // identically to Step-16, the main difference being the various smoothers + // defined above and the fact that we need different interface edge matrices + // for in and out since our problem is non-symetric. (In reality, for this + // tutorial these interface matrices are empty since we are only using global + // refinement, and thus have no refinement edges. However, we have still + // included both here since if one made the simple switch to an adaptively + // refined method, the program would still run correctly.) + + // The last thing to note is that since our problem is non-symetric, we must + // use an appropriate Krylov subspace method. We choose here to + // use GMRES since it offers the guarentee of residual reduction in each + // iteration. The major disatvantage to GMRES is that, for each iteration, we + // must store an additional temporary vector as well as compute an additional + // scalar product. However, the goal of this tutorial is to have very low + // iteration counts by using a powerful GMG preconditioner, so this should not + // be a factor. If the user is interested, another sutaible method offered in + // deal.II would be BiCGStab. + template void AdvectionProblem::solve() { @@ -909,7 +1054,8 @@ namespace Step63 std::cout << " Solving with GMRES to tol " << solve_tolerance << "..." << std::endl; - SolverGMRES<> solver(solver_control); + SolverGMRES<> solver(solver_control, + SolverGMRES<>::AdditionalData(50, true)); Timer time; time.start(); @@ -926,11 +1072,14 @@ namespace Step63 } + // @sect4{AdvectionProblem::output_results} + + // Here we output the solution and cell ordering in a .vtu format. template void AdvectionProblem::output_results(const unsigned int cycle) const { - // Here we generate an index for each cell to visualize the ordering used + // 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 @@ -988,6 +1137,12 @@ namespace Step63 } + // @sect4{AdvectionProblem::run} + + // As in most tutorials, this function creates/refines the mesh and calls + // the various functions defined above to setup, assemble, solve, and output + // the results. + template void AdvectionProblem::run() { @@ -998,8 +1153,13 @@ namespace Step63 if (cycle == 0) { - GridGenerator::hyper_cube_with_cylindrical_hole( - triangulation, 0.3, 1.0, 0.5, 1, false); + // We are solving on the square [-1,1]^dim with a hole + // of radius 3/10 units centered at the origin. + GridGenerator::hyper_cube_with_cylindrical_hole(triangulation, + 0.3, + 1.0); + + // Set manifold for the inner (curved) boundary. static const SphericalManifold manifold_description( Point(0, 0)); triangulation.set_manifold(1, manifold_description); @@ -1028,6 +1188,14 @@ namespace Step63 } // namespace Step63 +// @sect4{The main function} + +// Here the main function is like most tutorials. The only interesting bit +// is that we require the user to pass a .prm file as a sole command line +// argument (see Step-19 for a complete discussion of parameter files). If no +// parameter file is given, the program will output the contents of a sample +// parameter file with all default values to the screen that the user can then +// copy and paste into their own .prm file. int main(int argc, char *argv[]) {