Importing a solution from file in InitialConditions class

I/O
1

Hi everyone,
I am solving a time dependent problem, and I am building the InitialConditions class. My problem is a bit complicated so I want to use the steady state solution that I computed previously as an initial condition. Given the structure of my problem, I have to evaluate the initial conditions inside the InitialConditions class, but I cannot find the right way to import the solution there. I am using a mixed finite element formulation.

For reference, my spaces are

V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
S = Finiteelement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q, S]))

Here is an example of what I am trying to do (where w0.h5) contains my steady state solution.

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def eval(self, values, x):
        fFile = HDF5File(MPI.comm_world, "w0.h5", "r")
        fFile.read(values, "/f")
        fFile.close()
    def value_shape(self):
        return (4,)

which gives the error ‘numpy.ndarray’ object has no attribute ‘_cpp_object’.

I have also tried to do this

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.uin = 2
        self.Qfun = 2.3710
        self.x2 = 0.9
        self.x1 = 0.7

    def eval(self, values, x):
        w0 = Function(W)
        fFile = HDF5File(MPI.comm_world, "w0.h5", "r")
        fFile.read(values, "/f")
        fFile.close()
        u0, p0, T0 = split(w0)
        values[0] = u0[0]
        values[1] = u0[1]
        values[2] = p0
        values[3] = T0
    def value_shape(self):
        return (4,)

which then gives the error " values[0] = u0[0] , setting an array element with a sequence"

I can just use fFile outside of the class, when I am starting the time-stepping. This works only for one time-step, and then my solution stops getting updated properly (like I mentioned here), so I think I need to import the file somehow inside the class.

Does anyone know how to do it?

2

u0[0], u0[1], p0, and T0 are all sequence-like objects, because they are scalar functions defined on your entire mesh (as opposed to the value of your function at a single point in the mesh). I believe what you need to do is evaluate the functions at the point x, e.g. u0[0](x), u0[1](x):

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.uin = 2
        self.Qfun = 2.3710
        self.x2 = 0.9
        self.x1 = 0.7
        self.w0 = Function(W)
        fFile = HDF5File(MPI.comm_world, "w0.h5", "r")
        fFile.read(self.w0, "/f")  # I HAVE REMOVED values AND
                                   # REPLACED WITH self.w0
        fFile.close()
        self.u0, self.p0, self.T0 = split(self.w0)

    def eval(self, values, x):
        values[0] = self.u0[0](x)
        values[1] = self.u0[1](x)
        values[2] = self.p0(x)
        values[3] = self.T0(x)
    def value_shape(self):
        return (4,)

Note that I have edited the line with fFile.read in your code, since I think this is what you intended.

Furthermore, your code appears to be quite inefficient, because you read the file "w0.h5" every time the eval function is called, and it will be called repeatedly for many points in your mesh. Therefore, I have moved the fFile.read function into the __init__ method and made w0 a class variable self.w0.

2 Likes
3

(I know we’re only supposed to reply if we have any additional details relevant to the solution, but man I was really struggling, and your explanation was clear, informative, and 100% helpful. Thank you so much. I will read more about Python in general, as I think there are way too many things I haven’t grasped, and I can’t progress in FEniCS without it. )

1 Like