Thanks a lot @dokken. It is strange that I never had this warning before though. And this line: opts[f"{option_prefix}ksp_max_it"]= 1000 can also be found in the dolfinx tutorial, I guess it returns this warning now?
Anyway, by modifying the max_it parameter using your suggestion, I could push a little bit further the resolution of my problem. But it’s still far from satisfying. Krylov solver fails now, with the message:
Traceback (most recent call last):
File "/root/shared/Nye_example.py", line 132, in <module>
n, converged = solver.solve(TempVolt)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py", line 46, in solve
n, converged = super().solve(u.vector)
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 76, Error in external library
I am seriously considering to translate my code into using a linear solver instead of the non linear one. Do you know about possible drawbacks of using the NonlinearSolver as opposed to the Linear one?
By the way, my full code is:
import numpy as np
from dolfinx import log
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative)
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/rectangles.msh", MPI.COMM_WORLD, gdim=2)
# Reminder of the mesh:
#Physical Surface("material", 11) = {2, 1, 3};
#Physical Curve("bottom_left", 12) = {10};
#Physical Curve("top_right", 13) = {6};
inj_current_curve = 12
out_current_curve = 13 # This corresponds to the curve the current leaves the material.
reading_voltage_surface_0 = 12
reading_voltage_surface_1 = 13
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), 2)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
rho = 10
sigma = 1.0 / rho
kappa = 144
S_xx = 100e-6
S_xx = 0.0
S_yy = 3 * S_xx
Seebeck_tensor = as_tensor([[S_xx, 0], [0, S_yy]])
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs = locate_dofs_topological(
ME.sub(1), mesh.topology.dim-1, left_facets)
left_dofs_temp = locate_dofs_topological(
ME.sub(0), mesh.topology.dim-1, left_facets)
right_dofs_temp = locate_dofs_topological(
ME.sub(0), mesh.topology.dim-1, right_facets)
T_cold = 300.0
bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold+10.0), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
the_current = 0.21 # Current, in amperes.
J = the_current / \
assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-kappa * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
F_T = Fourier_term + Joule_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J * ds(out_current_curve) + v * J * ds(inj_current_curve)
weak_form = F_T + F_V
# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
Jac = derivative(weak_form,TempVolt,dTV)
print(f''' ------- Pre-processing --------
Current density: {J}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
This should correspond to the current injected: {assemble_scalar(form(J*ds(out_current_curve)))}
''')
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-9
solver.report = True
solver.max_it = 10000
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"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
I cannot provide the mesh file here due to the limitation of Discourse.
Edit: I fixed the error using the suggestion given in PETSc error 76 dolfinx - #4 by hawkspar. In short, using ksp instead of lu fixed the problem.