Time-dependent mesh and nonlinear problem (PETSC Error after migrating from FEniCSX 0.7.2 to 0.8.0)

General
1

Hello,
I am solving a time-dependent Cahn-Hilliard equation, where the mesh is re-defined in each time step.
A new mesh is created at each time step, and so is the problem definition.
I have a Minimal Working Example as follows, which works well using FEniCSX 0.7.2.

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
from math import tanh

import ufl
from ufl import dx, grad, inner, dot
from basix.ufl import element, mixed_element

from dolfinx import log, plot, fem
from dolfinx.fem import Function, functionspace
from dolfinx.io import VTKFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

log.set_output_file("log.out")

# Define parameters

lmbda = 1.0e-4  # non-local energy parameter
alpha = 1.0e+0 # local energy parameter
xi = (lmbda/(2.0*alpha))**0.5  # diffuse-interface thickness
mcoef = 1.0e-2 # mobility
dt = 5.0e-03  # time step
T = 20 * dt # terminal time
t = 0.0 # current time

# Create mesh

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

# define finite element 

P2 = element("Lagrange", msh.basix_cell(), 2)

# define mixed finite element and function space for c and mu

ME = mixed_element([P2, P2])
W = functionspace(msh, ME)

# define Test functions

(q, v) = ufl.TestFunctions(W)

# define Solutions

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

# Split mixed functions

c, mu = ufl.split(u)
c0, _ = ufl.split(u0)

# Interpolate initial condition for c

def ellips(x):
    radi = 0.3
    x0 = 0.5
    y0 = 0.5
    r = ( - ( (x[0]-x0)**2 + 4.0*(x[1]-y0)**2 )**0.5 + radi ) / (2.0 * xi)
    c = 0.5*( np.tanh(r) + 1.0 )
    return c

u0.sub(0).interpolate(ellips)
u0.x.scatter_forward()

# define nonlinear problem

dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
F0 += inner(c, q) * dx - inner(c0, q) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1

problem = NonlinearProblem(F, u)

# define Newton solver

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

# Customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options

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

# write vtk files

ch = u0.sub(0).collapse()
ch.name="c"

muh = u0.sub(1).collapse()
muh.name="mu"

vtk = VTKFile(MPI.COMM_WORLD, "output1/cmu.pvd", "w")
vtk.write_function([ch,muh], 0)

# time stepping

step = 0

while (t < T):
    
    t += dt
    step += 1

    # solve

    r = solver.solve(u)

    print(f"Step {step}: num iterations: {r[0]}")

    # new mesh msh2

    msh2 = create_unit_square(MPI.COMM_WORLD, 10+step, 10+step, CellType.triangle)

    # define Functions in msh2

    P2_2 = element("Lagrange", msh2.basix_cell(), 2)
    ME_2 = mixed_element([P2_2, P2_2])
    W_2 = functionspace(msh2, ME_2)
    X_2 = functionspace(msh2, P2_2)

    c_2 = Function(X_2)
    c0_2 = Function(X_2)
    mu_2 = Function(X_2)
    mu0_2 = Function(X_2)

    # project solution u defined in msh to Functions in msh2
    
    c_interpolation_data = fem.create_nonmatching_meshes_interpolation_data(
        c_2.function_space.mesh._cpp_object,
        c_2.function_space.element,
        u.sub(0).function_space.mesh._cpp_object,
        padding=1.e-10)
    
    mu_interpolation_data = fem.create_nonmatching_meshes_interpolation_data(
        mu_2.function_space.mesh._cpp_object,
        mu_2.function_space.element,
        u.sub(1).function_space.mesh._cpp_object,
        padding=1.e-10)
    
    c0_2.interpolate(u.sub(0), nmm_interpolation_data=c_interpolation_data)
    mu0_2.interpolate(u.sub(1), nmm_interpolation_data=mu_interpolation_data)

    # create mesh and function

    msh = msh2
    P2 = P2_2
    ME = ME_2
    W = W_2
    
    u = Function(W_2)
    u0 = Function(W_2)
    
    # interpolate previous solution u0

    u0.sub(0).interpolate(c0_2)
    u0.sub(1).interpolate(mu0_2)

    u0.x.scatter_forward()

    # Test functions
    # q for c and v for mu

    (q, v) = ufl.TestFunctions(W)

    # split u into c and mu

    c, mu = ufl.split(u)
    c0, _ = ufl.split(u0)
    
    # define nonlinear problem

    dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

    F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
    F0 += inner(c, q) * dx - inner(c0, q) * dx
    F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
    F = F0 + F1

    problem = NonlinearProblem(F, u)
    
    # define newton solver
    
    solver = NewtonSolver(MPI.COMM_WORLD, problem) ## PETSC ERROR
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6

    # Customize the linear solver used inside the NewtonSolver by
    # modifying the PETSc options

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

    # define c and mu from u to write output files

    ch = u0.sub(0).collapse()
    ch.name="c"

    muh = u0.sub(1).collapse()
    muh.name="mu"

    # write files
    
    vtk.write_function([ch, muh], t)

# close files

vtk.close()

As I migrate to FEniCSX 0.8.0, I have made only small updates in the code as follows.
However, it returns an PETSC Error in the middle of the time step, in a random way.

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
from math import tanh

import ufl
from ufl import dx, grad, inner, dot
from basix.ufl import element, mixed_element

from dolfinx import log, plot, fem
from dolfinx.fem import Function, functionspace
from dolfinx.io import VTKFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

log.set_output_file("log.out")

# Define parameters

lmbda = 1.0e-4  # non-local energy parameter
alpha = 1.0e+0 # local energy parameter
xi = (lmbda/(2.0*alpha))**0.5  # diffuse-interface thickness
mcoef = 1.0e-2 # mobility
dt = 5.0e-03  # time step
T = 20 * dt # terminal time
t = 0.0 # current time

# Create mesh

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

# define finite element 

P2 = element("Lagrange", msh.basix_cell(), 2)

# define mixed finite element and function space for c and mu

ME = mixed_element([P2, P2])
W = functionspace(msh, ME)

# define Test functions

(q, v) = ufl.TestFunctions(W)

# define Solutions

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

# Split mixed functions

c, mu = ufl.split(u)
c0, _ = ufl.split(u0)

# Interpolate initial condition for c

def ellips(x):
    radi = 0.3
    x0 = 0.5
    y0 = 0.5
    r = ( - ( (x[0]-x0)**2 + 4.0*(x[1]-y0)**2 )**0.5 + radi ) / (2.0 * xi)
    c = 0.5*( np.tanh(r) + 1.0 )
    return c

u0.sub(0).interpolate(ellips)
u0.x.scatter_forward()

# define nonlinear problem

dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
F0 += inner(c, q) * dx - inner(c0, q) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1

problem = NonlinearProblem(F, u)

# define Newton solver

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

# Customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options

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

# write vtk files

ch = u0.sub(0).collapse()
ch.name="c"

muh = u0.sub(1).collapse()
muh.name="mu"

vtk = VTKFile(MPI.COMM_WORLD, "output1/cmu.pvd", "w")
vtk.write_function([ch,muh], 0)

# time stepping

step = 0

while (t < T):
    
    t += dt
    step += 1

    # solve

    r = solver.solve(u)

    print(f"Step {step}: num iterations: {r[0]}")

    # new mesh msh2

    msh2 = create_unit_square(MPI.COMM_WORLD, 10+step, 10+step, CellType.triangle)

    # define Functions in msh2

    P2_2 = element("Lagrange", msh2.basix_cell(), 2)
    ME_2 = mixed_element([P2_2, P2_2])
    W_2 = functionspace(msh2, ME_2)
    X_2 = functionspace(msh2, P2_2)

    c_2 = Function(X_2)
    c0_2 = Function(X_2)
    mu_2 = Function(X_2)
    mu0_2 = Function(X_2)

    # project solution u defined in msh to Functions in msh2
    
    c_interpolation_data = fem.create_nonmatching_meshes_interpolation_data(
        c_2.function_space.mesh,
        c_2.function_space.element,
        u.sub(0).function_space.mesh,
        padding=1.e-10)
    
    mu_interpolation_data = fem.create_nonmatching_meshes_interpolation_data(
        mu_2.function_space.mesh,
        mu_2.function_space.element,
        u.sub(1).function_space.mesh,
        padding=1.e-10)
    
    c0_2.interpolate(u.sub(0), nmm_interpolation_data=c_interpolation_data)
    mu0_2.interpolate(u.sub(1), nmm_interpolation_data=mu_interpolation_data)

    # create mesh and function

    msh = msh2
    P2 = P2_2
    ME = ME_2
    W = W_2
    
    u = Function(W_2)
    u0 = Function(W_2)
    
    # interpolate previous solution u0

    u0.sub(0).interpolate(c0_2)
    u0.sub(1).interpolate(mu0_2)

    u0.x.scatter_forward()

    # Test functions
    # q for c and v for mu

    (q, v) = ufl.TestFunctions(W)

    # split u into c and mu

    c, mu = ufl.split(u)
    c0, _ = ufl.split(u0)
    
    # define nonlinear problem

    dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

    F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
    F0 += inner(c, q) * dx - inner(c0, q) * dx
    F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
    F = F0 + F1

    problem = NonlinearProblem(F, u)
    
    # define newton solver
    
    solver = NewtonSolver(MPI.COMM_WORLD, problem) ## PETSC ERROR
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6

    # Customize the linear solver used inside the NewtonSolver by
    # modifying the PETSc options

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

    # define c and mu from u to write output files

    ch = u0.sub(0).collapse()
    ch.name="c"

    muh = u0.sub(1).collapse()
    muh.name="mu"

    # write files
    
    vtk.write_function([ch, muh], t)

# close files

vtk.close()

Can anyone help with this issue?
The original one works well with FEniCSX 0.7.2, but the updated one for FEniCSX 0.8.0 produces an PETSC Error in the middle of time step in a random way.

2

Could you try v0.9.0?

3

Thanks Dokken!, it is working fine on v0.9.0.

The updated code for v0.9.0 is as follows, for anyone who might find it interesting.

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
from math import tanh

import ufl
from ufl import dx, grad, inner, dot
from basix.ufl import element, mixed_element

from dolfinx import log, plot, fem
from dolfinx.fem import Function, functionspace
from dolfinx.io import VTKFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

log.set_output_file("log.out")

# Define parameters

lmbda = 1.0e-4  # non-local energy parameter
alpha = 1.0e+0 # local energy parameter
xi = (lmbda/(2.0*alpha))**0.5  # diffuse-interface thickness
mcoef = 1.0e-2 # mobility
dt = 5.0e-03  # time step
T = 20 * dt # terminal time
t = 0.0 # current time

# Create mesh

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

# define finite element 

P2 = element("Lagrange", msh.basix_cell(), 2)

# define mixed finite element and function space for c and mu

ME = mixed_element([P2, P2])
W = functionspace(msh, ME)

# define Test functions

(q, v) = ufl.TestFunctions(W)

# define Solutions

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

# Split mixed functions

c, mu = ufl.split(u)
c0, _ = ufl.split(u0)

# Interpolate initial condition for c

def ellips(x):
    radi = 0.3
    x0 = 0.5
    y0 = 0.5
    r = ( - ( (x[0]-x0)**2 + 4.0*(x[1]-y0)**2 )**0.5 + radi ) / (2.0 * xi)
    c = 0.5*( np.tanh(r) + 1.0 )
    return c

u0.sub(0).interpolate(ellips)
u0.x.scatter_forward()

# define nonlinear problem

dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
F0 += inner(c, q) * dx - inner(c0, q) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1

problem = NonlinearProblem(F, u)

# define Newton solver

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

# Customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options

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

# write vtk files

ch = u0.sub(0).collapse()
ch.name="c"

muh = u0.sub(1).collapse()
muh.name="mu"

vtk = VTKFile(MPI.COMM_WORLD, "output1/cmu.pvd", "w")
vtk.write_function([ch,muh], 0)

# time stepping

step = 0

while (t < T):
    
    t += dt
    step += 1

    # solve

    r = solver.solve(u)

    print(f"Step {step}: num iterations: {r[0]}")

    # new mesh msh2

    msh2 = create_unit_square(MPI.COMM_WORLD, 10+step, 10+step, CellType.triangle)
    cmap = msh2.topology.index_map(msh2.topology.dim)
    num_cells = cmap.size_local + cmap.num_ghosts
    all_cells = np.arange(num_cells, dtype=np.int32)

    # define Functions in msh2

    P2_2 = element("Lagrange", msh2.basix_cell(), 2)
    ME_2 = mixed_element([P2_2, P2_2])
    W_2 = functionspace(msh2, ME_2)
    X_2 = functionspace(msh2, P2_2)

    c_2 = Function(X_2)
    c0_2 = Function(X_2)
    mu_2 = Function(X_2)
    mu0_2 = Function(X_2)

    # project solution u defined in msh to Functions in msh2
    
    c_interpolation_data = fem.create_interpolation_data(
        c_2.function_space,
        u.sub(0).function_space,
        all_cells,
        padding=1.e-10)
    
    mu_interpolation_data = fem.create_interpolation_data(
        mu_2.function_space,
        u.sub(1).function_space,
        all_cells,
        padding=1.e-10)
    
    c0_2.interpolate_nonmatching(u.sub(0), all_cells, interpolation_data=c_interpolation_data)
    mu0_2.interpolate_nonmatching(u.sub(1), all_cells, interpolation_data=mu_interpolation_data)

    # create mesh and function

    msh = msh2
    P2 = P2_2
    ME = ME_2
    W = W_2
    
    u = Function(W_2)
    u0 = Function(W_2)
    
    # interpolate previous solution u0

    u0.sub(0).interpolate(c0_2)
    u0.sub(1).interpolate(mu0_2)

    u0.x.scatter_forward()

    # Test functions
    # q for c and v for mu

    (q, v) = ufl.TestFunctions(W)

    # split u into c and mu

    c, mu = ufl.split(u)
    c0, _ = ufl.split(u0)
    
    # define nonlinear problem

    dfdc = 2.0 * alpha * c * ( 1.0 - c ) * ( 1.0 - 2.0 * c )

    F0 = dt * mcoef * inner(grad(mu), grad(q)) * dx
    F0 += inner(c, q) * dx - inner(c0, q) * dx
    F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
    F = F0 + F1

    problem = NonlinearProblem(F, u)
    
    # define newton solver
    
    solver = NewtonSolver(MPI.COMM_WORLD, problem) ## PETSC ERROR
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6

    # Customize the linear solver used inside the NewtonSolver by
    # modifying the PETSc options

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

    # define c and mu from u to write output files

    ch = u0.sub(0).collapse()
    ch.name="c"

    muh = u0.sub(1).collapse()
    muh.name="mu"

    # write files
    
    vtk.write_function([ch, muh], t)

# close files

vtk.close()