# Stokes PETSC solver but in a doubled-size space
**URL:** https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389
**Category:** General Discussion
**Created:** [April 11, 2023, 4:44pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389 "2023-04-11T16:44:16Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 11, 2023, 4:44pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/1 "2023-04-11T16:44:16Z")
Dear all
I learnt from [https://github.com/mbarzegary/navier-stokes-solver an efficient](https://github.com/mbarzegary/navier-stokes-solver) @prj @mojtaba.barzegari an efficient PETSC solver for Stokes solution
[u, v, w, p] in FEM space [P2, P2, P2, P1]:
> set(A, sparams = "-ksp\_monitor -ksp\_type fgmres -ksp\_converged\_reason "
> + "-pc\_type fieldsplit -pc\_fieldsplit\_type schur "
> + "-fieldsplit\_velocity\_pc\_type gamg -fieldsplit\_velocity\_pc\_gamg\_sym\_graph true -fieldsplit\_pressure\_ksp\_max\_it 5 "
> + “-fieldsplit\_pressure\_pc\_type jacobi -fieldsplit\_velocity\_ksp\_type preonly -pc\_fieldsplit\_schur\_fact\_type full”,
> fields = ux, names = names);
with
> string[int] names(2);
> names[0] = “velocity”;
> names[1] = “pressure”;
Now, I have **two velocity fields** and **two pressure fields** : [u1, v1, w1, u2, v2, w2, p1, p2] in a FEM space [P2,P2,P2,P2,P2,P2, P1,P1], **is the parameter set-up still the same?**
BTW, the second velocity or pressure field is from the adjoint space.
Many thanks in advance!
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 11, 2023, 5:45pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/2 "2023-04-11T17:45:17Z")
It all depends on how you want to precondition your linear system. But I do think it’s a bad idea to group everything together.
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 11, 2023, 5:54pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/3 "2023-04-11T17:54:51Z")
Thank you @prj .
The reason I solve the primal and adjoint euqations (which are **fully coupled** together) in a monolithic manner is because I want to **avoid iterations** from one another – convergence is also painful.
Look at the linear system **[[K, B’], [B, 0]],** which is quite similar to the Stokes system, but just both the sizes of K and B are doubled. So the preconditioner can also be similarly built from the Schur complement. Do you think the above parameter set will solve this in its most efficient manner?
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 11, 2023, 6:06pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/4 "2023-04-11T18:06:57Z")
I doubt it will be as efficient as solving this in a staggered way, but I’d love to be proven wrong!
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [April 11, 2023, 9:48pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/5 "2023-04-11T21:48:52Z")
Is it necessary for you to write a separate problem for your adjoint equations? I agree that the monolithic approach is probably not optimal.
Once you create your linear operator for the primal problem (to be solved with `KSPSolve()`), the discrete adjoint operator is given by its (conjugate) transpose and can be solved using `KPSSolveConjugateTranspose()` [complex] or `KSPSolveTranspose()` [real] without having to reconstruct the matrix or refactorize the preconditioner (though you will need to perform the linear iterations again to solve the adjoint system). Note that the same preconditioner may or may not work efficiently for the adjoint problem.
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 12, 2023, 8:42am UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/6 "2023-04-12T08:42:21Z")
Thank you @prj @cmd both. Will compare the performance of the monolithic and decoupled formulations in details…
BTW, my problem is a large-scale 3D FSI interaction, and the mesh is reconstructed in advance for every time frames. Will introduce this work in the CFC conference 2023 25-28 April at Cannes, hope to talk with you guys there if you will attend as well by any chance.
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 12, 2023, 1:34pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/7 "2023-04-12T13:34:05Z")
Dear @cmd @cmd : I implemented the decoupled version, however, the linear solver for the adjoint equations struggles to converge (although the primal equation solver converges in 10 to 20 iterations)… do you expect this?
[test.zip](https://community.freefem.org/uploads/short-url/jLyjfDdSl8yEhxB1RyBjFH9EBIk.zip) (1.4 MB)
The test is a pure fluid problem, not FSI. The adjoint equation couples the primal equations through a volume integral, while the primal equation couples the adjoint ones through a boundary condition.
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [April 12, 2023, 2:33pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/8 "2023-04-12T14:33:01Z")
Hi there. It is hard to say, but since you are using different BC’s for your primal and adjoint equations, I do understand why the convergence behavior could be quite different. Are you sure that your variational forms for the system are correct, including the bilinear concomitant in the adjoint system?
From what I can see, the equations are, in general, only one-way coupled: the primal system is independent of the adjoint one. Also, in the particular case you’ve sent of g=0, the primal system is homogeneous, and the solution is trivial. In this case, the adjoint system is also homogeneous and has a trivial solution.
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 12, 2023, 2:57pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/9 "2023-04-12T14:57:40Z")
The primal - adjoint system is derived from an optimal control problem: minimising the viscous dissipation by a boundary velocity profile, constrained by the PDEs (NS equations in this case). Using the Lagrange multiplier (adjoint variables) method to cancel out the constraints.
Imagine you take the variation with respect to the primal variables (which gives you the adjoint equations), the coupled-volume integral in the adjoint equations is from the variation of the objective function (viscous dissipation in this case).
The variation of the control variable (boundary velocity in this case) provides the optimality condition which relates the control variable to the adjoin variables, this would go back to the primal equations as a boundary condtion.
There was a previous paper, [A monolithic one-velocity-field optimal control formulation for fluid–structure interaction problems with large solid deformation - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/S0889974622000470)
using a distributed control, but is similar to a boundary control in the current code.
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [April 12, 2023, 3:41pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/10 "2023-04-12T15:41:06Z")
In the script you’ve sent, I see the following definitions:
```auto
varf primal([ux,uy,uz,p],[uhx,uhy,uhz,ph]) =
int3d(Th)( nuf/2*(DD(u):DD(uh))
-div(uh)*p-div(u)*ph - penal*p*ph
)
+ int3d(Th)(g*uhy)
+ on(1,2,3, ux=ubx, uy=uby, uz=ubz)
+ on(4,5,6,7,8,9,ux=0,uy=0,uz=0);
varf adjoint([vx,vy,vz,q],[vhx,vhy,vhz,qh]) =
int3d(Th)( nuf/2*(DD(v):DD(vh))
-div(vh)*q-div(v)*qh - penal*q*qh
)
+ int2d(Th,1,2,3)(vect(vh)'*vect(v))
+ int3d(Th)(- nuf/4*(DD(u):DD(vh)) )
+ on(4,5,6,7,8,9,vx=0,vy=0,vz=0);
```
where `g=0.0`, `ubx=uby=ubz=0.0`, `penal=0.0`, `nuf=1.0`. Hence, we could simplify the forms to:
```auto
varf primal([ux,uy,uz,p],[uhx,uhy,uhz,ph]) =
int3d(Th)( nuf/2*(DD(u):DD(uh))
-div(uh)*p-div(u)*ph
)
+ on(1,2,3,4,5,6,7,8,9,ux=0,uy=0,uz=0);
varf adjoint([vx,vy,vz,q],[vhx,vhy,vhz,qh]) =
int3d(Th)( nuf/2*(DD(v):DD(vh))
-div(vh)*q-div(v)*qh
)
+ int2d(Th,1,2,3)(vect(vh)'*vect(v))
+ int3d(Th)(- nuf/4*(DD(u):DD(vh)) )
+ on(4,5,6,7,8,9,vx=0,vy=0,vz=0);
```
Since the solution to the primal problem is trivial, there will be no forcing on the adjoint problem and that solution will also be trivial. If my understanding here is correct, maybe this isn’t a good test case?
Are you sure that the BC written for the adjoint problem on surfaces 1, 2, and 3 is well-posed?
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 12, 2023, 4:13pm UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/11 "2023-04-12T16:13:57Z")
Yes, your understanding is correct, except that [ubx, uby, ubz] is not zero, but
> ```
> ubx=ul2x+1./alpha*r(rl2)'*r(rl2)*vx - 1./alpha*r(rl2)'*vect(v)*(x-rl2x);
> uby=ul2y+1./alpha*r(rl2)'*r(rl2)*vy - 1./alpha*r(rl2)'*vect(v)*(y-rl2y);
> ubz=ul2z+1./alpha*r(rl2)'*r(rl2)*vz - 1./alpha*r(rl2)'*vect(v)*(z-rl2z);
>
> ```
where [ul2x, ul2y, ul2z] is given which is not zero – linear velocity measure from laboratory. The rest part is from the angular velocity, which is related to the ajoint solution, and only zero only in the first iteration step.
The boundary integral for the adjoint equation is from variation of L with respect to **u:**
with
is the solid domain (the worm), we can forget this at the moment and only consider a fluid problem in \Omeg - \Omega\_s.
The well possedness is a good question, it is difficult to analyse. But, compare with the classical boundary control problems, I just decomposed the boundary velocity to its linear and angular components. What do you think about the well posedness?
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [April 13, 2023, 11:09am UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/12 "2023-04-13T11:09:19Z")
Ah, I see. Thanks for clarifying. Unfortunately, I do not know about the well-posedness. You are enforcing some type of Robin-like condition resulting from the sum of the implicit Neumann conditions from integrating the viscous and pressure terms by parts, plus the `vect(vh)'*vect(v)` term. I’m not sure how this behaves.
---
### Author:  [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [April 13, 2023, 11:20am UTC](https://community.freefem.org/t/stokes-petsc-solver-but-in-a-doubled-size-space/2389/13 "2023-04-13T11:20:40Z")
Yes, I will implement a regularisation method (boundary force or Lagrange multiplier to enforce the Dirichlet:
[Dirichlet control.pdf](https://community.freefem.org/uploads/short-url/9QYPvawuZwKnZULsEkp0FfrqoKL.pdf) (2.8 MB)
)… in this case, the Adjoint equation would have a homogeneous boundary condition (but we have to solve a surface Laplace equation for the control), but it would still be interacting with the Primal equations through the volume integral, which is the reason to slow down the linear solver I think – I guess the numerical influence of the Robin boundary is trivial (although it may not be quite right mathematically).