MueLu factories for transfer operators

For this example, we reuse the Recirc2D example as introduced in Test example. The resulting linear systems are (slightly) non-symmetric and classical smoothed aggregation methods may be able to solve the problem but are not optimal in sense of convergence.

Multigrid setup phase - algorithmic design

Smoothed aggregation based algebraic multigrid methods originally have not been designed for non-symmetric linear systems. Inappropriately smoothed transfer operators may significantly deteriorate the convergence rate or even break convergence completely.

Unsmoothed transfer operators

Warning

Insert and link missing figures

Before we introduce smoothed aggregation methods for non-symmetric linear systems, we first go back one step and demonstrate how to use non-smoothed transfer operators which are eligible for non-symmetric linear systems. Figure XML Interface gives a simplified example how to build the coarse level matrix \(A_c\) using the fine level matrix \(A\) only. First, we “somehow” build aggregates using the information of the fine level matrix \(A\). The aggregates are then used to build the tentative non-smoothed prolongation operator. The restrictor is just the transpose of the (tentative) prolongator and finally the coarse level matrix \(A_c\) is calculated by the triple product \(A_c=RAP\).

In Figure XML Interface, the SaPFactory has been added after the TentativePFactory. Therein the non-smoothed transfer operator from the TentativePFactory is smoothed using information of the fine level matrix \(A\). This transfer operator design is used per default when the user does not specify its own transfer operator design. The default settings are optimal for symmetric positive definite systems. However, they might be problematic for our non-symmetric problem.

Smoothed transfer operators for non-symmetric systems

Warning

Insert missing figures

In case of non-symmetric linear systems it is \(A\neq A^T\). Therefore it is a bad idea just to use the transposed of the smoothed prolongation operator for the restrictor. Let \(\widehat{P}\) be the non-smoothed tentative prolongation operator. Then the smoothed prolongation operator \(P\) is built using

\(P = \bigl(I-\omega A\bigr) \widehat{P}\)

with some reasonable smoothing parameter \(\omega>0\). The standard restrictor is

\(R = P^T = \widehat{P}^T - \omega \widehat{P}^T A^T = \widehat{P}^T\bigl(I-\omega A^T\bigr).\)

That is, the restrictor would be smoothed using the information of \(A^T\). However, for non-symmetric systems we want to use the information of matrix \(A\) for smoothing the restriction operator, too. The restriction operator shall we built by the formula

\(R = P^T = \widehat{P}^T - \omega \widehat{P}^T A.\)

This corresponds to apply the same smoothing strategy to the non-smoothed restriction operator \(\widehat{R}=\widehat{P}^T\), which is applied to the (tentative) prolongation operator with using \(A^T\) as input instead of matrix \(A\). Figure muelu_factories_for_transfer_operators/figure_simpledesignpamg shows the changed factory design. The dashed line denotes, that the same smoothing strategy is used than for the prolongation operator. The concept is known as Petrov-Galerkin smoothed aggregation approach in the literature. A more advanced transfer operator smoothing strategy for non-symmetric linear systems that is based on the Petrov-Galerkin approach is described in [1]. Another approach based on SchurComplement approximations can be found in [2].

Insert missing figures here

XML Interface

Unsmoothed transfer operators

To construct a multigrid hierarchy with unsmoothed transfer operators one can use the following XML file (stored in ../../../test/tutorial/s3a.xml)

<ParameterList name="MueLu">

   <!-- Factory collection -->
  <ParameterList name="Factories">

    <ParameterList name="UncoupledAggregationFact">
      <Parameter name="factory"              type="string" value="UncoupledAggregationFactory"/>
      <Parameter name="aggregation: ordering"                    type="string" value="natural"/>
      <Parameter name="aggregation: max selected neighbors"            type="int"    value="0"/>
      <Parameter name="aggregation: min agg size"           type="int"    value="4"/>
    </ParameterList>

    <ParameterList name="myTentativePFact">
      <Parameter name="factory"                        type="string" value="TentativePFactory"/>
    </ParameterList>
    
  </ParameterList>

  <!-- Definition of the multigrid preconditioner -->
  <ParameterList name="Hierarchy">

    <Parameter name="use kokkos refactor"                   type="bool"     value="true"/>
    <Parameter name="max levels"                            type="int"      value="10"/>
    <Parameter name="coarse: max size"                      type="int"      value="10"/>
    <Parameter name="verbosity"                             type="string"   value="High"/>

    <ParameterList name="All">
      <Parameter name="Aggregates"           type="string"   value="UncoupledAggregationFact"/>
      <Parameter name="Nullspace"                    type="string"   value="myTentativePFact"/>
      <Parameter name="P"                            type="string"   value="myTentativePFact"/>
      <Parameter name="CoarseSolver"                     type="string"   value="DirectSolver"/>
    </ParameterList>

  </ParameterList>
</ParameterList>

Beside the TentiativePFactory which is responsible to generate the unsmoothed transfer operators we also introduce the UncoupledAggregationFactory with this example. In the Factories section of the XML file, you find both an entry for the aggregation factory and the prolongation operator factory with its parameters. In the Hierarchy section the defined factories are just put in into the multigrid setup algorithm. That is, the factory with the name UncoupledAggregationFact is used to generate the Aggregates and the myTentativePFact is responsible for generating both the (unsmoothed) prolongation operator P and the (coarse) near null space vectors Nullspace.

Note

For some more details about the (hidden) NullspaceFactory, which is internally used to handle the null space information and the dependencies, the reader might refer to Setup of block transfer operators.

Note

One can also use the Ptent variable for registering a TentativePFactory. This makes the TentativePFactory somewhat special in its central role for generating an aggregation based multigrid hierarchy. MueLu is smart enough to understand that you want to use the near null space vectors generated by the factory registered as Ptent for setting up the transfer operators.

That is, the following code would explicitly use the TentativePFactory object that is created as myTentativePFact. Since no factory is specified for the prolongation operator P, MueLu would decide to use a smoothed aggregation prolongation operator (represented by the SaPFactory), which correctly uses the factory for Ptent for the unsmoothed transfers with all its dependencies.

<ParameterList name="MueLu">
<ParameterList name="Factories">
        <ParameterList name="myTentativePFact">
        <Parameter name="factory" type="string" value="TentativePFactory"/>
        </ParameterList>
</ParameterList>

<ParameterList name="Hierarchy">
        <ParameterList name="Levels">
        <Parameter name="Ptent" type="string" value="myTentativePFact"/>
        </ParameterList>
</ParameterList>
</ParameterList>

Exercise 1

Create a sublist in the Factories part of the XML file for the restriction operator factory. Use a TransPFactory which builds the transposed of P to generate the restriction operator R. Register your restriction factory in the Hierarchy section to generate the variable R.

Smoothed aggregation for non-symmetric problems

Next, let’s try smoothed transfer operators for the non-symmetric linear system and compare the results of the transfer operator designs. Take a look at the XML file (in ../../../test/tutorial/s3b.xml).

<ParameterList name="MueLu">

  <!-- Factory collection -->
  <ParameterList name="Factories">

    <!-- Note that ParameterLists must be defined prior to being used -->

    <ParameterList name="UncoupledAggregationFact">
      <Parameter name="factory"                             type="string" value="UncoupledAggregationFactory"/>
      <Parameter name="aggregation: ordering"               type="string" value="natural"/>
      <Parameter name="aggregation: max selected neighbors" type="int"    value="0"/>
      <Parameter name="aggregation: min agg size"           type="int"    value="4"/>
    </ParameterList>

    <ParameterList name="myTentativePFact">
      <Parameter name="factory"                             type="string" value="TentativePFactory"/>
    </ParameterList>
    <ParameterList name="myProlongatorFact">
      <Parameter name="factory"                             type="string" value="SaPFactory"/>
      <Parameter name="P"                                   type="string" value="myTentativePFact"/>
      <Parameter name="sa: damping factor"                  type="double" value="1.0"/>
    </ParameterList>
    <ParameterList name="myTentRestrictorFact">
      <Parameter name="factory"                             type="string" value="TransPFactory"/>
      <Parameter name="P"                                   type="string" value="myTentativePFact"/>
    </ParameterList>
    <ParameterList name="mySymRestrictorFact">
      <Parameter name="factory"                             type="string" value="TransPFactory"/>
      <Parameter name="P"                                   type="string" value="myProlongatorFact"/>
    </ParameterList>
    <ParameterList name="myNonsymRestrictorFact">
      <Parameter name="factory"                             type="string" value="GenericRFactory"/>
      <Parameter name="P"                                   type="string" value="myProlongatorFact"/>
    </ParameterList>

    
    <ParameterList name="SymGaussSeidel">
      <Parameter name="factory"                             type="string" value="TrilinosSmoother"/>
      <Parameter name="type"                                type="string" value="RELAXATION"/>
      <ParameterList name="ParameterList">
        <Parameter name="relaxation: type"                  type="string" value="Symmetric Gauss-Seidel"/>
        <Parameter name="relaxation: sweeps"                type="int"    value="10"/>
        <Parameter name="relaxation: damping factor"        type="double" value="0.8"/>
      </ParameterList>
    </ParameterList>

  </ParameterList>

  <!-- Definition of the multigrid preconditioner -->
  <ParameterList name="Hierarchy">

    <Parameter name="use kokkos refactor"                   type="bool"     value="true"/>
    <Parameter name="max levels"                            type="int"      value="10"/>
    <Parameter name="coarse: max size"                      type="int"      value="10"/>
    <Parameter name="verbosity"                             type="string"   value="High"/>

    <ParameterList name="All">
      <Parameter name="Smoother"                    type="string"   value="SymGaussSeidel"/>
      <Parameter name="Aggregates"                  type="string"   value="UncoupledAggregationFact"/>
      <Parameter name="Nullspace"                   type="string"   value="myTentativePFact"/>
      <Parameter name="P"                           type="string"   value="myProlongatorFact"/>
      <Parameter name="R"                           type="string"   value="mySymRestrictorFact"/>
      <Parameter name="CoarseSolver"                type="string"   value="DirectSolver"/>
    </ParameterList>

  </ParameterList>
</ParameterList>

The interesting part is the Factories section where several different factories for the restriction operator are defined

Description

• myTentRestrictorFact: just uses the transposed of the unsmoothed prolongator for restriction.

• mySymRestrictorFact: uses the transposed of the smoothed prolongator for restriction.

• myNonsymRestrictorFact: uses the special non-symmetric smoothing for the restriction operator (based on the SaPFactory smoothing factory).

Note

The MueLu framework is very flexible and allows for arbitrary combinations of factories. However, be aware that the TentativePFactory cannot be used as input for the GenericRFactory. That is no problem since this combination not really makes sense. If you are using the TentativePFactory as your final prolongation operator you always have to use the TransPFactory for generating the restriction operators.

Exercise 2

Run the Recirc2D example with the different restriction operator strategies and compare the results for the iterative solver. What do you observe? What is the best choice for the transfer operators in the non-symmetric case?

Exercise 3

Change the myProlongatorFact from type SaPFactory to PgPFactory, which uses automatically calculated local damping factors instead of a global damping factor (with some user parameter sa: damping factor). Note that the PgPFactory might not accept the sa: damping factor parameter such that you have to comment it out (using <!– … –>).

Exercise 4

Try to set up a multigrid hierarchy with unsmoothed transfer operators for the transition from the finest level to level 1 and then use smoothed aggregation for the coarser levels (starting from level 1).

Footnotes

[1]

Sala, M. and Tuminaro, R. S., A new Petrov-Galerkin Smoothed Aggregation Preconditionerfor nonsymmetric Linear Systems, SIAM J. Sci. Comput., 2008, 31, p. 143–166

[2]

Wiesner, T. A., Tuminaro, R. S., Wall, W. A. and Gee, M. W., Multigrid transfers fornonsymmetric systems based on Schur complements and Galerkin projections., Numer.Linear Algebra Appl., 2013, doi: 10.1002/nla.1889