Dear Heitor,
>
> 1. Do you mean that I need only to modify step-31
> assemble_stokes_preconditioner in order to build a
> stokes_preconditioner_matrix S in accordance with your suggestions,
> i.e., step-31 will do the rest of the work once provided with this new
> S matrix?
Basically yes. It depends a bit on how you implement the preconditioner,
see my answer to question 3 below. If you follow my previous suggestion,
you just construct one additional matrix and fill it with the Laplacian
terms.
> 2. I have already a pressure mass matrix M_p ready as block(1,1) of my
> global mass matrix M; so I need to compute a laplacian pressure matrix
> block which I understand to be something like (phi_grad_p[i],
> phi_grad_p[j]) * stokes_fe_values.JxW(q) inserted in the same loop as
> in step-31 assemble_stokes_preconditioner; is that correct?
Yes, that would result in a pressure Laplace matrix.
> 3. How do I build S from M_p and L_p? You gave me the recipe for S^-1
> \approx L_p^-1 + dt * eta * (M_p)^-1, but as I understand this should
> not be read too literally: inverting L_p and M_p is forbidden! Besides
> the assemble_stokes_preconditioner expects S not S^-1; I am missing
> something ...
Here it gets interesting: You build the matrix S in
assemble_stokes_preconditioner. When you use a preconditioner (like ILU,
SSOR, etc), you are forming an approximation of the inverse. This means
that a multiplication with the preconditioner acts as S^{-1}.
What I suggested was to form a Laplace matrix on the pressure space and
a Laplace matrix on the pressure space, to create preconditioners for
both of them and then adding the action of the preconditioners, which is
doing something like
L_p^-1 + dt * eta * (M_p)^-1.
(You actually implement a vmult operation that first does the vmult of
the first preconditioner, then the vmult of the second one, and then
adds the vectors together).
> 4. As you hinted L_p is singular unless some constraint in the
> pressure space is imposed; in step-31 there are no boundary conditions
> for the pressure so I am left wondering what to do here ...
You can start by trying to form a pressure matrix as
(phi_grad_p[i] * phi_grad_p[j] + dt * eta * phi_p[i]*phi_p[j]) *
stokes_fe_values.JxW(q)
I'm not sure if that is really to give you a good approximation of the
Schur complement, but at least the matrix has full rank and you do not
need to set boundary conditions explicitly (you implicitly assume
Neumann boundary conditions if you do like this).
If we return to the question on how to treat L_p: I don't know a precise
answer. If I remember the experiments of some of my colleagues
correctly, you set homogeneous Neumann boundary conditions on L_p on all
boundaries where you have no-slip conditions on the velocity, and you
set homogeneous Dirichlet conditions on all boundaries where you have
natural conditions in the velocity system. However, I don't know any
serious theory behind that. It's probably best to try and look at the
results.
Best,
Martin
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
