Skip to content

Commit 09bf038

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

File tree

18 files changed

+274
-504
lines changed

18 files changed

+274
-504
lines changed

tests/fullscale/poroelasticity/mandel/Makefile.am

Lines changed: 7 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -13,34 +13,21 @@ TESTS = test_pylith.py
1313
dist_check_SCRIPTS = test_pylith.py
1414

1515
dist_noinst_PYTHON = \
16+
generate_gmsh.py \
1617
meshes.py \
1718
TestMandel.py \
1819
mandel_soln.py \
19-
mandel_gendb.py \
20-
TestMandelCompaction.py \
21-
mandel_compaction_soln.py \
22-
mandel_compaction_gendb.py
20+
mandel_gendb.py
2321

2422
dist_noinst_DATA = \
25-
geometry.jou \
26-
bc.jou \
27-
mesh_tri.jou \
28-
mesh_tri.exo \
29-
mesh_quad.jou \
30-
mesh_quad.exo \
31-
mandel.cfg \
23+
pylithapp.cfg \
24+
mesh_tri.msh \
25+
mesh_quad.msh \
3226
mandel_tri.cfg \
33-
mandel_quad.cfg \
34-
mandel_compaction.cfg \
35-
mandel_compaction_tri.cfg \
36-
mandel_compaction_quad.cfg
27+
mandel_quad.cfg
3728

3829
noinst_TMP = \
39-
mandel_bc.spatialdb \
40-
mandel_ic.spatialdb \
41-
mandel_compaction_bc.spatialdb \
42-
mandel_compaction_ic.spatialdb
43-
30+
mandel_disp.timedb
4431

4532

4633
export_datadir = $(abs_builddir)

tests/fullscale/poroelasticity/mandel/TestMandel.py

Lines changed: 7 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717

1818
import unittest
1919

20-
from pylith.testing.FullTestApp import (FullTestCase, Check, check_data)
20+
from pylith.testing.FullTestApp import (FullTestCase, Check)
2121

2222
import meshes
2323
import mandel_soln
@@ -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", "trace_strain"],
39+
final_time_only=True,
3940
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,
41+
tolerance=0.01,
4742
),
4843
Check(
4944
mesh_entities=["poroelastic"],
@@ -62,35 +57,11 @@ 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", "x_pos", "y_neg", "y_pos"],
60+
mesh_entities=["bc_xneg", "bc_xpos", "bc_yneg"],
7861
filename="output/{name}-{mesh_entity}_info.h5",
7962
vertex_fields=["initial_amplitude"],
8063
defaults=defaults,
8164
),
82-
Check(
83-
mesh_entities=["x_neg", "x_pos", "y_neg", "y_pos"],
84-
vertex_fields=["displacement"],
85-
defaults=defaults,
86-
tolerance=0.5,
87-
),
88-
Check(
89-
mesh_entities=["x_neg", "x_pos", "y_neg", "y_pos"],
90-
vertex_fields=["pressure"],
91-
defaults=defaults,
92-
scale=1.0e+6,
93-
),
9465
]
9566

9667
def run_pylith(self, testName, args):
@@ -105,7 +76,7 @@ def setUp(self):
10576
self.mesh = meshes.Quad()
10677
super().setUp()
10778

108-
TestCase.run_pylith(self, self.name, ["mandel.cfg", "mandel_quad.cfg"])
79+
TestCase.run_pylith(self, self.name, ["mandel_quad.cfg"])
10980
return
11081

11182

@@ -117,7 +88,7 @@ def setUp(self):
11788
self.mesh = meshes.Tri()
11889
super().setUp()
11990

120-
TestCase.run_pylith(self, self.name, ["mandel.cfg", "mandel_tri.cfg"])
91+
TestCase.run_pylith(self, self.name, ["mandel_tri.cfg"])
12192
return
12293

12394

tests/fullscale/poroelasticity/mandel/bc.jou

Lines changed: 0 additions & 34 deletions
This file was deleted.
Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
#!/usr/bin/env nemesis
2+
3+
import gmsh
4+
from pylith.meshio.gmsh_utils import BoundaryGroup, MaterialGroup, GenerateMesh
5+
6+
7+
class App(GenerateMesh):
8+
"""
9+
Block is DOMAIN_X by DOMAIN_Y with discretization size DX.
10+
11+
p4------------p3
12+
| |
13+
| |
14+
| |
15+
| |
16+
| |
17+
p1------------p2
18+
"""
19+
20+
DOMAIN_X = 8.0e+3
21+
DOMAIN_Y = 1.0e+3
22+
DX = 0.2e+3
23+
24+
def __init__(self):
25+
self.cell_choices = {
26+
"required": True,
27+
"choices": ["tri", "quad"],
28+
}
29+
self.filename = "mesh.msh"
30+
31+
def create_geometry(self):
32+
"""Create geometry."""
33+
lx = self.DOMAIN_X
34+
ly = self.DOMAIN_Y
35+
x0 = 0.0
36+
y0 = 0.0
37+
38+
p1 = gmsh.model.geo.add_point(x0, y0, 0.0)
39+
p2 = gmsh.model.geo.add_point(x0 + lx, y0, 0.0)
40+
p3 = gmsh.model.geo.add_point(x0 + lx, y0 + ly, 0.0)
41+
p4 = gmsh.model.geo.add_point(x0, y0 + ly, 0.0)
42+
43+
self.l_yneg = gmsh.model.geo.add_line(p1, p2)
44+
self.l_xpos = gmsh.model.geo.add_line(p2, p3)
45+
self.l_ypos = gmsh.model.geo.add_line(p3, p4)
46+
self.l_xneg = gmsh.model.geo.add_line(p4, p1)
47+
48+
c1 = gmsh.model.geo.add_curve_loop(
49+
[self.l_yneg, self.l_xpos, self.l_ypos, self.l_xneg]
50+
)
51+
self.s_domain = gmsh.model.geo.add_plane_surface([c1])
52+
53+
gmsh.model.geo.synchronize()
54+
55+
def mark(self):
56+
"""Mark geometry for materials, boundary conditions, faults, etc."""
57+
materials = (MaterialGroup(tag=1, entities=[self.s_domain]),)
58+
for material in materials:
59+
material.create_physical_group()
60+
61+
face_groups = (
62+
BoundaryGroup(
63+
name="boundary_xneg",
64+
tag=10,
65+
dim=1,
66+
entities=[self.l_xneg],
67+
),
68+
BoundaryGroup(
69+
name="boundary_xpos",
70+
tag=11,
71+
dim=1,
72+
entities=[self.l_xpos],
73+
),
74+
BoundaryGroup(
75+
name="boundary_yneg",
76+
tag=12,
77+
dim=1,
78+
entities=[self.l_yneg],
79+
),
80+
BoundaryGroup(
81+
name="boundary_ypos",
82+
tag=13,
83+
dim=1,
84+
entities=[self.l_ypos],
85+
),
86+
)
87+
for group in face_groups:
88+
group.create_physical_group()
89+
90+
def generate_mesh(self, cell):
91+
"""Generate the mesh. Should also include optimizing the mesh quality."""
92+
gmsh.option.setNumber("Mesh.MeshSizeMin", self.DX)
93+
gmsh.option.setNumber("Mesh.MeshSizeMax", self.DX)
94+
if cell == "quad":
95+
# Generate a tri mesh and then recombine cells to form quadrilaterals.
96+
# We use the Frontal-Delaunay for Quads algorithm.
97+
gmsh.option.setNumber("Mesh.Algorithm", 8)
98+
gmsh.model.mesh.generate(2)
99+
gmsh.model.mesh.recombine()
100+
else:
101+
gmsh.model.mesh.generate(2)
102+
gmsh.model.mesh.optimize("Laplace2D")
103+
104+
105+
if __name__ == "__main__":
106+
App().main()

tests/fullscale/poroelasticity/mandel/geometry.jou

Lines changed: 0 additions & 14 deletions
This file was deleted.

tests/fullscale/poroelasticity/mandel/mandel_gendb.py

Lines changed: 13 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -24,71 +24,19 @@ class GenerateDB(object):
2424
def run(self):
2525
"""Generate the database.
2626
"""
27-
# Domain
28-
x1 = numpy.arange(-0.1, 10.1, 0.1)
29-
y1 = numpy.arange(-0.1, 1.01, 0.1)
30-
x, y = numpy.meshgrid(x1, y1)
31-
32-
xy = numpy.zeros((len(x1) * len(y1), 2), dtype=numpy.float64)
33-
xy[:, 0] = x.ravel()
34-
xy[:, 1] = y.ravel()
35-
36-
from mandel_soln import AnalyticalSoln
37-
soln = AnalyticalSoln()
38-
disp = soln.initial_displacement(xy)
39-
pres = soln.initial_pressure(xy)
40-
trace_strain = soln.initial_trace_strain(xy)
41-
42-
from spatialdata.geocoords.CSCart import CSCart
43-
cs = CSCart()
44-
cs.inventory.spaceDim = 2
45-
cs._configure()
46-
data = {
47-
'x': x1,
48-
'y': y1,
49-
'points': xy,
50-
'coordsys': cs,
51-
'data_dim': 2,
52-
'values': [{'name': "initial_amplitude_x",
53-
'units': "m",
54-
'data': numpy.ravel(disp[0, :, 0])},
55-
{'name': "initial_amplitude_y",
56-
'units': "m",
57-
'data': numpy.ravel(disp[0, :, 1])},
58-
{'name': "initial_pressure",
59-
'units': "Pa",
60-
'data': numpy.ravel(pres[0, :])},
61-
{'name': "initial_trace_strain",
62-
'units': "none",
63-
'data': numpy.ravel(trace_strain[0, :])}]}
64-
65-
from spatialdata.spatialdb.SimpleGridAscii import SimpleGridAscii
66-
io = SimpleGridAscii()
67-
io.inventory.filename = "mandel_bc.spatialdb"
68-
io._configure()
69-
io.write(data)
70-
data["values"] = [
71-
{
72-
'name': "displacement_x",
73-
'units': "m",
74-
'data': numpy.ravel(disp[0, :, 0])
75-
}, {
76-
'name': "displacement_y",
77-
'units': "m",
78-
'data': numpy.ravel(disp[0, :, 1])
79-
}, {
80-
'name': "pressure",
81-
'units': "Pa",
82-
'data': numpy.ravel(pres[0, :])
83-
}, {
84-
'name': "trace_strain",
85-
'units': "none",
86-
'data': numpy.ravel(trace_strain[0, :])
87-
}]
88-
io.inventory.filename = "mandel_ic.spatialdb"
89-
io._configure()
90-
io.write(data)
91-
return
27+
import mandel_soln
28+
soln = mandel_soln.AnalyticalSoln()
29+
locs = numpy.array([[0, mandel_soln.y_max]])
30+
t_steps = numpy.arange(0.0, 100.0e+5, 1.0e+4)
31+
32+
# Need to restore tsteps in mandel_soln for checking fields
33+
soln_tsteps = mandel_soln.tsteps
34+
mandel_soln.tsteps = t_steps
35+
displacement = soln.displacement(locs)
36+
37+
from spatialdata.spatialdb.TimeHistoryIO import write
38+
write(t_steps, displacement[:,0,1], units="m", filename="mandel_disp.timedb")
39+
mandel_soln.tsteps = soln_tsteps
9240

9341

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

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

1515

1616
# End of file

0 commit comments

Comments
 (0)