Bug with low constants.

Question

Hi Fenics persons,

I have been used fenics for transient. When i use low value for constants for diffusive constant and low value for domain L. The grath of solution stop to work in some time. If i increase the Dirichlet conditions value it works for high value.

I tried to use python double precision numpy constant but it is doesnt work too. np.dtype('Float64') or np.dtype('d')

If i use u01 = Constant(1.0) it stop to work at around 0.2 sec
If i use u01 = Constant(10.0) it stop to work at around 5.5 sec
If i use u01 = Constant(100.0) it works well.

I am afraid about to use low constants. But i dont know if it is the problem. Is it?

You can confirm it running a simple code:

from dolfin import *
import numpy as np
# defining mesh
n=100

L=1E-6
mesh = IntervalMesh(n,0,L)



arq_saida="flux"
o1 = open(arq_saida,'w')

# definig Function space on this mesh using Lagrange polynoimals of degree 2.
V = FunctionSpace(mesh, "CG", 2)
boundaries = FacetFunction('size_t', mesh)
left   = AutoSubDomain(lambda x: near(x[0], 0.0))
right  = AutoSubDomain(lambda x: near(x[0], L))

left  .mark(boundaries, 1)
right .mark(boundaries, 2)
# definign boundary values

dsN = Measure("ds", subdomain_id=2, subdomain_data=boundaries)

def boundary0(x):
    return x[0] < DOLFIN_EPS
def boundary1(x):
    return x[0] > L - DOLFIN_EPS

dt = Constant(0.003536)
t = float(dt)
Tf = 20.
D=1.4*10**(-14)
g = Constant(0.0)

# Define boundary condition
u01 = Constant(10.0)
u02 = Constant(0.0)
bc0 = DirichletBC(V, u01, boundary0)
bc1 = DirichletBC(V, u02, boundary1)
bcs=[bc0,bc1]
u0 = interpolate(g, V)
#Variational problem at each time
u= Function(V)
c= Function(V)
v= TestFunction(V)
f= Constant(0.0)
F= (u*v*dx + dt*D*inner(grad(u), grad(v))*dx-u0*v*dx + dt*f*v*dx)
n = FacetNormal(mesh)

while t <= Tf:
    g.t = t
    solve(F == 0, u, bcs)
    flux = -dot(grad(u), n)*dsN(2)
    flux = assemble(flux)
    o1.write("%s"%t)
    o1.write(" %s\n" %flux)
    u0.assign(u)
    t += float(dt)
    print ("%s %f" %(t , flux))
    plot(u, title='Solution at t = %g' % t)
o1.close()
interactive()

How can i solve this problem?

Kind regards,
LeoCosta

Answer 1

It's not a bug (it's a convergence issue!). The problem is the scale of your solution (the solution is fine considering the default tolerances of the newton solver). You need to rescale your problem (changing the units of the lengths and physical parameters, for example) or change the tolerances of the nonlinear solver. Try this:

solve(F == 0, u, bcs, solver_parameters={"newton_solver": {"relative_tolerance": 1e-18, "absolute_tolerance": 1e-20}})