Energy depends on normalised vector field - what if magnitude goes to 0?

General
1

Dear all,
I need to solve a nonlinear problem for a vector field \mathbf{p} that arises from minimising an energy that has terms that depend on the direction \hat{\mathbf{p}}= \mathbf{p}/|\mathbf{p}|. However, the solver never finishes (but also doesn’t complain of anything). The problem is coming from the normalisation, as I get a solution if I use \mathbf{p} instead of \hat{\mathbf{p}} in the energy. I guess there are points where |{\mathbf{p}}|=0. Any ideas on how to tackle this problem?

Here’s a MWE:

import numpy as np
import ufl
import gmsh
from dolfinx import fem, nls
from dolfinx.io import gmshio
from petsc4py.PETSc import ScalarType
from mpi4py import MPI

model = gmsh.model; geom = gmsh.model.geo
mesh_comm = MPI.COMM_WORLD

def generateSphere3D(RR):
    """Create a Gmsh model of a sphere."""
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    model.add("model_acorn")
    sphere = model.occ.addSphere(0, 0, 0, RR, tag=1)
    model.occ.synchronize()
    model.add_physical_group(dim=3, tags=[sphere])    
    geom.synchronize()
    model.mesh.generate(3)
    return gmshio.model_to_mesh(model, mesh_comm, 0, gdim=3)

domain, _, _ = generateSphere3D(1)

#### define the space of funtions and the functions that we will use
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
V = fem.FunctionSpace(domain, element) # space of vector functions with 3 components
m = fem.Function(V, name='m') # 3-dimensional vector OP
dm = ufl.TrialFunction(V) # variation of the OP
phi = ufl.TestFunction(V) # test function

#### derive the weak formulation directly from the energy functional
def ET(u):
    # example energy
    normu = ufl.sqrt(ufl.dot(u,u))
    u1 = u/normu
    fT = normu**2 * (ufl.div(u1))**2 + ufl.dot(ufl.grad(normu),ufl.grad(normu))
    return fT*ufl.dx
F = ET(m)
dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy

#### define and solve the problem
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)

def solveProblem(prob):
    solver = nls.petsc.NewtonSolver(mesh_comm, prob)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6 
    solver.max_it = 100
    solver.report = True
    minit = [1.,0,0] # initial guess for the unknown function
    minit = ufl.as_vector([fem.Constant(domain, ScalarType(minit[0])),
                           fem.Constant(domain, ScalarType(minit[1])),
                           fem.Constant(domain, ScalarType(minit[2]))])
    expr = fem.Expression(minit, V.element.interpolation_points())
    m.interpolate(expr)
    solver.nonzero_initial_guess = True
    n, converged = solver.solve(m)
    assert(converged)

solveProblem(problem)
2

Please confirm this guess by computing the norm of p (or exporting p to an output format, and the visualize the result with paraview or similar). You probably don’t want to wait 100 iterations, so change that value to something small.