add function setSchurComplementTolerance

cpppetsc/src/CMakeLists.txt

@@ -65,6 +65,7 @@ add_library(${PROJECT_NAME}-cpppetsc
             ./nestVectors.cc
             ./readFromViewer.cc
             ./setName.cc
+            ./setSchurComplementTolerance.cc
             ./vertexDataOffsets.cc
             ./writeToViewer.cc
 )

cpppetsc/src/include/ae108/cpppetsc/setSchurComplementTolerance.h

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "ae108/cpppetsc/Matrix.h"
+#include "ae108/cpppetsc/ParallelComputePolicy_fwd.h"
+#include "ae108/cpppetsc/SequentialComputePolicy_fwd.h"
+
+namespace ae108 {
+namespace cpppetsc {
+
+/**
+ * @brief When computing the effect of the Schur complement when multiplied with
+ * a vector, PETSc uses a linear solver. This function configures the tolerance
+ * used by this linear solver.
+ *
+ * @param matrix Valid nonzero pointer to a matrix created by
+ * `asSchurComplement`.
+ */
+template <class Policy>
+void setSchurComplementTolerance(
+    Matrix<Policy> *const matrix,
+    const double relativeTolerance = PETSC_DEFAULT,
+    const double absoluteTolerance = PETSC_DEFAULT);
+
+extern template void
+setSchurComplementTolerance(Matrix<SequentialComputePolicy> *, const double,
+                            const double);
+extern template void
+setSchurComplementTolerance(Matrix<ParallelComputePolicy> *, const double,
+                            const double);
+
+} // namespace cpppetsc
+} // namespace ae108
+
+#include "ae108/cpppetsc/UniqueEntity.h"
+#include <petscksp.h>
+
+namespace ae108 {
+namespace cpppetsc {
+
+template <class Policy>
+void setSchurComplementTolerance(Matrix<Policy> *const matrix,
+                                 const double relativeTolerance,
+                                 const double absoluteTolerance) {
+  assert(matrix);
+
+  const auto ksp = [&]() {
+    auto ksp = KSP{};
+    Policy::handleError(MatSchurComplementGetKSP(matrix->data(), &ksp));
+    return ksp;
+  }();
+
+  Policy::handleError(KSPSetTolerances(
+      ksp, relativeTolerance, absoluteTolerance, PETSC_DEFAULT, PETSC_DEFAULT));
+}
+
+} // namespace cpppetsc
+} // namespace ae108
\ No newline at end of file

cpppetsc/src/setSchurComplementTolerance.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/cpppetsc/setSchurComplementTolerance.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/SequentialComputePolicy.h"
+
+namespace ae108 {
+namespace cpppetsc {
+
+template void setSchurComplementTolerance(Matrix<SequentialComputePolicy> *,
+                                          const double, const double);
+template void setSchurComplementTolerance(Matrix<ParallelComputePolicy> *,
+                                          const double, const double);
+
+} // namespace cpppetsc
+} // namespace ae108
\ No newline at end of file

cpppetsc/test/CMakeLists.txt

@@ -47,6 +47,7 @@ add_executable(${PROJECT_NAME}-CppPetscTests
                ./nestMatrices_Test.cc
                ./nestVectors_Test.cc
                ./setName_Test.cc
+               ./setSchurComplementTolerance_Test.cc
                ./vertexDataOffsets_Test.cc
 )
 

cpppetsc/test/setSchurComplementTolerance_Test.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/cpppetsc/Matrix.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/SequentialComputePolicy.h"
+#include "ae108/cpppetsc/asSchurComplement.h"
+#include "ae108/cpppetsc/computeElementsOfMatrix.h"
+#include "ae108/cpppetsc/setSchurComplementTolerance.h"
+#include "ae108/cppptest/Matchers.h"
+#include <gmock/gmock.h>
+
+using ae108::cppptest::ScalarEqIfLocal;
+using testing::Not;
+using testing::Test;
+using testing::Types;
+
+namespace ae108 {
+namespace cpppetsc {
+namespace {
+
+template <class Policy> struct setSchurComplementTolerance_Test : Test {
+  using matrix_type = Matrix<Policy>;
+
+  const matrix_type mat_00 = matrix_type::fromList({
+      {10., 0.},
+      {0., 20.},
+  });
+  const matrix_type mat_01 = matrix_type::fromList({
+      {3., 0.},
+      {0., 4.},
+  });
+  const matrix_type mat_10 = matrix_type::fromList({
+      {5., 0.},
+      {0., 6.},
+  });
+  const matrix_type mat_11 = matrix_type::fromList({
+      {7., 0.},
+      {0., 8.},
+  });
+
+  const matrix_type mat = matrix_type::fromList({
+      {10., 0., 3., 0.},
+      {0., 20., 0., 4.},
+      {5., 0., 7., 0.},
+      {0., 6., 0., 8.},
+  });
+};
+
+using Policies = Types<SequentialComputePolicy, ParallelComputePolicy>;
+
+TYPED_TEST_CASE(setSchurComplementTolerance_Test, Policies);
+
+TYPED_TEST(
+    setSchurComplementTolerance_Test,
+    schur_complement_created_using_submatrices_has_correct_elements_with_low_tolerance) {
+  auto schur = asSchurComplement(&this->mat_00, &this->mat_01, &this->mat_10,
+                                 &this->mat_11);
+  setSchurComplementTolerance(&schur, 1e-17, 1e-17);
+  const auto matrix = computeElementsOfMatrix(schur);
+
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 0, 7. - 3. * 5. / 10.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 0, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 1, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 1, 8. - 4. * 6. / 20.));
+}
+
+TYPED_TEST(
+    setSchurComplementTolerance_Test,
+    schur_complement_created_using_submatrices_has_incorrect_elements_with_high_tolerance) {
+  auto schur = asSchurComplement(&this->mat_00, &this->mat_01, &this->mat_10,
+                                 &this->mat_11);
+  setSchurComplementTolerance(&schur, .9, .9);
+  const auto matrix = computeElementsOfMatrix(schur);
+
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 0, 7.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 0, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 1, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 1, 8.));
+}
+
+TYPED_TEST(
+    setSchurComplementTolerance_Test,
+    schur_complement_created_using_indices_has_correct_elements_with_low_tolerance) {
+  auto schur = asSchurComplement(&this->mat, {2, 3});
+  setSchurComplementTolerance(&schur, 1e-17, 1e-17);
+  const auto matrix = computeElementsOfMatrix(schur);
+
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 0, 7. - 3. * 5. / 10.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 0, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 1, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 1, 8. - 4. * 6. / 20.));
+}
+
+TYPED_TEST(
+    setSchurComplementTolerance_Test,
+    schur_complement_created_using_indices_has_incorrect_elements_with_high_tolerance) {
+  auto schur = asSchurComplement(&this->mat, {2, 3});
+  setSchurComplementTolerance(&schur, .9, .9);
+  const auto matrix = computeElementsOfMatrix(schur);
+
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 0, 7.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 0, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(0, 1, 0.));
+  EXPECT_THAT(matrix, ScalarEqIfLocal(1, 1, 8.));
+}
+
+} // namespace
+} // namespace cpppetsc
+} // namespace ae108
\ No newline at end of file