I am having trouble outputting multiple fields to an XDMF file for my coupled problem. To demonstrate this error, I reproduced it in the Cahn-Hilliard equation demo.
I replace the solution loop of the demo with the following:
# Output file
xdmf_file = XDMFFile(mesh.mpi_comm(), "cahn_hilliard.xdmf")
# Step in time
t = 0.0
T = 20*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
c_out = u.split()[0]
mu_out = u.split()[1]
c_out.rename("c", "tmp")
mu_out.rename("mu", "tmp")
xdmf_file.write(c_out, t)
xdmf_file.write(mu_out, t)
This executes without error, but when I try to open the XDMF file in ParaView, I get the following error (extracted from ParaView's console)
Warning: In /home/buildslave/dashboards/buildbot/paraview-pvbinsdash-linux-shared-release_opengl2_qt4_superbuild/source-paraview/VTK/IO/Xdmf2/vtkXdmfReader.cxx, line 433
vtkXdmfReader (0x4669270): Data type generated (vtkMultiBlockDataSet) does not match data type expected (vtkUnstructuredGrid). Reader may not produce valid data.
Outputting multiple fields works when I don't use mixed elements, for example when I calculate the stresses and output them alongside the displacements. I am not acquainted with the format or how FEniCS outputs to it, so I don't know how to solve this problem.
One hint that I have is dimensionality of the function space. I had asked this question related to the same model, and the two problems might be connected.
Edit: full code
import random
from dolfin import *
# Class representing the intial conditions
class InitialConditions(Expression):
def __init__(self, **kwargs):
random.seed(2 + MPI.rank(mpi_comm_world()))
def eval(self, values, x):
values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
# Model parameters
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
# Create mesh and build function space
mesh = UnitSquareMesh(96, 96)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
xdmf_file = XDMFFile(mesh.mpi_comm(), "cahn_hilliard.xdmf")
# Step in time
t = 0.0
T = 20*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
c_out = u.split()[0]
mu_out = u.split()[1]
c_out.rename("c", "tmp")
mu_out.rename("mu", "tmp")
xdmf_file.write(c_out, t)
xdmf_file.write(mu_out, t)
