# eikonal equation

Hi,

I'd like to solve the pde |grad(u)|**2 = n**2. To start , I tried the following example but it does not work. What would be the correct variational form for this?

Thanks. Chaffra

mesh = UnitInterval(51, )

class RefractiveIndex

def eval(self, values,x):

if x[0]<0.5:

else:

Q = FunctionSpace(

u = TrialFunction(Q)

v = TestFunction(Q)

n = RefractiveIndex()

a = inner(grad(u), grad(u))*v*dx

L = n**2*v*dx

problem = VariationalProb

problem.solve()

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Mikael Mortensen

- Solved:

- Last query:

- Last reply:

## This question was reopened

Hi,

There is an Eikonal solver in the CBC.RANS project on

https:/

Basically, you need to declare a nonlinear problem. Try something like

y = Function(Q)

eps = Constant(0.01)

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

grad(v))*dx

J = derivative(F, y, u)

bc = DirichletBC(Q, Constant(0), DomainBoundary())

problem = VariationalProb

y = problem.solve()

where the last term in F is for stabilization. There are other possibilities

for

stabilizing as well, but this is the one that works for me.

Hope it works. Look up the nonlinear demos if it doesn't. (Never tried it in

1D before)

Best regards

Mikael

On 9 May 2011 21:56, Chaffra <email address hidden> wrote:

> New question #156790 on DOLFIN:

> https:/

>

> Hi,

>

> I'd like to solve the pde |grad(u)|**2 = n**2. To start , I tried the

> following example but it does not work. What would be the correct

> variational form for this?

>

> Thanks. Chaffra

>

>

> mesh = UnitInterval(51, )

>

> class RefractiveIndex

> def eval(self, values,x):

> if x[0]<0.5:

> values[0] = 2.0

> else:

> values[0] = 3.0

>

> Q = FunctionSpace(

> u = TrialFunction(Q)

> v = TestFunction(Q)

>

> n = RefractiveIndex()

> a = inner(grad(u), grad(u))*v*dx

> L = n**2*v*dx

>

> problem = VariationalProb

> problem.solve()

>

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

>

> _______

> Mailing list: https:/

> Post to : <email address hidden>

> Unsubscribe : https:/

> More help : https:/

>

Chaffra Affouda (chaffra) said : | #2 |

Thanks Mikael Mortensen, that solved my question.

Chaffra Affouda (chaffra) said : | #3 |

I do have a related question:

How do I specify the following boudary condition for my problem. I have a initial propagation vector let's say ko = [1,0,0] at x[0] ==0 and k = (1/n)*grad(y). Do I write something like:

F = inner(grad(y), grad(y))*v*dx - n**2*v*dx + eps*inner(grad(u) + inner(ko,

I am not sure how to write this kind of boundaries...

Thanks,

Chaffra

Do you have Neuman on Eikonal's equation? Is that even possible? Anyway,

I'm not sure how to do that since there are no ds terms naturally occuring

in this equation (except the stabilization). Perhaps if you rewrite using

the chain rule? Sorry I can't be of more assistanse.

Best regards

Mikael

On 10 May 2011 14:41, Chaffra <email address hidden> wrote:

> Question #156790 on DOLFIN changed:

> https:/

>

> Status: Solved => Open

>

> Chaffra is still having a problem:

> I do have a related question:

>

> How do I specify the following boudary condition for my problem. I have

> a initial propagation vector let's say ko = [1,0,0] at x[0] ==0 and k =

> (1/n)*grad(y). Do I write something like:

>

> F = inner(grad(y), grad(y))*v*dx - n**2*v*dx + eps*inner(grad(u) +

> inner(ko,

>

> I am not sure how to write this kind of boundaries...

>

> Thanks,

> Chaffra

>

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

>

> _______

> Mailing list: https:/

> Post to : <email address hidden>

> Unsubscribe : https:/

> More help : https:/

>

Chaffra Affouda (chaffra) said : | #5 |

Thanks for your answers Mikael. I would be interested in knowing more about the chain rule implementation if you have more info/resources on how to do it. Also using you code above I can solve the eikonal equation for my case but the solution is zero everywhere. I assume this has to do with the boundary conditions but is there a working demo example I could work from?

Many thanks,

Chaffra

--new code--

mesh = UnitInterval(51, )

class RefractiveIndex

def eval(self, values,x):

if x[0]<0.5:

else:

def light_exposed_

return near(x[0],0.0) and on_boundary

light_boundary = AutoSubDomain(

#k_markers = MeshFunction(

#k_markers.

#light_

Q = FunctionSpace(

u = Function(Q)

du = TrialFunction(Q)

v = TestFunction(Q)

eps = Constant(0.01)

n = RefractiveIndex()

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

J = derivative(F,u,du)

bc = DirichletBC(

problem = VariationalProb

problem.solve(u)

plot(u)

Hi Chaffra,

There were a few errors in the posted script. The following works for me

mesh = UnitInterval(51, )

class RefractiveIndex

def eval(self, values,x):

if x[0]<0.5:

else:

def light_exposed_

return near(x[0],0.0) and on_boundary

light_boundary = AutoSubDomain(

#k_markers = MeshFunction(

#k_markers.

#light_

Q = FunctionSpace(

u = Function(Q)

du = TrialFunction(Q)

v = TestFunction(Q)

eps = Constant(0.01)

n = RefractiveIndex()

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(

problem = VariationalProb

problem.solve(u)

plot(u)

Best regards

Mikael

On 10 May 2011 15:50, Chaffra <email address hidden> wrote:

> Question #156790 on DOLFIN changed:

> https:/

>

> Status: Answered => Open

>

> Chaffra is still having a problem:

> Thanks for your answers Mikael. I would be interested in knowing more

> about the chain rule implementation if you have more info/resources on

> how to do it. Also using you code above I can solve the eikonal equation

> for my case but the solution is zero everywhere. I assume this has to do

> with the boundary conditions but is there a working demo example I could

> work from?

>

> Many thanks,

> Chaffra

>

> --new code--

>

> mesh = UnitInterval(51, )

>

> class RefractiveIndex

> def eval(self, values,x):

> if x[0]<0.5:

> values[0] = 1.0

> else:

> values[0] = 100.0

>

> def light_exposed_

> return near(x[0],0.0) and on_boundary

> light_boundary = AutoSubDomain(

>

> #k_markers = MeshFunction(

> #k_markers.

> #light_

>

> Q = FunctionSpace(

> u = Function(Q)

> du = TrialFunction(Q)

> v = TestFunction(Q)

> eps = Constant(0.01)

>

> n = RefractiveIndex()

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

> J = derivative(F,u,du)

>

> bc = DirichletBC(

>

> problem = VariationalProb

> problem.solve(u)

> plot(u)

>

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

>

> _______

> Mailing list: https:/

> Post to : <email address hidden>

> Unsubscribe : https:/

> More help : https:/

>

Chaffra Affouda (chaffra) said : | #7 |

Thanks Mikael Mortensen, that solved my question.