proper fix for algebraic constraints

diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h
index dd4317eeac6e7e833253fc0db2c8ab9797159edb..9b2a068d50533bd880b5aacb04a68b84390ad2c6 100644 (file)
--- 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<VectorType &>(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));
       }
   }