Muelu: Incorrect results with pressure oscillations for Stokes flow equations

[dgmaths9]

Hi @jhux2,

I am trying to solve the Stokes flow equations accounting change in density due to pressure. The application is the magma flow in the underground reservoir with initial pressure around 30 MPa. For discretization I use stabilized finite element method with linear elements (P1 or Q1). In my code, the DOFs are arranged as p, v1, v2 for each node i.e.

p, v1, v2, p,v1, v2, ...

The global matrix is unsymmetric.

I first tried with GMRES from Belos and RILUK preconditioner from IFPACK2 that produced right results although taking too many iterations and thus more computational time. As an improvement I wanted to use MUELU for faster computations, however, with MUELU the results are not correct. I see that pressure largely oscilates (see attached figure).

So far I have tried these options

1. multigrid algorithm: unsmoothed and smoother: ILU
2. multigrid algorithm: pg and smoother: ILU

In smoother I also tried the level of fill 0 and 1. WIth both options the results (with minor difference) were incorrect.

I also tried Relaxation: Jacobi with sweeps 1 and damping 0.9, but the simulations crashed immediately.

Here is the MueLu precondition configuration I used for the unsmoothed case:

<ParameterList name="MueLu">

  <!-- ===========  GENERAL ================ -->
  <Parameter name="verbosity" type="string" value="none"/>
  <Parameter name="max levels" type="int" value="10"/>
  <Parameter name="coarse: max size" type="int" value="2000"/>
  <Parameter name="multigrid algorithm" type="string" value="unsmoothed"/>
  <Parameter name="problem: symmetric" type="bool" value="false"/>
  <Parameter name="number of equations" type="int" value="3"/>

  <!-- ===========  SMOOTHING  =========== -->
  <Parameter name="smoother: type" type="string" value="ILU"/>

  <!-- ===========  AGGREGATION  =========== -->
  <Parameter name="aggregation: type" type="string" value="uncoupled"/>
  <Parameter name="aggregation: min agg size" type="int" value="1"/>
  <Parameter name="aggregation: max agg size" type="int" value="3"/>

 <!-- ===========  REPARTITIONING  =========== -->      
    <Parameter        name="repartition: enable"                  type="bool"     value="true"/>
    <Parameter        name="repartition: partitioner"             type="string"   value="zoltan"/>
    <Parameter        name="repartition: start level"             type="int"      value="2"/>
    <Parameter        name="repartition: min rows per proc"       type="int"      value="800"/>
    <Parameter        name="repartition: max imbalance"           type="double"   value="1.1"/>
    <Parameter        name="repartition: remap parts"             type="bool"     value="false"/>

I would appreciate If you can suggest what to try to get the correct results.

Many thanks,
Deepak

[github-actions[bot]]

Automatic mention of the @trilinos/muelu team

[cgcgcg]

When you say that the result is not correct: does MueLu fail to converge? I'm assuming you're using MueLu as a preconditioner in Belos. Does the solve achieve the desired tolerance?

EDIT1: Also, what do you mean by "the simulations crashed immediately"? MueLu should not crash.

EDIT2: Oh, this is Stokes. So you're effectively mushing together all degrees of freedom without regard for pressure or velocities. If the (2,2) block has a zero diagonal then a crash is not that surprising.

[dgmaths9]

Hi Christian,

Yes, I'm using MueLu as a preconditioner within Belos. The GMRES solver converges well with MueLu — the relative residuals reach around 1.e-9. However, when I visualize the results, I notice a strange behavior in the pressure field (as shown in the figure). In contrast, using the ILU preconditioner from IFPACK produces a much smoother pressure profile, without any noticeable overshoots or undershoots. So I suspect the issue might be related to the MueLu preconditioner.

1. By "simulations crashed," I meant that the pressure overshoots become extremely large and unphysical — to the point where the simulation logic breaks down.

2. Exactly — I was hoping to retain my current DOF ordering, as changing it would require a major restructuring of the code.

Just to add, I'm using SUPG stabilization, so the matrix diagonal is not zero.

[cgcgcg]

A priori this doesn't sound like a problem with the solver. If the linear system is solved to tolerance the error might be elsewhere. However, I would suggest to double check the residual in case Belos only prints out an implicit residual. It might also be worth checking whether you're using left or right preconditioning.

[dgmaths9]

Dear Christian,

Thank you very much for your support.

I am using a right preconditioner with Flexible GMRES, and I compute the residual error explicitly as follows:

        typedef Belos::MultiVecTraits<ST, MV>     MVT;
        typedef Belos::OperatorTraits<ST, MV, OP>  OPT;

        std::vector<double> rhs_norm(1);        
        MVT::MvNorm( *b, rhs_norm );
        Epetra_MultiVector resid(A->RowMap(), 1);
        OPT::Apply( *A, *x, resid );
        MVT::MvAddMv( -1.0, resid, 1.0, *b, resid );
        std::vector<double> res_error(1);        
        MVT::MvNorm( resid, res_error);
        abs_res_err_ = res_error[0];
        rel_res_err_ = res_error[0]/rhs_norm[0];     

In my simulation, the pressure is on the order of $10^7$, while the velocity is on the order of $10{-4}$. Could this large disparity in magnitudes be the reason for the inaccurate solution?

I have used the same MueLu preconditioner on problems where both pressure and velocity are small and comparable in scale, and in those cases, the solution is accurate.

[cgcgcg]

So, to recap:

• You computed the residual norm outside of Belos, using the monolithic system. The norm has similar magnitude no matter what solver option you are using, yet you see significant error in your simulations when using MueLu.
• The conjecture is that the sub-blocks of the system are scaled quite differently, and hence we do not actually converge to small enough tolerance in one of the solution components (velocity).

Before rescaling your system, you could reduce the solver tolerance and check if that indeed restores correct behavior. Otherwise we'd need to look at MueLu parameters to see if any of the dropping thresholds need to be adjusted.

[dgmaths9]

• I was using the solver tolerance = $10^{-9}$. Do you suggest to reduce it further?

• As my dofs are arranged as p, v1, v2 per node, my linear system does not have a clear block structure for pressure and velocity.

[cgcgcg]

Yes, I would reduce that further for testing the conjecture. If that is indeed successful you can look into rescaling the velocity or using a norm that is more appropriate.

[dgmaths9]

Hi Christian, following your suggestion, we tried to reduce the solver tolerance to $10^{-14}$, but even that did not produce good results (the orange line in the attached fig).

Do you have suggestion on what can be tried next?

[cgcgcg]

Have you tried other preconditioners? I'm wondering whether it's something specific in MueLu that fails due to the scaling. You could also try to equilibrate the system, see e.g.

Trilinos/packages/ifpack2/example/RelaxationWithEquilibration.cpp

Line 1762 in af33497

void equilibrateMatrix ()

Trilinos/packages/ifpack2/example/RelaxationWithEquilibration.cpp

Line 1801 in af33497

preScaleRightHandSides (Tpetra::MultiVector<scalar_type, local_ordinal_type,

[dgmaths9]

The red line corresponds to GMRES with an ILU preconditioner from IFPACK, which produced correct results. However, the drawback was that the linear solver required a large number of iterations, resulting in increased computation time. To address this, we attempted to accelerate the solution process using MueLu, but we are currently encountering issues with it.

[jhux2]

@dgmaths9 I wonder if the ILU settings are different between the GMRES/RILUK and GMRES/MueLu results. Please post the multigrid summary (section labeled "Summary") from the GMRES/MueLu/ILU run. If you have setup output from the successful GMRES/RILUK output, please post that as well, thanks.

[dgmaths9]

Hi @jhux2,

thanks for chiming in, here is the multigrid summary section from the GMRES/MueLu/ILU run:

---

--- Multigrid Summary ---

Scalar = double
Number of levels = 6
Operator complexity = 2.15
Smoother complexity = 0.00
Cycle type = V
level rows nnz nnz/row c ratio procs
0 567918 11867040 20.90 72
1 172617 6099011 35.33 3.29 72
2 28581 1144923 40.06 6.04 57
3 22617 2563156 113.33 1.26 45
4 14214 2288730 161.02 1.59 28
5 7347 1571238 213.86 1.93 14
Smoother (level 0) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)' ***** Condition number estimate = 44472.2 Smoother (level 1) both : Ifpack_AdditiveSchwarz, ov = 0, local solver = ***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 70724.2
Smoother (level 2) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)' ***** Condition number estimate = 78412.9 Smoother (level 3) both : Ifpack_AdditiveSchwarz, ov = 0, local solver = ***** IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 73875.8
Smoother (level 4) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 70069.1
Smoother (level 5) pre : MueLu::AmesosSmoother{type = Klu}
Smoother (level 5) post : no smoother

and here is GMRES/RILUK setup output:

Block Size = 1 [unused]
Maximum Restarts = 50 [unused]
Maximum Iterations = 3000 [unused]
Convergence Tolerance = 1e-14 [unused]
Num Blocks = 500 [unused]

[jhux2]

And what does MueLu + RILUK from Ifpack2 look like?

[dgmaths9]

verbosity = high
coarse: max size = 10000
multigrid algorithm = pg
problem: symmetric = 0
number of equations = 3
smoother: type = RILUK
aggregation: type = uncoupled
aggregation: min agg size = 1
aggregation: max agg size = 1
repartition: enable = 1
repartition: partitioner = zoltan
repartition: start level = 2
repartition: min rows per proc = 500
repartition: max imbalance = 1.1
repartition: remap parts = 0
max levels = 10 [default]
rap: algorithm = galerkin [default]
smoother: params ->
fact: level-of-fill = 0 [unused]
relaxation: type = Jacobi [unused]
relaxation: sweeps = 1 [unused]
relaxation: damping factor = 0.9 [unused]

Clearing old data (if any)
Replacing coordinate map
Using default factory (AmalgamationFactory[234]) for building 'UnAmalgamationInfo'.

Level 0
Setup Smoother (MueLu::IfpackSmoother{type = ILU})
Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 43746.5
Using default factory (AmalgamationFactory[243]) for building 'UnAmalgamationInfo'.

Level 1
Build (MueLu::RebalanceTransferFactory)
Prolongator smoothing (PG-AMG) (MueLu::PgPFactory)
Build (MueLu::TentativePFactory)
Build (MueLu::UncoupledAggregationFactory)
Build (MueLu::CoalesceDropFactory)
Build (MueLu::AmalgamationFactory)
...

---

--- Multigrid Summary ---

Scalar = double
Number of levels = 6
Operator complexity = 2.20
Smoother complexity = 0.00
Cycle type = V

level rows nnz nnz/row c ratio procs
0 567918 11867040 20.90 72
1 172617 6098542 35.33 3.29 72
2 28593 1146114 40.08 6.04 57
3 23217 2814849 121.24 1.23 46
4 14508 2636374 181.72 1.60 29
5 6732 1551132 230.41 2.16 13

Smoother (level 0) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 43746.5

Smoother (level 1) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 70568.4

Smoother (level 2) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 77857.7

Smoother (level 3) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 96758.2

Smoother (level 4) both : Ifpack_AdditiveSchwarz, ov = 0, local solver =
***** `IFPACK ILU (fill=0, relax=0.000000, athr=0.000000, rthr=1.000000)'
***** Condition number estimate = 89727.1

Smoother (level 5) pre : MueLu::AmesosSmoother{type = Klu}
Smoother (level 5) post : no smoother

[jhux2]

Hmmm, that's not very revealing. You could check against the converging GMRES/RILUK (no AMG) to see if there's an obvious difference in the settings.