Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

Hello everyone,

I am currently solving a GLS stabilized form of the Navier-Stokes equation using DEALII.

The residual of the system looks similar to the regular Incompressible Navier-Stokes, except that a stiffness matrix that is dependent on the element size is added to the P-P block.

I have a great success / scaling when solving Manufactured Solutions and I have been able to reach relatively decent P1-P! meshes (say 2056x2056 or 256x256x256) on a workstation computer.

I am currently using the GMRES Trilinos Wrapper with ILU preconditioning with a relatively high fill-in level (4).

One of the issue of my system is that the pressure is defined up to a constant. On coarse mesh this does not affect the GMRES solver. However, on finer mesh, it seems that the GMRES Solver is greatly affected by this near-singularity of the matrix system.

I have often read in the literature that for stabilized method, the best way was to remove the "mode" associated to a constant pressure. I believe this implies some sort of projection of the residual in a space without a pressure constant?

I know this is a very broad question, but how would I go and implement such a thing in my solver?

Are there any examples in DEALII where a Poisson problem is solved but with strictly Neumann boundaries? That would be a very similar problem that could guide me in the right direction.

Thank you very much for your help. The support of this forum is greatly appreciated.

[Bruno Turcksin]

Bruno,

On Tuesday, March 5, 2019 at 9:00:13 AM UTC-5, Bruno Blais wrote:

Are there any examples in DEALII where a Poisson problem is solved but with strictly Neumann boundaries? That would be a very similar problem that could guide me in the right direction.

Yes, take a look at step-11. I think one of the tutorial solving the Stokes equation also shows how to deal with that problem. If you can use it, method 1, i.e. just fix one dof, is the easiest.

Best,

Bruno

[Wolfgang Bangerth]

> One of the issue of my system is that the pressure is defined up to a
> constant. On coarse mesh this does not affect the GMRES solver. However, on
> finer mesh, it seems that the GMRES Solver is greatly affected by this
> near-singularity of the matrix system.
> I have often read in the literature that for stabilized method, the best way
> was to remove the "mode" associated to a constant pressure. I believe this
> implies some sort of projection of the residual in a space without a pressure
> constant?

There are two parts of this problem:

1/ A deep theorem in linear algebra states that because the constant pressures
are in the kernel of the matrix (i.e., Ay=0 for all vectors y that correspond
to a constant pressure and zero velocity), that the *range* of the matrix A
has dimension at most n-1. As a consequence, if you have a linear system
A x = f
then it is only possible to find a vector x with
|| A x - f || = 0
if f is in the range of the matrix A. That will not generally be the case, due
to discretization and integration errors. If f is not in the range of A, there
is no way for GMRES to reduce the residual below a certain threshold, no
matter how many iterations one runs.

The way to solve this is to solve
A x = f - Pf
instead where Pf is the projection onto the subspace not reachable by A. This
is easy enough if you have a uniform mesh, but it's a bit complicated if
that's not the case. It's also complicated if one uses an enriched pressure
space. Here is ASPECT's implementation of this step:
https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
You will be able to copy this and simplify it for your case.


2/ There is now a null space and GMRES will find one of the infinitely many
solutions of
A x = f
This may not be the solution you are looking for, so you probably want to
update the additive constant in the pressure after solution. How you want to
do this depends on the application. Here again is the implementation in ASPECT:
https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L753

I hope this helps! Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Bruno Blais]

Dear Bruno and Wolfgang,

Thank you for your answers.
 I believe Wolfgang's answer is exactly what i had in mind (but said in clear words...). I will look at the Aspect code and try to port that to mine.

Thank you for the very detailed answer.

Best

Bruno

[Bruno Blais]

Dear Wolfgang,

A quick question. I think I understand what is done in

https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041

Weirdfully, I found the "pickaxe" version easier to understand (Love the comments by the way, this is awesome)

A question I have is related to the pressure_shape_function_integrals member.

My understanding is that this is the integration over the element cellI of the shape function associated with pressure.

However, I have looked over the ASPECT documentation and I have not been able to find where this is implemented. Is it a vector that is manually filled or is there a helper function that can be used to automatically generate it?

Thanks for everything

Bruno

[Wolfgang Bangerth]

> A quick question. I think I understand what is done in
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
> Weirdfully, I found the "pickaxe" version easier to understand (Love the
> comments by the way, this is awesome)

:-) I suspect that Timo gets credit for that one!


> A question I have is related to the _pressure_shape_function_integrals_ member.

> My understanding is that this is the integration over the element cellI of the
> shape function associated with pressure.

I believe it's actually the integral over all cells the shape function is
associated with.

> However, I have looked over the ASPECT documentation and I have not been able
> to find where this is implemented. Is it a vector that is manually filled or
> is there a helper function that can be used to automatically generate it?

Yes, that required a bit of searching. This is filled here:
https://github.com/geodynamics/aspect/blob/master/source/simulator/assemblers/stokes.cc#L601

[Bruno Blais]

Dear Wolfgang,

Sorry for an additional question.

I have tried to implement by myself something similar to what is done in ASPECT (I guess?, it is a bit harder for me to grasp ASPECT concepts because of the introspection and etc.).

However, I seem to have failed miserably at one point.

My understanding is that I can get

   A x = f- Pf

By:

- First I calculate my projection operator. It is the integral of the pressure shape function over ALL elements associated with a DOF. I calculate it by looping over the elements and integrating with a quadrature (see code below)

- Then I assemble my system and my right-hand side as usual

- Before I solve my linear system, I modify my Right-hand-side so that it does not include the constant mode. I litteraly calculate f= f- Pf

- I solve my linear system for this new RHS.

However, my GMRES stops very quickly after a certain number of newton iteration with the following :     AztecOO::Iterate error code -4: GMRES Hessenberg ill-conditioned

Clearly I am doing something wrong / I understood something wrong regarding this, because I know the rest of my code works.

In my system there is no block matrices and everything is lumped in the same matrix. I havent re-ordered by components either.

Calculation of the projection:

template <int dim>
void
GLSNavierStokesSolver<dim>::initializePressureRHSCorrection()
{
  TimerOutput::Scope t(computing_timer, "pressure_RHS_projector");
  pressure_shape_function_integrals=0;
  QGauss<dim>   quadrature_formula(degreeQuadrature_);
  FEValues<dim> fe_values (fe,
                           quadrature_formula,
                           update_values |
                           update_quadrature_points |
                           update_JxW_values );

  const unsigned int                    dofs_per_cell = fe.dofs_per_cell;
  const unsigned int                    n_q_points    = quadrature_formula.size();
  const FEValuesExtractors::Scalar      pressure (dim);
  Vector<double>                        local_pressure_shape_function_integrals(dofs_per_cell);
  std::vector<types::global_dof_index>  local_dof_indices (dofs_per_cell);


  typename DoFHandler<dim>::active_cell_iterator
  cell = dof_handler.begin_active(),
  endc = dof_handler.end();
  for (; cell!=endc; ++cell)
    {
      if (cell->is_locally_owned())
        {
          fe_values.reinit(cell);
          local_pressure_shape_function_integrals=0;
          for (unsigned int i=0; i<dofs_per_cell;++i)
            {
              for (unsigned int q=0; q<n_q_points; ++q)
                {
                  const double phi_p = fe_values[pressure].value(i, q);
                  local_pressure_shape_function_integrals(i) += phi_p * fe_values.JxW(q);
                }
            }
          cell->get_dof_indices (local_dof_indices);
          zero_constraints.distribute_local_to_global(local_pressure_shape_function_integrals,
                                                      local_dof_indices,
                                                      pressure_shape_function_integrals);
        }
    }
  pressure_shape_function_integrals.compress (VectorOperation::add);
}

Projection of the System RHS onto it's space without the constant:

template <int dim>
void GLSNavierStokesSolver<dim>::correctRHSClosedSystem()
{


  // calculate projection of system_rhs on pressure_shape_function_integrals
  for (unsigned int dof_i=0 ; dof_i < pressure_shape_function_integrals.size() ; ++dof_i)
    pressure_projection[dof_i] = pressure_shape_function_integrals[dof_i]*system_rhs[dof_i];

  // Add it to RHS
  system_rhs.add(-1.,pressure_projection);
}

Thank you again for everything,

This is greatly appreciated.

Bruno

[Wolfgang Bangerth]

On 3/7/19 7:45 AM, Bruno Blais wrote:
>
> However, my GMRES stops very quickly after a certain number of newton
> iteration with the following :     AztecOO::Iterate error code -4: GMRES
> Hessenberg ill-conditioned

This is surprising and suggests that something is wrong with the matrix, not
the right hand side. At least that's what I think this probably means. Did you
have the same error when you did not modify the right hand side? What happens
if you just artificially create a right hand side that you know is in the
range, for example by picking a random vector and multiplying it by the matrix?

[Bruno Blais]

On Thursday, 7 March 2019 23:35:15 UTC-5, Wolfgang Bangerth wrote:

On 3/7/19 7:45 AM, Bruno Blais wrote:
>
> However, my GMRES stops very quickly after a certain number of newton
> iteration with the following :     AztecOO::Iterate error code -4: GMRES
> Hessenberg ill-conditioned

This is surprising and suggests that something is wrong with the matrix, not
the right hand side. At least that's what I think this probably means. Did you
have the same error when you did not modify the right hand side? What happens
if you just artificially create a right hand side that you know is in the
range, for example by picking a random vector and multiplying it by the matrix?

I did not have the same error when I did not modify the RHS. Which I also find very surprising.

I also found that I did not get the same error if I modify the RHS, but that I also modify the Newton correction (if I project my newton correction such that corr = corr - P*corr)

I will look at what you suggested, that's an extremely good idea, and get back at you.

Thanks!