dolfin team mailing list archive

[Alexander Jarosch <alexanj@xxxxx>]

Anders Logg wrote:

> On Tue, May 23, 2006 at 02:46:32PM +0000, Alexander Jarosch wrote:
>  
>
> > Hi,
> >
> > did anyone play around with the stokes solver on more complex 
> > geometries? I only seem to get something senseful using a stabalized 
> > stokes like:
> >
> > scalar = FiniteElement("Lagrange", "triangle", 1)
> > vector = FiniteElement("Vector Lagrange", "triangle", 2)
> > system = vector + scalar
> >
> > (v, q) = TestFunctions(system)
> > (u, p) = TrialFunctions(system)
> >
> > f = Function(vector)
> > h = Function(scalar)
> > nu = Function(scalar)
> >
> > beta  = 0.2
> > delta = beta*h*h
> >
> > a = (nu*dot(grad(v), grad(u)) - div(v)*p + q*div(u) + delta*dot(grad(q), 
> > grad(p)))*dx
> > L = dot(v + mult(delta, grad(q)), f)*dx
> >
> > which is fine, but when I compare the solution to another FEM package, 
> > they do not quite match. Was there already some benchmarking done?
> >
> > cheers,
> >
> > Alex
> >    
>
> Do you really need to stabilize here? You're already using Taylor-Hood
> elements anyway. Have you tried a stabilized equal-order method?
> (Change 2 to 1 in the example above.)
>
> I have only tested the Stokes solver on simple geometries and I
> haven't compared to any other solvers.
>
> /Anders
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>
>  
Woops sorry, I posted the wrong ffc form: I used a stabilized 
equal-order method. But I have to set the pressure on the fluid surface 
to get at least the same velocity patterns. I also tried a form like that:

scalar = FiniteElement("Lagrange", "triangle", 1)
vector = FiniteElement("Vector Lagrange", "triangle", 1)
system = vector + scalar

(v, q) = TestFunctions(system)
(u, p) = TrialFunctions(system)

f = Function(vector)
h = Function(scalar)
nu = Function(scalar)

beta  = 0.2
delta = beta*h*h

a = (nu*dot(grad(v), grad(u)) - div(v)*p + nu*q*div(u) + 
delta*dot(grad(q), grad(p)))*dx
L = dot(v + mult(delta, grad(q)), f)*dx

mind the nu*q*div(u) term. This form seems to produce the velocity 
patterns with just fixing the pressure to one point on the surface but 
has constant offsets in the velocity values. The Taylor - Hood U2/P1 is 
not stable on a complex geometry, the pressure fluctuates.

Alex

--
Alexander H. Jarosch

Jarðvísindastofnun Háskólans
Institute of Earth Sciences, University of Iceland
Náttúrufræðahús, Askja
Building of Natural Sciences, Askja
Sturlugata 7
IS - 101 Reykjavík
Iceland

Tel.: +354 525 4906
http://raunvis.hi.is/~jarosch/