PREALLOCATE and ALLOW_NEW_NONZEROS

seamplex · Jul 31, 2023 · e712bf1 · e712bf1
1 parent 3129adc
commit e712bf1
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 5 deletions.
diff --git a/src/feenox.h b/src/feenox.h
@@ -1817,6 +1817,7 @@ struct feenox_t {
 
 #ifdef HAVE_PETSC    
     char *petsc_options;
+    PetscBool pre_allocate;        // preallocate? onyl works for petsc >= 3.19
     PetscBool allow_new_nonzeros;  // flag to set MAT_NEW_NONZERO_ALLOCATION_ERR to false, needed in some rare cases
     PetscBool petscinit_called;    // flag
 

diff --git a/src/parser/parser.c b/src/parser/parser.c
@@ -2985,6 +2985,13 @@ int feenox_parse_problem(void) {
 ///kw_pde+PROBLEM+detail The latter should have at least one symmetry BC whilst the former does not.
     } else if (strcasecmp(token, "ALLOW_UNRESOLVED_BCS") == 0) {
       feenox.pde.unresolved_bcs = unresolved_bcs_allow;
+
+///kw_pde+PROBLEM+usage [ PREALLOCATE ]
+    } else if (strcasecmp(token, "PREALLOCATE") == 0) {
+      feenox.pde.pre_allocate = PETSC_TRUE;
+///kw_pde+PROBLEM+usage [ ALLOW_NEW_NONZEROS ]
+    } else if (strcasecmp(token, "ALLOW_NEW_NONZEROS") == 0) {
+      feenox.pde.allow_new_nonzeros = PETSC_TRUE;
 
 ///kw_pde+PROBLEM+detail If the special variable `end_time` is zero, FeenoX solves a static
 ///kw_pde+PROBLEM+detail  problem---although the variable `static_steps` is still honored.

diff --git a/src/pdes/init.c b/src/pdes/init.c
@@ -598,11 +598,11 @@ Mat feenox_problem_create_matrix(const char *name) {
   Mat A = NULL;
   petsc_call_null(MatCreate(PETSC_COMM_WORLD, &A));
   petsc_call_null(MatSetSizes(A, feenox.pde.size_local, feenox.pde.size_local, feenox.pde.size_global, feenox.pde.size_global));
+  if (feenox.pde.pre_allocate == PETSC_TRUE || PETSC_VERSION_LT(3,19,0)) {
+    petsc_call_null(MatMPIAIJSetPreallocation(A, feenox.pde.width, PETSC_NULLPTR, feenox.pde.width, PETSC_NULLPTR));
+    petsc_call_null(MatSeqAIJSetPreallocation(A, feenox.pde.width, PETSC_NULLPTR));
+  }
   petsc_call_null(MatSetFromOptions(A));
-//#if PETSC_VERSION_LT(3,19,0)
-  petsc_call_null(MatMPIAIJSetPreallocation(A, feenox.pde.width, PETSC_NULLPTR, feenox.pde.width, PETSC_NULLPTR));
-  petsc_call_null(MatSeqAIJSetPreallocation(A, feenox.pde.width, PETSC_NULLPTR));
-//#endif
 
   // this flag needs the matrix type to be set, and we had just set it with setfromoptions   
   petsc_call_null(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
@@ -611,7 +611,7 @@ Mat feenox_problem_create_matrix(const char *name) {
     petsc_call_null(PetscObjectSetName((PetscObject)(A), name));
   }
 
-  if (feenox.pde.allow_new_nonzeros) {
+  if (feenox.pde.allow_new_nonzeros == PETSC_TRUE) {
     petsc_call_null(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
   }