How to construct functions of calculated solutions?

1

Hey guys,

so I wrote a script calculating two components of a timedependent PDE system m and c. What I would like to do is whenever I write the solution at time t to a file that I calculate a function f(m(x,t),c(x,t)) that I also want to access through a pvd file in the end.

My first attempt was including the the function into my PDE system but that runs into trouble because the boundary conditions aren’t fulfilled.

In particular: f1 = scaling*c*(m+1)**2/(2*m)+c-1 where scaling = 1 at the moment.

The important snippet of the code:

vtkfile_u1 << (u1_,t)
vtkfile_u2 << (u2_,t)


for n in range(steps):
    
    t+=dt
    
    solve(F==0, u)
	
    u_n.assign(u)
    
    if n%stampspace==0:
        u1_, u2_ = u_n.split()
        vtkfile_u1 << (u1_,t)
        vtkfile_u2 << (u2_,t)

I would like to include the computation in the if statement, this way it wouldn’t be too bad if it were costly.

Ideas?

Thank you!

Sarah

2

In this case, f is relatively simple, and the easiest approach might be to just use a Calculator filter in Paraview. If f is more complicated, or needed elsewhere in the FEniCS script, you might consider a few options, as demonstrated below:

from dolfin import *
mesh = UnitIntervalMesh(10)
el = FiniteElement("CG",mesh.ufl_cell(),1)
el_x_el = MixedElement(2*[el,])
V = FunctionSpace(mesh,el_x_el)
u = project(Constant((1,2)),V)
m,c = u.split()
scaling = 1.0

# Output using Expression.  This directly interpolates $f$ into a CG
# FunctionSpace of degree 1.
f_expr = Expression("scaling*c*pow(m+1.0,2)/(2*m)+c-1.0",
                    scaling=scaling,m=m,c=c,degree=1)
V_out = FunctionSpace(mesh,"CG",1)
f_out = Function(V_out)
f_out.interpolate(f_expr)
File("f_expr.pvd") << f_out

# Output using $L^2$ projection.  This is more flexible (if $f$ involves,
# e.g., derivatives, conditionals, etc.), but requires solving a linear system.
f = scaling*c*pow(m+1.0,2)/(2*m)+c-1.0
File("f_proj.pvd") << project(f, V_out)

# Lumped-mass projection.  This remains flexible and avoids solving a linear
# system, but is only an approximation of the $L^2$ projection, and requires
# an extra user-defined function.  This is also more robust than full $L^2$
# projection if $f$ is discontinuous (e.g., from a UFL conditional).
def lumped_project(f):
    v = TestFunction(V_out)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V_out)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u
File("f_lump.pvd") << lumped_project(f)
3

Thanks, it took me some while to figure out how to save the expressions and I also checked out the paraview filter, which is awesome! Really helpful here :slight_smile: