Simulation stops working if one constant is changed from 1E-15 to 1E-16

Question

For my chemistry bachelors degree I wrote a simulation of drying polymer films with FEniCS. It all seems to work fine till I change the permeability from 1E-15 to 1E-16 it stops working, after the first time-step the Newton and the linear solver always finished in 0 iterations. The problem also occurs if I approach 1E-15 asymptotic (lets say 0.999999E-15).

Cause it's bound to a specific value I think it's a problem of some setting but changing the tolerance of the solver so far has no effect.

Here is the code, the problematic constant is kappa0

from dolfin import *
import numpy, sys, time

#variables to create mesh and function space
T = 0.001
dt_out = T/10.
nx = 1.4E-5
nz = 1.4E-5
cnx = 40 #cell devision x
cnz = 40 #cell devision z
degree = 1 #degree of function
kappa0 = 0.999999E-15 # permeability in m2 
myu = 1E-3 # viscosity of the fluid in Pas
K = 1.0E8 #compression modul in Pa
u0 = 0.01 #start condition

# Create mesh and define function space
mesh = RectangleMesh(Point(0, 0), Point(nx, nz), cnx, cnz, "right")
V = FunctionSpace(mesh, 'Lagrange', degree)

# define boundary conditions: Dirichlet Boundary on top
def GammaD(x, on_boundary):
    return on_boundary and near(x[1], nz) 

# Choice of nonlinear coefficient
m = 2./3.

def q(u):
    return kappa0*u**m

# Boundary condition
g = 0.0025
bc = DirichletBC(V, g, GammaD)

# Initial condition
u = Function(V) 
u_1 = Function(V)
u = interpolate(Expression("u0",u0=u0),V)
u_1 = interpolate(Expression("u0",u0=u0),V)

dt = 0.0001      # time step

f = Expression("0")

# Define variational problem
v = TestFunction(V)
du = TrialFunction(V)
F = dt*inner((q(u)*K)/myu*nabla_grad(u), nabla_grad(v))*dx+u*v*dx-(u_1 + dt*f)*v*dx
mass = (u-g)*dx
uni = interpolate(Expression("1"),V)
volume = uni*dx


J = derivative(F, u, du)

problem = NonlinearVariationalProblem(F, u, bc, J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-21
prm['newton_solver']['relative_tolerance'] = 1E-20
prm['newton_solver']['maximum_iterations'] = 100

# Compute solution
t = dt
t_out = dt_out
viz = plot(u,
           scale=0.0,
           title='Drying Polymere Film',
           axes=False
           )
viz.elevate(65)
viz.set_min_max(0,u0)
viz.plot(u_1)
interactive()
while t <= T:
    solver.solve()
    M = assemble(mass)
    V = assemble(volume)

    print 'time =', t, 'M =', M, 'V =', V

    # plot
    if (t>=t_out) :
      viz.plot(u)
      interactive()
      t_out += dt_out

    t += dt
    u_1.assign(u)

Comments

Just a comment about the solver tolerance: Your line prm['newton_solver']['relative_tolerance'] = 1E-20 does not have any function i guess, it is unreachable since it is way below DOLFIN_EPS. So obviously if you change this tolerance it has no effect.

Answer 1

You should really nondimensionalize your equations before solving. Directly working with permeabilities around machine precision makes no sense... Even in the cases where the solver actually computes a solution I am highly suspicious about the accuracy, since it is likely that truncation errors elsewhere dominate.

After you have obtained a solution, you can easily recover dimensions for reporting your results.