Spatial Derivatives in non linear coefficient

variational formulation
8

Take a look at some tips and tricks with Newton Solvers.

There are a few things to consider here:

  • You should carefully choose your initial guess
  • The exponent H**5 is horrific… think about it in the context of your boundary condition that H = 100 on the left and right of the domain, that value is going to be 10^{10}. I couldn’t get this to work so added the parameter H**n where 1 < n < 2 seems to work.
  • Further to the previous point, the gradient changes (\nabla H)^2 could easily blow out numerical precision.

Perhaps you could consider some kind of non-dimensionalisation?

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

#Constants
A = 10e-16 # temperature constant
rho = 911 # density
g = 9.81 # gravity

# Horizontal range
Lx = 1000 # 
# nodes
Nx = 51 # Nodes along x

interval = IntervalMesh(Nx-1, 0, Lx)

# function space
V = FunctionSpace(interval, "Lagrange", 1)

#Here we set the boundary condition
H_D = Expression('100', degree=0)

#Here we set up our boundary definition
def boundary(x, on_boundary):
   return on_boundary

#Here we apply our boundary definition, Dirichilet because we're specifying value
bc = DirichletBC(V, H_D, boundary)


B=(2/5)*A*(rho*g)**3
H = Function(V)
H_x= grad(H)

# Some non-trivial intial guess
H.interpolate(Expression("100.0 + x[0]", degree=1))

v = TestFunction(V)

n = Constant(1.0)
F = (B*H**n)*dot(H_x, H_x)*dot(H_x, grad(v))*dx

# Setup the problem and solve
solve(F == 0, H, bc)

plot(H)
plt.show()
1 Like