Error creating a vizualizaton .vtk of a UFL function

I/O
1

Hi,
I have reworked the Poisson solver from the FEniCS tutorial, to solve an optimal control problem constrained by a Poisson equation. I am having trouble creating a .vtk file for my function u, which describes the control of the problem. I’ve included the entirety of the code below, since the function u is involved in several steps. I get the following error:

RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_530/1586109424.py in <module>
     77 vtkfile = File('poisson/optimalcontrol.pvd')
---> 78 vtkfile << u
     79 vtkfile = File('CondGDplots/optimalstate.pvd')
     80 vtkfile << y
~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/io/__init__.py in __lshift__(self, u)
     19             self.write(u[0], u[1])
     20     else:
---> 21         self.write(u)
     22 
     23 
RuntimeError: Unable to cast Python instance to C++ type (compile in debug mode for details)
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(16, 16)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Constant(0)
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Init problem data
x0, x1 = MeshCoordinates(mesh)
yd = conditional(x0>0.5,1,-1)
beta = Constant(10)
lmbd = 0.01

# Initialize u
u = Constant(0)

for i in range(10):
    #Solve state y_n
    y = TrialFunction(V)
    v = TestFunction(V)
    f = beta*u
    a = dot(grad(y), grad(v))*dx
    L = f*v*dx
    # Compute solution
    y = Function(V)
    solve(a == L, y, bc)
    #vtkfile = File('CondGDplots/state.pvd')
    #vtkfile << y

    #Solve adjoint state p_n
    p = TrialFunction(V)
    v = TestFunction(V)
    f = y-yd
    a = dot(grad(p), grad(v))*dx
    F = f*v*dx
    # Compute solution
    p = Function(V)
    solve(a == F, p, bc)
    #vtkfile = File('CondGDplots/adjointstate.pvd')
    #vtkfile << p
    #break

    #Choose step direction
    df = lmbd*u+beta*p
    w = conditional(df >= 0 ,10,-10)
    

    #Calculate step state y(v)
    yw = TrialFunction(V)
    v = TestFunction(V)
    f = beta*w
    a = dot(grad(yw), grad(v))*dx
    L = f*v*dx
    # Compute solution
    yw = Function(V)
    solve(a == L, yw, bc)

    #Choose step length
    dyd = y-yd
    dys = yw-y
    dus = w-u
    g1 = assemble((dyd*dys+lmbd*u*dus)*dx)
    g2 = assemble(0.5*(dys*dys+lmbd*dus*dus)*dx)
    gminarg = -0.5*g1/g2
    alpha = min(max(0,gminarg),1)
    print(alpha)
    if alpha < 1e-8:
        break
    u = u+alpha*dus

vtkfile = File('poisson/optimalcontrol.pvd')
vtkfile << u
vtkfile = File('CondGDplots/optimalstate.pvd')
vtkfile << y

Thank you in advance!

2

Please note that you need to project u and y into a suitable function space for visualization.
Also note that to visualize functions from other spaces than “CG” 1, it is recommended to use XDMFFile.write_checkpoint, as it makes it possible to visualize arbitrary order CG and DG spaces.

1 Like
4

Thank you for the help! In the initial tests, I did not need to project y before visualization. Which operations did I apply to u that hindered this?

5

You can only save functions to file.
When you add values to your

You change the type of u.

6

I appreciate you bear with me. Once again, thanks!