Splitting mixed function space in dolfinx

dolfinx
1

Hello everybody,

I want to assign the result ‘up’ based on mixed functions to a previous timestep. Knowing from former fenics, I naively wanted to adopt

(_ua, _pa) = up.split(deepcopy=True)

But using

_ua, _pa = up.split()

within dolfinx, I get unexpected behavior. To confirm this, an MWE is

from mpi4py import MPI
from dolfinx.generation import UnitSquareMesh
from dolfinx.cpp.mesh import CellType
from dolfinx import fem, FunctionSpace 
from ufl import TestFunctions, TrialFunctions, VectorElement, FiniteElement, MixedElement, dot, grad, div, dx, lhs, rhs
import numpy as np

mesh = UnitSquareMesh(MPI.COMM_WORLD, 4, 4, CellType.triangle)    

U = VectorElement("Lagrange", mesh.ufl_cell(), 1)
P = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([U, P]))

(ua_n2, pa_n2) = TrialFunctions(W) 
(qu, qp) = TestFunctions(W)
up = fem.Function(W)

ua_n1 = fem.Function(FunctionSpace(mesh, U))
pa_n1 = fem.Function(FunctionSpace(mesh, P))



F = (\
        dot(\
            ua_n2-ua_n1
            +grad(pa_n2) 
            , qu)+
        dot(
            pa_n2-pa_n1
            +div(ua_n2)
            , qp)
         )*dx
         
problem = fem.LinearProblem(lhs(F),  rhs(F) , u=up)
problem.solve()
_ua, _pa = up.split()

print(f'ua_n1 pre: {np.shape(ua_n1.x.array)}')
print(f'pa_n1 pre: {np.shape(pa_n1.x.array)}')

print(f'u_a post: {np.shape(_ua.x.array)}')
print(f'_ua post: {np.shape(_ua.compute_point_values()[:])}')    

ua_n1.x.array[:] = _ua.x.array
ua_n1.x.scatter_forward()

As you can see, scattering _ua to ua_n1 doesn’t work caused by a compatibility mismatch. Does anybody have an idea, how I definitely get the values of _ua in order to update ua_n1? It seems that the missing deepcopy-keyword is the reason for this mismatch.

Many thanks!

2

Consider the minimal code example, which illustrates how to copy data from a sub space to the parent space:

from dolfinx.fem import Function, FunctionSpace
from dolfinx.generation import UnitSquareMesh
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType
from mpi4py import MPI
from ufl import FiniteElement, MixedElement, VectorElement

mesh = UnitSquareMesh(MPI.COMM_WORLD, 4, 4, CellType.triangle)

u_el = VectorElement("Lagrange", mesh.ufl_cell(), 1)
p_el = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([u_el, p_el]))

up = Function(W)

U, U_to_W = W.sub(0).collapse(True)

ua_n1 = Function(U)


def f(x):
    return (x[0], x[1])


ua_n1.interpolate(f)

up.x.array[U_to_W] = ua_n1.x.array

with XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(up.sub(0).collapse())

The key is to base the sub function space on the parent mixed space, using W.sub(i).collapse(True) which will return you the sub space, and the map from the sub space to the parent space.

2 Likes
3

Great, that solves my problem. For those interested in the minimal working example, please consider

from mpi4py import MPI
from dolfinx.generation import UnitSquareMesh
from dolfinx.cpp.mesh import CellType
from dolfinx import fem, FunctionSpace 
from ufl import TestFunctions, TrialFunctions, VectorElement, FiniteElement, MixedElement, dot, grad, div, dx, lhs, rhs
import numpy as np

mesh = UnitSquareMesh(MPI.COMM_WORLD, 1, 1, CellType.triangle)    

U_el = VectorElement("Lagrange", mesh.ufl_cell(), 1)
P_el = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([U_el, P_el]))

(ua_n2, pa_n2) = TrialFunctions(W) 
(qu, qp) = TestFunctions(W)
up = fem.Function(W)

U, U_to_W = W.sub(0).collapse(True)
P, P_to_W = W.sub(1).collapse(True)

ua_n1 = fem.Function(U)
pa_n1 = fem.Function(P)

F = (\
        dot(\
            ua_n2-ua_n1
            +grad(pa_n2) 
            , qu)+
        dot(
            pa_n2-pa_n1
            +div(ua_n2)
            , qp)
         )*dx
         

         
problem = fem.LinearProblem(lhs(F),  rhs(F) , u=up)
problem.solve()
_ua, _pa = up.split()


ua_n1.x.array[:] = up.x.array[U_to_W]
ua_n1.x.scatter_forward()    

pa_n1.x.array[:] = up.x.array[P_to_W]
pa_n1.x.scatter_forward()
2 Likes
4

Hi I used DOLFINx 0.6.0 and made minor changes on your code. But it returned me 0 for the solution.

from mpi4py import MPI
from dolfinx import fem, mesh
from ufl import (TestFunctions, TrialFunctions, VectorElement,
                 FiniteElement, MixedElement, dot, grad, div, dx, lhs, rhs)
import numpy as np

nx, ny = 10, 10
msh = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)  

U_el = VectorElement("Lagrange", msh.ufl_cell(), 1)
P_el = FiniteElement("Lagrange", msh.ufl_cell(), 1)
W = fem.FunctionSpace(msh, MixedElement([U_el, P_el]))

(ua_n2, pa_n2) = TrialFunctions(W) 
(qu, qp) = TestFunctions(W)
up = fem.Function(W)

U, U_to_W = W.sub(0).collapse()
P, P_to_W = W.sub(1).collapse()

ua_n1 = fem.Function(U)
pa_n1 = fem.Function(P)

F = (dot(ua_n2 - ua_n1 + grad(pa_n2), qu) + dot(pa_n2 - pa_n1 + div(ua_n2), qp)) * dx
problem = fem.petsc.LinearProblem(lhs(F),  rhs(F) , u=up)
problem.solve()

_ua, _pa = up.split()
ua_n1.x.array[:] = up.x.array[U_to_W]
ua_n1.x.scatter_forward()    

pa_n1.x.array[:] = up.x.array[P_to_W]
pa_n1.x.scatter_forward()
1 Like