Hello! I am trying to run the sample code in http://www.diva-portal.org/smash/get/diva2:1368307/FULLTEXT01.pdf and it is giving me some errors. Towards the end of the code in line 222 it says that the object does not have the property to do extrapolation. I am wondering if the methods above have an error in them. If anyone can help that would be much appreciated.
from fenics import *
from mshr import *
import numpy as np
# Flow Parameters (km, yr, Mpa)
n_i = 3.0 # glen's flow law exponent
A = 100.0 # glen's flow law parameter
rho = 9.037e-13 # density
g = 9.760e+12 # gravitational acceleration'
# Horizontal Extent
Lx = 1000 # (km)
# Number of Nodes
Nx = 101 # nodes along x
Nz = 11 # nodes along z
# Define EISMINT Domain
def EISMINT():
x = np.linspace(0, Lx, Nx)
tops = 0.1*np.ones(len(x))
bots = np.zeros(len(x))
return x, bots, tops
xvec, z_bot, z_top = EISMINT()
# Map the Geometry Enclosed by bottomSurf and topSurf to Structured Mesh
def createMesh(xvec, bottomsurf, topsurf, Nx, Nz):
hl = Nx - 1 # number of horizontal layers
vl = Nz - 1 # number of vertical layers
# generate mesh on unitsquare
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 1.0), hl, vl)
# Extract x and y coordinates from mesh
x = mesh.coordinates()[:, 0]
y = mesh.coordinates()[:, 1]
x0 = min(xvec)
x1 = max(xvec)
# Map coordinates on unit square to the computational domain
xnew = x0 + x*(x1 - x0)
zs = np.interp(xnew, xvec, topsurf)
zb = np.interp(xnew, xvec, bottomsurf)
ynew = zb + y*(zs - zb)
xy_new_coor = np.array([xnew, ynew]).transpose()
mesh.coordinates()[:] = xy_new_coor
return mesh
# Class describing the bottom boundary of the domain
class bottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool(on_boundary and near(x[1], 0.0) )
# Class describing the left boundary of the domain
class leftBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool(on_boundary and near(x[0], 0.0) )
# Class describing the right boundary of the domain
class rightBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool(on_boundary and near(x[0], Lx) )
# Function extracting values along surface z_top
def boundaryValues(u, xvec, z_top):
# Define projection space
V = FunctionSpace(u.function_space().mesh(), "Lagrange", 2)
# Project velocity onto V
ux = project(u[0], V)
uz = project(u[1], V)
# Allow extrapolation, in case coordinates are slightly outside the domain
ux.set_allow_extrapolation(True)
uz.set_allow_extrapolation(True)
# Preallocate array for velocity components
ux_b = np.zeros(len(xvec))
uz_b = np.zeros(len(xvec))
# Extract velocity along z_top
for i in range(len(xvec)):
ux_b[i] = ux(xvec[i], z_top[i])
uz_b[i] = uz(xvec[i], z_top[i])
return ux_b, uz_b
# Solve Stoke's equation on doman described by mesh
def solveStokes(mesh):
# Define Taylor−Hood element
V = VectorElement("Lagrange", triangle, 2)
Q = FiniteElement("Lagrange", triangle, 1)
TH = V*Q
W = FunctionSpace(mesh, TH)
# Test function
y = TestFunction(W)
(v, q) = split(y)
# Trial function
w = TrialFunction(W)
(u, p) = split(w)
# Initial guess for nonlinear solver
w0 = Function(W)
(u0, p0) = split(w0)
# Create boundary objects
bottom = bottomBoundary()
left = leftBoundary()
right = rightBoundary()
bc_bot = DirichletBC(W.sub(0), ((0.0, 0.0)), bottom)
bc_left = DirichletBC(W.sub(0), ((0.0, 0.0)), left)
bc_right = DirichletBC(W.sub(0), ((0.0, 0.0)), right)
bcs = [bc_bot, bc_left, bc_right]
# Glen's flow law
def viscosity(u):
strain_rate = sym(grad(u))
eff_strain_rate = 0.5*tr(strain_rate*strain_rate)
return 0.5*A**(-1.0/n_i)*(eff_strain_rate + 1e-10)**((1.0-n_i)/(2.0*n_i))
# Force term
f = Constant((0.0, -rho*g))
# Define bilinear form
a = ( 2*viscosity(u0)*inner(sym(grad(u)), grad(v)) - div(v)*p - div(u)*q ) * dx
# Define linear form
L = inner(f, v) * dx
# Define function for holding the solution
w = Function(W)
(u, p) = split(w)
(ux, u_z) = split(u)
# Set up fixed point iteration algorithm for solving nonlinear problem
i = 0
max_iter = 100
tol = 1e-6
for i in range(max_iter):
# Define linearized variational problem of nonlinear PDE
pde = LinearVariationalProblem(a, L, w, bcs)
solver = LinearVariationalSolver(pde)
# Solve variational problem
solver.solve()
# Estimate the error
diff = w.vector().get_local() - w0.vector().get_local()
if (np.linalg.norm(diff, ord=np.Inf) < tol):
break
# Update the initial guess
w0.assign(w)
if (i+1 >= max_iter):
print("Warning: maximum number of iterations reached")
else:
print("Nonlinear solver finished in", i+1, "iterations")
# Solve surface equation
def solveSurface(xvec, z_top, u, ddx, dt, eps):
# Define 1D mesh
interval = IntervalMesh(Nx-1, 0, Lx)
# Define 1D function space
V = FunctionSpace(interval, "Lagrange", 1)
# Extract velocity on boundary
ux_b, uz_b = boundaryValues(u, xvec, z_top)
ux = Function(V)
uz = Function(V)
ux.vector().set_local(ux_b)
uz.vector().set_local(uz_b)
# Define trial function on V
v = TrialFunction(V)
# Define test function on V
s = TestFunction(V)
# Function for holding the old solution
s0 = Function(V)
s0.vector().set_local(np.copy(z_top))
# Accumulation function (km/yr)
a_s_expr = "fmax(0, fmin(0.0005, s*( R - fabs(x[0] - 0.5*L) ) ))"
a_s = Expression(a_s_expr, s=1e-5, R=200, L=Lx, degree=1)
# Define bilinear and linear forms
a = Constant(1.0/dt)*s*v*dx + ux*s.dx(0)*v*dx + Constant(eps)*s.dx(0)*v.dx(0)*dx
L = Constant(1.0/dt)*s0*v*dx + (a_s + uz)*v*dx
# Function for holding the new solution
h = Function(V)
# Define linear variational problem
pde = LinearVariationalProblem(a, L, h, [])
solver = LinearVariationalSolver(pde)
# Solve variational problem
solver.solve()
return h.vector().get_local()
t = 0 # initial time (yr)
T = 2000 # simulation time (yr)
dt = 25 # time stepping (yrs)
eps = 0.2 # artificial viscosity
ddx = Lx/(Nx-1) # horizontal grid resolution (km)
while(t < T):
# Generate mesh on domain bounded by curves z bot and z top
mesh = createMesh(xvec, z_bot, z_top, Nx, Nz)
# Solve stokes
w = solveStokes(mesh)
w.set_allow_extrapolation(True)
(u, p) = w.split()
# Solve surface equation
ux_b, uz_b, h = solveSurface(z_top, u, ddx, dt, eps)
z_top = h.vector().get_local()
# Update time
t += float(dt)