How to store solution to four coupled equations

Question

I have a system of four coupled nonlinear partial differential equations. My solution function is defined as:

fullspace = dolf.MixedFunctionSpace([space1, space2, space3, space4])
fullsol = dolf.Function(fullspace)

where space1,2,3,4 are created as:

space1 = dolf.FunctionSpace(mesh,"Lagrange",1)

I solve the equations and want to store the solution to a file. I try:

fil = dolf.File("data/testdata.pvd")
fil << fullsol

but I get an error:

*** Error:   Unable to write data to VTK file.
*** Reason:  Don't know how to handle vector function with dimension other than 2 or 3.
*** Where:   This error was encountered inside VTKFile.cpp.

It seems that VTK file format doesn't support solution functions with four components. I also tried using pickle, but that gives me:

TypeError: can't pickle SwigPyObject objects

What other option do I have for storing my solutions in such a way that Fenics can understand them later?

EDIT: I also tried .raw file format. For that, saving into file seems to work, but reading doesn't. I get an error:

*** Error:   Unable to read object from file.
*** Reason:  Cannot read objects of type Function from RAW files.
*** Where:   This error was encountered inside GenericFile.cpp.

Answer 1

Use the HDF5File interface (run as python foo.py read/write):

from dolfin import *
from numpy.random import seed, random
import sys

action = sys.argv[1]
assert action in ["write", "read"]

seed(0)

# Mesh
mesh = UnitSquareMesh(6,6)
V = FunctionSpace(mesh, "CG", 3)
Q = TensorFunctionSpace(mesh, "DG", 0)
C = VectorFunctionSpace(mesh, "CR", 1)
R = FunctionSpace(mesh, "R", 0)

W = V*Q*C*R

u = Function(W)

if action == "write":
    u.vector()[:] = random(W.dim())
    hdf5file = HDF5File("foo.hdf5", 'w')
    hdf5file.write(u, "u")
    del hdf5file
elif action == "read":
    hdf5file = HDF5File("foo.hdf5", 'r')
    hdf5file.read(u, "u")
    del hdf5file

print norm(u)

It's both a lot faster and less storage intensive than xml. Alternatively, you can use .xml.gz instead of .xml.

Comments

Thanks, I will try this.

Answer 2

Okay, I'm going to answer my own question. I found a list of supported file formats and tried them all. It seems that .xml format works. This doesn't feel like the optimal solution though, because .xml file for a 4-component solution for a 2D grid of size 512x128 take 25MB space. Since I will be calculating the solution for a large number of different parameters, the storage space required is huge (tens of gigabytes).