dolfin team mailing list archive

["Johan Hoffman" <jhoffman@xxxxxxxxxx>]

Ok, very good!

/Johan

> I just increased the delta2 for the stabilization of the
> incompressibility and used the form as written on the very end of this
> email. beta2 is now 40 to 50 times smaller than the viscosity and it
> reproduces the velocities similar to if I use a nu*q*div(u) term in my
> form. So it seems the problem is due to stability of incompressibility
> as Johan pointed out. And now if I set the surface pressure = 0 as a
> boundary condition, and not only on one point, everything looks quite
> okay.
>
> Thanks allot for the input
>
> /Alex
>
> Alexander Jarosch wrote:
>
>> Johan Hoffman wrote:
>>
>>>> On Tue, May 30, 2006 at 10:07:47AM +0000, Alexander Jarosch wrote:
>>>>
>>>>
>>>>> Hi everybody.
>>>>>
>>>>> There was a short discussion about Stokes flow on complex
>>>>> geometries. I
>>>>> would like to report that a ffc form like the one below works well on
>>>>> complex geometries, I still only have a factor of 1.2 offset in
>>>>> magnitude if I compare it with another fem model.  I would like to
>>>>> get
>>>>> some feedback about this form, whether it is okay or not from the
>>>>> developers point of view.
>>>>>
>>>>> much appreciated,
>>>>>
>>>>> Alex
>>>>>
>>>>> ffc form -------------------
>>>>>
>>>>> 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
>>>>>
>>>>
>>>> I think it looks a little strange. The difference from what you had
>>>> before seems to be that you have included the viscosity in the
>>>> incompressibility term:
>>>>
>>>>    nu*q*div(u)
>>>>
>>>> Is this what makes the difference? If so, I don't know why this
>>>> helps. It's certainly ok to put it there since you're solving for
>>>> div(u) = 0 anyway, so nu*div(u) = 0 is also correct.
>>>>
>>>> In some way, it corresponds to a smaller penalty for compressibility
>>>> when the viscosity is small which might be reasonable.
>>>>
>>>> /Anders
>>>>
>>>
>>>
>>> I think this term is somewhat strange, and I haven't seen it before.
>>> I do
>>> not see why the viscosity should multiply the incompressibility.
>>>
>>> If the problem is connected to the stability of the incompressibility
>>> constraint, I propose to instead add to the left hand side of the form
>>> a
>>> stabilization term for the incompressibility, of the form:
>>>
>>> delta2*div(v)*div(u)
>>>
>>> with delta2 = C2*h*h, and C2 of unit order.
>>>
>>> Also, what is the value of your parameter beta in the definition of
>>> your
>>> delta? It might be worthhile to try to experiment with higher values of
>>> beta.
>>>
>>> /Johan
>>>
>>>
>>>
>>>
>>>
>>>> _______________________________________________
>>>> DOLFIN-dev mailing list
>>>> DOLFIN-dev@xxxxxxxxxx
>>>> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>>>>
>>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> DOLFIN-dev mailing list
>>> DOLFIN-dev@xxxxxxxxxx
>>> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>>>
>>>
>>>
>> Thanks for all the input. The idea of multiplying the viscosity with
>> the incompressibility was also strange to me, but it seems to make the
>> difference.
>> I tried a form proposed by Johan like this (is that what you had in
>> mind?):
>>
>> beta  = 0.4
>> beta2 = 0.4
>> delta = beta*h*h
>> delta2 = beta2*h*h
>>
>> a = (nu*dot(grad(v), grad(u)) - div(v)*p + q*div(u) +
>> delta*dot(grad(q), grad(p)) + delta2*div(v)*div(u))*dx
>> L = dot(v + mult(delta, grad(q)), f)*dx
>>
>> and also played around with beta and beta2 values. It did not fix the
>> problem, which is: I work at very high viscosities, 1e13 Pa s and my
>> surface is curved. Now if I use the form above or a normal stokes
>> form, the pressure on the surface, where it should be 0 is not 0, the
>> pressure is 0 on the endpoints, but not along the curved surface. The
>> pressure seems not to realize that the surface is curved. Now if I
>> multiply the viscosity into the incompressibility term like
>>
>> nu*q*div(u)
>>
>> suddenly the pressure is correct on the surface boundary, and the
>> velocities as well, as mentioned.
>>
>> Any idea why this could be?
>>
>> Thanks in advance
>>
>> /Alex
>>
>> _______________________________________________
>> DOLFIN-dev mailing list
>> DOLFIN-dev@xxxxxxxxxx
>> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>>
>
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>