Leave boundary dofs out of matrix construction

Asked by Nico Schlömer on 2013-05-05

I'm implementing a Navier--Stokes solver for axisymmetric geometries, and in the weak forumulation, the term

   u[0] * v[0] / r^2

appears (cf. <https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Cylindrical_coordinates>). If the domain contains r=0 -- which is perfectly legal --, e.g. the unit square, the 1/r^2 looks a little frightening at first. This isn't critical if the trial space and test space only contain functions at least proportional to r, corresponding to homogenous Dirichlet conditions for u[0] at r=0. This is quite easy to enforce on the test space (namely by applying said Dirichlet boundary conditions).
In Dolfin, the resulting matrix would look like

   [I 0 ] [u_r] = [0]
   [b A] [u ] [f]

where b corresponds to the terms where u[0] is nonzero at the boundary, and thus b would contain NaNs. They don't play a particular role since u_r is 0 there, but I don't if every backend picks that up. Also, applying some sort of LU solver, in fact any method that does something with the matrix entries (e.g., AMG), is asking for all sorts of trouble.

Question:
Is is possible to restrict the trial space as well to directly get an equation system of the form

  [A][u] = [f]

?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Garth Wells
Solved:
2013-05-08
Last query:
2013-05-08
Last reply:
2013-05-05
Jan Blechta (blechta) said : #1

On Sun, 05 May 2013 02:21:17 -0000
Nico Schlömer <email address hidden> wrote:
> New question #228170 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/228170
>
> I'm implementing a Navier--Stokes solver for axisymmetric geometries,
> and in the weak forumulation, the term
>
> u[0] * v[0] / r^2

You should simplify bad terms by multiplication by volume element
2*pi*r*dx - UFL does not reduce fractions except very special cases.

>
> appears (cf.
> <https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Cylindrical_coordinates>).

Note this is a Laplace form - stress form is more suitable for some
problems.

> If the domain contains r=0 -- which is perfectly legal --, e.g. the
> unit square, the 1/r^2 looks a little frightening at first. This
> isn't critical if the trial space and test space only contain
> functions at least proportional to r, corresponding to homogenous
> Dirichlet conditions for u[0] at r=0. This is quite easy to enforce
> on the test space (namely by applying said Dirichlet boundary
> conditions). In Dolfin, the resulting matrix would look like
>
> [I 0 ] [u_r] = [0]
> [b A] [u ] [f]
>
> where b corresponds to the terms where u[0] is nonzero at the
> boundary, and thus b would contain NaNs. They don't play a particular

I think there are no NaNs as Gauss qaudrature is used therefore 1/r
term in not evaluated at r=0. When Dirichlet conditions for u_r are to
be enforced on r=0 bad terms are not evaluated at r=0, only
corresponding dofs are found and manipulated in linear system's matrix
and vector.

> role since u_r is 0 there, but I don't if every backend picks that
> up. Also, applying some sort of LU solver, in fact any method that
> does something with the matrix entries (e.g., AMG), is asking for all
> sorts of trouble.
>
> Question:
> Is is possible to restrict the trial space as well to directly get an
> equation system of the form
>
> [A][u] = [f]
>
> ?
>

Nico Schlömer (nschloe) said : #2

> You should simplify bad terms by multiplication by volume element
> 2*pi*r*dx - UFL does not reduce fractions except very special cases.

Makes sense. -- After all, looking at the problem from a 3D angle, one would really like to weight the equations (and with it the residual) with the volume of the rotational 3D element.

> I think there are no NaNs as Gauss qaudrature is used therefore 1/r term in not evaluated at r=0.

Ah, okay. In face, the PDE isn't defined in r=0 either, so the only thing one has to take care of is that the functional doesn't grow indefinitely towards r=0. Is there any documentation about the Gauss points?

Best Garth Wells (garth-wells) said : #3

On 5 May 2013 21:11, Nico Schlömer <email address hidden> wrote:
> Question #228170 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/228170
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
>> You should simplify bad terms by multiplication by volume element
>> 2*pi*r*dx - UFL does not reduce fractions except very special cases.
>
> Makes sense. -- After all, looking at the problem from a 3D angle, one
> would really like to weight the equations (and with it the residual)
> with the volume of the rotational 3D element.
>
>> I think there are no NaNs as Gauss qaudrature is used therefore 1/r
> term in not evaluated at r=0.
>
> Ah, okay. In face, the PDE isn't defined in r=0 either, so the only
> thing one has to take care of is that the functional doesn't grow
> indefinitely towards r=0. Is there any documentation about the Gauss
> points?
>

Take a look at the FFC code:

    ffc/quadrature_schemes.py

Garth

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Nico Schlömer (nschloe) said : #4

Thanks Garth Wells, that solved my question.