Skip to content
Commits on Source (5)
......@@ -31,7 +31,7 @@ import pyopencl.tools as cl_tools
from grudge.array_context import MPIPyOpenCLArrayContext
from grudge.shortcuts import set_up_rk4
from grudge import DiscretizationCollection
from grudge import make_discretization_collection
from mpi4py import MPI
......@@ -47,7 +47,7 @@ class WaveTag:
pass
def main(ctx_factory, dim=2, order=4, visualize=False):
def main(dim=2, order=4, visualize=True):
comm = MPI.COMM_WORLD
num_parts = comm.size
......@@ -83,7 +83,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
else:
local_mesh = comm.scatter(None)
dcoll = DiscretizationCollection(actx, local_mesh, order=order)
dcoll = make_discretization_collection(actx, local_mesh, order=order)
def source_f(actx, dcoll, t=0):
source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
......@@ -196,7 +196,7 @@ if __name__ == "__main__":
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
main(
dim=args.dim,
order=args.order,
visualize=args.visualize)
......@@ -34,7 +34,8 @@ THE SOFTWARE.
"""
from typing import Mapping, Optional, Union, TYPE_CHECKING, Any
from meshmode.discretization.poly_element import ModalGroupFactory
from meshmode.discretization.poly_element import (
InterpolatoryEdgeClusteredGroupFactory, ModalGroupFactory)
from pytools import memoize_method, single_valued
......@@ -79,9 +80,6 @@ def _normalize_discr_tag_to_group_factory(
Mapping[DiscretizationTag, ElementGroupFactory]],
order: Optional[int]
) -> Mapping[DiscretizationTag, ElementGroupFactory]:
from meshmode.discretization.poly_element import \
default_simplex_group_factory
if discr_tag_to_group_factory is None:
if order is None:
raise TypeError(
......@@ -89,8 +87,7 @@ def _normalize_discr_tag_to_group_factory(
)
discr_tag_to_group_factory = {
DISCR_TAG_BASE: default_simplex_group_factory(
base_dim=dim, order=order)}
DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order=order)}
else:
discr_tag_to_group_factory = dict(discr_tag_to_group_factory)
......@@ -102,7 +99,7 @@ def _normalize_discr_tag_to_group_factory(
)
discr_tag_to_group_factory[DISCR_TAG_BASE] = \
default_simplex_group_factory(base_dim=dim, order=order)
InterpolatoryEdgeClusteredGroupFactory(order)
assert discr_tag_to_group_factory is not None
......
......@@ -76,7 +76,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import sys
from warnings import warn
from typing import Hashable, Union, Type, Optional, Any, Tuple
from dataclasses import dataclass, replace
......@@ -491,11 +490,6 @@ def __getattr__(name):
raise AttributeError(f"module {__name__} has no attribute {name}")
if sys.version_info < (3, 7):
for name in _deprecated_name_to_new_name:
globals()[name] = globals()[_deprecated_name_to_new_name[name]]
# }}}
......
......@@ -28,6 +28,7 @@ THE SOFTWARE.
import numpy as np
from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc
from grudge.models import HyperbolicOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
......@@ -113,6 +114,8 @@ class WeakWaveOperator(HyperbolicOperator):
v = w[1:]
actx = u.array_context
base_dd = as_dofdesc("vol", DISCR_TAG_BASE)
# boundary conditions -------------------------------------------------
# dirichlet BCs -------------------------------------------------------
......@@ -160,9 +163,12 @@ class WeakWaveOperator(HyperbolicOperator):
dcoll,
sum(flux(tpair) for tpair in op.interior_trace_pairs(
dcoll, w, comm_tag=self.comm_tag))
+ flux(op.bv_trace_pair(dcoll, self.dirichlet_tag, w, dir_bc))
+ flux(op.bv_trace_pair(dcoll, self.neumann_tag, w, neu_bc))
+ flux(op.bv_trace_pair(dcoll, self.radiation_tag, w, rad_bc))
+ flux(op.bv_trace_pair(
dcoll, base_dd.trace(self.dirichlet_tag), w, dir_bc))
+ flux(op.bv_trace_pair(
dcoll, base_dd.trace(self.neumann_tag), w, neu_bc))
+ flux(op.bv_trace_pair(
dcoll, base_dd.trace(self.radiation_tag), w, rad_bc))
)
)
)
......
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
8
2 3 "inflow"
2 4 "outflow"
2 5 "flow"
2 6 "isothermal_wall"
2 7 "wall_interface"
2 8 "wall_farfield"
3 1 "fluid"
3 2 "wall_insert"
$EndPhysicalNames
$Nodes
117
1 0 -0.01 0
2 0.1 -0.01 0
3 0.1 0.01 0
4 0 0.01 0
5 0.12 -0.01 0
6 0.12 0.01 0
7 0 -0.01 0.02
8 0 0.01 0.02
9 0.1 0.01 0.02
10 0.1 -0.01 0.02
11 0.12 0.01 0.02
12 0.12 -0.01 0.02
13 0.009999999999982485 -0.01 0
14 0.01999999999995601 -0.01 0
15 0.02999999999992731 -0.01 0
16 0.0399999999998959 -0.01 0
17 0.04999999999986855 -0.01 0
18 0.0599999999998943 -0.01 0
19 0.06999999999991978 -0.01 0
20 0.07999999999994742 -0.01 0
21 0.08999999999997288 -0.01 0
22 0.1 -2.65987232239695e-14 0
23 0.0900000000000228 0.01 0
24 0.08000000000004992 0.01 0
25 0.07000000000008301 0.01 0
26 0.06000000000011067 0.01 0
27 0.05000000000013697 0.01 0
28 0.04000000000011013 0.01 0
29 0.0300000000000834 0.01 0
30 0.02000000000005492 0.01 0
31 0.01000000000002775 0.01 0
32 0 2.65987232239695e-14 0
33 0.1099999999999756 0.01 0
34 0.1099999999999756 -0.01 0
35 0.12 -2.65987232239695e-14 0
36 0 -2.65987232239695e-14 0.02
37 0.009999999999982485 0.01 0.02
38 0.01999999999995601 0.01 0.02
39 0.02999999999992731 0.01 0.02
40 0.0399999999998959 0.01 0.02
41 0.04999999999986855 0.01 0.02
42 0.0599999999998943 0.01 0.02
43 0.06999999999991978 0.01 0.02
44 0.07999999999994742 0.01 0.02
45 0.08999999999997288 0.01 0.02
46 0.1 2.65987232239695e-14 0.02
47 0.0900000000000228 -0.01 0.02
48 0.08000000000004992 -0.01 0.02
49 0.07000000000008301 -0.01 0.02
50 0.06000000000011067 -0.01 0.02
51 0.05000000000013697 -0.01 0.02
52 0.04000000000011013 -0.01 0.02
53 0.0300000000000834 -0.01 0.02
54 0.02000000000005492 -0.01 0.02
55 0.01000000000002775 -0.01 0.02
56 0 -0.01 0.009999999999973929
57 0 0.01 0.009999999999973929
58 0.1 0.01 0.009999999999973929
59 0.1 -0.01 0.009999999999973929
60 0.1099999999999756 0.01 0.02
61 0.12 2.65987232239695e-14 0.02
62 0.1100000000000244 -0.01 0.02
63 0.12 0.01 0.009999999999973929
64 0.12 -0.01 0.009999999999973929
65 0.01000000000000377 7.12693062460239e-15 0
66 0.02000000000000486 1.908999274446981e-15 0
67 0.03000000000000474 5.090664731859039e-16 0
68 0.04000000000000338 1.272666182964841e-16 0
69 0.05000000000000278 -1.762611085103303e-29 0
70 0.0600000000000022 -1.272666182965546e-16 0
71 0.07000000000000106 -5.090664731861571e-16 0
72 0.07999999999999929 -1.908999274447778e-15 0
73 0.08999999999999873 -7.126930624604066e-15 0
74 0.1099999999999878 -1.329936161198474e-14 0
75 0 0 0.009999999999986964
76 0.09000000000000011 0.01 0.009999999999979494
77 0.08000000000000039 0.01 0.009999999999983416
78 0.07000000000000083 0.01 0.009999999999986017
79 0.06000000000000119 0.01 0.009999999999987496
80 0.05000000000000147 0.01 0.009999999999987975
81 0.04000000000000164 0.01 0.009999999999987496
82 0.03000000000000169 0.01 0.009999999999986017
83 0.02000000000000146 0.01 0.009999999999983416
84 0.0100000000000009 0.01 0.009999999999979494
85 0.1 0 0.009999999999986964
86 0.01000000000000092 -0.01 0.009999999999979494
87 0.02000000000000149 -0.01 0.009999999999983416
88 0.03000000000000173 -0.01 0.009999999999986017
89 0.04000000000000167 -0.01 0.009999999999987496
90 0.05000000000000149 -0.01 0.009999999999987975
91 0.06000000000000121 -0.01 0.009999999999987496
92 0.07000000000000083 -0.01 0.009999999999986017
93 0.0800000000000004 -0.01 0.009999999999983416
94 0.09000000000000012 -0.01 0.009999999999979494
95 0.01000000000000377 -7.12693062460239e-15 0.02
96 0.02000000000000486 -1.908999274446981e-15 0.02
97 0.03000000000000474 -5.090664731859039e-16 0.02
98 0.04000000000000338 -1.272666182964841e-16 0.02
99 0.05000000000000278 1.762611085103303e-29 0.02
100 0.0600000000000022 1.272666182965546e-16 0.02
101 0.07000000000000106 5.090664731861571e-16 0.02
102 0.07999999999999929 1.908999274447778e-15 0.02
103 0.08999999999999873 7.126930624604066e-15 0.02
104 0.1099999999999878 0.01 0.009999999999986964
105 0.12 0 0.009999999999986964
106 0.11 -0.01 0.009999999999986964
107 0.11 1.329936161198474e-14 0.02
108 0.00999999999999956 -1.42247325030107e-19 0.009999999999992529
109 0.02000000000000087 6.21898366137663e-19 0.009999999999996447
110 0.0300000000000011 7.528699885738973e-19 0.009999999999999053
111 0.04000000000000203 1.942890293093493e-19 0.01000000000000053
112 0.05000000000000147 1.734723475976137e-19 0.01000000000000102
113 0.06000000000000093 3.608224830031225e-19 0.01000000000000053
114 0.07000000000000051 7.476658181459641e-19 0.009999999999999053
115 0.08000000000000101 3.365363543394754e-19 0.009999999999996451
116 0.09000000000000097 6.019490461639405e-19 0.00999999999999253
117 0.1100000000000001 -4.770489558935599e-19 0.01
$EndNodes
$Elements
168
1 3 2 6 1 4 31 65 32
2 3 2 6 1 31 30 66 65
3 3 2 6 1 30 29 67 66
4 3 2 6 1 29 28 68 67
5 3 2 6 1 28 27 69 68
6 3 2 6 1 27 26 70 69
7 3 2 6 1 26 25 71 70
8 3 2 6 1 25 24 72 71
9 3 2 6 1 24 23 73 72
10 3 2 6 1 23 3 22 73
11 3 2 6 1 32 65 13 1
12 3 2 6 1 65 66 14 13
13 3 2 6 1 66 67 15 14
14 3 2 6 1 67 68 16 15
15 3 2 6 1 68 69 17 16
16 3 2 6 1 69 70 18 17
17 3 2 6 1 70 71 19 18
18 3 2 6 1 71 72 20 19
19 3 2 6 1 72 73 21 20
20 3 2 6 1 73 22 2 21
21 3 2 8 2 2 22 74 34
22 3 2 8 2 34 74 35 5
23 3 2 8 2 22 3 33 74
24 3 2 8 2 74 33 6 35
25 3 2 3 16 4 57 75 32
26 3 2 5 16 4 57 75 32
27 3 2 3 16 57 8 36 75
28 3 2 5 16 57 8 36 75
29 3 2 3 16 32 75 56 1
30 3 2 5 16 32 75 56 1
31 3 2 3 16 75 36 7 56
32 3 2 5 16 75 36 7 56
33 3 2 6 20 3 58 76 23
34 3 2 6 20 58 9 45 76
35 3 2 6 20 23 76 77 24
36 3 2 6 20 76 45 44 77
37 3 2 6 20 24 77 78 25
38 3 2 6 20 77 44 43 78
39 3 2 6 20 25 78 79 26
40 3 2 6 20 78 43 42 79
41 3 2 6 20 26 79 80 27
42 3 2 6 20 79 42 41 80
43 3 2 6 20 27 80 81 28
44 3 2 6 20 80 41 40 81
45 3 2 6 20 28 81 82 29
46 3 2 6 20 81 40 39 82
47 3 2 6 20 29 82 83 30
48 3 2 6 20 82 39 38 83
49 3 2 6 20 30 83 84 31
50 3 2 6 20 83 38 37 84
51 3 2 6 20 31 84 57 4
52 3 2 6 20 84 37 8 57
53 3 2 4 24 2 59 85 22
54 3 2 5 24 2 59 85 22
55 3 2 7 24 2 59 85 22
56 3 2 4 24 59 10 46 85
57 3 2 5 24 59 10 46 85
58 3 2 7 24 59 10 46 85
59 3 2 4 24 22 85 58 3
60 3 2 5 24 22 85 58 3
61 3 2 7 24 22 85 58 3
62 3 2 4 24 85 46 9 58
63 3 2 5 24 85 46 9 58
64 3 2 7 24 85 46 9 58
65 3 2 6 28 1 56 86 13
66 3 2 6 28 56 7 55 86
67 3 2 6 28 13 86 87 14
68 3 2 6 28 86 55 54 87
69 3 2 6 28 14 87 88 15
70 3 2 6 28 87 54 53 88
71 3 2 6 28 15 88 89 16
72 3 2 6 28 88 53 52 89
73 3 2 6 28 16 89 90 17
74 3 2 6 28 89 52 51 90
75 3 2 6 28 17 90 91 18
76 3 2 6 28 90 51 50 91
77 3 2 6 28 18 91 92 19
78 3 2 6 28 91 50 49 92
79 3 2 6 28 19 92 93 20
80 3 2 6 28 92 49 48 93
81 3 2 6 28 20 93 94 21
82 3 2 6 28 93 48 47 94
83 3 2 6 28 21 94 59 2
84 3 2 6 28 94 47 10 59
85 3 2 6 29 7 36 95 55
86 3 2 6 29 55 95 96 54
87 3 2 6 29 54 96 97 53
88 3 2 6 29 53 97 98 52
89 3 2 6 29 52 98 99 51
90 3 2 6 29 51 99 100 50
91 3 2 6 29 50 100 101 49
92 3 2 6 29 49 101 102 48
93 3 2 6 29 48 102 103 47
94 3 2 6 29 47 103 46 10
95 3 2 6 29 36 8 37 95
96 3 2 6 29 95 37 38 96
97 3 2 6 29 96 38 39 97
98 3 2 6 29 97 39 40 98
99 3 2 6 29 98 40 41 99
100 3 2 6 29 99 41 42 100
101 3 2 6 29 100 42 43 101
102 3 2 6 29 101 43 44 102
103 3 2 6 29 102 44 45 103
104 3 2 6 29 103 45 9 46
105 3 2 8 42 3 33 104 58
106 3 2 8 42 58 104 60 9
107 3 2 8 42 33 6 63 104
108 3 2 8 42 104 63 11 60
109 3 2 8 46 5 64 105 35
110 3 2 8 46 64 12 61 105
111 3 2 8 46 35 105 63 6
112 3 2 8 46 105 61 11 63
113 3 2 8 50 2 59 106 34
114 3 2 8 50 59 10 62 106
115 3 2 8 50 34 106 64 5
116 3 2 8 50 106 62 12 64
117 3 2 8 51 9 60 107 46
118 3 2 8 51 60 11 61 107
119 3 2 8 51 46 107 62 10
120 3 2 8 51 107 61 12 62
121 5 2 1 1 4 32 65 31 57 75 108 84
122 5 2 1 1 57 75 108 84 8 36 95 37
123 5 2 1 1 31 65 66 30 84 108 109 83
124 5 2 1 1 84 108 109 83 37 95 96 38
125 5 2 1 1 30 66 67 29 83 109 110 82
126 5 2 1 1 83 109 110 82 38 96 97 39
127 5 2 1 1 29 67 68 28 82 110 111 81
128 5 2 1 1 82 110 111 81 39 97 98 40
129 5 2 1 1 28 68 69 27 81 111 112 80
130 5 2 1 1 81 111 112 80 40 98 99 41
131 5 2 1 1 27 69 70 26 80 112 113 79
132 5 2 1 1 80 112 113 79 41 99 100 42
133 5 2 1 1 26 70 71 25 79 113 114 78
134 5 2 1 1 79 113 114 78 42 100 101 43
135 5 2 1 1 25 71 72 24 78 114 115 77
136 5 2 1 1 78 114 115 77 43 101 102 44
137 5 2 1 1 24 72 73 23 77 115 116 76
138 5 2 1 1 77 115 116 76 44 102 103 45
139 5 2 1 1 23 73 22 3 76 116 85 58
140 5 2 1 1 76 116 85 58 45 103 46 9
141 5 2 1 1 32 1 13 65 75 56 86 108
142 5 2 1 1 75 56 86 108 36 7 55 95
143 5 2 1 1 65 13 14 66 108 86 87 109
144 5 2 1 1 108 86 87 109 95 55 54 96
145 5 2 1 1 66 14 15 67 109 87 88 110
146 5 2 1 1 109 87 88 110 96 54 53 97
147 5 2 1 1 67 15 16 68 110 88 89 111
148 5 2 1 1 110 88 89 111 97 53 52 98
149 5 2 1 1 68 16 17 69 111 89 90 112
150 5 2 1 1 111 89 90 112 98 52 51 99
151 5 2 1 1 69 17 18 70 112 90 91 113
152 5 2 1 1 112 90 91 113 99 51 50 100
153 5 2 1 1 70 18 19 71 113 91 92 114
154 5 2 1 1 113 91 92 114 100 50 49 101
155 5 2 1 1 71 19 20 72 114 92 93 115
156 5 2 1 1 114 92 93 115 101 49 48 102
157 5 2 1 1 72 20 21 73 115 93 94 116
158 5 2 1 1 115 93 94 116 102 48 47 103
159 5 2 1 1 73 21 2 22 116 94 59 85
160 5 2 1 1 116 94 59 85 103 47 10 46
161 5 2 2 2 74 22 2 34 117 85 59 106
162 5 2 2 2 117 85 59 106 107 46 10 62
163 5 2 2 2 35 74 34 5 105 117 106 64
164 5 2 2 2 105 117 106 64 61 107 62 12
165 5 2 2 2 33 3 22 74 104 58 85 117
166 5 2 2 2 104 58 85 117 60 9 46 107
167 5 2 2 2 6 33 74 35 63 104 117 105
168 5 2 2 2 63 104 117 105 11 60 107 61
$EndElements
\ No newline at end of file
......@@ -23,11 +23,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
from meshmode.mesh import TensorProductElementGroup
import numpy as np
import numpy.linalg as la
from grudge.array_context import PytestPyOpenCLArrayContextFactory
from arraycontext import pytest_generate_tests_for_array_contexts
pytest_generate_tests = pytest_generate_tests_for_array_contexts(
[PytestPyOpenCLArrayContextFactory])
......@@ -428,44 +430,86 @@ def test_tri_diff_mat(actx_factory, dim, order=4):
# {{{ divergence theorem
def test_2d_gauss_theorem(actx_factory):
@pytest.mark.parametrize("case", ["circle", "tp_box2", "tp_box3", "gh-403"])
def test_gauss_theorem(actx_factory, case, visualize=False):
"""Verify Gauss's theorem explicitly on a mesh"""
pytest.importorskip("meshpy")
from meshpy.geometry import make_circle, GeometryBuilder
from meshpy.triangle import MeshInfo, build
if case == "circle":
from meshpy.geometry import make_circle, GeometryBuilder
from meshpy.triangle import MeshInfo, build
geob = GeometryBuilder()
geob.add_geometry(*make_circle(1))
mesh_info = MeshInfo()
geob.set(mesh_info)
geob = GeometryBuilder()
geob.add_geometry(*make_circle(1))
mesh_info = MeshInfo()
geob.set(mesh_info)
mesh_info = build(mesh_info)
mesh_info = build(mesh_info)
from meshmode.mesh.io import from_meshpy
from meshmode.mesh import BTAG_ALL
from meshmode.mesh.io import from_meshpy
mesh = from_meshpy(mesh_info, order=1)
elif case == "gh-403":
# https://github.com/inducer/meshmode/issues/403
from meshmode.mesh.io import read_gmsh
mesh = read_gmsh("gh-403.msh")
elif case.startswith("tp_box"):
dim = int(case[-1])
mesh = mgen.generate_regular_rect_mesh(
a=(-0.5,)*dim,
b=(0.5,)*dim,
nelements_per_axis=(4,)*dim,
group_cls=TensorProductElementGroup)
else:
raise ValueError(f"unknown case: {case}")
mesh = from_meshpy(mesh_info, order=1)
from meshmode.mesh import BTAG_ALL
actx = actx_factory()
dcoll = DiscretizationCollection(actx, mesh, order=2)
volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
dcoll = make_discretization_collection(actx, mesh, order=2)
volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
x_volm = actx.thaw(volm_disc.nodes())
def f(x):
if len(x) == 3:
x0, x1, x2 = x
elif len(x) == 2:
x0, x1 = x
x2 = 0
else:
raise ValueError("unsupported dimensionality")
return flat_obj_array(
actx.np.sin(3*x[0]) + actx.np.cos(3*x[1]),
actx.np.sin(2*x[0]) + actx.np.cos(x[1])
)
actx.np.sin(3*x0) + actx.np.cos(3*x1) + 2*actx.np.cos(2*x2),
actx.np.sin(2*x0) + actx.np.cos(x1) + 4*actx.np.cos(0.5*x2),
actx.np.sin(1*x0) + actx.np.cos(2*x1) + 3*actx.np.cos(0.8*x2),
)[:dcoll.ambient_dim]
f_volm = f(x_volm)
int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_volm))
div_f = op.local_div(dcoll, f_volm)
int_1 = op.integral(dcoll, "vol", div_f)
prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
normal = geo.normal(actx, dcoll, BTAG_ALL)
int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal))
f_dot_n = prj_f.dot(normal)
int_2 = op.integral(dcoll, BTAG_ALL, f_dot_n)
if visualize:
from grudge.shortcuts import make_visualizer, make_boundary_visualizer
vis = make_visualizer(dcoll)
bvis = make_boundary_visualizer(dcoll)
vis.write_vtk_file(
f"gauss-thm-{case}-vol.vtu", [("div_f", div_f),])
bvis.write_vtk_file(
f"gauss-thm-{case}-bdry.vtu", [
("f_dot_n", f_dot_n),
("normal", normal),
])
assert abs(int_1 - int_2) < 1e-13
......