]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Some versions of PETSc do not apply Shift if zero!
authorStefano Zampini <stefano.zampini@gmail.com>
Sat, 24 Jun 2023 22:15:26 +0000 (00:15 +0200)
committerStefano Zampini <stefano.zampini@gmail.com>
Sun, 25 Jun 2023 00:00:43 +0000 (03:00 +0300)
include/deal.II/lac/petsc_snes.templates.h
include/deal.II/lac/petsc_ts.templates.h

index bd5516e410507c406b8a8892a0a372099cf41fca..6feeda5b5a33d0c80cfa8e890ffdf3f70d4bc0a0 100644 (file)
@@ -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)
index 03c5b16a7fc08c23f9c50283cfb77b451a960989..0d72db1db09da873f2af1451e7678b63bddf283a 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.