From: Stefano Zampini <stefano.zampini@gmail.com>
Date: Sat, 24 Jun 2023 22:15:26 +0000 (+0200)
Subject: Some versions of PETSc do not apply Shift if zero!
X-Git-Tag: v9.5.0-rc1~21^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f494d501321661f9391f182f1c995f6d9b996812;p=dealii.git

Some versions of PETSc do not apply Shift if zero!
---

diff --git a/include/deal.II/lac/petsc_snes.templates.h b/include/deal.II/lac/petsc_snes.templates.h
index bd5516e410..6feeda5b5a 100644
--- a/include/deal.II/lac/petsc_snes.templates.h
+++ b/include/deal.II/lac/petsc_snes.templates.h
@@ -447,6 +447,7 @@ namespace PETScWrappers
           PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
           PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
         }
+      user->need_dummy_assemble = false;
 
       // Handle the Jacobian-free case
       // This call allows to resample the linearization point
@@ -531,6 +532,7 @@ namespace PETScWrappers
     if (energy)
       AssertPETSc(SNESSetObjective(snes, snes_objective, this));
 
+    this->need_dummy_assemble = false;
     if (setup_jacobian)
       {
         AssertPETSc(SNESSetJacobian(snes,
@@ -543,7 +545,6 @@ namespace PETScWrappers
 
 
         // Do not waste memory by creating a dummy AIJ matrix inside PETSc.
-        this->need_dummy_assemble = false;
         if (!P)
           {
 #  if DEAL_II_PETSC_VERSION_GTE(3, 13, 0)
diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h
index 03c5b16a7f..0d72db1db0 100644
--- a/include/deal.II/lac/petsc_ts.templates.h
+++ b/include/deal.II/lac/petsc_ts.templates.h
@@ -568,20 +568,22 @@ namespace PETScWrappers
           PETSC_ERR_LIB,
           PETSC_ERROR_INITIAL,
           "Failure in ts_ijacobian_with_setup from dealii::PETScWrappers::TimeStepper");
-
       petsc_increment_state_counter(P);
 
       // Handle older versions of PETSc for which we cannot pass a MATSHELL
       // matrix to DMSetMatType. This has been fixed from 3.13 on.
       // What we need to do for older version of PETSc is instead to have
       // a zero matrix with all diagonal entries present.
+      // We use PETSC_MACHINE_EPSILON because some versions of PETSc skip
+      // MatShift if we pass 0.0
       if (user->need_dummy_assemble)
         {
           PetscCall(MatZeroEntries(P));
           PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
           PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
-          PetscCall(MatShift(P, 0.0));
+          PetscCall(MatShift(P, PETSC_MACHINE_EPSILON));
         }
+      user->need_dummy_assemble = false;
 
       // Handle the Jacobian-free case
       // This call allow to resample the linearization point
@@ -668,19 +670,8 @@ namespace PETScWrappers
           PETSC_ERR_LIB,
           PETSC_ERROR_INITIAL,
           "Failure in ts_rhsjacobian from dealii::PETScWrappers::TimeStepper");
-
       petsc_increment_state_counter(P);
 
-      // Handle older versions of PETSc for which we cannot pass a MATSHELL
-      // matrix to DMSetMatType
-      if (user->need_dummy_assemble)
-        {
-          PetscCall(MatZeroEntries(P));
-          PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
-          PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
-          PetscCall(MatShift(P, 0.0));
-        }
-
       // Handle the Jacobian-free case
       // This call allow to resample the linearization point
       // of the MFFD tangent operator