eikonal equation

Hi,

There is an Eikonal solver in the CBC.RANS project on
https://launchpad.net/cbc.rans.

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 = VariationalProblem(J, -F, bcs)
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://answers.launchpad.net/dolfin/+question/156790
>
> 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(Expression):
> def eval(self, values,x):
> if x[0]<0.5:
> values[0] = 2.0
> else:
> values[0] = 3.0
>
> Q = FunctionSpace(mesh,'CG',1)
> u = TrialFunction(Q)
> v = TestFunction(Q)
>
> n = RefractiveIndex()
> a = inner(grad(u), grad(u))*v*dx
> L = n**2*v*dx
>
> problem = VariationalProblem(a, L, bcs=[])
> 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://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Thanks Mikael Mortensen, that solved my question.

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,grad(y))*ds ?

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://answers.launchpad.net/dolfin/+question/156790
>
> 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,grad(y))*ds ?
>
> 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://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

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(Expression):
    def eval(self, values,x):
        if x[0]<0.5:
            values[0] = 1.0
        else:
            values[0] = 100.0

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

#k_markers = MeshFunction("uint", mesh, mesh.topology().dim() - 1)
#k_markers.set_all(0)
#light_boundary.mark(k_markers,1)

Q = FunctionSpace(mesh,'CG',1)
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(grad(du),grad(v))*dx
J = derivative(F,u,du)

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

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

Thanks Mikael Mortensen, that solved my question.