Skip to content

Commit 7de0c41

Browse files
thomasgibsonthomasgibson
committed
Clean up quadrature projection routine
1 parent 58681b4 commit 7de0c41

File tree

3 files changed

+58
-66
lines changed

3 files changed

+58
-66
lines changed

‎grudge/flux_differencing.py‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -61,22 +61,20 @@ def get_skew_symetric_hybridized_diff_mats(
6161
base_grp, quad_vol_grp, face_quad_grp):
6262
from meshmode.discretization.poly_element import diff_matrices
6363
from modepy import faces_for_shape, face_normal
64-
from grudge.projection import volume_quadrature_l2_projection_matrix
6564
from grudge.interpolation import (
6665
volume_quadrature_interpolation_matrix,
6766
surface_quadrature_interpolation_matrix
6867
)
68+
from grudge.op import reference_inverse_mass_matrix
6969

7070
# {{{ Volume operators
7171

7272
weights = quad_vol_grp.quadrature_rule().weights
7373
vdm_q = actx.to_numpy(
74-
volume_quadrature_interpolation_matrix(actx, base_grp, quad_vol_grp)
75-
)
76-
mass_mat = weights * vdm_q.T @ vdm_q
77-
p_mat = actx.to_numpy(
78-
volume_quadrature_l2_projection_matrix(actx, base_grp, quad_vol_grp)
79-
)
74+
volume_quadrature_interpolation_matrix(actx, base_grp, quad_vol_grp))
75+
inv_mass_mat = actx.to_numpy(
76+
reference_inverse_mass_matrix(actx, base_grp))
77+
p_mat = inv_mass_mat @ (vdm_q.T * weights)
8078

8179
# }}}
8280

@@ -112,7 +110,7 @@ def get_skew_symetric_hybridized_diff_mats(
112110

113111
# {{{ Hybridized (volume + surface) operators
114112

115-
q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat
113+
q_mats = [p_mat.T @ (weights * vdm_q.T @ vdm_q) @ diff_mat @ p_mat
116114
for diff_mat in diff_matrices(base_grp)]
117115
e_mat = vf_mat @ p_mat
118116
q_skew_hybridized = np.asarray(

‎grudge/projection.py‎

Lines changed: 45 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
-----------
66
77
.. autofunction:: project
8+
.. autofunction:: volume_quadrature_project
89
"""
910

1011
__copyright__ = """
@@ -31,10 +32,11 @@
3132
THE SOFTWARE.
3233
"""
3334

35+
import numpy as np
3436

3537
from functools import partial
3638

37-
from arraycontext import ArrayContext, map_array_container
39+
from arraycontext import map_array_container
3840
from arraycontext.container import ArrayOrContainerT
3941

4042
from grudge.discretization import DiscretizationCollection
@@ -75,65 +77,59 @@ def project(
7577
return dcoll.connection_from_dds(src, tgt)(vec)
7678

7779

78-
# {{{ Projection matrices
80+
def volume_quadrature_project(
81+
dcoll: DiscretizationCollection, dd_q, vec) -> ArrayOrContainerT:
82+
"""Projects a field on the quadrature discreization, described by *dd_q*,
83+
into the polynomial space described by the volume discretization.
7984
80-
def volume_quadrature_l2_projection_matrix(
81-
actx: ArrayContext, base_element_group, vol_quad_element_group):
82-
"""todo.
83-
"""
84-
@keyed_memoize_in(
85-
actx, volume_quadrature_l2_projection_matrix,
86-
lambda base_grp, vol_quad_grp: (base_grp.discretization_key(),
87-
vol_quad_grp.discretization_key()))
88-
def get_ref_l2_proj_mat(base_grp, vol_quad_grp):
89-
from grudge.interpolation import volume_quadrature_interpolation_matrix
90-
from grudge.op import reference_inverse_mass_matrix
91-
92-
vdm_q = actx.to_numpy(
93-
volume_quadrature_interpolation_matrix(
94-
actx, base_grp, vol_quad_grp
95-
)
96-
)
97-
weights = vol_quad_grp.quadrature_rule().weights
98-
inv_mass_mat = actx.to_numpy(reference_inverse_mass_matrix(actx, base_grp))
99-
return actx.freeze(actx.from_numpy(inv_mass_mat @ (vdm_q.T * weights)))
100-
101-
return get_ref_l2_proj_mat(base_element_group, vol_quad_element_group)
102-
103-
# }}}
104-
105-
106-
def volume_quadrature_project(dcoll: DiscretizationCollection, dd_q, vec):
107-
"""todo.
85+
:arg dd_q: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
86+
:arg vec: a :class:`~meshmode.dof_array.DOFArray` or an
87+
:class:`~arraycontext.container.ArrayContainer` of them.
88+
:returns: a :class:`~meshmode.dof_array.DOFArray` or an
89+
:class:`~arraycontext.container.ArrayContainer` like *vec*.
10890
"""
10991
if not isinstance(vec, DOFArray):
11092
return map_array_container(
11193
partial(volume_quadrature_project, dcoll, dd_q), vec
11294
)
11395

96+
from grudge.geometry import area_element
97+
from grudge.interpolation import volume_quadrature_interpolation_matrix
98+
from grudge.op import inverse_mass
99+
114100
actx = vec.array_context
115101
discr = dcoll.discr_from_dd("vol")
116102
quad_discr = dcoll.discr_from_dd(dd_q)
103+
jacobians = area_element(
104+
actx, dcoll, dd=dd_q,
105+
_use_geoderiv_connection=actx.supports_nonscalar_broadcasting)
117106

118-
return DOFArray(
119-
actx,
120-
data=tuple(
121-
actx.einsum("ij,ej->ei",
122-
volume_quadrature_l2_projection_matrix(
123-
actx,
124-
base_element_group=bgrp,
125-
vol_quad_element_group=qgrp
126-
),
127-
vec_i,
128-
arg_names=("Pq_mat", "vec"),
129-
tagged=(FirstAxisIsElementsTag(),))
130-
131-
for bgrp, qgrp, vec_i in zip(
132-
discr.groups,
133-
quad_discr.groups,
134-
vec
107+
@keyed_memoize_in(
108+
actx, volume_quadrature_project,
109+
lambda base_grp, vol_quad_grp: (base_grp.discretization_key(),
110+
vol_quad_grp.discretization_key()))
111+
def get_mat(base_grp, vol_quad_grp):
112+
vdm_q = actx.to_numpy(
113+
volume_quadrature_interpolation_matrix(
114+
actx, base_grp, vol_quad_grp
115+
)
116+
)
117+
weights = np.diag(vol_quad_grp.quadrature_rule().weights)
118+
return actx.freeze(actx.from_numpy(vdm_q.T @ weights))
119+
120+
return inverse_mass(
121+
dcoll,
122+
DOFArray(
123+
actx,
124+
data=tuple(
125+
actx.einsum("ij,ej,ej->ei",
126+
get_mat(bgrp, qgrp),
127+
jac_i,
128+
vec_i,
129+
arg_names=("vqw_t", "jac", "vec"),
130+
tagged=(FirstAxisIsElementsTag(),))
131+
for bgrp, qgrp, vec_i, jac_i in zip(
132+
discr.groups, quad_discr.groups, vec, jacobians)
135133
)
136134
)
137135
)
138-
139-
# vim: foldmethod=marker

‎test/test_sbp_ops.py‎

Lines changed: 7 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -74,11 +74,11 @@ def test_reference_element_sbp_operators(actx_factory, dim, order):
7474

7575
from meshmode.discretization.poly_element import diff_matrices
7676
from modepy import faces_for_shape, face_normal
77-
from grudge.projection import volume_quadrature_l2_projection_matrix
7877
from grudge.interpolation import (
7978
volume_quadrature_interpolation_matrix,
8079
surface_quadrature_interpolation_matrix
8180
)
81+
from grudge.op import reference_inverse_mass_matrix
8282

8383
for vgrp, qgrp, qfgrp in zip(volm_discr.groups,
8484
quad_discr.groups,
@@ -93,18 +93,16 @@ def test_reference_element_sbp_operators(actx_factory, dim, order):
9393

9494
weights = qgrp.quadrature_rule().weights
9595
vdm_q = actx.to_numpy(
96-
volume_quadrature_interpolation_matrix(actx, vgrp, qgrp)
97-
)
98-
mass_mat = weights * vdm_q.T @ vdm_q
99-
p_mat = actx.to_numpy(
100-
volume_quadrature_l2_projection_matrix(actx, vgrp, qgrp)
101-
)
96+
volume_quadrature_interpolation_matrix(actx, vgrp, qgrp))
97+
inv_mass_mat = actx.to_numpy(
98+
reference_inverse_mass_matrix(actx, vgrp))
99+
p_mat = inv_mass_mat @ (vdm_q.T * weights)
102100

103101
# }}}
104102

105103
# Checks Pq @ Vq = Minv @ Vq.T @ W @ Vq = I
106104
assert np.allclose(p_mat @ vdm_q,
107-
np.identity(len(mass_mat)), rtol=1e-15)
105+
np.identity(len(inv_mass_mat)), rtol=1e-15)
108106

109107
# {{{ Surface operators
110108

@@ -130,7 +128,7 @@ def test_reference_element_sbp_operators(actx_factory, dim, order):
130128

131129
# {{{ Hybridized (volume + surface) operators
132130

133-
q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat
131+
q_mats = [p_mat.T @ (weights * vdm_q.T @ vdm_q) @ diff_mat @ p_mat
134132
for diff_mat in diff_matrices(vgrp)]
135133
e_mat = vf_mat @ p_mat
136134
qtilde_mats = 0.5 * np.asarray(

0 commit comments

Comments
 (0)