I am working on solving a set of two non-linear differential equations that are also time dependent. I am using the version fenics-dolfinx 0.9.0 and fenics-ufl 2024.2.0 and I installed them using conda-forge on my macOS Ventura13.1 (chip M1). My MWE is the following:
import pygmsh
import meshio
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import fem, io
from ufl import TestFunction, TrialFunction, inner, dx, grad, SpatialCoordinate
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
# %% Create mesh
# Generate the spherical surface mesh with pygmsh
with pygmsh.geo.Geometry() as geom:
radius = 1.0 # Radius of the sphere
sphere = geom.add_ball([0, 0, 0], radius, mesh_size=0.05)
mesh = geom.generate_mesh(dim=2)
# Filter out mixed cells and retain only triangles
triangle_cells = mesh.cells_dict["triangle"]
mesh = meshio.Mesh(
points=mesh.points,
cells=[("triangle", triangle_cells)],
)
# Save the filtered mesh in XDMF format
meshio.write("SphericalSurfaceMesh.xdmf", mesh)
# %% Load the mesh and integrate the system
# Load the filtered mesh into FEniCSx
with io.XDMFFile(MPI.COMM_WORLD, "SphericalSurfaceMesh.xdmf", "r") as infile:
fenics_mesh = infile.read_mesh(name="Grid")
Du = fem.Constant(fenics_mesh, 1.0)
# Define the function space for the scalar fields u1 and u2
V = fem.functionspace(fenics_mesh, ("CG", 1))
# Initial conditions
x = SpatialCoordinate(fenics_mesh)
u1 = fem.Function(V) # Current solution for u1
u2 = fem.Function(V) # Current solution for u2
u1_minus = fem.Function(V) # Previous solution for u1
u2_minus = fem.Function(V) # Previous solution for u2
def initial_conditions_u1(x):
return 1.0 + x[0]
def initial_conditions_u2(x):
return 1.0 + x[1]
u1.interpolate(initial_conditions_u1)
u2.interpolate(initial_conditions_u2)
u1_minus.x.array[:] = u1.x.array[:]
u2_minus.x.array[:] = u2.x.array[:]
# Time-stepping parameters
T = 1.0
dt = 0.1
t = dt
# Define test functions
v1 = TestFunction(V)
v2 = TestFunction(V)
# Define nonlinear functions for FHN
def f(u1, u2):
return -u1**3 + u1 - u2
def g(u1, u2):
return 0.1 * (u1 - 0.5 * u2 + 0.1)
# Residuals for u1 and u2
F1 = (inner(u1, v1) - inner(u1_minus, v1) + dt*Du*inner(grad(u1), grad(v1)) - dt*inner(f(u1, u2), v1)) * dx
F2 = (inner(u2, v2) - inner(u2_minus, v2) + dt*Du*inner(grad(u2), grad(v2)) - dt*inner(g(u1, u2), v2)) * dx
# Combine the residuals into a nonlinear problem
F_total = F1 + F2
# Define the nonlinear problem
problem = NonlinearProblem(F_total, fem.Function(V))
solver = NewtonSolver(MPI.COMM_WORLD, problem)
# Configure the solver
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True
# PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
opts["ksp_type"] = "gmres"
opts["pc_type"] = "hypre"
opts["pc_hypre_type"] = "boomeramg"
opts["ksp_rtol"] = 1e-8
ksp.setFromOptions()
# Time-stepping loop
while t <= T:
# Solve the nonlinear problem
n, converged = solver.solve(u1) # Use u1 as the overall solution
if not converged:
raise RuntimeError(f"Newton solver did not converge at time {t}")
# Update previous solutions
u1_minus.x.array[:] = u1.x.array[:]
u2_minus.x.array[:] = u2.x.array[:]
t += dt
# Extract the mesh coordinates for visualization
coordinates = fenics_mesh.geometry.x # Nodal coordinates
x, y, z = coordinates[:, 0], coordinates[:, 1], coordinates[:, 2]
# Plot the solution
u1_values = u1.x.array
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(projection="3d")
ax.scatter(x, y, z, c=u1_values, cmap="viridis", marker="o", s=10)
plt.show()
And the error is the following:
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
File ~/miniconda3/envs/fenicsx-env/lib/python3.13/site-packages/spyder_kernels/customize/utils.py:209, in exec_encapsulate_locals(code_ast, globals, locals, exec_fun, filename)
207 if filename is None:
208 filename = "<stdin>"
--> 209 exec_fun(compile(code_ast, filename, "exec"), globals, None)
210 finally:
211 if use_locals_hack:
212 # Cleanup code
File ~/Desktop/FEMSphere_v3.py:106
103 # Time-stepping loop
104 while t <= T:
105 # Solve the nonlinear problem
--> 106 n, converged = solver.solve(u1) # Use u1 as the overall solution
107 if not converged:
108 raise RuntimeError(f"Newton solver did not converge at time {t}")
File ~/miniconda3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/nls/petsc.py:51, in NewtonSolver.solve(self, u)
48 def solve(self, u: fem.Function):
49 """Solve non-linear problem into function u. Returns the number
50 of iterations and if the solver converged."""
---> 51 n, converged = super().solve(u.x.petsc_vec)
52 u.x.scatter_forward()
53 return n, converged
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 73, Object is in wrong state
It is probably due to how I am defining the variables used by the solver, but I have not found a different way to do it. I have seen similar questions in the forum but the solutions were dependent in the ‘ufl.split’ function, which I do not use in here, or in redefinition of other variables. For what I have read it might also be that the way that I update the variable at t-1 is not the correct way of doing it, but is the only replacement that I found for the old function .assign.