Derivatives of normal vectors; dolfinx Newton solver doesn't converge for highly nonlinear problems

General
1

This is a post for solutions of two issues encountered in my problem.

My differential system has two PDEs both in the second order. Basically I have two unknown functions in a 2D region in the second order. So I have four boundary conditions. Two are Dirichlet, which are imposed strongly. The other two are Neumann, which are imposed through the variational form. And this system of equations is highly nonlinear. To type them done explicitly, I fully utilized the package Unified Form Language.

The first issue is that, in my problem, when I impose the Neumann bc, there is a term to take spatial derivatives on the normal vector to the domain boundary. Initially I used the UFL function: define n = ufl.FacetNormal and take the derivative to its component, say (n[0]).dx(0). And it turned out to be zero, which is understandable since this facet normal is totally a geometrical quantity defined only on the boundary. To solve this issue, I have to define my own normal vector using UFL expressions:

x = ufl.SpatialCoordinate(mesh)
n = ufl.as_vector([x[0] / ufl.sqrt((x[0])**2 + (x[1])**2), x[1] / ufl.sqrt((x[0])**2 + (x[1])**2)])

for a circular mesh as an example. Then I can safely take the derivative for this vector.

The second issue is that, my problem is very much nonlinear. I was using the integrated Newton solver from dolfinx, say,

problem = dolfinx.fem.petsc.NonlinearProblem(F, uxuy, bcs=[bc_)

solver = dolfinx.nls.petsc.NewtonSolver(mpi4py.MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-10
solver.max_it = 10
solver.report = True

ksp = solver.krylov_solver
opts = petsc4py.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()

uxuy.interpolate(lambda x: [x[0] * ((x[0])**2 + (x[1])**2)**0.18, x[1] * ((x[0])**2 + (x[1])**2)**0.18])

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

But it turned that it was extremely difficult to find a proper initial guess to make it converge. So I resorted to custom Newton Solvers. And it gave me a significantly better experience. Though sometimes I still need to try for a valid initial.

I feel I have stated clearly the issues and the solutions so I don’t copy and paste my long code from several blocks to save some time and space. Please reply if it’s not clear enough.

2

In summary:

  1. To take the derivative to the normal vector, one shouldn’t use the normal vector defined in UFL: ufl.FacetNormal. Instead, one should defined his own normal vector using ufl.SpatialCoodinate() depending on the boundary shape of his mesh to get differentiated properly.

  2. For a highly nonlinear problem, or just a complicated problem, it is highly recommended to use Custom Newton Solvers in the tutorial by Dokken.

1 Like
3

Thanks for sharing your solutions. I would love to study your code and how you constructed the custom solver. If it’s not something confidential please feel free to share here or message me directly.

4

Hello, I didnt share it just bcz its too long and not necessary. I just copy and paste the custom newton solver from the tutorial, not by myself. I believe there is an alternative using the tutorial from multiphenicsx which I didnt choose to use.