From 44b23d9c4e6d6ea124fc7b161549507ddb363513 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sun, 20 Sep 2026 20:56:27 +0000 Subject: [PATCH] Put the material mode index last The multilinear forms yield the material basis as the leading axis of kDivM, kDivMT, rDivM and V3mTo2nTWDivM. yateto assigns stride one to the leading index and aligns it to the vector width, which would give the material mode index the smallest stride and pad it to the vector width. The fused simulation index is the one that wants that end, so move the material mode index to the trailing, slowest varying axis instead. Adds tests asserting that a material order of one reproduces the matrices for constant material for all four families, and that the material index ends up last. AI-generated. Model: Claude Opus 5 --- src/seissol_matrices/base.py | 13 +++++ src/seissol_matrices/dg_matrices.py | 55 +++++++++++---------- src/seissol_matrices/dr_matrices.py | 7 ++- tests/test_material_order.py | 74 +++++++++++++++++++++++++++++ 4 files changed, 121 insertions(+), 28 deletions(-) create mode 100644 tests/test_material_order.py diff --git a/src/seissol_matrices/base.py b/src/seissol_matrices/base.py index 2137683..c22c079 100644 --- a/src/seissol_matrices/base.py +++ b/src/seissol_matrices/base.py @@ -1,6 +1,19 @@ import numpy as np +def material_index_last(tensor): + """ + Move the material mode index of a tensor from the leading to the + trailing axis. + + The multilinear forms yield the material basis as the leading axis. + yateto assigns stride one to the leading index and aligns it to the + vector width, so the material index belongs at the trailing, slowest + varying end instead. + """ + return np.moveaxis(tensor, 0, -1) + + def collocate(basis, points): # takes a basis object, and a points array # of the form: npoints × dim diff --git a/src/seissol_matrices/dg_matrices.py b/src/seissol_matrices/dg_matrices.py index 114d6a2..a976bdd 100644 --- a/src/seissol_matrices/dg_matrices.py +++ b/src/seissol_matrices/dg_matrices.py @@ -165,27 +165,29 @@ def kDivM(self, dim, matorder=None): mass = self.mass_matrix() if matorder is None: stiffness = self.stiffness_matrix(dim) - else: - stiffness = self.multilinear_form( - False, [None, dim, None], order=[matorder, self.order, self.order] - ) - return np.linalg.solve(mass, stiffness) + return np.linalg.solve(mass, stiffness) + + stiffness = self.multilinear_form( + False, [None, dim, None], order=[matorder, self.order, self.order] + ) + return base.material_index_last(np.linalg.solve(mass, stiffness)) def kDivMT(self, dim, matorder=None): mass = self.mass_matrix() if matorder is None: stiffness = self.stiffness_matrix(dim) - sT = stiffness.T - else: - correct1 = self.multilinear_form( - False, [dim, None, None], order=[matorder, self.order, self.order] - ) - correct2 = self.multilinear_form( - False, [None, None, dim], order=[matorder, self.order, self.order] - ) - sT = correct1 + correct2 + return np.linalg.solve(mass, stiffness.T) - return np.linalg.solve(mass, sT) + # the material sits inside the derivative, so the weak form transpose + # picks up the term in which the material basis is differentiated + differentiated_material = self.multilinear_form( + False, [dim, None, None], order=[matorder, self.order, self.order] + ) + differentiated_basis = self.multilinear_form( + False, [None, None, dim], order=[matorder, self.order, self.order] + ) + sT = differentiated_material + differentiated_basis + return base.material_index_last(np.linalg.solve(mass, sT)) def face_to_face_parametrisation(self, x, side): # implement Dumbser & Käser, 2006 Table 2 b) @@ -252,17 +254,18 @@ def rDivM(self, side, matorder=None): mass = self.mass_matrix() if matorder is None: matrix = self.rT(side).T - else: - prematrix = self.multilinear_form( - True, - [None, None, None], - order=[matorder, self.order, self.order], - side=[None, None, side], - dim=[3, 2, 3], - ) - facemass = self.face_generator.mass_matrix() - matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1)) - return np.linalg.solve(mass, matrix) + return np.linalg.solve(mass, matrix) + + prematrix = self.multilinear_form( + True, + [None, None, None], + order=[matorder, self.order, self.order], + side=[None, None, side], + dim=[3, 2, 3], + ) + facemass = self.face_generator.mass_matrix() + matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1)) + return base.material_index_last(np.linalg.solve(mass, matrix)) def collocate_volume(self, points): return base.collocate(self.generator, points) diff --git a/src/seissol_matrices/dr_matrices.py b/src/seissol_matrices/dr_matrices.py index a22ed1b..064ebda 100755 --- a/src/seissol_matrices/dr_matrices.py +++ b/src/seissol_matrices/dr_matrices.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import numpy as np +from seissol_matrices import base from seissol_matrices import basis_functions from seissol_matrices import dg_matrices from seissol_matrices import quad_points @@ -79,9 +80,11 @@ def V3mTo2nTWDivM(self, a, b, matorder=None): for i in range(n): W[i, i] = weights[i] - matrixT = matrix.T if matorder is None else matrix.transpose((0, 2, 1)) + if matorder is None: + return np.linalg.solve(mass, np.dot(matrix.T, W)) - return np.linalg.solve(mass, np.dot(matrixT, W)) + matrixT = matrix.transpose((0, 2, 1)) + return base.material_index_last(np.linalg.solve(mass, np.dot(matrixT, W))) def quadpoints(self): points = self.quadrule.points() diff --git a/tests/test_material_order.py b/tests/test_material_order.py new file mode 100644 index 0000000..20a781c --- /dev/null +++ b/tests/test_material_order.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 + +import unittest +import numpy as np + +from seissol_matrices import dg_matrices +from seissol_matrices import dr_matrices +from seissol_matrices import quad_points + + +class abstract_tester(object): + def compare(self, a, b): + self.assertEqual(a.shape, b.shape) + np.testing.assert_allclose(a, b, rtol=0.0, atol=1.2e-11) + + def test_kDivM_material_index_last(self): + for dim in range(3): + matrix = self.generator.kDivM(dim, 2) + self.assertEqual(matrix.shape, (self.nbf, self.nbf, 4)) + + def test_kDivM_constant_material(self): + for dim in range(3): + self.compare( + self.generator.kDivM(dim, 1)[:, :, 0], self.generator.kDivM(dim) + ) + + def test_kDivMT_constant_material(self): + for dim in range(3): + self.compare( + self.generator.kDivMT(dim, 1)[:, :, 0], self.generator.kDivMT(dim) + ) + + def test_rDivM_constant_material(self): + for side in range(4): + self.compare( + self.generator.rDivM(side, 1)[:, :, 0], self.generator.rDivM(side) + ) + + def test_V3mTo2nTWDivM_constant_material(self): + for a in range(4): + for b in range(4): + self.compare( + self.dr_generator.V3mTo2nTWDivM(a, b, 1)[:, :, 0], + self.dr_generator.V3mTo2nTWDivM(a, b), + ) + + +def setUpClassFromOrder(cls, order): + cls.order = order + cls.nbf = order * (order + 1) * (order + 2) // 6 + cls.generator = dg_matrices.dg_generator(order, 3) + cls.dr_generator = dr_matrices.dr_generator(order, quad_points.stroud(order + 1)) + + +class test_material_order_2(abstract_tester, unittest.TestCase): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 2) + + +class test_material_order_3(abstract_tester, unittest.TestCase): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 3) + + +class test_material_order_4(abstract_tester, unittest.TestCase): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 4) + + +if __name__ == "__main__": + unittest.main()