Eikonal equation

[Karlosjusus]

Hi everyone!
I am trying to calculate the distance from each point to the upper and lower walls using the Eikonal equation.

My code:

import dolfinx
import ufl
import numpy as np
from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import mesh, fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.FunctionSpace(domain, ("CG", 1))

def walls(x):
    a = np.isclose(x[1],0)
    b = np.isclose(x[1],1)
    return np.logical_or(a,b)

fdim = domain.topology.dim -1
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
wall_f = dolfinx.fem.locate_dofs_topological(V, fdim, wall_facets)

ubc = fem.Function(V)

bc = fem.dirichletbc(ubc, wall_f)


d = fem.Function(V)
v = ufl.TestFunction(V)

F = (ufl.dot(ufl.grad(d), ufl.grad(d)) - 1) * v * ufl.dx

problem = NonlinearProblem(F, d, bcs=[bc])
solver = NewtonSolver(domain.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-8
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()

# Resolver el problema
num_iterations, converged = solver.solve(d)

But I am getting a null solution…

Does anyone have any ideas??
Thank you so much!

[dokken]

Your variatonal form is not correct, as you are missing a sqrt.

Solving the unsmoothed version of the eikonal equations are notoriously difficult. Usually, one smooths the equation by adding a Laplace operator.
See:

Euclidean distance map of mesh (distance to surface)

Hi, one option (which works in parallel out of the box) is to get the distance by solving the Eikonal equation, see Johan Hake’s answer here. However, if the boundary is complex there are some difficulties in getting the Newton solver to converge on the nonlinear problem that is the Eikonal equation. You can try to make it more robust following the discussion in the linked post. An alternative presented below is to decompose the boundary, solve simpler problems and then get the final distance as…

and
https://fenicsproject.org/qa/1446/distance-to-boundary/

Me and collaborators are working on a new mixed method/variational form for solving the Eikonal equation. I will post it here once we have a preprint out.

[Karlosjusus]

I’m trying to adding a Lapalce operator and get:

that it is not the distance to the walls…

This is my new code:

import dolfinx
import ufl
import numpy as np
from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import mesh, fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver


domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
V = fem.FunctionSpace(domain, ("CG", 1))

def walls(x):
    a = np.isclose(x[1],0)
    b = np.isclose(x[1],1)
    return np.logical_or(a,b)

fdim = domain.topology.dim -1
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
wall_f = dolfinx.fem.locate_dofs_topological(V, fdim, wall_facets)

ubc = fem.Function(V)

bc = fem.dirichletbc(ubc, wall_f)



d0 = ufl.TrialFunction(V)
v = ufl.TestFunction(V)



f = fem.Constant(domain, PETSc.ScalarType(1.0))


a = ufl.inner(ufl.grad(d0), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx


problem = fem.petsc.LinearProblem(a, L, bcs=[bc])
d0 = problem.solve()



hmax = 1e-6

d = fem.Function(V)
d.x.array[:] = d0.x.array[:]

eps = fem.Constant(domain, PETSc.ScalarType(hmax / 100))

F = ufl.sqrt(ufl.dot(ufl.grad(d), ufl.grad(d)) - f) * v * ufl.dx + eps*ufl.inner(ufl.grad(d), ufl.grad(v))*ufl.dx



problem = NonlinearProblem(F, d, bcs=[bc])
solver = NewtonSolver(domain.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-8
opts[f"{option_prefix}pc_type"] = "hypre"
opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
ksp.setFromOptions()


num_iterations, converged = solver.solve(d)

[dokken]

Please be more careful when writing up your variational form, as there is yet another typo:

@Karlosjusus

F = ufl.sqrt(ufl.dot(ufl.grad(d), ufl.grad(d)) - f) * v * ufl.dx + eps*ufl.inner(ufl.grad(d), ufl.grad(v))*ufl.dx

should be
F = ufl.sqrt(ufl.dot(ufl.grad(d), ufl.grad(d)))*v*ufl.dx - f * v * ufl.dx + eps*ufl.inner(ufl.grad(d), ufl.grad(v))*ufl.dx
which yields

[Karlosjusus]

I just tried your formula and I get my same result as above…

[dokken]

Please check that your solver actually converges.
Here is the complete code I ran:

import dolfinx
import ufl
import numpy as np
from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import mesh, fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver


domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)
V = fem.functionspace(domain, ("CG", 1))

def walls(x):
    a = np.isclose(x[1],0)
    b = np.isclose(x[1],1)
    return np.logical_or(a,b)

fdim = domain.topology.dim -1
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
wall_f = dolfinx.fem.locate_dofs_topological(V, fdim, wall_facets)

ubc = fem.Function(V)

bc = fem.dirichletbc(ubc, wall_f)



d0 = ufl.TrialFunction(V)
v = ufl.TestFunction(V)



f = fem.Constant(domain, PETSc.ScalarType(1.0))


a = ufl.inner(ufl.grad(d0), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx


problem = fem.petsc.LinearProblem(a, L, bcs=[bc])
d0 = problem.solve()



hmax = 1e-6

d = fem.Function(V)
d.x.array[:] = d0.x.array[:]

eps = fem.Constant(domain, PETSc.ScalarType(hmax / 100))

F = ufl.sqrt(ufl.dot(ufl.grad(d), ufl.grad(d)))*v*ufl.dx - f * v * ufl.dx + eps*ufl.inner(ufl.grad(d), ufl.grad(v))*ufl.dx



problem = NonlinearProblem(F, d, bcs=[bc])
solver = NewtonSolver(domain.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-8
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
num_iterations, converged = solver.solve(d)

with io.XDMFFile(domain.comm, "eikonal.xdmf","w") as writer:
    writer.write_mesh(domain)
    writer.write_function(d)