ValueError: vector::reserve for writing a xdmf file

Errors
1

Hey guys,
first of all: this is my very first Topic so if there are any mistakes/ comments from your side to match this forums rules, please let me know:
I am currently working on my bachelors thesis which involves solving the NSE with fenics. Since I am not that experiences with numeric methods i am trying to get to know this software better by solving random (mostly easy) equations and working my way up to NSE.
Now i got the following code which gives me :

ERROR

Traceback (most recent call last):
File “RAW.py”, line 107, in
xdmf_file_u.write(u, t_start)
ValueError: vector::reserve

CODE

from future import print_function
import fenics
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mshr import *
from fenics import XDMFFile
from fenics import derivative

Umgebung definieren

l1 = 1
l2 = 1
domain = Rectangle(Point(0, 0), Point(2*l1, l2))

Mesh definieren

mesh = generate_mesh(domain, 128)

Element definieren

P1 = FiniteElement(“P”, triangle, 1)

Funktionsraum erstellen

V = FunctionSpace(mesh, P1)

Zeitdiskretisierung

f = Constant(5)
t_start = 0.0
TZeit = 1/f # Anzahl der Oszillationen
num_steps = 25
dt = TZeit/ (num_steps)

Funktionen definieren

u = TrialFunction(V)
u_n = Function(V)
w = TestFunction(V)

Anfangsbedingungen + Projektion der Anfangslösung

u_init = Constant(19.0)
u_n = project(u_init, V) # source term

Randbedingungen (Schallwelle)

amplitudeu = Constant(1)
Schallwelle = Expression((‘u_init + amplitudeu * sin(2 * pi * f * t_start)’), amplitudeu=amplitudeu, f=f, t_start=t_start, u_init=u_init, degree=1, domain=mesh)

tol = 1E-14

def boundary_Dl(x, on_boundary):
if on_boundary:
if near(x[0], 0, tol):
return True
else:
return False
else:
return False

bcl = DirichletBC(V, Schallwelle, boundary_Dl)

Gleichung in Variationalform

F = inner(nabla_grad(u) - nabla_grad(u_n), nabla_grad(w)) * dx

a = lhs(F)
L = rhs(F)

pvd files erstellen

xdmf_file_u = XDMFFile(‘RAWtest/u.xdmf’)

Lösungsprozess

for n in range(num_steps):

u = Function(V)
solve(a == L , u, bcl)


# Speicherprzozess
xdmf_file_u.write(u, t_start)


# nächster Schritt

t_start += dt
Schallwelle.t_start = t_start
bcl = DirichletBC(V, Schallwelle, boundary_Dl)
u_n.assign(u)

I think the problem is related to solving the eq. F for u and putting those results in a xdmf before assigning u to u_n for the next time step. There are examples for creating another function for example u_ which acts as a place holder just for saving purposes but so far I have not realised how to correctly impliment that and if that even is a solution for my problem. Thanks for everyone reading this (:

2

Please use 3x` encapsulation of all code, as it is currently not well formatted, i.e.

```python
# add code here
```

Secondly, could you use a RectangleMesh from dolfin instead of the one from mshr? Does this change the behavior of the error?

Thirdly, add a full traceback of the error message

3 
from __future__ import print_function
import fenics
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mshr import *
from fenics import XDMFFile
from fenics import derivative

l1 = 1
l2 = 1
nx = ny = 30
mesh = RectangleMesh(Point(0.0), Point(2*l1, l2), nx, ny)

P1 = FiniteElement("P", triangle, 1)

V = FunctionSpace(mesh, P1)

f = Constant(5)
t_start = 0.0
TZeit = 1/f                                                                           # Anzahl der Oszillationen
num_steps = 25
dt = TZeit/ (num_steps)

u = TrialFunction(V)
u_n = Function(V)
w = TestFunction(V)

u_init = Constant(19.0)
u_n = project(u_init, V)       

amplitudeu = Constant(1)
Schallwelle = Expression(('u_init + amplitudeu * sin(2 * pi * f * t_start)'), amplitudeu=amplitudeu,  f=f, t_start=t_start, u_init=u_init, degree=1, domain=mesh)

tol = 1E-14

def boundary_Dl(x, on_boundary):
    if on_boundary:
        if near(x[0], 0, tol):
            return True
        else:
            return False
    else:
        return False

bcl = DirichletBC(V, Schallwelle, boundary_Dl)

F = inner(nabla_grad(u) - nabla_grad(u_n), nabla_grad(w)) * dx

a = lhs(F)
L = rhs(F)

u_solution = Function(V)
u = Function(V)

xdmf_file_u = XDMFFile('RAWtest/u.xdmf')

for n in range(num_steps):

    # Solve for u
    solve(a == L, u, bcl)

    # Assign values to the storage function

    u_solution.assign(u)

    # Write to XDMF file

    xdmf_file_u.write(u_solution, t_start) -------------------------- ValueError: vector::reserve

    # Update for the next time step

    t_start += dt
    Schallwelle.t_start = t_start
    bcl = DirichletBC(V, Schallwelle, boundary_Dl)
    u_n.assign(u)

I changed the mesh but I get the same error. The full error message is: Traceback (most recent call last):
File “RAW.py”, line 106, in
xdmf_file_u.write(u_solution, t_start)
ValueError: vector::reserve
I marked the relevant line in the code

4

Your code works for me if I use

t_start += float(dt)

instead of

t_start += dt

Since dt depends on f = Constant(5), it will result in a ufl type rather than a float. Converting it to float ensures that t_start remains a real number rather than a ufl expression.

5

Hey,
thank you very much. I am now able to save my data in a xdmf file! The results lack physical meaning though so I might need your help again when I am able to show the actual problem.