LinearProblem error with mixed elements

dolfinx
1

Hello, everybody,

I am trying to solve a system of PDEs with a weak formulation like this:

image
image876×143 15.7 KB

and

image
image877×195 23 KB

where I should find X_{h}, H_{h} and a_{h} each time.

My MWE is as follows:

# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression, form
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem,assemble_matrix,assemble_vector, apply_lifting
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-04  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.2)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)


    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1],gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

X = Function(ME2)  # current solution
X0 = Function(ME2)  # solution from previous converged step

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution
(Xh0, H0) = ufl.TrialFunctions(ME2)  # solution from previous converged step

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

X.sub(0).interpolate(x_exp)
X.sub(1).interpolate(H_exp)
# -

u0.x.array[:] = u.x.array
X0.x.array[:] = X.x.array

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5, 0])

Xh0

# Weak statement of the equations

k = dt
d = 8.8676
gamma = 379.2100
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 3

F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
   -(gamma*(u1*u1*u2-u1+a))*q*dx \
  + ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
  -(gamma*(-u1*u1*u2+b))*v*dx 

#F_x = inner(((Xh - Xh0) / k),n_oriented)*phi*dx + ks*H*phi*dx \
#  +kb*inner(grad(H),grad(phi))*dx - dot(kp,u_vec)*phi*dx \
#    + inner(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx

ax = inner((Xh  / k),
            n_oriented)*phi*dx + ks*H*phi*dx \
    + kb*inner(grad(H),grad(phi))*dx  \
    + inner(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx

Lx = inner((Xh0  / k),n_oriented)*phi*dx + dot(kp,u_vec)*phi*dx

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

problemX = LinearProblem(ax, Lx)

...

The nonlinear problem works well, but I cannot define the Linear Problem for the second system. I have already tried to implement the problem with different formulations but am not able to find what produces the error. I think I’m confusing the use of ‘Function’ and ‘TrialFunction’ because I’m not sure if I’ve defined the initial conditions correctly either.

The error I obtain is something like this:

Traceback (most recent call last):
  File "/home/sdcardozof/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 68, in get_cached_module
    open(c_filename, "x")
FileExistsError: [Errno 17] File exists: '/home/sdcardozof/.cache/fenics/libffcx_forms_04fbc351847d98607131abeed21bfe8eee7e567f.c'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/mnt/d/Universidad/UBC/CodesFenics/SchnakenbergFx2D_circle_moving_V2.py", line 196, in <module>
    problemX = LinearProblem(ax, Lx)
               ^^^^^^^^^^^^^^^^^^^^^
...

  File "/home/sdcardozof/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 627, in __del__
    self._solver.destroy()
    ^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'

Thanks in advance

2

Hi.

The problem is in the line:

Lx = inner((Xh0 / k), n_oriented)*phi*dx + dot(kp, u_vec)*phi*dx

because Xh0 is a TrialFunction of ME2, so it can’t be on the RHS. If (Xh0, H0) are going to be updated, then you can define them as

tempsol = Function(ME2)
(Xh0, H0) = ufl.split(tempsol)  # solution from previous converged step

This solves the error. For the initial boundary conditions, I think the MWE is quite incomplete and has many unused variables in comparison with the posted formulation, so it is very hard to someone to help you.

Cheers.

2 Likes
3

Hi Jesus,

I believe it works properly; nevertheless, I have some concerns regarding the initial conditions. Basically, I am implementing a method that solves for a new mesh at time n+1. The variable X0 represents the initial mesh coordinates. I’ve implemented the initial conditions as follows:

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

X0.sub(0).interpolate(x_exp)

But when I run my code, I obtain strange results, and I guess that could be due to the initial conditions. I tried to write the initial condition of the variable X0 in a VTK file, but it seems to have just one component in the x-direction, and it should be outward (at least in the initial condition).

This is my MWE:

# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression, form
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem,assemble_matrix,assemble_vector, apply_lifting
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-05  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)


    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# X = Function(ME2)  # current solution
X0 = Function(ME2)  # solution from previous converged step
X0.name = "X0"

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution
#(Xh0, H0) = ufl.TrialFunctions(ME2)  # solution from previous converged step

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

# # Split mixed functions
# Xh, H = ufl.split(X)
Xh0, H0 = ufl.split(X0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

X0.x.scatter_forward()

print (f"X array  {X0.sub(0).x.array[:]}")

# -

#X0.x.array[:] = uh.x.array

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])

# Weak statement of the equations

k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 2

F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
   -(gamma*(u1*u1*u2-u1+a))*q*dx \
  + ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
  -(gamma*(-u1*u1*u2+b))*v*dx 

ax = dot(((Xh) / k),n_oriented)*phi*dx + kb*inner(grad(H),grad(phi))*dx  - 1/2 *kb* H0**2 * H *phi*dx + ks*H*phi*dx\
    + dot(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx
    #+ kb*inner(grad(H),grad(phi))*dx  + ks*H*phi*dx\
Lx =  dot(((Xh0) / k),n_oriented)*phi*dx + dot(kp,u_vec)*phi*dx

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

X_h = Function(ME2)
X_h.name = "X_h"

problemX = LinearProblem(ax, Lx, u = X_h)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 3 * dt
else:
    T = 5 * dt

# Get the sub-space for c and the corresponding dofs 

normal_vec.interpolate(normal_expr)

V0, dofs = ME.sub(0).collapse()
W02, dofs2 = ME2.sub(0).collapse()

u1 = u.sub(0)
u2 = u.sub(1)

Xh0 = X0.sub(0)
H0 = X0.sub(1)

file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu1.xdmf", "w")
# file2 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu2.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/X_coord.xdmf", "w")

file1.write_mesh(mesh)
#file2.write_mesh(mesh)
file3.write_mesh(mesh)
file4.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

norm_stop = 0.1e-6

fileL2 = open('L2.dat', 'w') 

z = 0
inc = 0 
impresion = 1

#print (f"MESH array {mesh.geometry.x[:,]}")

while (t < T):
    t += dt

    problemX.solve()
    r = solver.solve(u)
    inc += 1

    print(f"Step {int(t/dt)}: num iterations: {r[0]}")

    l2_norm = np.linalg.norm(u.x.array[dofs]-u0.x.array[dofs])/dt
    #print (f"L2_norm = {l2_norm}")

    fileL2.write(f"{t}\t{l2_norm}\n")
    u0.x.array[:] = u.x.array
    X0.x.array[:] = X_h.x.array

    num=inc/impresion-1 

    # print (f"MESH array shape {mesh.geometry.x[:,].shape}")
    # print (f"X array shape {X_h.x.array[dofs2].shape}")
    # print (f"MESH array {mesh.geometry.x[:,]}")
    # print (f"X array  {X_h.x.array[dofs2]}")

    #print (f"X array  {X0.sub(0).x.array[:]}")

    X_coord = X_h.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)
    DeltaX = mesh.geometry.x[:,[0,1]]-X_coord

    #print (f"Delta X array  {DeltaX}")
    mesh.geometry.x[:,0] = X_coord[:,0]
    mesh.geometry.x[:,1] = X_coord[:,1]

    normal_vec.interpolate(normal_expr)

    if (int(num)==z) or (int(num-z)==0):
        us, vs = ufl.split(u)

        # Save solution to file (VTK)

        name = "mesh_at_t"+str(t)
        mesh.name = name

        file1.write_mesh(mesh)
        file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        file3.write_mesh(mesh)
        file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        file4.write_mesh(mesh)
        file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        z += 1
    

    # if l2_norm < norm_stop:
    #     print("l2_norm convergence")
    #     break

file1.close()

fileL2.close()

print(l2_norm)
#print(H.x.array)

How can I define the initial conditions properly?

Thank you.

4

Let’s start by reducing the example to only contain initial conditions, as that is what you are asking about:


from mpi4py import MPI

import ufl
from basix.ufl import element, mixed_element
from dolfinx import log
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.io import gmshio
from ufl import div

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-05  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys

    print("This demo requires gmsh to be installed")
    sys.exit(0)


def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle."""
    model.add(name)
    model.setCurrent(name)

    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model


gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

mesh, ct, ft = gmshio.model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

mesh.name = "test"
ct.name = f"{mesh.name}_cells"
ft.name = f"{mesh.name}_facets"


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)


P1 = element("Lagrange", mesh.basix_cell(), 2)
ME = functionspace(mesh, mixed_element([P1, P1]))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x / ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign * n

V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
V2 = element("Lagrange", mesh.basix_cell(), 2)

ME2 = functionspace(mesh, mixed_element([V1, V2]))

# Trial and test functions of the space `ME` are now defined:

X0 = Function(ME2)  # solution from previous converged step
X0.name = "X0"


V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign * n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

X0.x.scatter_forward()

X00 = X0.sub(0).collapse()
X00.name = "X00"
X01 = X0.sub(1).collapse()
X01.name = "X01"
from dolfinx.io import VTXWriter

with VTXWriter(mesh.comm, "X0.bp", [X00], engine="BP4") as bp:
    bp.write(0.0)


with VTXWriter(mesh.comm, "X1.bp", [X01], engine="BP4") as bp:
    bp.write(0.0)

yielding

image
image903×732 46.6 KB