From f494d501321661f9391f182f1c995f6d9b996812 Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Sun, 25 Jun 2023 00:15:26 +0200 Subject: [PATCH] Some versions of PETSc do not apply Shift if zero! --- include/deal.II/lac/petsc_snes.templates.h | 3 ++- include/deal.II/lac/petsc_ts.templates.h | 17 ++++------------- 2 files changed, 6 insertions(+), 14 deletions(-) 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 -- 2.39.5