Skip to content

Commit 087c622

Browse files
author
Brad Aagaard
committed
Update Cryer full-scale test to use meaningful material properties. Fix checks.
1 parent b199519 commit 087c622

File tree

9 files changed

+216
-188
lines changed

9 files changed

+216
-188
lines changed

‎tests/fullscale/poroelasticity/cryer/Makefile.am‎

Lines changed: 3 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -13,24 +13,17 @@ TESTS = test_pylith.py
1313
dist_check_SCRIPTS = test_pylith.py
1414

1515
dist_noinst_PYTHON = \
16+
generate_gmsh.py \
1617
meshes.py \
1718
TestCryer.py \
1819
cryer_soln.py
1920

2021
dist_noinst_DATA = \
21-
geometry.jou \
22-
bc.jou \
23-
mesh_tet.jou \
24-
mesh_tet.exo \
25-
mesh_hex.jou \
26-
mesh_hex.exo \
27-
cryer.cfg \
22+
pylithapp.cfg \
23+
mesh_tet.msh \
2824
cryer_tet.cfg \
2925
cryer_hex.cfg
3026

31-
noinst_TMP = \
32-
cryer_bc.spatialdb \
33-
cryer_ic.spatialdb
3427

3528
export_datadir = $(abs_builddir)
3629
include $(top_srcdir)/tests/data.am

‎tests/fullscale/poroelasticity/cryer/TestCryer.py‎

Lines changed: 10 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -35,15 +35,10 @@ def setUp(self):
3535
self.checks = [
3636
Check(
3737
mesh_entities=["domain"],
38-
vertex_fields=["displacement"],
38+
vertex_fields=["displacement", "pressure"],
39+
final_time_only=True,
40+
tolerance=0.09,
3941
defaults=defaults,
40-
tolerance=0.5,
41-
),
42-
Check(
43-
mesh_entities=["domain"],
44-
vertex_fields=["pressure"],
45-
defaults=defaults,
46-
scale=1.0e+6,
4742
),
4843
Check(
4944
mesh_entities=["poroelastic"],
@@ -62,34 +57,16 @@ def setUp(self):
6257
defaults=defaults,
6358
),
6459
Check(
65-
mesh_entities=["poroelastic"],
66-
vertex_fields = ["displacement"],
67-
defaults=defaults,
68-
tolerance=0.5,
69-
),
70-
Check(
71-
mesh_entities=["poroelastic"],
72-
vertex_fields = ["pressure"],
73-
defaults=defaults,
74-
scale=1.0e+6,
75-
),
76-
Check(
77-
mesh_entities=["x_neg", "y_neg", "z_neg", "surface_pressure"],
60+
mesh_entities=["bc_xneg", "bc_yneg", "bc_zneg", "bc_shell_pressure"],
7861
filename="output/{name}-{mesh_entity}_info.h5",
7962
vertex_fields=["initial_amplitude"],
8063
defaults=defaults,
8164
),
8265
Check(
83-
mesh_entities=["x_neg", "y_neg", "z_neg", "surface_pressure"],
84-
vertex_fields=["displacement"],
85-
defaults=defaults,
86-
tolerance=0.5,
87-
),
88-
Check(
89-
mesh_entities=["x_neg", "y_neg", "z_neg", "surface_pressure"],
90-
vertex_fields=["pressure"],
66+
mesh_entities=["bc_shell_traction"],
67+
filename="output/{name}-{mesh_entity}_info.h5",
68+
cell_fields=["initial_amplitude"],
9169
defaults=defaults,
92-
scale=1.0e+6,
9370
),
9471
]
9572

@@ -105,7 +82,7 @@ def setUp(self):
10582
self.mesh = meshes.Hex()
10683
super().setUp()
10784

108-
TestCase.run_pylith(self, self.name, ["cryer.cfg", "cryer_hex.cfg"])
85+
TestCase.run_pylith(self, self.name, ["cryer_hex.cfg"])
10986
return
11087

11188

@@ -117,14 +94,14 @@ def setUp(self):
11794
self.mesh = meshes.Tet()
11895
super().setUp()
11996

120-
TestCase.run_pylith(self, self.name, ["cryer.cfg", "cryer_tet.cfg"])
97+
TestCase.run_pylith(self, self.name, ["cryer_tet.cfg"])
12198
return
12299

123100

124101
# -------------------------------------------------------------------------------------------------
125102
def test_cases():
126103
return [
127-
TestHex,
104+
#TestHex,
128105
TestTet,
129106
]
130107

‎tests/fullscale/poroelasticity/cryer/cryer_hex.cfg‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
[pylithapp.metadata]
22
base = [cryer.cfg]
33
keywords = [hexahedral cells]
4-
arguments = [cryer.cfg, cryer_hex.cfg]
4+
arguments = [cryer_hex.cfg]
55

66
[pylithapp.problem]
77
defaults.name = cryer_hex
@@ -10,7 +10,7 @@ defaults.name = cryer_hex
1010
# mesh_generator
1111
# ----------------------------------------------------------------------
1212
[pylithapp.mesh_generator.reader]
13-
filename = mesh_hex.exo
13+
filename = mesh_hex.msh
1414

1515

1616
# End of file

‎tests/fullscale/poroelasticity/cryer/cryer_soln.py‎

Lines changed: 56 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -13,59 +13,54 @@
1313
# Owing to the symmetry of the problem, we only need consider the quarter
1414
# domain case.
1515
#
16-
# -F
17-
# ----------
18-
# | |
19-
# Ux=0 | | P=0
20-
# | |
21-
# | |
22-
# ----------
23-
# Uy=0
16+
# This is based on Cheng (Poroelasticity, Section 7.10) and its implementation in PETSc tutorial
17+
# src/ts/tutorials/ex53.c.
18+
#
19+
# Notes:
20+
# - The traction loading is impulsive, so we use a custom PETSc TS (pylith/utils/TSAdaptImpulse),
21+
# which uses a very small time step for the first step before using the user-specified time step.
22+
# - The accuracy of the solution is poor in the first few time steps due to the impulsive loading
23+
# so we only check the last time step in the simulation, which is more accurate and select a
24+
# a tolerance appropriate for the discretization size.
2425
#
2526
# Dirichlet boundary conditions
26-
# Ux(0,y) = 0
27-
# Uy(x,0) = 0
27+
# Ux(0,y,z) = 0
28+
# Uy(x,0,z) = 0
29+
# Uy(x,y,0) = 0
2830
# Neumann boundary conditions
29-
# \tau_normal(x,ymax) = -1*Pa
31+
# \tau_normal(rmax) = -1*Pa
3032

3133
import numpy
3234

3335
# Physical properties
34-
G = 3.0
36+
G = 30.0e+9
3537
rho_s = 2500
3638
rho_f = 1000
37-
K_fl = 8.0
38-
K_sg = 10.0
39-
K_d = 4.0
39+
K_fl = 80.0e+9
40+
K_sg = 100.0e+9
41+
K_d = 40.0e+9
4042
alpha = 0.6
4143
phi = 0.1
42-
k = 1.5
43-
mu_f = 1.0
44-
P_0 = 1.0
45-
R_0 = 1.0
46-
ndim = 3
44+
k = 1.5e-13
45+
mu_f = 0.001
46+
P_0 = 10.0e+3
47+
R_0 = 10.0e+3
4748

48-
M = 1.0 / ( phi / K_fl + (alpha - phi) /K_sg)
49+
M = 1.0 / ( phi / K_fl + (alpha - phi) /K_sg)
4950
kappa = k/mu_f
5051
K_u = K_d + alpha*alpha*M
5152
S = (3*K_u + 4*G) / (M*(3*K_d + 4*G)) #(1/M) + ( (3*alpha*alpha) / (3*K_d + 4*G) )#
5253
c = kappa / S
5354
nu = (3*K_d - 2*G) / (2*(3*K_d + G))
5455
nu_u = (3*K_u - 2*G) / (2*(3*K_u + G))
55-
U_R_inf = -1.*(P_0*R_0*(1.-2.*nu))/(2.*G*(1.+nu))
56+
U_R_inf = -(P_0*R_0*(1-2*nu))/(2*G*(1+nu))
5657
eta = (alpha*(1-2*nu))/(2*(1-nu))
5758

58-
xmin = 0.0 # m
59-
xmax = 1.0 # m
60-
ymin = 0.0 # m
61-
ymax = 1.0 # m
62-
zmin = 0.0 # m
63-
zmax = 1.0 # m
64-
65-
# Time steps
66-
ts = 0.0028666667 # sec
67-
nts = 2
68-
tsteps = numpy.arange(0.0, ts * nts, ts) + ts # sec
59+
dt0 = 1.0e-6 * 1.0e+8
60+
dt = 5.0e+4
61+
t_end = 100.0e+4
62+
time = numpy.concatenate((numpy.array([dt0]), dt0+numpy.arange(dt, t_end+0.5*dt, dt)))
63+
tsteps = numpy.array([time[-1]])
6964

7065
# ----------------------------------------------------------------------
7166
class AnalyticalSoln(object):
@@ -74,30 +69,27 @@ class AnalyticalSoln(object):
7469
"""
7570
SPACE_DIM = 3
7671
TENSOR_SIZE = 4
77-
ITERATIONS = 50
78-
EPS = 1e-25
72+
ITERATIONS = 100
7973

8074
def __init__(self):
8175
self.fields = {
8276
"displacement": self.displacement,
8377
"pressure": self.pressure,
84-
#"trace_strain": self.trace_strain,
8578
"porosity": self.porosity,
8679
"solid_density": self.solid_density,
8780
"fluid_density": self.fluid_density,
8881
"fluid_viscosity": self.fluid_viscosity,
8982
"shear_modulus": self.shear_modulus,
90-
"undrained_bulk_modulus": self.undrained_bulk_modulus,
9183
"drained_bulk_modulus": self.drained_bulk_modulus,
9284
"biot_coefficient": self.biot_coefficient,
9385
"biot_modulus": self.biot_modulus,
9486
"isotropic_permeability": self.isotropic_permeability,
9587
"initial_amplitude": {
96-
"x_neg": self.zero_vector,
97-
"y_neg": self.zero_vector,
98-
"z_neg": self.zero_vector,
99-
"surface_traction": self.surface_traction,
100-
"surface_pressure": self.zero_scalar
88+
"bc_xneg": self.zero_vector,
89+
"bc_yneg": self.zero_vector,
90+
"bc_zneg": self.zero_vector,
91+
"bc_shell_traction": self.surface_traction,
92+
"bc_shell_pressure": self.zero_scalar
10193
}
10294
}
10395
return
@@ -157,14 +149,6 @@ def fluid_viscosity(self, locs):
157149
fluid_viscosity = mu_f * numpy.ones((1, npts, 1), dtype=numpy.float64)
158150
return fluid_viscosity
159151

160-
def undrained_bulk_modulus(self, locs):
161-
"""
162-
Compute undrained bulk modulus field at locations.
163-
"""
164-
(npts, dim) = locs.shape
165-
undrained_bulk_modulus = K_u * numpy.ones((1, npts, 1), dtype=numpy.float64)
166-
return undrained_bulk_modulus
167-
168152
def drained_bulk_modulus(self, locs):
169153
"""
170154
Compute undrained bulk modulus field at locations.
@@ -206,30 +190,26 @@ def displacement(self, locs):
206190
displacement = numpy.zeros((ntpts, npts, dim), dtype=numpy.float64)
207191

208192
x_n = self.cryerZeros()
209-
center = numpy.where(~locs.any(axis=1))[0]
210193
R = numpy.sqrt(locs[:,0]*locs[:,0] + locs[:,1]*locs[:,1] + locs[:,2]*locs[:,2])
211194
theta = numpy.nan_to_num( numpy.arctan( numpy.nan_to_num( numpy.sqrt(locs[:,0]**2 + locs[:,1]**2) / locs[:,2] ) ) )
212-
phi = numpy.nan_to_num( numpy.arctan( numpy.nan_to_num( locs[:,1] / locs[:,0] ) ) )
213-
R_star = R.reshape([R.size,1]) / R_0
195+
R_star = R.reshape((R.size,1)) / R_0
214196

215-
x_n.reshape([1,x_n.size])
216-
217-
E = numpy.square(1-nu)*numpy.square(1+nu_u)*x_n - 18*(1+nu)*(nu_u-nu)*(1-nu_u)
197+
x_n.reshape((1,x_n.size))
198+
sqrt_xn = numpy.sqrt(x_n)
218199

219-
t_track = 0
200+
E = (1-nu)**2*(1+nu_u)**2*x_n - 18*(1+nu)*(nu_u-nu)*(1-nu_u)
220201

221-
for t in tsteps:
202+
for i_t, t in enumerate(tsteps):
222203
t_star = (c*t)/(R_0**2)
223204
r_exact_N = R_star.ravel() - numpy.nan_to_num(numpy.sum(((12*(1 + nu)*(nu_u - nu)) / \
224-
((1 - 2*nu)*E*R_star*R_star*x_n*numpy.sin(numpy.sqrt(x_n))) ) * \
225-
(3*(nu_u - nu) * (numpy.sin(R_star*numpy.sqrt(x_n)) - R_star*numpy.sqrt(x_n)*numpy.cos(R_star*numpy.sqrt(x_n))) + \
226-
(1 - nu)*(1 - 2*nu)*R_star*R_star*R_star*x_n*numpy.sin(numpy.sqrt(x_n))) * \
205+
((1 - 2*nu)*E*R_star**2*x_n*numpy.sin(sqrt_xn)) ) * \
206+
(3*(nu_u - nu) * (numpy.sin(R_star*sqrt_xn) - R_star*sqrt_xn*numpy.cos(R_star*sqrt_xn)) + \
207+
(1 - nu)*(1 - 2*nu)*R_star**3*x_n*numpy.sin(sqrt_xn)) * \
227208
numpy.exp(-x_n*t_star),axis=1))
228209

229-
displacement[t_track, :, 0] = (r_exact_N*U_R_inf)*numpy.cos(phi)*numpy.sin(theta)
230-
displacement[t_track, :, 1] = (r_exact_N*U_R_inf)*numpy.sin(phi)*numpy.sin(theta)
231-
displacement[t_track, :, 2] = (r_exact_N*U_R_inf)*numpy.cos(theta)
232-
t_track += 1
210+
displacement[i_t, :, 0] = (r_exact_N*U_R_inf)*numpy.cos(phi)*numpy.sin(theta)
211+
displacement[i_t, :, 1] = (r_exact_N*U_R_inf)*numpy.sin(phi)*numpy.sin(theta)
212+
displacement[i_t, :, 2] = (r_exact_N*U_R_inf)*numpy.cos(theta)
233213

234214
return displacement
235215

@@ -244,37 +224,32 @@ def pressure(self, locs):
244224
center = numpy.where(~locs.any(axis=1))[0]
245225

246226
x_n = self.cryerZeros()
227+
sqrt_xn = numpy.sqrt(x_n)
247228

248229
R = numpy.sqrt(locs[:,0]*locs[:,0] + locs[:,1]*locs[:,1] + locs[:,2]*locs[:,2])
249230
R_star = R.reshape([R.size,1]) / R_0
250231
x_n.reshape([1,x_n.size])
251232

252233
E = (1-nu)**2 * (1+nu_u)**2 * x_n - 18*(1+nu)*(nu_u-nu)*(1-nu_u)
253234

254-
t_track = 0
255-
256-
for t in tsteps:
235+
for i_t, t in enumerate(tsteps):
257236

258237
t_star = (c*t)/(R_0**2)
259-
pressure[t_track,:,0] = numpy.sum( ( (18*numpy.square(nu_u-nu) ) / (eta*E) ) * \
260-
( (numpy.sin(R_star*numpy.sqrt(x_n))) / (R_star*numpy.sin(numpy.sqrt(x_n)) ) - 1 ) * \
238+
pressure[i_t,:,0] = P_0*numpy.sum( ( (18*numpy.square(nu_u-nu) ) / (eta*E) ) * \
239+
( (numpy.sin(R_star*sqrt_xn)) / (R_star*numpy.sin(sqrt_xn) ) - 1 ) * \
261240
numpy.exp(-x_n*t_star) , axis=1)
262241

263-
# Account for center value
264-
#pressure[t_track,center] = numpy.sum( (8*eta*(numpy.sqrt(x_n) - numpy.sin(numpy.sqrt(x_n)))) / ( (x_n - 12*eta + 16*eta*eta)*numpy.sin(numpy.sqrt(x_n)) ) * numpy.exp(-x_n * t_star) )
265-
pressure[t_track,center,0] = numpy.sum( ( (18*numpy.square(nu_u-nu) ) / (eta*E) ) * \
266-
( (numpy.sqrt(x_n)) / (numpy.sin(numpy.sqrt(x_n)) ) - 1 ) * \
242+
pressure[i_t,center,0] = P_0*numpy.sum( ( (18*numpy.square(nu_u-nu) ) / (eta*E) ) * \
243+
( sqrt_xn / (numpy.sin(sqrt_xn) ) - 1 ) * \
267244
numpy.exp(-x_n*t_star))
268-
t_track += 1
269-
270245
return pressure
271246

272-
# Series functions
273247

248+
# Series functions
274249

275250
def cryerZeros(self):
276251

277-
f = lambda x: numpy.tan(numpy.sqrt(x)) - (6*(nu_u - nu)*numpy.sqrt(x))/(6*(nu_u - nu) - (1 - nu)*(1 + nu_u)*x) # Compressible Constituents
252+
f = lambda x: numpy.tan(numpy.sqrt(x)) - (6*(nu_u - nu)*numpy.sqrt(x))/(6*(nu_u - nu) - (1 - nu)*(1 + nu_u)*x) # Compressible Constituents
278253

279254
n_series = self.ITERATIONS
280255
a_n = numpy.zeros(n_series) # initializing roots array
@@ -326,7 +301,9 @@ def cryerZeros(self):
326301
def surface_traction(self, locs):
327302
(npts, dim) = locs.shape
328303
traction = numpy.zeros((1, npts, self.SPACE_DIM), dtype=numpy.float64)
329-
traction[:,:,-1] = -1.0
304+
traction[:,:,-1] = -P_0
330305

331306
return traction
307+
308+
332309
# End of file

‎tests/fullscale/poroelasticity/cryer/cryer_tet.cfg‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
[pylithapp.metadata]
22
base = [cryer.cfg]
33
keywords = [tetrahedral cells]
4-
arguments = [cryer.cfg, cryer_tet.cfg]
4+
arguments = [cryer_tet.cfg]
55

66
[pylithapp.problem]
77
defaults.name = cryer_tet
@@ -10,7 +10,7 @@ defaults.name = cryer_tet
1010
# mesh_generator
1111
# ----------------------------------------------------------------------
1212
[pylithapp.mesh_generator.reader]
13-
filename = mesh_tet.exo
13+
filename = mesh_tet.msh
1414

1515

1616
# End of file

0 commit comments

Comments
 (0)