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()