Hi all,
First time posting and related to a similar question posted a few days ago (i.e., PETSc error code is: 73, Object is in wrong state for simple hyperelastic problem).
After following through the tutorials provided for FEniCSx, I tried to solve a 2D nonlinear Reaction-Diffusion PDE (Lagrangian Formulation) over a circular domain containing two reaction species with No flux Neumann Boundary Conditions. Attempting to formulate the problem and solve the nonlinear time-discretised system results in the following Error code 73 in PETSc:
line 158, in <module>
num_its, converged = solver.solve(c)
^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx/lib/python3.12/site-packages/dolfinx/nls/petsc.py", line 47, in solve
n, converged = super().solve(u.x.petsc_vec)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 73, Object is in wrong state
As a side note, for the various terms when formulating the problem there is an additional quirk in the solver (which may or may not be related but results in the same error code). I have left in all of term1-4 for both Func_A and Func_I as there is an additional issue related to the solver.convergence_criterion setting. If solver.convergence_criterion = “incremental”, then all of the individual terms result in Error 73. However, if solver.convergence_criterion = “residual”, only term3_A and term3_I result in Error 73. Not sure if this is related or not hence why it was kept in the MWE.
Moreover, I’m not sure if I’ve set up my initial conditions properly. These should be the constants 0 and 1. Is this how you would set up constant initial conditions? Not sure if this is what’s causing the problem.
Additional Information:
OS: macOS Sonoma 14.5
FEniCSx Version: 0.8.0
Installation Procedure (FEniCSx/dolfinx): micromamba (conda-forge)
Installation Procedure (gmsh): micromamba (pip)
The MWE code is as follows:
#-------------
# IMPORTS
#-------------
import ufl
import numpy as np
import gmsh
from mpi4py import MPI
from dolfinx import fem, mesh
from dolfinx.io import gmshio
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import basix
#-------------
# CREATE MESH
#-------------
# Initialise the GMSH module
gmsh.initialize()
# Create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures
membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1) # x, y, z, x-radius, y-radius
gmsh.model.occ.synchronize()
# Make membrane a physical surface for GMSH to recognise when generating the mesh
gdim = 2 # 2D Geometric Dimension of the Mesh
gmsh.model.addPhysicalGroup(gdim, [membrane], 0) # Dimension, Entity tag, Physical tag
# Generate 2D Mesh with uniform mesh size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.25)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.25)
gmsh.model.mesh.generate(gdim)
# Create a mesh from the GMSH model
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
# Finalize GMSH
gmsh.finalize()
#---------------------
# CREATE FEM FUNCTIONS
#---------------------
# Time Discretisation
t = 0 # Initial Time
T = 200 # Final Time
Nsteps = 2000 # Number of Time Steps
dt = T/Nsteps # Time Step Size
# Global Spatial Coordinates as a Vector-Valued Expression i.e. x[0] = x, x[1] = y
x = ufl.SpatialCoordinate(domain) # This may not be needed
P1 = basix.ufl.element(family="Lagrange", cell=domain.topology.cell_name(), degree=1)
element = basix.ufl.mixed_element([P1, P1])
C = fem.functionspace(mesh=domain, element=element)
# Create function space for concentrations
c = fem.Function(C)
c.name = "c"
c_n = fem.Function(C)
c_n.name = "c_n"
c_A, c_I = c.split()
c_A_n, c_I_n = c_n.split()
w_A, w_I = ufl.TestFunction(C)
#-------------------------------------------------------
# CONSTANT INITIAL CONDITIONS AND NO BOUNDARY CONDITIONS
#-------------------------------------------------------
# Initial Conditions of Species I=1 and Species A=0
c_A_0 = fem.Constant(domain, 0.0)
c_I_0 = fem.Constant(domain, 1.0)
c_A.interpolate(lambda x: np.full((x.shape[1],), c_A_0))
c_I.interpolate(lambda x: np.full((x.shape[1],), c_I_0))
c_A_n.interpolate(lambda x: np.full((x.shape[1],), c_A_0))
c_I_n.interpolate(lambda x: np.full((x.shape[1],), c_I_0))
#---------------------------
# REACTION-DIFFUSION PROBLEM
#---------------------------
# Reaction Function
def reaction_function(cA, cI):
delta = fem.Constant(domain, 3.0) # Deactivation Rate
b_0 = fem.Constant(domain, 2.0) # Basal Activation Rate
c_hill = fem.Constant(domain, 1.0) # Hill Function Magnitude
cA_0 = fem.Constant(domain, 0.5)
return (b_0 + c_hill * (cA**6/(cA_0**6 + cA**6))) * cI - delta * cA
# Setting up tensors
spa_dim = len(x) # Spatial Dimension
I = ufl.Identity(spa_dim) # Identity Tensor
F = ufl.variable(ufl.Identity(spa_dim)) # Deformation Gradient Tensor (e.g. Identity Tensor)
F_inv = ufl.inv(F)
F_inv_tra = ufl.transpose(F_inv)
J = fem.Constant(domain, 1.0) # Jacobian Function (e.g. Constant or Pre-Defined Function)
J_n = J
# Examine Isotropic Diffusion where D_A=D*I (Diffusivity Constant multiplied by the Identity Tensor)
D_A = 2 * I
D_A_eq = J * F_inv * D_A * F_inv_tra
D_I = 4 * I
D_I_eq = J * F_inv * D_I * F_inv_tra
# Measure Set Up
metadata = {"quadrature_degree": 3}
dx = ufl.Measure('dx', domain=domain, metadata=metadata)
# Weak Form Set Up
term1_A = J * (c_A - c_A_n)/dt * w_A * dx
term2_A = ufl.inner(D_A_eq * ufl.grad(c_A), ufl.grad(w_A)) * dx
term3_A = -J * reaction_function(c_A, c_I) * w_A * dx
term4_A = c_A * (J - J_n)/dt * w_A * dx
Func_A = term1_A + term2_A + term3_A + term4_A
term1_I = J * (c_I - c_I_n)/dt * w_I * dx
term2_I = ufl.inner(D_I_eq * ufl.grad(c_I), ufl.grad(w_I)) * dx
term3_I = J * reaction_function(c_A, c_I) * w_I * dx
term4_I = c_I * (J - J_n)/dt * w_I * dx
Func_I = term1_I + term2_I + term3_I + term4_I
Func = Func_A + Func_I
# Create Nonlinear Problem and Solver
problem = NonlinearProblem(Func, c, bcs=[])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental" # Incremental causes State issues but Residual only causes for term3_A and term3_I
solver.report = True
# Update Function for each time step of the reaction-diffusion problem
for i in range(Nsteps):
t += dt
# Solve PDE for the current time step
num_its, converged = solver.solve(c)
assert converged
c.x.scatter_forward
# Update solution (c) at previous time step (c_n)
c_A_n.x.array[:] = c_A.x.array
c_I_n.x.array[:] = c_I.x.array
Cheers,
Volkan
…
).