format/ordering of Stokes matrix

Discussion about coding and new developments
Post Reply
jonasthies
Posts: 7
Joined: 02 Sep 2011, 14:17
Antispam: Yes

format/ordering of Stokes matrix

Post by jonasthies »

Hi!

I'm trying to put together a linear system solver for the Stokes problem that rolls out of Elmer,
in the simplest case just with constant viscosity in 2D. I observe the following in the stiffness matrix:

* P2/P1 triangles - the matrix has 3x (num nodes) rows/cols - I would expect about 2.5* (num nodes)
as the pressure uses only 3 of the 6 nodes in a P2 element? I guess this is to keep the matrix structure conform
with the stabilized case?

* P2P1, the P/P part (in Matlab notation K(3:3:end,3:3:end) is not initialized to 0 but contains entries of order 10^-300,
a bug?

* There are substantial nonzero couplings between u and v velocities, even though Convect=False. (so K(1:3:end,2:3:end) is not empty).
I *think* this shouldn't be so, at least for regular P2/P1 FEM? Am I mistaken to assume 'interleaved' ordering like this?

Jonas
Juha
Site Admin
Posts: 357
Joined: 21 Aug 2009, 15:11

Re: format/ordering of Stokes matrix

Post by Juha »

Hi Jonas,

the P2/P1 element implementation using the nodal elements constraints the edge node
pressures such that the pressure is linear over the element. The constraint is added to
the local system, and no effort to eliminate (condensate) the corresponding degrees of
freedom is made. This, among other things, makes post processing output easier although
also the linear system somewhat larger.

The P/P part is nullified (and the constraint added) explicitely in NavierStokes.src at code
block starting at about line 980 (?)

BTW. Instead of the Equation section keyword "NS Convect=False" i'd recommend using
the Solver section Keyword "Flow Model=Stokes".

The code discretizes the operator div(2\mu*\epsilon), not the operator div(mu*grad(u)),
hence the coupling terms between the velocity components. This has an effect on the
natural boundary conditions, even if these otherwise would be equivalent. I could add
an option to allow using the latter form instead of the former...

I would also recommend looking at the fem/tests/StokesPFEM test case and the
(Navier-)Stokes solver there for easier testing of different options. It can do the
P2/P1 stabilization using P-elements, where it's a lot more natural than using the
nodal elements. That is only a very basic N-S implementation lacking most of the
options within the FlowSolver though.

-Juha
jonasthies
Posts: 7
Joined: 02 Sep 2011, 14:17
Antispam: Yes

Re: format/ordering of Stokes matrix

Post by jonasthies »

Hello Juha,

thanks for the reply, that's very illuminating. I still get uninitialized entries in the global matrix
but don't see why. I'll get back to you if I have more questionmarks.

Jonas
Post Reply