Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add generic nodal/modal bilinear form routines #103

Merged
merged 33 commits into from
Feb 5, 2025
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
88ac3b5
add nodal and modal bilinear form routines
a-alveyblanc Dec 5, 2024
ee62ab0
use resampling routine; make ruff and mypy happy
a-alveyblanc Dec 5, 2024
26e1258
clean it up
a-alveyblanc Dec 5, 2024
c99dc72
test_derivatives -> test_weak_d_dr for clarity
a-alveyblanc Dec 5, 2024
19b1cda
update moda/nodal quadrature mass routines to use new routines
a-alveyblanc Dec 5, 2024
6cbcfd5
fix failing docs
a-alveyblanc Dec 5, 2024
8ea7aed
fix failing ruff test
a-alveyblanc Dec 6, 2024
e40277e
separate quad domain and no quad domain cases; update tests
a-alveyblanc Dec 12, 2024
3881f15
minor typo fix
a-alveyblanc Dec 13, 2024
1de0b7f
Merge branch 'main' of https://github.com/inducer/modepy into nodal-a…
a-alveyblanc Dec 14, 2024
72175ad
undo changes to nodal_/modal_quad_mass_matrix
a-alveyblanc Dec 14, 2024
c026e74
clean things up, make new routine for full bilinear form vs half (?) …
a-alveyblanc Dec 14, 2024
fb7fb1f
fix doc failures
a-alveyblanc Dec 14, 2024
d3d5ac6
add order 0 test back to test_bilinear_forms
a-alveyblanc Dec 14, 2024
c063d23
typo fix, add an additional assert in nodal_bilinear_forms
a-alveyblanc Dec 14, 2024
d3f9e36
*_quadrature_operator -> *_quad_operator
a-alveyblanc Dec 14, 2024
6adb3e7
clarify docs
a-alveyblanc Dec 14, 2024
f914a0a
more doc clarification
a-alveyblanc Dec 14, 2024
6f86301
Merge branch 'main' of https://github.com/inducer/modepy into nodal-a…
a-alveyblanc Dec 15, 2024
7b82ab5
deprecate old functions, delete tests, beef up nodal bilinear form
a-alveyblanc Dec 23, 2024
e5c496c
restore nodal/modal quad mass (keep deprecated)
a-alveyblanc Dec 23, 2024
6b2f053
update docs
a-alveyblanc Dec 23, 2024
e1b5ae1
move facemass to use general routines
a-alveyblanc Dec 31, 2024
6157efc
update docs, remove modal routines, deprecate mass, face mass routines
a-alveyblanc Jan 29, 2025
5cb7c28
Merge branch 'main' of https://github.com/inducer/modepy into nodal-a…
a-alveyblanc Jan 29, 2025
e117c24
fix typos, add unisolvency check
a-alveyblanc Jan 29, 2025
d026f65
Update modepy/matrices.py
a-alveyblanc Jan 30, 2025
b82eea9
update according to feedback
a-alveyblanc Jan 30, 2025
4617088
add back quad_mass_matrix test
a-alveyblanc Jan 30, 2025
17cebb3
Merge branch 'nodal-and-modal-bilinear-forms' of github.com:a-alveybl…
a-alveyblanc Jan 30, 2025
d6d53fa
fix docs
a-alveyblanc Jan 30, 2025
98c8c81
Merge branch 'main' of https://github.com/inducer/modepy into nodal-a…
a-alveyblanc Feb 5, 2025
19206fc
improve arg names and docs
a-alveyblanc Feb 5, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions contrib/weights-from-festa-sommariva.py
Original file line number Diff line number Diff line change
Expand Up @@ -76,17 +76,17 @@ def generate_festa_sommariva_quadrature_rules(outfile):
#
# DEGREE: <degree>

_re_degree = re.compile(r"DEGREE:\s+(\d+)")
_re_xw = re.compile(r"xw\s?=\s?\[(.+?)\];", re.DOTALL)
re_degree = re.compile(r"DEGREE:\s+(\d+)")
re_xw = re.compile(r"xw\s?=\s?\[(.+?)\];", re.DOTALL)

# NOTE: some degrees have multiple quadrature rules with a repeated
# header, so we just take the unique ones here
degrees = np.unique(
np.fromiter(_re_degree.findall(mfile), dtype=np.int64)
np.fromiter(re_degree.findall(mfile), dtype=np.int64)
)

rules = {}
for imatch, match in enumerate(_re_xw.findall(mfile)):
for imatch, match in enumerate(re_xw.findall(mfile)):
d = degrees[imatch]
assert d == imatch + 1, (d, imatch + 1)

Expand Down
4 changes: 4 additions & 0 deletions modepy/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,12 @@
inverse_mass_matrix,
mass_matrix,
modal_mass_matrix_for_face,
modal_quad_bilinear_form,
modal_quad_mass_matrix,
modal_quad_mass_matrix_for_face,
multi_vandermonde,
nodal_mass_matrix_for_face,
nodal_quad_bilinear_form,
nodal_quad_mass_matrix,
nodal_quad_mass_matrix_for_face,
resampling_matrix,
Expand Down Expand Up @@ -158,11 +160,13 @@
"legendre_gauss_tensor_product_nodes",
"mass_matrix",
"modal_mass_matrix_for_face",
"modal_quad_bilinear_form",
"modal_quad_mass_matrix",
"modal_quad_mass_matrix_for_face",
"monomial_basis_for_space",
"multi_vandermonde",
"nodal_mass_matrix_for_face",
"nodal_quad_bilinear_form",
"nodal_quad_mass_matrix",
"nodal_quad_mass_matrix_for_face",
"node_tuples_for_space",
Expand Down
85 changes: 71 additions & 14 deletions modepy/matrices.py
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,8 @@
.. autofunction:: modal_quad_mass_matrix_for_face
.. autofunction:: nodal_quad_mass_matrix_for_face
.. autofunction:: spectral_diag_nodal_mass_matrix
.. autofunction:: modal_quad_bilinear_form
.. autofunction:: nodal_quad_bilinear_form

Differentiation is also convenient to express by using :math:`V^{-1}` to
obtain modal values and then using a Vandermonde matrix for the derivatives
Expand Down Expand Up @@ -351,6 +353,68 @@ def mass_matrix(
return la.inv(inverse_mass_matrix(basis, nodes))


def modal_quad_bilinear_form(
quadrature: Quadrature,
test_functions: Sequence[Callable[[np.ndarray], np.ndarray]]
) -> np.ndarray:
r"""Using the *quadrature*, provide a matrix :math:`A` that
satisfies:

.. math::

\displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j,

where :math:`\phi_i` are the *test_functions* at the nodes
:math:`r_j` of *quadrature*, with corresponding weights :math:`w_j`.
"""
modal_operator = np.empty((len(test_functions), len(quadrature.weights)))

for i, test_f in enumerate(test_functions):
modal_operator[i] = test_f(quadrature.nodes) * quadrature.weights

return modal_operator


def nodal_quad_bilinear_form(
test_functions: Sequence[Callable[[np.ndarray], np.ndarray]],
trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]],
quadrature: Quadrature,
nodes: np.ndarray,
uses_quadrature_domain: bool = False
) -> np.ndarray:
r"""Using *quadrature*, provide a matrix :math:`A` that satisfies:

.. math::

\displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j,

where :math:`\phi_i` are the Lagrange basis functions obtained from
*test_functions* at *nodes*, :math:`w_j` and :math:`r_j` are the weights
and nodes from *quadrature*, and :math:`u_j` are point values of the trial
function at either *nodes* or the *quadrature* nodes depending on whether
a quadrature domain is used (as signified by *uses_quadrature_domain*).

If *uses_quadrature_domain* is set to False, then an interpolation operator
is used to make the operator :math:`N \times N` where :math:`N` is the
number of nodes in *nodes*. Otherwise, the operator has shape :math:`N
\times N_q` where :math:`N_q` is the number of nodes associated with
*quadrature*. Default behavior is to assume a quadrature domain is not used.
"""
if len(test_functions) != nodes.shape[1]:
raise ValueError("volume_nodes not unisolvent with trial_functions")

vdm_out = vandermonde(trial_functions, nodes)

nodal_operator = la.solve(
vdm_out.T, modal_quad_bilinear_form(quadrature, test_functions))

if uses_quadrature_domain:
return nodal_operator
else:
return nodal_operator @ resampling_matrix(
trial_functions, quadrature.nodes, nodes)


def modal_quad_mass_matrix(
quadrature: Quadrature,
test_functions: Sequence[Callable[[np.ndarray], np.ndarray]],
Expand All @@ -367,12 +431,7 @@ def modal_quad_mass_matrix(

.. versionadded :: 2024.2
"""
modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights)))

for i, test_f in enumerate(test_functions):
modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights

return modal_mass_matrix
return modal_quad_bilinear_form(quadrature, test_functions)


def nodal_quad_mass_matrix(
Expand All @@ -399,13 +458,12 @@ def nodal_quad_mass_matrix(
if nodes is None:
nodes = quadrature.nodes

if len(test_functions) != nodes.shape[1]:
raise ValueError("volume_nodes not unisolvent with test_functions")

vdm = vandermonde(test_functions, nodes)

return la.solve(vdm.T,
modal_quad_mass_matrix(quadrature, test_functions))
return nodal_quad_bilinear_form(
test_functions,
test_functions,
quadrature,
nodes,
uses_quadrature_domain=True)


def spectral_diag_nodal_mass_matrix(
Expand Down Expand Up @@ -549,5 +607,4 @@ def nodal_quad_mass_matrix_for_face(
return la.solve(vol_vdm.T,
modal_quad_mass_matrix_for_face(face, face_quad, test_functions))


# vim: foldmethod=marker
64 changes: 64 additions & 0 deletions modepy/test/test_matrices.py
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,70 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None:
assert np.max(err) < 1e-14


@pytest.mark.parametrize("shape_cls", [mp.Hypercube, mp.Simplex])
@pytest.mark.parametrize("dim", [1, 2, 3])
@pytest.mark.parametrize("order", [0, 1, 2, 4])
@pytest.mark.parametrize("nodes_on_bdry", [False, True])
@pytest.mark.parametrize("test_weak_d_dr", [False, True])
def test_bilinear_forms(
shape_cls: type[mp.Shape],
dim: int,
order: int,
nodes_on_bdry: bool,
test_weak_d_dr: bool
) -> None:
shape = shape_cls(dim)
space = mp.space_for_shape(shape, order)

quad_space = mp.space_for_shape(shape, 2*order)
quad = mp.quadrature_for_space(quad_space, shape)

if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1):
if nodes_on_bdry:
nodes = mp.legendre_gauss_lobatto_tensor_product_nodes(
shape.dim, order,
)
else:
nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order)
elif isinstance(shape, mp.Simplex):
if nodes_on_bdry:
nodes = mp.warp_and_blend_nodes(shape.dim, order)
else:
nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes
else:
raise AssertionError()

basis = mp.orthonormal_basis_for_space(space, shape)

if test_weak_d_dr and order not in [0, 1]:
mass_inv = mp.inverse_mass_matrix(basis, nodes)

for ax in range(dim):
f = 1 - nodes[ax]**2
fp = -2*nodes[ax]

weak_operator = mp.nodal_quad_bilinear_form(
basis.derivatives(ax),
basis.functions,
quad,
nodes,
uses_quadrature_domain=False)

err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp)
assert err <= 1e-12
else:
quad_mass_mat = mp.nodal_quad_bilinear_form(
basis.functions,
basis.functions,
quad,
nodes,
uses_quadrature_domain=False)

vdm_mass_mat = mp.mass_matrix(basis, nodes)
err = la.norm(quad_mass_mat - vdm_mass_mat) / la.norm(vdm_mass_mat)
assert err < 1e-14


# You can test individual routines by typing
# $ python test_modes.py 'test_routine()'

Expand Down
Loading