Nonlinear solver converging to wrong solution

dolfinx
1

Hi,

I have to solve steady-state heat flow problem with radiation boundary conditions. I begin by considering a simple 2D MWE in which there is a heat source with left and right facets of the domain imposed with radiation condition. However, the nonlinear solver converges to a uniform heat map which obviously does not satisfy the radiation condition -dT(x=0)/dx = dT(x=1)/dx = T^4. I would appreciate if experts here could help me out.

Below is a MWE evaluated in FEniCSx 0.5.0 installed through Conda:

import ufl, dolfinx, mpi4py, petsc4py, pyvista
import numpy as np

# create a unit square geometry
mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 20, 20)
# tag mesh boundaries, https://jorgensd.github.io/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# define variational form
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))
uh2 = dolfinx.fem.Function(V)
v = ufl.TestFunction(V)

# source heat flux
Qs = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(1))

ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
F = ufl.inner(ufl.grad(uh2), ufl.grad(v))*ufl.dx - v*Qs*ufl.dx
# apply radiation boundary on left and right edges
F +=  (uh2**4 * v)*ds(1) + (uh2**4 * v)*ds(2)

# create nonlinear solver
# https://jorgensd.github.io/dolfinx-tutorial/chapter2/nonlinpoisson_code.html
problem = dolfinx.fem.petsc.NonlinearProblem(F, uh2)
solver = dolfinx.nls.petsc.NewtonSolver(mpi4py.MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = False

ksp = solver.krylov_solver
opts = petsc4py.PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "pardiso"
ksp.setFromOptions()

# dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
# dolfinx.log.set_log_level()
n, converged = solver.solve(uh2)
assert(converged)
print(f"Number of interations: {n:d}")

# visualize result
u_topology, u_cell_types, u_geometry = dolfinx.plot.create_vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh2.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

Thanks!

2

Are you running dolfinx in real mode ? The Newton solver can be actually quite sensitive to mode change.

Also if you are running in complex mode, you can no longer simply do uh2**4*v, you have to use inner (there’s a complex conjugate in there).

3

I did not explicitly switch into complex mode so that means it should be real, right? Is there a way to make sure that I am working in real mode?

4

Check the type of ScalarType

from petsc4py import PETSc
print(PETSc.ScalarType)
5

I get this in the output:

<class ‘numpy.float64’>

6

The problem is solved.

The solver was unable to converge with a zero initial guess. Doing the following made it to converge:

uh2 = dolfinx.fem.Function(V)
uh2.x.array[:] = 1