1d supg scheme

Asked by Derek Monahan

I am trying to reproduce the 1d SUPG advec-diff example in ch. 2 of "Finite element for flow problems" by Donea & Huerta, however I can't recover their solution. The problem is to solve the equation

a*u_x - nu*u_xx = s, where s = 5.0*exp(-100*(x-1/8)**2) - 5.0*exp(-100*(x-1/4)**2)

I am assuming order 1 CG elements and so the second derivative term in the stabalization is dropped and the element by element summation of the stablizing integral is replaced by an integral over the entire domain. My FEniCS implementation is as follows (compare to fig 2.13 in Donea & Huerta):

    a = 1.0 # vel coef (1d a vector = scalar)
    nu = 0.01 # diff. coef => Pe = 5.0

    # domain setup
    mesh = UnitInterval(10)
    V = FunctionSpace(mesh, "CG", 1)
    def u0_boundary(x, on_boundary):
        return on_boundary
    bc = DirichletBC(V, 0.0, u0_boundary)
    u = TrialFunction(V)
    w = TestFunction(V)

    # stabilization coeficient, tau
    h = mesh.hmin() # uniform mesh
    Pe = a*h/2.0/nu; beta = 1.0/numpy.tanh(Pe) - 1/Pe
    nu2 = beta*a*h/2.0; tau = nu2/(a*a)

    # variational forms
    s = Expression('5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) - 5*exp(-100*(x[0]-1./4)*(x[0]-1./4))')
    a_ = w*dot(a, nabla_grad(u))*dx + inner(nabla_grad(w), nu*nabla_grad(u))*dx
    L_ = w*s*dx

    # stabilization term
    Lu_ = dot(a, nabla_grad(u)) # linear elem => div(nabla_grad(u)) == 0
    Pw_ = dot(a, nabla_grad(w)) # stabal. operator

    # stabilized forms
    a_ += tau*Pw_*Lu_*dx
    L_ += tau*Pw_*s*dx

    # solve
    u = Function(V)
    solve(a_ == L_, u, bc)
    plot(u).interactive()

The value I calculate on the 0.2 node is 0.44, where as in the implentation in Donea & Huerta is ~ 0.5. I've clearly mis-understood somthing here.

I've searched for SUPG on the form and got no results and I've gone through the Ca2+ example in the FEniCS book, but it seems to be a bit over kill for this simple problem. Any help/suggestions would be very much appricated.

Derek

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Martin Sandve Alnæs
Solved:
Last query:
Last reply:
Revision history for this message
Martin Sandve Alnæs (martinal) said :
#1

One possible culprit is the Expression s, which will be interpolated into a
piecewise polynomial space. Try setting

x = mesh.ufl_cell().x

And change
s = Expression('...')
into just
s = ...

This will probably give you too high quadrature degree, but try it first
and see.

Martin
Den 28. feb. 2013 18:16 skrev "Derek Monahan" <
<email address hidden>> følgende:

> New question #223092 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/223092
>
> I am trying to reproduce the 1d SUPG advec-diff example in ch. 2 of
> "Finite element for flow problems" by Donea & Huerta, however I can't
> recover their solution. The problem is to solve the equation
>
> a*u_x - nu*u_xx = s, where s = 5.0*exp(-100*(x-1/8)**2) -
> 5.0*exp(-100*(x-1/4)**2)
>
> I am assuming order 1 CG elements and so the second derivative term in the
> stabalization is dropped and the element by element summation of the
> stablizing integral is replaced by an integral over the entire domain. My
> FEniCS implementation is as follows (compare to fig 2.13 in Donea & Huerta):
>
> a = 1.0 # vel coef (1d a vector = scalar)
> nu = 0.01 # diff. coef => Pe = 5.0
>
> # domain setup
> mesh = UnitInterval(10)
> V = FunctionSpace(mesh, "CG", 1)
> def u0_boundary(x, on_boundary):
> return on_boundary
> bc = DirichletBC(V, 0.0, u0_boundary)
> u = TrialFunction(V)
> w = TestFunction(V)
>
> # stabilization coeficient, tau
> h = mesh.hmin() # uniform mesh
> Pe = a*h/2.0/nu; beta = 1.0/numpy.tanh(Pe) - 1/Pe
> nu2 = beta*a*h/2.0; tau = nu2/(a*a)
>
> # variational forms
> s = Expression('5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) -
> 5*exp(-100*(x[0]-1./4)*(x[0]-1./4))')
> a_ = w*dot(a, nabla_grad(u))*dx + inner(nabla_grad(w),
> nu*nabla_grad(u))*dx
> L_ = w*s*dx
>
> # stabilization term
> Lu_ = dot(a, nabla_grad(u)) # linear elem => div(nabla_grad(u)) == 0
> Pw_ = dot(a, nabla_grad(w)) # stabal. operator
>
> # stabilized forms
> a_ += tau*Pw_*Lu_*dx
> L_ += tau*Pw_*s*dx
>
> # solve
> u = Function(V)
> solve(a_ == L_, u, bc)
> plot(u).interactive()
>
> The value I calculate on the 0.2 node is 0.44, where as in the
> implentation in Donea & Huerta is ~ 0.5. I've clearly mis-understood
> somthing here.
>
> I've searched for SUPG on the form and got no results and I've gone
> through the Ca2+ example in the FEniCS book, but it seems to be a bit over
> kill for this simple problem. Any help/suggestions would be very much
> appricated.
>
> Derek
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

Revision history for this message
Derek Monahan (derek-d-monahan) said :
#2

Hi Martin,

Thanks for such a quick reply.

I've tried your suggestion, and it failed to run, giving me the following Traceback:

Invalid number of indices (1) for tensor expression of rank 0:
 SpatialCoordinate(Cell('interval', Space(1)))

Traceback (most recent call last):
  File "supgAbriv.py", line 39, in supg_advec_diff
    s = 5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) - 5*exp(-100*(x[0]-1./4)*(x[0]-1./4))
  File "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/exproperators.py", line 310, in _getitem
    a = Indexed(self, indices)
  File "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/indexed.py", line 45, in __init__
    % (len(indices), expression.rank(), expression))
  File "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/log.py", line 148, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Invalid number of indices (1) for tensor expression of rank 0:
 SpatialCoordinate(Cell('interval', Space(1)))

Revision history for this message
Best Martin Sandve Alnæs (martinal) said :
#3

Ok, that's changed in trunk. You can either change x[0] to x in your
expression, or set
x = [mesh.ufl_cell().x]

You can also set the polynomial degree to interpolate the Expression into
by writing
s = Expression('...', degree=3)

Martin

On 28 February 2013 20:45, Derek Monahan <
<email address hidden>> wrote:

> Question #223092 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/223092
>
> Derek Monahan posted a new comment:
> Hi Martin,
>
> Thanks for such a quick reply.
>
> I've tried your suggestion, and it failed to run, giving me the
> following Traceback:
>
> Invalid number of indices (1) for tensor expression of rank 0:
> SpatialCoordinate(Cell('interval', Space(1)))
>
> Traceback (most recent call last):
> File "supgAbriv.py", line 39, in supg_advec_diff
> s = 5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) -
> 5*exp(-100*(x[0]-1./4)*(x[0]-1./4))
> File
> "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/exproperators.py",
> line 310, in _getitem
> a = Indexed(self, indices)
> File
> "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/indexed.py",
> line 45, in __init__
> % (len(indices), expression.rank(), expression))
> File
> "/scratch/dmonahan/FEniCS_install/installed_dir/FEniCS/lib/python2.6/site-packages/ufl/log.py",
> line 148, in error
> raise self._exception_type(self._format_raw(*message))
> ufl.log.UFLException: Invalid number of indices (1) for tensor expression
> of rank 0:
> SpatialCoordinate(Cell('interval', Space(1)))
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

Revision history for this message
Roland Siegbert (roland-siegbert) said :
#4

Hi,

Just solved the same problem a while ago. Can you send a plot? Your value sounds quite close.

Did you check a simpler expression for beta? I remember an approximation being drawn in the diagram in the book.

Best, Roland

Derek Monahan <email address hidden> wrote:

>New question #223092 on FEniCS Project:
>https://answers.launchpad.net/fenics/+question/223092
>
>I am trying to reproduce the 1d SUPG advec-diff example in ch. 2 of
>"Finite element for flow problems" by Donea & Huerta, however I can't
>recover their solution. The problem is to solve the equation
>
>a*u_x - nu*u_xx = s, where s = 5.0*exp(-100*(x-1/8)**2) -
>5.0*exp(-100*(x-1/4)**2)
>
>I am assuming order 1 CG elements and so the second derivative term in
>the stabalization is dropped and the element by element summation of
>the stablizing integral is replaced by an integral over the entire
>domain. My FEniCS implementation is as follows (compare to fig 2.13 in
>Donea & Huerta):
>
> a = 1.0 # vel coef (1d a vector = scalar)
> nu = 0.01 # diff. coef => Pe = 5.0
>
> # domain setup
> mesh = UnitInterval(10)
> V = FunctionSpace(mesh, "CG", 1)
> def u0_boundary(x, on_boundary):
> return on_boundary
> bc = DirichletBC(V, 0.0, u0_boundary)
> u = TrialFunction(V)
> w = TestFunction(V)
>
> # stabilization coeficient, tau
> h = mesh.hmin() # uniform mesh
> Pe = a*h/2.0/nu; beta = 1.0/numpy.tanh(Pe) - 1/Pe
> nu2 = beta*a*h/2.0; tau = nu2/(a*a)
>
> # variational forms
>s = Expression('5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) -
>5*exp(-100*(x[0]-1./4)*(x[0]-1./4))')
>a_ = w*dot(a, nabla_grad(u))*dx + inner(nabla_grad(w),
>nu*nabla_grad(u))*dx
> L_ = w*s*dx
>
> # stabilization term
> Lu_ = dot(a, nabla_grad(u)) # linear elem => div(nabla_grad(u)) == 0
> Pw_ = dot(a, nabla_grad(w)) # stabal. operator
>
> # stabilized forms
> a_ += tau*Pw_*Lu_*dx
> L_ += tau*Pw_*s*dx
>
> # solve
> u = Function(V)
> solve(a_ == L_, u, bc)
> plot(u).interactive()
>
>The value I calculate on the 0.2 node is 0.44, where as in the
>implentation in Donea & Huerta is ~ 0.5. I've clearly mis-understood
>somthing here.
>
>I've searched for SUPG on the form and got no results and I've gone
>through the Ca2+ example in the FEniCS book, but it seems to be a bit
>over kill for this simple problem. Any help/suggestions would be very
>much appricated.
>
>Derek
>
>--
>You received this question notification because you are a member of
>FEniCS Team, which is an answer contact for FEniCS Project.

--
Sent from the future. Please excuse my brevity.

Revision history for this message
Derek Monahan (derek-d-monahan) said :
#5

Great, both of those solutions worked - that is, setting:

    s = Expression('5.0*exp(-100*(x[0]-...', degree=3)

or:

    x = [mesh.ufl_cell().x] |
    s = 5*exp(-100*(x[0]-1./8)...)

Thank you so much for your help!

If it's not to much, could I ask what exactly this correction is doing? Is it replacing the continum 'x' in the expression with the linear element?

Hi Roland, I'm afraid I am unable to attach any figures to this post but my (corrected) results are:

beta = coth(Pe)-Pe:
  u_SU = [ -9.18e-17, 3.308e-01, 2.57e-01, -5.78e-02,
                 -1.37e-01, -1.40e-01, -1.40e-01, -1.40e-01,
                 -1.40e-01, -1.40e-01, 1.25e-16]
  u_SUPG - [ -9.52e-17, 3.43e-01, 5.01e-01, 1.27e-01,
                     -2.27e-02, -3.23e-02, -3.24e-02, -3.24e-02,
                     -3.24e-02, -3.24e-02, 2.88e-17]

beta = 1.0
  u_SU = [ -8.11e-17, 2.92e-01, 1.98e-01, -9.43e-02,
                 -1.66e-01, -1.70e-01, -1.70e-01, -1.70e-01,
                 -1.68e-01, -1.54e-01, 1.37e-16 ]
  u_SUPG = [ -9.14e-17, 3.29e-01, 4.93e-01, 1.24e-01,
                      -2.80e-02, -3.81e-02, -3.83e-02, -3.82e-02,
                      -3.79e-02, -3.48e-02, 3.09e-17]

Derek

Revision history for this message
Derek Monahan (derek-d-monahan) said :
#6

Thanks Martin Sandve Alnæs, that solved my question.

Revision history for this message
Martin Sandve Alnæs (martinal) said :
#7

Degree=3 causes interpolation to a CG3 space prior to assembly, while the
other solution is evaluation at all quadrature points. So there are
differences in accuracy and efficiency depending on the expression.

Martin
Den 1. mars 2013 13:11 skrev "Derek Monahan" <
<email address hidden>> følgende:

> Question #223092 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/223092
>
> Derek Monahan posted a new comment:
>
> Great, both of those solutions worked - that is, setting:
>
> s = Expression('5.0*exp(-100*(x[0]-...', degree=3)
>
> or:
>
> x = [mesh.ufl_cell().x]
> |
> s = 5*exp(-100*(x[0]-1./8)...)
>
> Thank you so much for your help!
>
> If it's not to much, could I ask what exactly this correction is doing?
> Is it replacing the continum 'x' in the expression with the linear
> element?
>
> Hi Roland, I'm afraid I am unable to attach any figures to this post but
> my (corrected) results are:
>
> beta = coth(Pe)-Pe:
> u_SU = [ -9.18e-17, 3.308e-01, 2.57e-01, -5.78e-02,
> -1.37e-01, -1.40e-01, -1.40e-01, -1.40e-01,
> -1.40e-01, -1.40e-01, 1.25e-16]
> u_SUPG - [ -9.52e-17, 3.43e-01, 5.01e-01, 1.27e-01,
> -2.27e-02, -3.23e-02, -3.24e-02, -3.24e-02,
> -3.24e-02, -3.24e-02, 2.88e-17]
>
> beta = 1.0
> u_SU = [ -8.11e-17, 2.92e-01, 1.98e-01, -9.43e-02,
> -1.66e-01, -1.70e-01, -1.70e-01, -1.70e-01,
> -1.68e-01, -1.54e-01, 1.37e-16 ]
> u_SUPG = [ -9.14e-17, 3.29e-01, 4.93e-01, 1.24e-01,
> -2.80e-02, -3.81e-02, -3.83e-02, -3.82e-02,
> -3.79e-02, -3.48e-02, 3.09e-17]
>
> Derek
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>