solving Eikonal equation

Asked by Pary Parsian on 2013-05-01

I am a new Fenics user, and I want to solve eikonal equation, but I don't know what is is the problem with my code.
Thanks,
Pary

from dolfin import *
import numpy
import sys

mesh = UnitInterval(51, )

class RefractiveIndex(Expression):
   def eval(self, values,x):
       if x[0]<0.5:
           values[0] = 1.0
       else:
           values[0] = 1.0

def light_exposed_area(x, on_boundary):
   return near(x[0],0.0) and on_boundary
light_boundary = AutoSubDomain(light_exposed_area)

V = FunctionSpace(mesh,'CG',1)
u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
eps = Constant(0.01)
n = RefractiveIndex()

F1 = inner(grad(u), grad(v))*dx - n**2*v*dx
a1, L1 = lhs(F1), rhs(F1)
A1, b1 = assemble_system(a1, L1)

for bc in bcs: bc.apply(A1, b1)
solve(A1, u.vector(), b1)

F = inner(grad(u), grad(u))*v*dx - n**2*v*dx + eps*inner(grad(u),grad(v))*dx
J = derivative(F,u,du)

bc = DirichletBC(V,Constant(0.0),light_boundary)

problem = NonlinearVariationalProblem(-F, u, bcs=[bc,], J=J))
solver = NonlinearVariationalSolver(problem)
solver.solve()
plot(u)

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
2013-05-08
Last query:
2013-05-08
Last reply:
2013-05-05
Jan Blechta (blechta) said : #1

On Wed, 01 May 2013 20:16:21 -0000
Pary Parsian <email address hidden> wrote:
> New question #227969 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227969
>
> I am a new Fenics user, and I want to solve eikonal equation, but I
> don't know what is is the problem with my code: from dolfin import *
> import numpy
> import sys
>
> mesh = UnitInterval(51, )
>
> class RefractiveIndex(Expression):
> def eval(self, values,x):
> if x[0]<0.5:
> values[0] = 1.0
> else:
> values[0] = 1.0
>
> def light_exposed_area(x, on_boundary):
> return near(x[0],0.0) and on_boundary
> light_boundary = AutoSubDomain(light_exposed_area)
>
> V = FunctionSpace(mesh,'CG',1)
> u = Function(V)
> du = TrialFunction(V)
> v = TestFunction(V)
> eps = Constant(0.01)
> n = RefractiveIndex()
>
>
> F1 = inner(grad(u), grad(v))*dx - n**2*v*dx

This is only rank 1 form. It does not contain trial function so
lhs(F1) is empty. Use u = TrialFunction(V)

> a1, L1 = lhs(F1), rhs(F1)
> A1, b1 = assemble_system(a1, L1)
>
> for bc in bcs: bc.apply(A1, b1)

You did not define bcs.

> solve(A1, u.vector(), b1)
>
>
> F = inner(grad(u), grad(u))*v*dx - n**2*v*dx +
> eps*inner(grad(u),grad(v))*dx J = derivative(F,u,du)
>
> bc = DirichletBC(V,Constant(0.0),light_boundary)
>
> problem = NonlinearVariationalProblem(-F, u, bcs=[bc,], J=J))

You provide Jacobian of F but tell DOLFIN that it is Jacobian of -F. Do
it consistently.

> solver = NonlinearVariationalSolver(problem)
> solver.solve()
> plot(u)
>

Pary Parsian (pary-karimi) said : #2

Oh, I made some mistakes! I have boundary condition just in one point (0,1)!
My code is not converging even thought I am using initial guess!
Thanks

from dolfin import *
import numpy
import sys

# Define parameters
nx = 4
ny = 4

# Create mesh and define function space
mesh=Rectangle(-1, -1, 1, 1, nx, ny,"left/right")
#Definig boundary condition
V = FunctionSpace(mesh,'CG',1)
u0 = Constant(0.0)
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14 # tolerance for coordinate comparisons
        return on_boundary and \
               (abs(x[0]) < tol and abs(x[1] - 1) < tol)

u0_boundary = DirichletBoundary()
bc = DirichletBC(V, u0, u0_boundary)

u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
eps = Constant(0.01)

f = Constant(1.0)
# Initialization problem to get good initial guess for nonlinear problem:
F1 = inner(grad(du), grad(v))*dx - f*v*dx
a1, L1 = lhs(F1), rhs(F1)
A1, b1 = assemble_system(a1, L1)
# Apply boundary conditions: DirichletBC = 0 on the boundary
bc.apply(A1, b1)
solve(A1, u.vector(), b1)

F = inner(grad(u), grad(u))*v*dx - f*v*dx + eps*inner(grad(u),grad(v))*dx
J = derivative(F,u,du)

problem = NonlinearVariationalProblem(-F, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
plot(u)

Best Johan Hake (johan-hake) said : #3

Try:

from dolfin import *

nx = 10
ny = 10

# Create mesh and define function space
mesh=RectangleMesh(-1, -1, 1, 1, nx, ny,"left/right")
#Definig boundary condition
V = FunctionSpace(mesh,'CG',1)
u0 = Constant(0.0)
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], -1.0) and near(x[1], 1.0)

u0_boundary = DirichletBoundary()
bc = DirichletBC(V, u0, u0_boundary, "pointwise")

u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
eps = Constant(0.1)

f = Constant(1.0)

# Initialization problem to get good initial guess for nonlinear problem:
F1 = inner(grad(du), grad(v))*dx - f*v*dx
a1, L1 = lhs(F1), rhs(F1)
solve(a1==L1, u, bc)

plot(u, interactive=True)

F = sqrt(inner(grad(u), grad(u)))*v*dx - f*v*dx +
eps*inner(grad(u),grad(v))*dx
solve(F==0, u, bc)

plot(u, interactive=True)

On 05/04/2013 10:16 PM, Pary Parsian wrote:
> Question #227969 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/227969
>
> Status: Answered => Open
>
> Pary Parsian is still having a problem:
> Oh, I made some mistakes! I have boundary condition just in one point (0,1)!
> My code is not converging even thought I am using initial guess!
> Thanks
>
> from dolfin import *
> import numpy
> import sys
>
>
> # Define parameters
> nx = 4
> ny = 4
>
> # Create mesh and define function space
> mesh=Rectangle(-1, -1, 1, 1, nx, ny,"left/right")
> #Definig boundary condition
> V = FunctionSpace(mesh,'CG',1)
> u0 = Constant(0.0)
> class DirichletBoundary(SubDomain):
> def inside(self, x, on_boundary):
> tol = 1E-14 # tolerance for coordinate comparisons
> return on_boundary and \
> (abs(x[0]) < tol and abs(x[1] - 1) < tol)
>
> u0_boundary = DirichletBoundary()
> bc = DirichletBC(V, u0, u0_boundary)
>
>
> u = Function(V)
> du = TrialFunction(V)
> v = TestFunction(V)
> eps = Constant(0.01)
>
> f = Constant(1.0)
> # Initialization problem to get good initial guess for nonlinear problem:
> F1 = inner(grad(du), grad(v))*dx - f*v*dx
> a1, L1 = lhs(F1), rhs(F1)
> A1, b1 = assemble_system(a1, L1)
> # Apply boundary conditions: DirichletBC = 0 on the boundary
> bc.apply(A1, b1)
> solve(A1, u.vector(), b1)
>
>
> F = inner(grad(u), grad(u))*v*dx - f*v*dx + eps*inner(grad(u),grad(v))*dx
> J = derivative(F,u,du)
>
> problem = NonlinearVariationalProblem(-F, u, bc, J=J)
> solver = NonlinearVariationalSolver(problem)
> solver.solve()
> plot(u)
>

Pary Parsian (pary-karimi) said : #4

Thanks Johan Hake, that solved my question.

jfpxtal (jan-fp) said : #5

Dear all,
I used the implementation described here to solve the eikonal equation and it works just fine.
Does any of you has any (citable) reference to a paper in which this method is discribed (and maybe analyzed)?

The only reference I found is which is kind of close:
Differential equation-based wall distance computation for
DES and RANS by P.G. Tucker, Journal of Computational Physics 190 (2003) 229–248.

Thanks and best regards,
Jan