NonlinearPDE_SNESProblem solver for Cahn-Hilliard tutorial

dolfinx
1

Hello everyone. I am interested in solving Cahn-Hilliard tutorial with a modification of adding bounds to the concentration value. Most of the previous discussion suggests to make use of SNES solver. So as per the major instruction for using the SNES solver for nonlinear PDE, I visited the GitHub and make use of class NonlinearPDE_SNESProblem and attempted to solve. The modified code of Cahn-Hilliard tutorial is,

import os
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element, mixed_element
from dolfinx.io import XDMFFile
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import la, default_real_type
from dolfinx.fem import (Function,form, functionspace)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_matrix,set_bc)
from dolfinx.mesh import create_unit_square, CellType
from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):
        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
        apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(F, [self.bc], x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()

lmbda = 1.0e-02  # surface parameter
dt = 5.0e-06  # time step
theta = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 24, 24, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)

# Zero u
u.x.array[:] = 0.0

# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 100 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)

# mu_(n+theta)  
mu_mid = (1.0 - theta) * mu0 + theta * mu

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
 
F = F0 + F1   
bc =[]
problem = NonlinearPDE_SNESProblem(F,u,bc)

# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.25  + 0.01 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

b = la.create_petsc_vector(ME.dofmap.index_map, ME.dofmap.index_map_bs)
J = create_matrix(problem.a)

# Create Newton solver and solve
snes = PETSc.SNES().create()
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)

snes.setTolerances(rtol=1.0e-9, max_it=10)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-9)
snes.getKSP().getPC().setType("lu")

# Output file
file = XDMFFile(MPI.COMM_WORLD, "SNES.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
#  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 =  10 * dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array

while (t < T):
    t += dt
    r = snes.solve(None, u.vector)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    file.write_function(c, t)
file.close()

I wanted to add bounds, but before that I tried to solve with SNES solver, unfortunately it causes error like shown below,

python3: src/petsc4py.PETSc.c:352134: __Pyx_PyCFunction_FastCall: Assertion `!PyErr_Occurred()' failed.
[Legion:37114] *** Process received signal ***
[Legion:37114] Signal: Aborted (6)
[Legion:37114] Signal code:  (-6)
[Legion:37114] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7f17e5dc4520]
[Legion:37114] [ 1] /lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7f17e5e189fc]
[Legion:37114] [ 2] /lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7f17e5dc4476]
[Legion:37114] [ 3] /lib/x86_64-linux-gnu/libc.so.6(abort+0xd3)[0x7f17e5daa7f3]
[Legion:37114] [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x2871b)[0x7f17e5daa71b]
[Legion:37114] [ 5] /lib/x86_64-linux-gnu/libc.so.6(+0x39e96)[0x7f17e5dbbe96]

Can someone please tell that the way of proceeding is correct or not to satisfy my aim of solving a time dependent bound problem as I am not sure of the reason for this error. Also please give me some ideas or tips in solving the same. Any ideas for my problem is greatly appreciated. Thanks in advance for all your time and support.

2

There is no point in applying BCs if you have the empty list

Drop the following lines:

  • apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
  • set_bc(F, [self.bc], x, -1.0)

and change

  • assemble_matrix(J, self.a, bcs=[self.bc]) into assemble_matrix(J, self.a)

With those changes, the code runs for me, although I didn’t open the resulting solution in paraview.

1 Like
3

By @francesco-ballarin pinpointimg, one can spot that
[[self.bc]] in apply lifting should be [self.bc] and similarly
[self.bc] in set_bc should be self.bc, as self.bc should be a list of boundary conditions

1 Like
4

Thanks a lot for your responses Mr.Dokken and Mr.Ballarin, it was very much helpful. I have run the code without error and tried to add the bounds as per your instructions in old topics using setVariableBounds(). But the result looks not as expected. The norm were showing only for few timesteps nonzero and rest everything zero. My edited code with included bounds is

import os
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element, mixed_element
from dolfinx.io import XDMFFile
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import la, default_real_type
from dolfinx.fem import (Function,form, functionspace)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_matrix,set_bc)
from dolfinx.mesh import create_unit_square, CellType
from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):
        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u
    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a)
        #assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()

lmbda = 1.0e-02  # surface parameter
dt = 5.0e-06  # time step
theta = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 24, 24, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)

# Zero u
u.x.array[:] = 0.0

# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 100 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)

# mu_(n+theta)  
mu_mid = (1.0 - theta) * mu0 + theta * mu

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
 
F = F0 + F1   
bc =[]
problem = NonlinearPDE_SNESProblem(F,u,bc)

# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.25  + 0.01 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

b = la.create_petsc_vector(ME.dofmap.index_map, ME.dofmap.index_map_bs)
J = create_matrix(problem.a)

# Create Newton solver and solve
snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonrsls') #Changed even to vinewtonssls
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(rtol=1.0e-8, max_it=50)

snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-8)
snes.getKSP().getPC().setType("lu")

###Setting constraint for concentration
c_min, c_max = Function(ME), Function(ME)
with c_min.vector.localForm() as loc:
    loc.set(0.0)
with c_max.vector.localForm() as loc:
    loc.set(1.0)
snes.setVariableBounds(c_min.vector, c_max.vector)

# Output file
file = XDMFFile(MPI.COMM_WORLD, "SNES.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
#  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 =  80 * dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array

while (t < T):
    t += dt
    r = snes.solve(None, u.vector)
    #print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    file.write_function(c, t)
file.close()

I tried with some previous topic recommendations of changing snes.set Type for both vinewtonrsls and vinewtonssls.
Could you please clarify me where I am going wrong?. Is this approach of SNES for constraining the concentration parameter between 0 and 1 for a problem like this cahn-hilliard is correct?.Because when I plot, the results are not good, the interface becomes invisible and not like distinct phase as a result in Newton method. Any tips or advice will be appreciated.

Thanks for your time and continuous support.

5

Sorry for the trouble. I missed a logic on applying the bound to only one of the dof in the mixed space. Since I have both concentration and chemical potential in my mixed space where I need to bound only the concentration. I have tried two possibilities but both were showing error.

Fun1,dof1 = ME.sub(0).collapse()
c_min = Function(ME)
c_max = Function(ME)
with c_min.sub(0).vector.localForm() as loc:
    loc.set(0.0)
with c_max.sub(0).vector.localForm() as loc:
    loc.set(1.0)
snes.setVariableBounds(c_min.vector, c_max.vector)

and the another approach is

Fun1,dof1 = ME.sub(0).collapse()
c_min = Function(Fun1)
c_max = Function(Fun1)
with c_min.vector.localForm() as loc:
    loc.set(0.0)
with c_max.vector.localForm() as loc:
    loc.set(1.0)
snes.setVariableBounds(c_min.vector, c_max.vector)

This approach yields mismatch between solution vector and bound vector. Could anyone please tell whether is it able to apply these bound to a mixed problem like mine?

Additionally, Is this approach of SNES for constraining the concentration parameter between 0 and 1 for a problem like this cahn-hilliard is correct?.Because when I plot, the results are not good, the interface becomes invisible and not like distinct phase as a result in Newton method. Any tips or advice will be appreciated.

Thanks for your time and continuous support.

6

Dear all, can anyone please help on the above regard?. I want to assign bound to only one dof in my mixed function space. When I tried with above two possibilities, it resulted in error. So please suggest me if something is wrong. Thanks in advance.

7

You have to set bounds for the second sub space as well.
I would set the bounds really high and really low, i.e. something along the lines of:

import os
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element, mixed_element
from dolfinx.io import XDMFFile
import numpy as np
import ufl
from dolfinx import la
from dolfinx.fem import (Function, form, functionspace)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
                               create_matrix)
from dolfinx.mesh import create_unit_square, CellType
from ufl import TrialFunction, derivative, dx, grad, inner


class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):
        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                      mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(
            addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a)
        # assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()


lmbda = 1.0e-02  # surface parameter
dt = 5.0e-06  # time step
theta = 0.5  # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 24, 24, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)

# Zero u
u.x.array[:] = 0.0

# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 100 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)

# mu_(n+theta)
mu_mid = (1.0 - theta) * mu0 + theta * mu

# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + \
    dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - \
    lmbda * inner(grad(c), grad(v)) * dx

F = F0 + F1
bc = []
problem = NonlinearPDE_SNESProblem(F, u, bc)

# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.25 + 0.01 *
                     (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()

b = la.create_petsc_vector(ME.dofmap.index_map, ME.dofmap.index_map_bs)
J = create_matrix(problem.a)

# Create Newton solver and solve
snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonrsls')  # Changed even to vinewtonssls
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(rtol=1.0e-8, max_it=50)

snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-8)
snes.getKSP().getPC().setType("lu")

# Setting constraint for concentration
Fun1, dof1 = ME.sub(0).collapse()
c_min = Function(ME)
c_max = Function(ME)
c_min.sub(0).interpolate(lambda x: np.zeros(x.shape[1], dtype=np.float64))
c_max.sub(0).interpolate(lambda x: np.ones(x.shape[1], dtype=np.float64))
c_min.sub(1).interpolate(lambda x: np.full(x.shape[1], -1e7, dtype=np.float64))
c_max.sub(1).interpolate(lambda x: np.full(x.shape[1], 1e7, dtype=np.float64))
snes.setVariableBounds(c_min.vector, c_max.vector)
# Output file
file = XDMFFile(MPI.COMM_WORLD, "SNES.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
#  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 = 80 * dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
u0.x.array[:] = u.x.array

while (t < T):
    t += dt
    r = snes.solve(None, u.vector)
    # print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    file.write_function(u.sub(0), t)
file.close()
1 Like
8

Thanks a lot Mr.Dokken. The information is very helpful in implementing bound in mixed function space. Unfortunately for problem like Cahn-hilliard is not producing better result with vinewtonrsls or vinewtonssls. All the interface gets invisible.

Hence these type of inequality constrain implementation like my requirement (Fischer-Burmeister constraint), Petsc TAO solver has the options to solve. But unfortunately I am not finding any good examples or solved code in TAO solver in dolfinx. Could you please recommend some ideas or examples on the same. It will be very much helpful.

Thank you very much as always for your time and considerations. I hope I can find few examples such that I can get some ideas for solving the same.

9

I’ve never used PETSc-TAO, so I cannot really comment on its usage.

10

So other than that, any other possibilities to solve a constrain problem with inequality constraint. I have tried already SNES with the class NonlinearPDE_SNESProblem, but with it result is not as expected.

Any other options for proceeding to solve a constrained problem? In older topics I am not succeeded in finding solver other than SNES. So could you please suggest any other ideas.

Thanks for your support.

11

While using a bound constrained nonlinear solver such as TAO is one way, there are multiple ways of doing this. You can impose a suitable constraint as a penalty/regularization term in the weak form itself, and then use an appropriate nonlinear solver (newton, quasi-newton, steepest descent etc.). For e.g., for phase-field fracture you may want to take a look at (non-exhaustive): here and/or here.

You will then need to decide how much violation of the penalty can you afford for your problem for deciding the penalty/regularization parameter.

3 Likes
12

Sorry for continuing the thread. As I attempted to try the same exact procedure for a coupled problem, it again resulting in the same error like before. Even this time, I am solving with boundary conditions, so included lifting, assembly. But still faced the same error. The code is as follows,

import os
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx import fem, default_scalar_type,mesh
from dolfinx import la
from dolfinx.fem import (Function, form, functionspace)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
                               create_matrix, create_vector, set_bc,apply_lifting)
from dolfinx.mesh import create_unit_square, CellType
from ufl import TrialFunction, derivative, ds, dx, grad, inner

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):
        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
        apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(F, [self.bc], x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()

#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(384,64), cell_type=mesh.CellType.quadrilateral)

#Create function space
V = ufl.VectorElement("Lagrange", msh.ufl_cell(),1) #For deformation map
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(),1) 
ME = fem.FunctionSpace(msh, ufl.MixedElement([V, P1, P1]))

#Creating Test functions and function
a, q, v= ufl.TestFunctions(ME)

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

# Split mixed functions

d, p, mu= ufl.split(u)
d0, p0, mu0= ufl.split(u0)

#Applying boundary condition to beam
def left_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))  
ME0, _ = ME.sub(0).collapse()
bottom_value1 = fem.Function(ME0)
bottom_value1.x.array[:]= 0
left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME.sub(0))

def right_corner(x):
    return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],6.0))
ME0y, _ = ME.sub(0).sub(1).collapse()    
bottom_value2 = fem.Function(ME0y)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological((ME.sub(0).sub(1),ME0y),0,right_corner_vertex) 
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0).sub(1))

bcs = [bc_left_corner, bc_right_corner]
#Defining for applying point load

#Defining traction and body forces
Traction=fem.Constant(msh,default_scalar_type((0,0)))
B=fem.Constant(msh,default_scalar_type((0,0)))

#Spatial dimension
dim=len(d)
#Identity tensor
I=ufl.variable(ufl.Identity(dim))
#Deformation gradient
F=ufl.variable(I+grad(d))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E =default_scalar_type(74000)
nu = default_scalar_type(0.33)
G = fem.Constant(msh, E / (2 * (1 + nu)))
K = fem.Constant(msh, E / (3 * (1 - 2 * nu)))
lmbda = fem.Constant(msh, E * nu / ((1 + nu) * (1 - 2 * nu)))
void_c = default_scalar_type(0.001)
gamma = default_scalar_type(10000) #Additional Penalty boundary term

#Defining parameters
p = ufl.variable(p)
epsilon = default_scalar_type(1/64)  # interface thickness
areainteng = default_scalar_type(0.3) #Area specific interface energy
dt = 5e-06  # time step
k = default_scalar_type(0.05) #Dissipation term
vp = p
void_c = default_scalar_type(0.001)

# #Bulk energy expression
Bulk_psi = (0.001-0.001*p + p)*(((G/2) * (((J)**(-2/3) * Ic) - 3)) +  \
                                ((K/2) * ((J**2 /2) - (1/2) - ufl.ln(J)) ))
# Hyper-elasticity
P = ufl.diff(Bulk_psi, F)
Int_Psi = (6*(1/epsilon)* areainteng)*((p**2)*(1-p)**2) 
dintdp = ufl.diff (Int_Psi, p)
dbulkdp = ufl.diff (Bulk_psi, p)
dfdc = -dbulkdp + dintdp

#Weak form equations
F0 = inner(P, grad(a)) * dx - inner(B, a) * dx - inner(Traction, a) * ds
F1 = inner(mu,q)*dx - inner(dfdc,q)*dx - (3*(epsilon)*areainteng*inner(grad(p),grad(q)))*dx 
F2 = inner(p,v)*dx - inner(p0,v)*dx + (dt*k*inner(grad(mu),grad(v)))*dx 
NF = F0+F1+F2
#Solving Nonlinear problem
problem=NonlinearPDE_SNESProblem(NF,u,bcs)

# Interpolate initial condition for the concentration
const_val = 0.5
u.sub(1).interpolate(lambda x: np.full((x.shape[1],), const_val))
u.x.scatter_forward()
b = la.create_petsc_vector(ME.dofmap.index_map, ME.dofmap.index_map_bs)
J = create_matrix(problem.a)
#Solver configurations
snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonrsls')  # Changed even to vinewtonssls
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(rtol=1.0e-8, max_it=50)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-8)
snes.getKSP().getPC().setType("lu")

# Step in time
t = 0.0
#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 1 * dt
else:
    T = 5* dt
u0.x.array[:] = u.x.array
while (t < T):
    t += dt   
    r = snes.solve(None, u.vector)
    u0.x.array[:] = u.x.array

My apologies for extending the same problem. I am not understanding, why it is not working for this coupled problem. Could you please help, as the error is not understandable.
Thanks for your time and considerations.