# Eikonal equation

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)

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!

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:

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.

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))

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))

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)

``````

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

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

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

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))

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))

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)
``````