dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Tuesday 21 October 2008 12:43:27 Heitor Pina wrote:
> Hi all
>
> I am trying to solve (PyDolfin) a convection-diffusion problem with a
> stabilization term (SU-Streamline Upwind). I define a SU class and an
> instance tau.
>
> ...
> class SU(Function):        # SU stabilization
>     def __init__(self, element, mesh, sigma, K):
>         Function.__init__(self, element)
>         self.sigma = sigma
>         self.K = K
>     def eval(self, values, x):
>         mh = MeshSize("triangle",mesh)
>         mhx = mh.eval(x)
>         m = 1.0/3.0     # For linear elements only
>         nv = self.sigma * sqrt(numpy.inner(vel, vel))
>         Re = m*nv*mhx / (2.0*self.K)
>         xi = min(Re, 1.0)
>         values[0] = mhx * xi / (2.0*nv)

Hello!

I am also doing some SUPG stuff in python. 

You should aquire the local MeshSize by 

  mhx = self.cell().diameter()

instead of

  mh = MeshSize("triangle",mesh)
  mhx = mh.eval(x)

The latter will create a function each time your SU function get called. 
self.cell() returns the local cell when the function is used during assembly, 
and its diameter() will return the size of the cell.

What is vel? I suppose it is an array of a constant velocity? If you have a 
non constant velocity you might want to define that in a function and hand 
that to the constructor too, and call it with x inside eval.

When you get it going you might also try using the compile_function interface 
instead of defining a userdefined function in python, it is _a lot_ faster. I 
can help you out with this if you want to try it.

Cheers!

Johan

> tau = SU(scalar, mesh)   # SU stabilization
> a_K = K * dot(grad(v), grad(u))*dx
> a_H = h * v*u*ds
> a_C = sigma * v*dot(vel, grad(u))*dx
> a_S = (tau * dot(vel,grad(v)) * dot(vel,grad(u)))*dx
> a = a_K + a_H + a_C + a_S
> A = assemble(a, mesh)
> ...
>
> The following message appears:
>
> Assembling matrix over cells (finished).
> Plot active, press 'q' to continue.
> Assembling matrix over cells (finished).
> Traceback (most recent call last):
>   File "/home/hpina/Desktop/Solar-Dolfin/conv-dif.py", line 169, in
> <module> A = assemble(a, mesh)
>   File "/usr/lib/python2.5/site-packages/dolfin/assemble.py", line 107, in
> assemble
>     cell_domains, exterior_facet_domains, interior_facet_domains,
> reset_tensor)
>   File "/usr/lib/python2.5/site-packages/dolfin/dolfin.py", line 8084, in
> cpp_assemble
>     return _dolfin.cpp_assemble(*args)
> RuntimeError: *** Error: Function contains no data.
>
>
> When I drop the a_S term the program works fine. I am at a loss so any help
> will be appreciated.
>
> Best regards
>
>   Heitor Pina