### Discrete form of PSPG Oseen equation get a sign flipped?

50
views
0
23 days ago by
Hello all,
This might be a question more related to a Finite element problem rather than FEniCS, but since I implement in FEniCS and has nowhere else to ask questions, I could only ask it here. What is more, this may give a simple example for people who want to use block preconditioner.

What I am trying to do is that:
1. In general FEniCS, solve the Oseen equation with:
F = inner(grad(u),grad(v))*dx - div(v)*p*dx + inner(gread(u)*w,v)*dx
a,L = lhs(F),rhs(F)​
A Taylor-hood element was used and the body force f is 0.

2. Now I want to solve it with GMRES, thus in cbc.block, I implement:
a11 = inner(grad(u),grad(v))*dx
a12 = -div(v)*p*dx
a21 = -div(u)*q*dx
a22 = 0

A = block_assemble([[a11,  a12],
[a21, -a22]], bcs=bcs)​
Note here for a21 the sign flipped. Now my first question is why is it flipped?
I then used the preconditoner:
[[F, D],
[B, C]] = A
Fp = AMG(A)
Dp = AMG(collapse(B*InvDiag(F)*D-C))
#no approximate here for Schur component as it is a test​
prec = block_mat([[Fp,  D],
[0 , -Dp]).scheme('sgs')
Ainv = LGMRES(A, precond=prec, tolerance=1e-10, maxiter=50, show=2)
I got almost same results from these two different method.

Next step I want to use P1P1 elements with SUPG/PSPG, here just take PSPG for example.
F = inner(grad(u),grad(v))*dx - div(v)*p*dx + inner(gread(u)*w,v)*dx
a,L = lhs(F),rhs(F)​​

a11 = inner(grad(u),grad(v))*dx
a12 = -div(v)*p*dx

A = block_assemble([[a11,  a12],
[a21, -a22]], bcs=bcs)​

(The preconditioner part is same.)
Here as you can see, PSPG was divided into two parts and added to a21 and a22 respectively.
My second question is: as the sign of div(u)*q*dx in a21 is flipped, why is not the sign of PSPG flipped?

The equations are taken from chapter 4 of "John, V. (2016). Finite element methods for incompressible flow problems. Springer International Publishing."

Or you may find it at "https://www.wias-berlin.de/people/john/LEHRE/NUM_NSE_17_18/num_nse_1_3_4.pdf"

I also put the screenshot of his formula here:
weak form:

And discrete form:

Thanks,

Community: FEniCS Project