-
Notifications
You must be signed in to change notification settings - Fork 0
/
minimal_model2D_monolithic.py
147 lines (119 loc) · 4.84 KB
/
minimal_model2D_monolithic.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
# TODO
# run_solver(dim, N, Pu, Pp, dt)
#if __name__ == '__main__':
# elasticity 3D
# quadliteral mesh
"""
2D Minimal model (of 1D Terzaghi), hydromechanical (HM), staggered scheme (stress-split)
spatial FE discretization for both H and M (Taylor Hood)
temporal Euler-backward discretization
incompressible fluid, incompressible solid (but linear elastic bulk)
"""
from __future__ import print_function
import fenics as fe
import numpy as np
import time
import matplotlib.pyplot as plt
import minimal_model_parameters as mmp
import stress_and_strain as sas
# DECLARATION
model=mmp.MMP()
ZeroScalar = fe.Constant((0))
ZeroVector = fe.Constant((0,0))
p_ref, _, _, k, mu = model.get_physical_parameters()
Nx, dt, dt_prog, Nt, _, _, _ = model.get_fem_parameters()
Length, Width, Height, K, Lame1, Lame2, k_mu, cc = model.get_dependent_parameters()
p_ic, p_bc, p_load = model.get_icbc()
material=sas.SAS(Lame1, Lame2, K) # stress and strain
dT=fe.Constant(dt) # make time step mutable
fe.set_log_level(30) # control info/warning/error messages
vtkfile_p = fe.File('results/mini_mono_pressure.pvd')
vtkfile_u = fe.File('results/mini_mono_displacement.pvd')
vtkfile_s_total = fe.File('results/mini_mono_totalstress.pvd') # xdmf for multiple fields
## MESH (simplex elements in 2D=triangles)
mesh = fe.UnitSquareMesh(Nx, Nx)
#fe.plot(mesh)
#mesh = fe.RectangleMesh.create([fe.Point(0, 0), fe.Point(Width, Length)], [Nx,Ny], fe.CellType.Type.quadrilateral)
Pp = fe.FiniteElement('P', fe.triangle, 1)
Pu = fe.VectorElement('P', fe.triangle, 2)
element = fe.MixedElement([Pp, Pu])
V = fe.FunctionSpace(mesh, element)
Vsigma = fe.TensorFunctionSpace(mesh, "P", 1)
pu_ = fe.Function(V, name="pressure_displacement") # function solved for and written to file
sigma_ = fe.Function(Vsigma, name="total_stress") # function solved for and written to file
# INITIAL CONDITIONS: undeformed, at rest, same pressure everywhere
pu_ic = fe.Expression(
(
p_ic, # p
"0.0","0.0", # (ux, uy, uz)
), degree = 2)
pu_n = fe.interpolate(pu_ic, V) # current value in time-stepping
pu_.assign(fe.interpolate(pu_ic, V)) # previous value in time-stepping
# BOUNDARY CONDITIONS
# DirichletBC, assign geometry via functions
tol = 1E-14
def top(x, on_boundary):
return on_boundary and fe.near(x[1], Height, tol)
bc_top = fe.DirichletBC(V.sub(0), p_bc, top) # drainage on top
def bottom(x, on_boundary):
return on_boundary and fe.near(x[1], 0.0, tol)
bc_bottom = fe.DirichletBC(V.sub(1), ZeroVector, bottom) # fixed bottom
def sides(x, on_boundary):
#return True
return on_boundary and ( fe.near(x[0], 0.0, tol) or fe.near(x[0], Length, tol) )
bc_sides_x = fe.DirichletBC(V.sub(1).sub(0), ZeroScalar, sides) # rollers on side, i.e. fix in x-direction
bc = [bc_top, bc_bottom, bc_sides_x]
# Neumann BC, assign geometry via subdomains to have ds accessible in variational problem
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
topSD = fe.AutoSubDomain(lambda x: fe.near(x[1], Height, tol))
topSD.mark(boundaries, 1) # accessable via ds(1)
ds = fe.ds(subdomain_data=boundaries)
traction_topM = fe.Constant((0, -p_load))
# FEM SYSTEM
# Define variational HM problem a(v,p)=L(v)
vp, vu = fe.TestFunctions(V)
pu = fe.TrialFunction(V)
p, u = fe.split(pu)
#vp, vu =fe.split(vpvu)
p_n, u_n = fe.split(pu_n)
sigma_.assign(fe.project(material.sigma(p_n, u_n), Vsigma))
Fdx = ( (fe.div(u)-fe.div(u_n))*vp
+ dT*k_mu*fe.dot(fe.grad(p), fe.grad(vp)) # Euler backward
+ fe.inner(material.sigma(p,u), material.epsilon(vu)) )*fe.dx
Fds = - fe.dot(vu, traction_topM)*ds(1)
F=Fdx+Fds
a, L=fe.lhs(F), fe.rhs(F)
# no non-zero Neumann BC for H, i.e. no prescribed in-outflows (only flow via DirichletBC possible)
# no sources (H)
# no body forces (M)
# TIME-STEPPING
tic = time.perf_counter()
t = 0.0
p_, u_ = pu_.split()
vtkfile_p << (p_, t)
#vtkfile_u << (u_, t)
#vtkfile_s_total << (sigma_, t)
y_mono=np.linspace(tol, Height-tol, Nx+1)
p_mono=np.zeros((Nt, Nx+1))
points=[ (0.5*Length, y_) for y_ in y_mono ]
for n in range(Nt): # time steps
t += dt
print(n+1,".step t=",t)
fe.solve(a == L, pu_, bc, solver_parameters={'linear_solver' : 'mumps'}) # solver_parameters={'linear_solver': 'gmres','preconditioner': 'ilu'}
pu_n.assign(pu_)
p_n, u_n = fe.split(pu_n)
sigma_.assign(fe.project(material.sigma(p_n, u_n), Vsigma))
p_mono[n, :] = np.array([ p_n(point) for point in points])
color_code=[0.9*(1-(n+1)/Nt)]*3
# plt.plot(y_mono, p_mono[n,:], color=color_code)
# p_, u_ = pu_.split()
vtkfile_p << (p_, t)
# vtkfile_u << (u_, t)
# vtkfile_s_total << (sigma_, t)
dt*=dt_prog
dT.assign(dt)
toc = time.perf_counter()
run_time=toc-tic
print(f"Monolithic scheme run time: {run_time:0.4f} seconds")
plt.plot(p_mono[:,0])