From d3eb3d9c2ad91488d7812a43020b01e80768d6eb Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Tue, 22 Aug 2023 17:58:00 +0900 Subject: [PATCH] proper fix for algebraic constraints --- include/deal.II/lac/petsc_ts.templates.h | 25 ++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h index dd4317eeac..9b2a068d50 100644 --- a/include/deal.II/lac/petsc_ts.templates.h +++ b/include/deal.II/lac/petsc_ts.templates.h @@ -1088,13 +1088,28 @@ namespace PETScWrappers PetscReal atol, rtol; AssertPETSc(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr)); + // Fill up vectors with individual tolerances, using -1 for the + // algebraic variables. + // We need to input them with the same vector type of our solution + // but we are not sure that the user's VectorType supports operator[] + // We thus first create standard vectors that we know support the + // operator, and then copy the values into the user's type operator. Vec av, rv; Vec py = const_cast(y).petsc_vector(); AssertPETSc(VecDuplicate(py, &av)); AssertPETSc(VecDuplicate(py, &rv)); - VectorType avdealii(av); - VectorType rvdealii(rv); + // Standard vectors + Vec avT, rvT; + PetscInt n, N; + AssertPETSc(VecGetLocalSize(py, &n)); + AssertPETSc(VecGetSize(py, &N)); + AssertPETSc(VecCreateMPI(this->get_mpi_communicator(), n, N, &avT)); + AssertPETSc(VecDuplicate(avT, &rvT)); + + // Fill-up the vectors + VectorBase avdealii(avT); + VectorBase rvdealii(rvT); avdealii = atol; rvdealii = rtol; for (auto i : algebraic_components()) @@ -1104,9 +1119,15 @@ namespace PETScWrappers } avdealii.compress(VectorOperation::insert); rvdealii.compress(VectorOperation::insert); + + // Copy, set and destroy + AssertPETSc(VecCopy(avT, av)); + AssertPETSc(VecCopy(rvT, rv)); AssertPETSc(TSSetTolerances(ts, atol, av, rtol, rv)); AssertPETSc(VecDestroy(&av)); AssertPETSc(VecDestroy(&rv)); + AssertPETSc(VecDestroy(&avT)); + AssertPETSc(VecDestroy(&rvT)); } } -- 2.39.5