Cannot copy function values with mixed finite elements

General
1

Hello,
I have an issue with copying the values of a function when using mixed finite elements.

I define my function spaced as

O = VectorFunctionSpace(mesh, 'P', 2, dim=2)
Q = FunctionSpace(mesh, 'P', 1)

and the mixed space as

P_v = VectorElement('P', triangle, 2)
P_w = FiniteElement('P', triangle, 1)
element = MixedElement([P_v, P_w])
V_vw = FunctionSpace(mesh, element)

I define the following functions

vw_ = Function(V_vw)
v_, w_ = split(vw_)
v_dummy = Function(O)
w_n_1 = Function(Q)
w_n_2 = Function(Q)

I in a time loop, I have the following : I obtain the solution of a mixed variational problem and stored it into vw_ with

solve(A1, vw_.vector(), b1, 'bicgstab', 'hypre_amg')

Then I want to

  • copy the numerical values of the function w_n_1 into w_n_2
  • copy the numerical values of the function w_ into w_n_1

If I do

w_n_2.assign(w_n_1)
v_temp, w_n_1 = vw_.split()

This works in the first step of the time loop, but in the second step I get the error message

Traceback (most recent call last):
  File "navier_stokes_membrane.py", line 292, in <module>
    w_n_2.assign(w_n_1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 408, in assign
    self._cpp_object._assign(rhs._cpp_object)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to initialize vector of degrees of freedom for function.
*** Reason:  Cannot re-initialize a non-empty vector. Consider creating a new function.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

What am I doing wrong ?

Thank you :grinning:

2

Please make a minimal reproducible code, i.e. make sure that the code is runnable, but remove all stuff relating to solving a PDE, i.e. defining a UnitSquareMesh, the spaces and functions mentioned above + the temporal loop where you split and assign should be sufficient to have a self contained example.

3

Here it is:

import math
from fenics import *
from mshr import *
import numpy as np
import meshio
import ufl as ufl


# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
domain = channel
mesh = generate_mesh(domain, 64)

O = VectorFunctionSpace(mesh, 'P', 2, dim=2)
Q = FunctionSpace(mesh, 'P', 1)

P_v = VectorElement('P', triangle, 2)
P_w = FiniteElement('P', triangle, 1)
element = MixedElement([P_v, P_w])
V_vw = FunctionSpace(mesh, element)


v_n = Function(O)

w_n = Function(Q)
w_n_1 = Function(Q)
w_n_2 = Function(Q)

vw_ = Function(V_vw)
v_, w_ = split(vw_)

vw = TrialFunction(V_vw)
v, w = split(vw)

v_temp = Function(O)


T = 1.0
num_steps = 2
dt = T / num_steps  


print("Starting time iteration ...", flush=True)
t = 0
for n in range(num_steps):

    t += dt
    
    #  this line causes the issue
    w_n_2.assign(w_n_1)
    
    v_temp, w_n_1 = vw_.split()    
    
    print("\t%.2f %%" % (100.0 * (t / T)), flush=True)

print("... done.", flush=True)

when I run it I get

$ python3 code.py 
Starting time iteration ...
	50.00 %
Traceback (most recent call last):
  File "code.py", line 103, in <module>
    w_n_2.assign(w_n_1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 408, in assign
    self._cpp_object._assign(rhs._cpp_object)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to initialize vector of degrees of freedom for function.
*** Reason:  Cannot re-initialize a non-empty vector. Consider creating a new function.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
4

Consider the following working example:

from fenics import *

mesh = UnitSquareMesh(10, 10)

P_v = VectorElement('P', triangle, 2)
P_w = FiniteElement('P', triangle, 1)
element = MixedElement([P_v, P_w])
V_vw = FunctionSpace(mesh, element)
O = V_vw.sub(0).collapse()
Q = V_vw.sub(1).collapse()

v_n = Function(O)

w_n = Function(Q)
w_n_1 = Function(Q)
w_n_2 = Function(Q)

vw_ = Function(V_vw)
v_, w_ = split(vw_)

vw = TrialFunction(V_vw)
v, w = split(vw)

v_temp = Function(O)


T = 1.0
num_steps = 2
dt = T / num_steps  


print("Starting time iteration ...", flush=True)
t = 0
for n in range(num_steps):

    t += dt
    
    #  this line causes the issue
    w_n_2.assign(w_n_1)
    
    v_temp, w_n_1 = vw_.split(deepcopy=True)    
    
    print("\t%.2f %%" % (100.0 * (t / T)), flush=True)

print("... done.", flush=True)
5

I see, may you please explain to me why it crashed before ?

6

Without the deepcopy, you are just making a reference to the sub function, but the underlying array is from the mixed space.

Please also note that I simplified the construction of sub spaces O and Q to ensure consistency

7

Oh I see, thank you.