How to evaluate variational formulation

1

Hi,

from SUPG I have the residual

res = u - u_n + E_dt*(dot(W, grad(u)) - a1U*div(grad(cpiU*u))).

If I try to save this via

vtkfileR = File('solutionResidual.pvd')
vtkfileR << res

I get:

Traceback (most recent call last):
  File "t-modell_02_SUPG.py", line 355, in <module>
    vtkfileR << res
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/io/__init__.py", line 21, in __lshift__
    self.write(u)
AttributeError: 'Sum' object has no attribute '_cpp_object'

How can I access/save the value of the residual?

2

You need to assemble the residual. Anything written with dot, grad, inner, div, dx etc are objects from the Unified form language (ufl). To obtain numerical values, such objects has to be assembled. To assemble a residual with boundary conditions, see for instance: Residual for Stokes equation not equal to zero - #3 by nate

3

Hmmm … I tried to assemble

res = assemble(u - u_n + E_dt*(dot(W, grad(u)) - a1U*div(grad(cpiU*u))))

but got this

Traceback (most recent call last):
  File "t-modell_02_SUPG.py", line 353, in <module>
    res  = assemble(u - u_n + E_dt*(dot(W, grad(u)) - a1U*div(grad(cpiU*u))))
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 64, in _create_dolfin_form
    raise TypeError("Invalid form type %s" % (type(form),))
TypeError: Invalid form type <class 'ufl.algebra.Sum'>
4

You need to multiply the residual by an integration measure; i.e. dx (if you want the residual over the cells), ds (exterior facets) etc.

5

I did this too:

Traceback (most recent call last):
  File "t-modell_02_SUPG.py", line 353, in <module>
    res  = assemble(u*dx - u_n*dx + E_dt*(dot(W, grad(u))*dx - a1U*div(grad(cpiU*u)))*dx)
  File "/usr/lib/python3/dist-packages/ufl/form.py", line 302, in __sub__
    return self + (-other)
TypeError: unsupported operand type(s) for +: 'Form' and 'Product'
6

Be careful with your brackets. You are multiplying the second to last term by dx twice.

7

Oh, now the error message becomes more unspecific:

Traceback (most recent call last):
  File "t-modell_02_SUPG.py", line 353, in <module>
    res  = assemble(u*dx - u_n*dx + E_dt*dot(W, grad(u))*dx - a1U*div(grad(cpiU*u))*dx)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
    return Form(form,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 85, in mpi_jit
    raise RuntimeError(error_msg)
RuntimeError: Compilation failed on root node.

The variational formulation one line later ist accepted (not in assemble but in solve) without any problem:

F    =   u*v*dx \
       - u_n*v*dx \
       + E_dt*dot(W,grad(u))*v*dx \
       + E_dt*a1U*dot(grad(cpiU*u),grad(v))*dx \
       + sum(integrals_R) + sum(integrals_N)
8

Please supply a minimal working code example. That means a complete code that with as few lines as possible can reproduce the error. A Good example should be no longer tha 25 lines, and have no boundary conditions or other parameters.

9 
from fenics import *
import numpy as np

p0   = Point(0.0,0.0)
p1   = Point(1.0,1.0)
mesh = RectangleMesh(p0,p1,10,10,'left')

Q  = FunctionSpace(mesh, 'CG', 1)
V  = VectorFunctionSpace(mesh, 'CG', 1)

u = TrialFunction(Q)
v = TestFunction(Q)

U = Function(Q)
        
u_n     = interpolate(Constant(1.0),Q)  

W = Function(V)
w = np.array(W.vector()).reshape((-1,2))
w[:,] = np.array([0.0,0.5])
W.vector()[:] = w.reshape((-1))

lambdaU = Function(Q)
rhoU    = interpolate(Constant(1),Q) 
cpiU    = Function(Q)
a1U     = lambdaU/rhoU 

E_dt  = Expression('dt',degree=0,dt=1.0)

res  = assemble(u*dx - u_n*dx + E_dt*dot(W, grad(u))*dx - a1U*div(grad(cpiU*u))*dx)

This is the condensed code. Running without mpirun (before it was with mpirun) returns:

...
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_1',).
10

u is a trial function, and has to be multiplied by a test function (and forms a matrix). The second argument u_n is a function, multiplied by dx will be assembled as a scalar value.

Please make sure that you are assembling something that makes sense. See for instance the stokes residual question above for how one forms a residual.

11

The initial code came basically from the advection SUPG example:

# Residual
r = u - u0 + dt*(dot(velocity, grad(u_mid)) - c*div(grad(u_mid)) - f)

# Galerkin variational problem
F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx \
                      + c*dot(grad(v), grad(u_mid))*dx)

# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx

# Create bilinear and linear forms
a = lhs(F)
L = rhs(F)

The idea now was to plot how the residual looks like.

How would I plot (save for ParaView) the residual over the domain?

12

Note that in r you have a trial function (u). This has to be replaced by your solution to the problem, using
For instance ufl.replace. See Derivative function has no attribute 'subs' - #2 by dokken on how to use replace. If there are no test or trial functions in the residual, the assembly (r*dx) simply becomes a scalar value.
If you want to look at the spatial variation of the residual, you should project it into a suitable function space and save the result to pvd or xdmf

13

Ok, works - thank you.