compiler: Edit argument to PetscInitialize

Author: ZoeLeibowitz

devito/finite_differences/finite_difference.py

@@ -9,7 +9,7 @@
 
 # Number of digits for FD coefficients to avoid roundup errors and non-deterministic
 # code generation
-_PRECISION = 9
+_PRECISION = 18
 
 
 @check_input

devito/petsc/initialize.py

@@ -12,7 +12,7 @@
 _petsc_initialized = False
 
 
-def PetscInitialize():
+def PetscInitialize(clargs=sys.argv):
     global _petsc_initialized
     if not _petsc_initialized:
         dummy = Symbol(name='d')
@@ -33,11 +33,12 @@ def PetscInitialize():
         # `argv_bytes` must be a list so the memory address persists
         # `os.fsencode` should be preferred over `string().encode('utf-8')`
         # in case there is some system specific encoding in use
-        argv_bytes = list(map(os.fsencode, sys.argv))
-        argv_pointer = (POINTER(c_char)*len(sys.argv))(
+        argv_bytes = list(map(os.fsencode, clargs))
+        argv_pointer = (POINTER(c_char)*len(clargs))(
             *map(lambda s: cast(s, POINTER(c_char)), argv_bytes)
         )
-        op_init.apply(argc=len(sys.argv), argv=argv_pointer)
+        # from IPython import embed; embed()
+        op_init.apply(argc=len(clargs), argv=argv_pointer)
 
         atexit.register(op_finalize.apply)
         _petsc_initialized = True

examples/petsc/solver_options.c

@@ -0,0 +1,105 @@
+
+static char help[] = "Modified snes tutorial 1\n\n";
+
+#include <petscsnes.h>
+
+extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
+extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
+
+int main(int argc, char **argv)
+{
+  SNES        snes;
+  KSP         ksp;
+  PC          pc;
+  Vec         x, r;
+  Mat         J;
+  PetscMPIInt size;
+  PetscBool   flg;
+  PetscFunctionBeginUser;
+
+  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
+  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
+
+  PetscCall(PetscOptionsSetValue(NULL, "-poisson_ksp_type", "cg"));
+  PetscCall(PetscOptionsInsert(NULL, &argc, &argv, NULL));
+  
+  PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
+  PetscCall(SNESSetOptionsPrefix(snes, "poisson_"));
+  PetscCall(SNESSetFromOptions(snes));
+
+  PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
+  PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
+  PetscCall(VecSetFromOptions(x));
+  PetscCall(VecDuplicate(x, &r));
+
+  PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
+  PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
+  PetscCall(MatSetFromOptions(J));
+  PetscCall(MatSetUp(J));
+
+  PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
+  PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
+
+  PetscCall(SNESGetKSP(snes, &ksp));
+  PetscCall(KSPGetPC(ksp, &pc));
+  PetscCall(PCSetType(pc, PCNONE));
+  PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
+
+
+//   PetscCall(VecSet(x, pfive));
+  PetscCall(SNESSolve(snes, NULL, x));
+
+  PetscCall(VecDestroy(&x));
+  PetscCall(VecDestroy(&r));
+  PetscCall(MatDestroy(&J));
+  PetscCall(SNESDestroy(&snes));
+  PetscCall(PetscFinalize());
+  return 0;
+}
+
+PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
+{
+  const PetscScalar *xx;
+  PetscScalar       *ff;
+
+  PetscFunctionBeginUser;
+
+  PetscCall(VecGetArrayRead(x, &xx));
+  PetscCall(VecGetArray(f, &ff));
+
+  /* Compute function */
+  ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
+  ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
+
+  /* Restore vectors */
+  PetscCall(VecRestoreArrayRead(x, &xx));
+  PetscCall(VecRestoreArray(f, &ff));
+  PetscFunctionReturn(PETSC_SUCCESS);
+}
+
+PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
+{
+  const PetscScalar *xx;
+  PetscScalar        A[4];
+  PetscInt           idx[2] = {0, 1};
+
+  PetscFunctionBeginUser;
+
+  PetscCall(VecGetArrayRead(x, &xx));
+
+  A[0] = 2.0 * xx[0] + xx[1];
+  A[1] = xx[0];
+  A[2] = xx[1];
+  A[3] = xx[0] + 2.0 * xx[1];
+  PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
+
+  PetscCall(VecRestoreArrayRead(x, &xx));
+
+  PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
+  PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
+  if (jac != B) {
+    PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
+    PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
+  }
+  PetscFunctionReturn(PETSC_SUCCESS);
+}

examples/petsc/solver_options.py

@@ -48,4 +48,6 @@
 import sys
 petsctools.options._commandline_options = sys.argv[1:]
 tmp = get_commandline_options()
-print("Command line options:", tmp)
+print("Command line options:", tmp)
+
+