Getting error 'Expected only linear combinations of Functions in the same FunctionSpaces' while adding solutions obtained from coupled system

General
1

Hi,
I am trying to solve time dependent coupled scheme and I want to sum the obtained solutions U_1, U_2 on each time step. However, I am getting error ‘Expected only linear combinations of Functions in the same FunctionSpaces’. And I searched the same error message post on the forum but not able to resolve the issue. So can anyone help me? The minimal code is given as:

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
set_log_level(LogLevel.ERROR)

T = 1.0 # final time 
for j in range(1,4):
    nx = 2**(j+1)
    T = 1
    num_steps = nx*nx
    dt = T / (num_steps)
    mesh = UnitSquareMesh(nx, nx)
    P1 = FiniteElement('P', triangle, 1)
    element = MixedElement([P1, P1])
    V = FunctionSpace(mesh, element)
    W = FunctionSpace(mesh,'CG', 1)
    R = Function(W)
    u_D = Expression('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', degree = 3, t=0)
    def boundary(x, on_boundary):
        return on_boundary
    bc1 = DirichletBC(V.sub(0), u_D, boundary)
    bc2 = DirichletBC(V.sub(1), u_D, boundary)
    bc = [bc1, bc2]
    #U = Function(V)
    U_1, U_2 = TrialFunctions(V) 
    u_g = Expression(('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', '0'), degree = 3, t=0)
    U_n = interpolate(u_g, V)
    U_n1, U_n2 = split(U_n)
    p, q = TestFunctions(V)
    e_time = 0
    for i in range(num_steps):
        t = (i+1)*dt
        u_g.t = t
        a = 3/4*U_1*p*dx + 1/4*U_2*p*dx + dt/2*dot(grad(U_1), grad(p))*dx - 9/4*U_1*q*dx + 5/4*U_2*q*dx + dt/2*dot(grad(U_2), grad(q))*dx - dt/2 *(2*pi*pi-1)*U_1*p*dx - dt/2*(2*pi*pi-1)*U_2*q*dx
        L = U_n1*p*dx - U_n2*q*dx 
        U = Function(V)
        solve(a==L, U, bc)
        U_1, U_2 = U.split(deepcopy=True)
        R.assign(U_1+ U_2)
    #errorin space using direct command
    R_in = interpolate(R, W)
    u_ex = interpolate(u_D, W)
    error_L2 = errornorm(u_ex, R_in, 'L2')
    error_H1 = errornorm(u_ex, R_in, 'H1')
    print ('dt=%.4f nx=%.2f error_L2=%8.4E error_H1=%8.4E' % (dt, nx, error_L2, error_H1))
plot(R)
plt.show()

complete error message is R.assign(U_1+ U_2)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 432, in assign
linear_comb = _check_and_contract_linear_comb(rhs, self, multi_index)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 161, in _check_and_contract_linear_comb
_assign_error()
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 34, in _assign_error
raise RuntimeError(“Expected only linear combinations of Functions in the same FunctionSpaces”)
RuntimeError: Expected only linear combinations of Functions in the same FunctionSpaces.

2

Here is a smaller example to illustrate how I would do this as your functions are in P1:

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

set_log_level(LogLevel.ERROR)
nx = 10
mesh = UnitSquareMesh(nx, nx)
P1 = FiniteElement("P", triangle, 1)
element = VectorElement("P", triangle, 1, dim=2)  # MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
W = V.sub(0).collapse()
R = Function(W)
u_D = Expression("exp(-t)*sin(pi*x[0])*sin(pi*x[1])", degree=3, t=0)


def boundary(x, on_boundary):
    return on_boundary


bc1 = DirichletBC(V.sub(0), u_D, boundary)
bc2 = DirichletBC(V.sub(1), u_D, boundary)
bc = [bc1, bc2]
# U = Function(V)
U_1, U_2 = TrialFunctions(V)
u_g = Expression(("exp(-t)*sin(pi*x[0])*sin(pi*x[1])", "0"), degree=3, t=0)
U_n = interpolate(u_g, V)
U_n1, U_n2 = split(U_n)
p, q = TestFunctions(V)
e_time = 0
dt = 0.1
i = 0
t = (i + 1) * dt
u_g.t = t
a = (
    3 / 4 * U_1 * p * dx
    + 1 / 4 * U_2 * p * dx
    + dt / 2 * dot(grad(U_1), grad(p)) * dx
    - 9 / 4 * U_1 * q * dx
    + 5 / 4 * U_2 * q * dx
    + dt / 2 * dot(grad(U_2), grad(q)) * dx
    - dt / 2 * (2 * pi * pi - 1) * U_1 * p * dx
    - dt / 2 * (2 * pi * pi - 1) * U_2 * q * dx
)
L = U_n1 * p * dx - U_n2 * q * dx
U = Function(V)
solve(a == L, U, bc)
U_1, U_2 = U.split(deepcopy=True)
U_1_to_vertex = dof_to_vertex_map(U_1.function_space())
U_2_to_vertex = dof_to_vertex_map(U_2.function_space())
vertex_to_R = vertex_to_dof_map(R.function_space())
R.vector().set_local(
    (U_1.vector().get_local()[U_1_to_vertex] + U_2.vector().get_local()[U_2_to_vertex])[
        vertex_to_R
    ]
)


R.rename("R", "R")
U_1.rename("U_1", "U_1")
U_2.rename("U_2", "U_2")

with XDMFFile("functions.xdmf") as xdmf:
    xdmf.write(R, 0.0)
    xdmf.write(U_1, 0.0)
    xdmf.write(U_2, 0.0)
3

Thanks Dokken for quick reply! It’s working now.
But I am not able to understand this syntax properly so it would be helpful to me if you can share any link regarding this.
And I am not getting the required error as I am using dg(1)cg(1) scheme(dg(1)-piece wise linear polynomials), and expected order of convergence for L2(0,T;H1(\Omega))(e_time1) and max_error are 2, but I am getting wrong order of convergence i.e., 1. However, for space variable it is working fine. I am solving u_t -\frac{1}{2\pi^2} \Delta u = 0, using the scheme

scheme
scheme1533×779 79.2 KB

with A= -\frac{1}{2\pi^2} \Delta , \quad k_n = dt, \quad f =0

The minimal code is given as

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

set_log_level(LogLevel.ERROR)
T = 1.0 # final time
e_time = 0 
for j in range(1,5):
    nx = 2**(j+1)
    T = 1
    num_steps = nx*nx
    dt = 4*T / (num_steps)
    mesh = UnitSquareMesh(nx, nx)
    P1 = FiniteElement('P', triangle, 1)
    element = VectorElement("P", triangle, 1, dim=2)
    V = FunctionSpace(mesh, element)
    W = V.sub(0).collapse()
    R = Function(W)
    t = 0
    u_D = Expression('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', degree = 3, t=t)

    def boundary(x, on_boundary):
        return on_boundary
    bc1 = DirichletBC(V.sub(0), u_D, boundary)
    bc2 = DirichletBC(V.sub(1), u_D, boundary)
    bc = [bc1, bc2]
    
    U_1, U_2 = TrialFunctions(V) 
    u_g = Expression(('exp(-t)*sin(pi*x[0])*sin(pi*x[1])', '0'), degree = 3, t=t)
    U_n = interpolate(u_g, V)
    U_n1, U_n2 = split(U_n)
    p, q = TestFunctions(V)
    e_time = 0
    for i in range(num_steps):
        t = (i+1)*dt
        
        a = (U_1 + U_2) * p * dx + dt/(2*pi*pi) * dot(grad(U_1), grad(p)) * dx + dt / (4*pi*pi) * dot(grad(U_2), grad(p)) * dx + 1 / 2 * U_2 * q * dx + 1 / (4*pi*pi) * dt * dot(grad(U_1), grad(q)) * dx + dt / (6*pi*pi) * dot(grad(U_2), grad(q)) * dx 
        L =  U_n1 * p * dx + U_n2*p*dx 
        U = Function(V)
        solve(a == L, U, bc)
    U_1, U_2 = U.split(deepcopy=True)
       
    U_1_to_vertex = dof_to_vertex_map(U_1.function_space())
    U_2_to_vertex = dof_to_vertex_map(U_2.function_space())
    vertex_to_R = vertex_to_dof_map(R.function_space())
    R.vector().set_local(
       (U_1.vector().get_local()[U_1_to_vertex] + U_2.vector().get_local()[U_2_to_vertex])[
           vertex_to_R
        ]
    )
    R.rename("R", "R")
    U_1.rename("U_1", "U_1")
    U_2.rename("U_2", "U_2")
    vertex_values_u_D = u_D.compute_vertex_values(mesh)
    vertex_values_R = R.compute_vertex_values(mesh)
    error_max = np.max(np.abs(vertex_values_u_D - vertex_values_R))  
    print(error_max)
    error_L2 = errornorm(u_D, R, 'L2')
    error_H1 = errornorm(u_D, R, 'H1')
    e_time = (e_time + dt*error_H1*error_H1)
    e_time1 = sqrt(e_time)
    print ('dt=%.4f nx=%.2f error_L2=%8.4E error_H1=%8.4E e_time1=%8.4E' % (dt, nx, error_L2, error_H1, e_time1 ))
plot(R)
plt.show()

Thank you in advance.

4

Your scheme is not using DG-1 anywhere, so Im not sure where that pops up in your code.
What I am doing in the code above is based on using first order continuous Lagrange elements in all components, as you did.

When doing this, I map each sub space to its vertex, so that we can map the sum back to the latter space using the vertex_to_dof_map.
FIrst I would start by saving U_1, U_2 and R to file, and inspect them in Paraview. Does it look like R is the sum of the to others? If so, something else in your code is not.
If they are not the same, please make a minimal reproducible example (similar to the one I made above).

5

Above defined heat equation can we written by using Space-time finite element method as

\begin{equation} \begin{aligned} (U_n, X) & + A(U, X) \, dt + \left(U^{+}_{n-1}, X^{+}\right)_{J_n} \\ & + (f, X) \, dt, \quad \forall X \in S_k^n, \quad 1 \leq n \leq N, \\ U_0 & = v. \end{aligned} \end{equation}

with the substitution of U(t) = U_1^{n} + \frac{(t-t_n)}{dt} U_2^n on J_n, for q = 2, i.e piecewise linear functions of t and obtained the system of equations in U_1, U_2 snipped attached in previous reply. ( applying the algorithm from the book entitled ’ Galerkin Finite Element Methods for Parabolic Problems’ on page no 206, https://www.math.hkust.edu.hk/~mamu/courses/531/parabolic.pdf ). That is why I am expecting 2 order of convergence.
As I did not find any code regarding space-time finite element method on Fenics, so my code may miss something but according to the plots of R and exact solutions it seems fine.
For my problem, the minimal code is the same as I provided in the previous response. It is a bit messy, and it is not feasible to elaborate on the entire algorithm here. The orders of convergence in L2 and H1 norms are accurate, but in the time error, I am observing first-order convergence. According to theory, it should be 2. I hope you understand the issue. Thanks!

6

Okay, I got the fact why I am getting only 1 order of convergence in time as I am calculating L^2(H^1), the error is dominated by H^1, but now I have different issue that
dof_to_vertex_map() Is just a utility function for simple ‘P1’ function space. Then for P2 elements how can I get R= U_1 + U_2 in this minimal code .

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

set_log_level(LogLevel.ERROR)
nx = 10
mesh = UnitSquareMesh(nx, nx)
P2 = FiniteElement("Lagrange", triangle, 2)
#element = VectorElement("P", triangle, 1, dim=2)  # MixedElement([P1, P1])
element = MixedElement([P2, P2]) 
V = FunctionSpace(mesh, element)
W = V.sub(0).collapse()
R = Function(W)
u_D = Expression("exp(-t)*sin(pi*x[0])*sin(pi*x[1])", degree=3, t=0)


def boundary(x, on_boundary):
    return on_boundary


bc1 = DirichletBC(V.sub(0), u_D, boundary)
bc2 = DirichletBC(V.sub(1), u_D, boundary)
bc = [bc1, bc2]
# U = Function(V)
U_1, U_2 = TrialFunctions(V)
u_g = Expression(("exp(-t)*sin(pi*x[0])*sin(pi*x[1])", "0"), degree=3, t=0)
U_n = interpolate(u_g, V)
U_n1, U_n2 = split(U_n)
p, q = TestFunctions(V)
e_time = 0
dt = 0.1
i = 0
t = (i + 1) * dt
u_g.t = t
a = (
    3 / 4 * U_1 * p * dx
    + 1 / 4 * U_2 * p * dx
    + dt / 2 * dot(grad(U_1), grad(p)) * dx
    - 9 / 4 * U_1 * q * dx
    + 5 / 4 * U_2 * q * dx
    + dt / 2 * dot(grad(U_2), grad(q)) * dx
    - dt / 2 * (2 * pi * pi - 1) * U_1 * p * dx
    - dt / 2 * (2 * pi * pi - 1) * U_2 * q * dx
)
L = U_n1 * p * dx - U_n2 * q * dx
U = Function(V)
solve(a == L, U, bc)
U_1, U_2 = U.split(deepcopy=True)
U_1_to_vertex = dof_to_vertex_map(U_1.function_space())
U_2_to_vertex = dof_to_vertex_map(U_2.function_space())
vertex_to_R = vertex_to_dof_map(R.function_space())
R.vector().set_local(
    (U_1.vector().get_local()[U_1_to_vertex] + U_2.vector().get_local()[U_2_to_vertex])[
        vertex_to_R
    ]
)


R.rename("R", "R")
U_1.rename("U_1", "U_1")
U_2.rename("U_2", "U_2")

with XDMFFile("functions.xdmf") as xdmf:
    xdmf.write(R, 0.0)
    xdmf.write(U_1, 0.0)
    xdmf.write(U_2, 0.0)

g
U_1_to_vertex = dof_to_vertex_map(U_1.function_space())
RuntimeError:
*** Error: Unable to tabulate dof to vertex map.
*** Reason: Can only tabulate dofs on vertices.
*** Where: This error was encountered inside DofMap.cpp.

7

As R is computed outside the temporal loop, I would argue that the simplest solution is to compute R = project(U_1+U_2, W), as it will be exact as U_1 and U_2 is in W1.

1 Like
8

If I use polynomial degree 2, then while using DG(1)CG(2), i.e., DG time stepping with standard Galerkin method in space, I am not getting the correct order of convergence in L2 and H1 norm. One thing here to mention is that while defining the solution U(t) = \tilde{U}_0 + \tilde{U}_1 (t-t_{n-1})/ \Delta t on J_n:= ( t_{n-1}, t_n], so I think the issue is in defining initial value because t_{n-1} is not in J_n, the minimal heat equation code for dg time stepping method for approximating polynomial is given below and here as I decrease the mesh size, error is increasing…

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym

t = sym.Symbol('t')
x, y = sym.symbols('x[0] x[1]')
u = t*x*(1-x)*y*(1-y)
#u = sym.exp(-t)*sym.sin(pi*x)*sym.sin(pi*y)
f = sym.diff(u,t)- sym.diff(sym.diff(u, x), x) - sym.diff(sym.diff(u, y), y) 
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

nx = 2**3
num_steps = nx
T = 1.0 
dt = T / num_steps # time step size
mesh = UnitSquareMesh(nx, nx)
P2 = FiniteElement("Lagrange", triangle, 2)
element = MixedElement([P2, P2])
V = FunctionSpace(mesh, element)
W = V.sub(0).collapse()
R = Function(W)

u_D = Expression(u_code, degree = 4, t=0)
def boundary(x, on_boundary):
    return on_boundary
bc1 = DirichletBC(V.sub(0), u_D, boundary)
bc2 = DirichletBC(V.sub(1), u_D, boundary)
bc = [bc1, bc2]

(U_1, U_2) = TrialFunctions(V)
(p, q) = TestFunctions(V)

u_g = Expression((u_code ,'0'), degree = 3, t=0)
U_n = interpolate(u_g, V)
#U_n = Function(V)
U_n1, U_n2 = split(U_n)

UF = [U_n1]
US = [U_n2]
e_time = 0
for i in range(num_steps):
    t = (i+1)*dt
    (U_1, U_2) = TrialFunctions(V)
    f = Expression(f_code, degree =4, t = (i+1)*dt)
    UF.append(U_1)
    US.append(U_2)
    u_D = Expression(u_code, degree = 4, t=(i+1)*dt)

    f4 = Expression(f_code, degree=4, t=(i+1)*dt) # final
    f0 = Expression(f_code, degree=4, t=i*dt)
    f1 = Expression(f_code, degree=4, t=(i+ 0.25)*dt)     # initial
    f2 = Expression(f_code, degree=4, t=(i+0.50)*dt)
    f3 = Expression(f_code, degree=4, t=(i+0.75)*dt)
    f = dt/ 12 *( f0 + 4*(f1+f3) + 2*f2 + f4) # using simpson's 1/3rd rule 
    g0 = -dt*f0 # for right hand side of second equation
    g1 = -0.75*f1
    g2 = -0.50*f2
    g3 = -0.25*f3
    g4 = 0
    g = dt/ 12 *( g0 + 4*(g1+g3) + 2*g2 + g4) # using simpson's 1/3rd rule
    # in append sense
    #F = (UF[i+1] + US[i+1]) * p * dx + dt * dot(grad(UF[i+1]), grad(p)) * dx + dt / 2 * dot(grad(US[i+1]), grad(p)) * dx + 1 / 2 * US[i+1] * q * dx + 1 / 2 * dt * dot(grad(UF[i+1]), grad(q)) * dx + dt / 3 * dot(grad(US[i+1]), grad(q)) * dx - UF[i]*p*dx - US[i]*p*dx  - f * p * dx - 1 / dt * g * q * dx 

    # in simple U_1 U_2 sene
    F = (U_1 + U_2) * p * dx + dt * dot(grad(U_1), grad(p)) * dx + dt / 2 * dot(grad(U_2), grad(p)) * dx + 1 / 2 * U_2 * q * dx + 1 / 2 * dt * dot(grad(U_1), grad(q)) * dx + dt / 3 * dot(grad(U_2), grad(q)) * dx - UF[i]*p*dx - US[i]*p*dx  - f * p * dx - 1 / dt * g * q * dx 
    # Update current time
    a, L = lhs(F), rhs(F)
    U = Function(V)
    solve(a == L, U, bc)
    U_1, U_2 = U.split()
    R = project(U_1 + U_2, W)
    
    #Compute error at vertices
    u_e = interpolate(u_D, W)
    error_L2 = errornorm(u_D, R, 'L2')
    error_H1 = errornorm(u_D, R, 'H1')
    e_time = (e_time + dt*error_H1*error_H1)
    e_time1 = sqrt(e_time)
    error_max = np.abs(u_e.vector().get_local() - R.vector().get_local()).max() 
print('nx = %.2f: error_max = %.3g, error_L2 = %.3g, error_H1 = %.3g, e_time1= %.3g' % (nx, error_max, error_L2, error_H1,e_time1))

scheme i follow while solving heat equation using dg time stepping method in time

scheme
scheme1533×779 79.2 KB

9

This question is nor related to the original post. Please open a separate issue.

I would Also suggest that you first try to minimize the amount of code, and include output from your code as part of any question.

10

Thank you for mentioning it. I have created a new post.