Newton solver not converging because of reacion term

Question

Hello, I am solving an advection-reaction equation:

F = mu*inner(nabla_grad(u), nabla_grad(v))*dx +R*inner(u,v)*dx
solve(F == 0, u, bcs)

With Dirichlet boundary conditions (u=1 on the left, u=0 on the right). U is a CG 1 function.

With the following reactive term, everything works:

R=Expression('cos((x[0]-0.5)*pi/0.2)*cos((x[1]-0.5)*pi/0.2)\
*(pow((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5),0.5)<=0.1)+0.01')

On the other hand, with this (similarly shaped, and smoother) reactive term, the newton solver does not converge.

R=Expression('exp(1.0-1.0/(1.0-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.1*0.1)))\
*(pow((x[0]-0.5)*(x[0]-0.5)/+(x[1]-0.5)*(x[1]-0.5),0.5)<=0.1)+1.0')

I noticed that, for coarse meshes (less than 100x100), if I add

R=interpolate(R,FunctionSpace(M,'DG',0))

Then the solver will converge. For a bit finer meshes the same is true if I interpolate R onto a CG1 space, but for very fine meshes the solver will never converge.

The first reactive term always works without any need of interpolating it.

Any idea?
Thanks.

Comments

If you have problems with quadrature of the expression, you can increase order of underlying element. In your case, you can do, for instance

R = Expression(..., degree=3)

---

Thanks for the reply Jan. It does not solve the problem though. I am quite lost here about what could be the issue. Is that solution working for you? Here's the complete script...

from dolfin import *

M=UnitSquareMesh(200,200)
V=FunctionSpace(M,'CG',1)

u=Function(V)
v=TestFunction(V)

class Left(SubDomain):
    def inside(self, x, on_boundary):
      return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
      return near(x[0], 1.0)

left = Left()
right = Right()

boundaries = FacetFunction("uint", M)
left.mark(boundaries, 1)
right.mark(boundaries, 3)

ds = Measure("ds")[boundaries]

bc1=DirichletBC(V,"0.0",boundaries,3)
bc0=DirichletBC(V,"1.0",boundaries,1)
bcs=[bc0,bc1]
mu=0.1

n=FacetNormal(M)

R=Expression('7.0*exp(1.0-1.0/(1.0-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.1*0.1)))\
*((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5)<=0.1*0.1)+1.0',degree=3)
#R=Expression('100*cos((x[0]-0.5)*pi/0.2)*cos((x[1]-0.5)*pi/0.2)\
#*(pow((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5),0.5)<=0.1)+0.01',degree=3)

F = mu*inner(nabla_grad(u), nabla_grad(v))*dx +R*inner(u,v)*dx

solve(F == 0, u, bcs,solver_parameters={"newton_solver": {"maximum_iterations": 3,'absolute_tolerance':1E-15}})

plot(u,title='u')
interactive()

edit: correct R

---

Works for me on FEniCS 1.3.0 (with PETSc backend) even with Expression(..., degree=1).

---

Sorry Jan, my bad, I commented the expression that is causing the problem for me, and left uncommented the one that works... Did you try with both the expressions for R? Sorry again

Answer 1

R takes infinite values

>>> R(0.5, 0.5)
nan

Comments

Oh god... I don't know why the hell copypasting I added a '/'... this is the correct form of R... Sorry again Jan

R=Expression('exp(1.0-1.0/(1.0-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.1*0.1)))\
  *((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5)<=0.1*0.1)+1.0',degree=1)

The one I wrote previously one had a / in the middle of the first line that came out of nowhere... This one is correctly valued (R(0.5,0.5)=2) and still the solver does not converge... Thanks for your patience with my stupid typos

I edited the previous comment with the correct R so now is ready to copypaste...

---

Still, it has some infinities

for c in cells(M):
    temp = R(c.midpoint())
    if temp != temp:
        print c.midpoint(), temp

---

Ok, thanks a lot for helping in the debug. It's still wierd IMO because that function is bounded .. I'll try to figure out why this is happening and post here any good result.

---

Divisor

(1.0-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))

vanishes on a set ${|x-(0.5, 0.5)|=1}.$

---

Here's how I solved it:

import sys
R=Expression('100.*exp(1.0-1.0/max((1.-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.01)),eps))\
      *((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5)<=0.01)+1.0',eps=sys.float_info.epsilon)

The problem was that

1.0/((1.-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.01))

Should tend to infinity close to the border circle centered in 0.5,0.5, but in some case it returned nan, and using

1.0/max((1.-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/(0.01)),eps))

solved the problem.

Thanks a lot jan