-
Notifications
You must be signed in to change notification settings - Fork 3
/
icetools_3d_demo.py
175 lines (139 loc) · 5.28 KB
/
icetools_3d_demo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
"""A 3D demo of the icetools code"""
__author__ = "Alexander H. Jarosch (research@alexj.at)"
__date__ = "2010-06-09 -- 2019-05-22"
__copyright__ = "Copyright (C) 2015 Alexander H. Jarosch"
__license__ = "GNU GPL Version 3"
__version__ = "1.2"
"""This file is part of icetools.
icetools is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
icetools is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with icetools. If not, see <http://www.gnu.org/licenses/>."""
from fenics import *
import ufl
parameters["form_compiler"]["quadrature_degree"] = 2
# define some general constants
g = -9.81 # gravitational constant
rho = 917.0 # fluid density
Aglen = 2.4e-24 # Glen flow parameter for temperate ice
nglen = 3.0 # Glen's n
glen_fact = 0.5 * Aglen**(-1.0/nglen)
# Define the body force f, i.e. gravity, driving the flow
alpha = 10.0 # inclination of glacier
alphar = alpha*pi/180
f_x0 = sin(alphar)*g*rho
f_x2 = cos(alphar)*g*rho
f = Constant((f_x0, 0, f_x2))
# generate a mesh
Lx = 200.
Ly = 200.
Lz = 100.
Nx = 20
Ny = 20
Nz = 10
mesh = BoxMesh(Point(0., 0., 0.), Point(Lx, Ly, Lz), Nx, Ny, Nz)
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two slave edges
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], Lx) and near(x[1], 0)) or
(near(x[0], 0) and near(x[1], Ly)))) and on_boundary)
def map(self, x, y):
if near(x[0], Lx) and near(x[1], Ly):
y[0] = x[0] - Lx
y[1] = x[1] - Ly
y[2] = x[2]
elif near(x[0], Lx):
y[0] = x[0] - Lx
y[1] = x[1]
y[2] = x[2]
elif near(x[1], Ly):
y[0] = x[0]
y[1] = x[1] - Ly
y[2] = x[2]
else:
y[0] = -1000
y[1] = -1000
y[2] = -1000
# Define the sub domains for the boundary conditions
class NoslipBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[2] < DOLFIN_EPS and on_boundary
# defining pbc as the periodic boundary
pbc = PeriodicBoundary()
# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH, constrained_domain=pbc)
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
# Apply a no-slip boundary condition for velocity
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, NoslipBoundary())
# Collect boundary conditions
bcs = [bc0]
# comment: to solve the linear problem: u & p needs to be Trialfunctions first
# before the variational form is defined.
nu = 8e13 # linear viscosity for for first linear guess
epsilon = sym(grad(u))
F = (2*nu*inner(epsilon, grad(v)) - div(u)*q - div(v)*p)*dx - inner(v, f)*dx
dF = derivative(F, w)
pde = NonlinearVariationalProblem(F, w, bcs, dF)
solver = NonlinearVariationalSolver(pde)
solver.parameters["symmetric"] = True
solver.parameters["newton_solver"]["maximum_iterations"] = 5
solver.parameters["newton_solver"]["linear_solver"] = 'mumps'
solver.solve()
# define the nonlinear stokes equation directly:
def visc(u):
# second invariant of strain
eps_dot = sqrt(0.5*inner(sym(grad(u)), sym(grad(u))))
nu_out = glen_fact * eps_dot**((1.0 - nglen)/nglen)
# return nu_out
return ufl.Min(nu_out, 2e15) # would introduce a viscosity limit
nu = visc(u)
epsilon = sym(grad(u))
F = (2*nu*inner(epsilon, grad(v)) - div(u)*q - div(v)*p)*dx - inner(v, f)*dx
dF = derivative(F, w)
pde = NonlinearVariationalProblem(F, w, bcs, dF)
solver = NonlinearVariationalSolver(pde)
solver.parameters["symmetric"] = True
solver.parameters["newton_solver"]["maximum_iterations"] = 80
solver.parameters["newton_solver"]["error_on_nonconvergence"] = False
solver.parameters["newton_solver"]["relaxation_parameter"] = 0.6
solver.parameters["newton_solver"]["relative_tolerance"] = 1E-4
solver.parameters["newton_solver"]["linear_solver"] = 'mumps'
solver.solve()
(u, p) = w.split(deepcopy=True)
# compare to analytical solution
H = Lz
Aterm1 = (2*Aglen)/(nglen+1)*(rho*g*sin(alphar))**nglen
ue_x0 = '%g * (pow(%g,%g+1) - pow(%g-x[2],%g+1))' % (Aterm1, H, nglen, H, nglen)
ue_x1 = '0.0'
ue_x2 = '0.0'
# Evaluate the equation on the mesh to create a 3D field of the solution
ue_E = Expression((ue_x0, ue_x1, ue_x0), degree=1)
Q = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, 'CG', 2)
u_e = interpolate(ue_E, V)
nu_out = project(nu, Q)
u_err = project(u - u_e, V)
# Save the final solution
uendfile_pvd = File("velocity_3D.pvd")
uendfile_pvd << u
u_eendfile_pvd = File("velocity_3D_analytical.pvd")
u_eendfile_pvd << u_e
u_errendfile_pvd = File("velocity_3D_err.pvd")
u_errendfile_pvd << u_err
pendfile_pvd = File("pressure_3D.pvd")
pendfile_pvd << p
nuendfile_pvd = File("nu_3D.pvd")
nuendfile_pvd << nu_out