Augmented Lagrangian Method for Asymetric Stokes Problem

[Félix Bunel]

Hi ! I'm Félix from France and I'm here to ask for help for a numerical project i'm working on.

I'm trying to solve the following system of equation in 2D. It's the so called Leslie-Eriksen problem. It's written here in it's simplified form and is in reality a bit more complex. But I would like to solve this in this form as a start. Fx and Fy are also fonction of phi. (more details in the attached PDF). Phi_lin can be quite large (from 1 to 100 ) and the parameter a can be from 1 to 4...

I'm planning to solve this problem by separting the equation on \phi from the stokes problem. My base example for this procedure is the tutorial step 31 which solve the boussinesq equation (Temperature T and Stokes problem).

So my plan is the following for each timestep :

- Solve equation 7 in phi, knowing u and v from the previous timestep

- Solve the stokes problem (8,9,10) in u, v, p, knowing phi

Solving equation (1) knowing u and v is not an issue. I can already solve it in its more complex form using a Trust Region Algorithm. In it's form written here it's linear and easy to solve but when I put the anysotropic complex terms it's get higly non-linear and needs at least a Newton method. You can go take a look at the PDF I put here if you want (sorry it's in french).

Now the problem I have is solving the stokes problem. As you can see, the problem is linear in u and v. I was planning in using Step 20, 21, 22 as an example in solving the stokes problem using a Schur Preconditioning. However, when I derived the weak formulation of my problem :

I realized that this was no longer a symmetrical problem because of the term in \philin. All the other terms are the usual stokes terms and are symmetric.

Fortunately for me, the dealII doc has an answer for everything and there is an example of using a Schur Preconditioning for an asymmetric matrix : step 57 which solves navier stokes equation using a newton method. It solves the problem using an augmented lagrangian based approach.

So here are my question. Question 1, 2 and 3 is just some confirmation to things I think I already know. Real question is number 4.

1 - Is there a simpler way to solve my stokes problem than using a Schur Preconditioner and the augmented lagrangian method ?

2 - Will the augmented Lagrangien Method also works for my terms ? I don't see an issue from the theory... In the tutorial and in the paper referenced it's used to solve the antisymetric terms that comes from the advectiv terms.

3 - In my case I have eta  = 1 in comparison to step 57 ?

4 - I don't understand how I should choose the gamma parameter and why its value does not change the results of the resolution. In step 57, it's value is set to 1 and does not change at all...

5 - Why is the Stokes FE is always set to :
    stokes_fe(FE_Q<2>(degree+1), 2,
              FE_Q<2>(degree), 1)

Can it work with  :

    stokes_fe(FE_Q<2>(degree), 2,
              FE_Q<2>(degree), 1)

I'm asking this because my main interest is in phi and with the degree+1, I have 8 times more dofs describing the flow than the ones describing phi.

Thanks in advance,

[Wolfgang Bangerth]

On 6/24/19 3:59 AM, Félix Bunel wrote:
>
> Fortunately for me, the dealII doc has an answer for everything

:-)

> 1 - Is there a simpler way to solve my stokes problem than using a Schur
> Preconditioner and the augmented lagrangian method ?

There are plenty of simpler methods, such as just using a direct solver for
the matrix, or using some kind of ILU decomposition as a preconditioner. The
problem is just that these are not very efficient (ILU) or don't scale to
large problems (direct solvers) or don't easily work in parallel (both). But
if all you want to do is play with the equations, then start with a direct
solver and replace it by something else when you go to large problems where
the direct solver is no longer sufficient. I often advise to use the direct
solver in the beginning because it's a 5-line piece of code that always work.

> 2 - Will the augmented Lagrangien Method also works for my terms ? I don't see
> an issue from the theory... In the tutorial and in the paper referenced it's
> used to solve the antisymetric terms that comes from the advectiv terms.

I don't know enough about this topic to say for sure, but my understanding is
that the penalty term step-57 is adding just improve the spectrum of the
matrix to make the (iterative) solver more efficient. This may or may not work
for your equations -- it's worth playing with it.

> 3 - In my case I have eta  = 1 in comparison to step 57 ?

Yes.

> 4 - I don't understand how I should choose the gamma parameter and why its
> value does not change the results of the resolution. In step 57, it's value is
> set to 1 and does not change at all...

I don't know enough about that as well. You might want to look into Timo
Heister's PhD thesis.

> 5 - Why is the Stokes FE is always set to :
>     stokes_fe(FE_Q<2>(degree+1), 2,
>               FE_Q<2>(degree), 1)
>
>
> Can it work with  :
>
>     stokes_fe(FE_Q<2>(degree), 2,
>               FE_Q<2>(degree), 1)

You will want to read up on the theory of saddle point problems (of which the
Stokes system is one). There, it is necessary to satisfy the so-called inf-sup
(or LBB, or Babuska-Brezzi) condition for the solution to be stable. Using the
first of the choice above, this condition is satisfied, but for the second
choice it is not.

> I'm asking this because my main interest is in phi and with the degree+1, I
> have 8 times more dofs describing the flow than the ones describing phi.

Then use degree+1 for phi as well :-) That's what we do in several of our codes.

Best
W.


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

[Félix Bunel]

Thank you for your answer !

I tried a direct solver and it sucessfully solved my problem. That way I can start working on the physical problem and I will see later for optimization.

As a thank you, here is a nice velocity field outputed by my code ! I don't think it's the correct solution because it has a weird symmetry (4 lobes), but that gives me something to work on. I believe I made a mistake in the writing of the weak formulation.

Thanks again for your answer and your lecture video which taught me a lot.

I hope I will be able to finish this project and publish something with the results. I will be sure to cite dealii if that happens again. (last project with dealii is already published  DOI:10.1103/PhysRevE.96.012705 )

[Wolfgang Bangerth]

On 6/28/19 6:36 AM, Félix Bunel wrote:
>
> I tried a direct solver and it sucessfully solved my problem. That way I can
> start working on the physical problem and I will see later for optimization.

Glad to hear that works now. I think using a direct solver for first
experiments is generally a good idea. Get things to work, and worry about
making it fast at a later time.

> As a thank you, here is a nice velocity field outputed by my code !

:-) Looks like a spiral galaxy :-)

> I hope I will be able to finish this project and publish something with the
> results. I will be sure to cite dealii if that happens again. (last project
> with dealii is already published  DOI:10.1103/PhysRevE.96.012705 )

Thanks in advance! (We already list your last project, I see.)