Error in trying to write Mixed Elements

I/O
1

I am solving a nonlinear coupled ODE, with 2 functions u and A and a lagrange multiplier r.
I am using Newton Rhapson; I have to control the step size (relaxation parameter) so that the residue doesn’t blow up. This means more no. of steps till convergence, so i am trying to write the output in a xdmf file which i can use as input for the next set of iterations. When I try to store the output using write_checkpoint, I get the following runtime error.

*** -------------------------------------------------------------------------
*** Error: Unable to find element family.
*** Reason: Element Mixed not yet supported.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0

*** DOLFIN version: 2018.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

My FunctionSpace is created using MixedElements, and I guess the error above remarks that I have to use another route (it’s an implementation issue). It is not clear to me if another route exists. Any method that exists should also be able to read in the written data.

A MWE is presented below:-

from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt
#Create mesh and define function space

Lx=10. #length of domain
mesh = fe.IntervalMesh(1000,0,Lx)
x = SpatialCoordinate(mesh)
VA = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for A
Vu = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for u
R = FiniteElement("Real", mesh.ufl_cell(), 0) #Element for Lagrange multiplier
V = FunctionSpace(mesh, MixedElement(VA, Vu, R)) #Creating functionSpace

#Newton Rhapson input
Aur = interpolate( Expression(("1","0.0", "1.5"), degree=2), V)#SC phase as initial cond.

#<NewtonRhapson Code taken out>

##Save solution in a .xdmf file
Aur_out = XDMFFile('test-Const.xdmf')
Aur_out.write_checkpoint(project(Aur,V), "Aur", 0, XDMFFile.Encoding.HDF5, False)
Aur_out.close()
2

Try replacing your I/O with

Aur_split = Aur.split()

##Save solution in a .xdmf file
for (component_index, component) in enumerate(Aur_split[:2]):
    Aur_out = XDMFFile(f"test-{component_index}.xdmf")
    Aur_out.write_checkpoint(component, "component-{component_index}", 0, XDMFFile.Encoding.HDF5, False)
    Aur_out.close()
with open(f"test-2.txt", "w") as file:
    print(float(Aur_split[2]), file=file)

Furthermore, note that project(Aur,V) in your I/O section is unnecessary, since Aur is already a Function on V.

3

So like i mentioned in my question, you would have to read in what you wrote. The initial condition for the NewtonRhapson Scheme is in Aur. How would you combine the separated files into Aur?

wrt “Project(Aur,V)”.Thanks for pointing that out. I meant to remove that from the MWE. I was part of the code from which I had originally obtained the “xdmf write” statement.

4

Have a look at dolfin.assign

1 Like
5

I think i closed the question too soon. I do think I am facing an issue when I try to read the file again. I decided to plot the files to show the difference. The value at the first node is being changed.

I have attached a MWE below. I should point out that I didn’t use your syntax since I wasn’t familiar with it, so I did it the more tedious way. I couldn’t read in the xdmf file when i used your command,

If you run the code, then the code will first print the initial condition to the Newton Rhapson scheme as input. It then stores it in a .xdmf file, reads it in, and plots it again. The value at the left end has changed. The plots are titled with “before” and “after”. It seems like the Lagrange multiplier is being used to modify the value at the left end.

from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt

#Create mesh and define function space
Lx=10.  #length of domain
mesh = fe.IntervalMesh(1000,0,Lx)
x = SpatialCoordinate(mesh)
VA = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for A
Vu = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for u
R = FiniteElement("Real", mesh.ufl_cell(), 0) #Element for Lagrange multiplier
V = FunctionSpace(mesh, MixedElement(VA, Vu, R)) #Creating MixedFunctionSpace
VFnSp = FunctionSpace(mesh, "Lagrange", 2)#FnSpace for u and A
RFnSp = FunctionSpace(mesh, "Real", 0)#Fn space for lagrange multiplier

#Newton Rhapson input
Aur = interpolate( Expression(("1","0.0", "3.5"), degree=2), V)#SC phase as initial cond.

Aur_split = Aur.split()

(A, u, r) = split(Aur)
plot(u)
plt.title(r"$u(x)-before$",fontsize=26)
plt.show()
plot(A)
plt.title(r"$A(x)e_2-before$",fontsize=26)
plt.show()


Aur_out = XDMFFile(f"test-0.xdmf")
Aur_out.write_checkpoint(Aur_split[0], "component-0", 0, XDMFFile.Encoding.HDF5, False)
Aur_out.close()
Aur_out = XDMFFile(f"test-1.xdmf")
Aur_out.write_checkpoint(Aur_split[1], "component-1", 0, XDMFFile.Encoding.HDF5, False)
Aur_out.close()
with open(f"test-2.txt", "w") as file:
    print(float(Aur_split[2]), file=file)



#Reading in the .xdmf file 
A = Function(VFnSp)
A_in =  XDMFFile("test-0.xdmf")
A_in.read_checkpoint(A,"component-0",0)
u = Function(VFnSp)
u_in =  XDMFFile("test-1.xdmf")
u_in.read_checkpoint(u,"component-1",0)
r = Function(RFnSp)
r = interpolate(Constant("1.5"), RFnSp) #Instead of reading it in, I'm just assigning it.

Aur2 = Function(V)

assign(Aur2, [A,u,r])

plot(u)
plt.title(r"$u(x)-after$",fontsize=26)
plt.show()
plot(A)
plt.title(r"$A(x)e_2-after$",fontsize=26)
plt.show()

Im sure i am making some mistake, but I cant figure out what it is. Any help is highly appreciated.

6

I missed the optional argument when suggesting you to use dolfin.Function.split.

Compare the two cases in the for loop below:

from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np


#Create mesh and define function space
Lx=10.  #length of domain
mesh = fe.IntervalMesh(1000,0,Lx)
x = SpatialCoordinate(mesh)
VA = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for A
Vu = FiniteElement("CG", mesh.ufl_cell(), 2) #Element for u
R = FiniteElement("Real", mesh.ufl_cell(), 0) #Element for Lagrange multiplier
V = FunctionSpace(mesh, MixedElement(VA, Vu, R)) #Creating MixedFunctionSpace
VFnSp = FunctionSpace(mesh, "Lagrange", 2)#FnSpace for u and A
RFnSp = FunctionSpace(mesh, "Real", 0)#Fn space for lagrange multiplier

#Newton Rhapson input
Aur = interpolate( Expression(("1","0.0", "3.5"), degree=2), V)#SC phase as initial cond.

def print_function_info(pre_text, function):
    local = function.vector().get_local()
    print(pre_text, np.min(local), np.max(local))

for deepcopy in (False, True):
    Aur_split = Aur.split(deepcopy=deepcopy)

    print_function_info(f"Component 0, deepcopy={deepcopy}, before", Aur_split[0])
    print_function_info(f"Component 1, deepcopy={deepcopy}, before", Aur_split[1])

    Aur_out = XDMFFile(f"test-0-{deepcopy}.xdmf")
    Aur_out.write_checkpoint(Aur_split[0], "component-0", 0, XDMFFile.Encoding.HDF5, False)
    Aur_out.close()
    Aur_out = XDMFFile(f"test-1-{deepcopy}.xdmf")
    Aur_out.write_checkpoint(Aur_split[1], "component-1", 0, XDMFFile.Encoding.HDF5, False)
    Aur_out.close()

    #Reading in the .xdmf file
    A = Function(VFnSp)
    A_in =  XDMFFile(f"test-0-{deepcopy}.xdmf")
    A_in.read_checkpoint(A,"component-0",0)
    u = Function(VFnSp)
    u_in =  XDMFFile(f"test-1-{deepcopy}.xdmf")
    u_in.read_checkpoint(u,"component-1",0)

    print_function_info(f"Component 1, deepcopy={deepcopy}, after", A)
    print_function_info(f"Component 1, deepcopy={deepcopy}, fater", u)