clean

fluid/mumps_solver.c

@@ -15,57 +15,6 @@ void mumps_set_options(DMUMPS_STRUC_C *id, int *options, int n){
   }
 }
 
-DMUMPS_STRUC_C *mumps_new(int n, int nnz, int *row, int *col) {
-  DMUMPS_STRUC_C *id = malloc(sizeof(DMUMPS_STRUC_C));
-  int myid, ierr;
-  ierr = MPI_Init(NULL, NULL);
-  //ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-  /*Initialize a MUMPS instance. Use MPI_COMM_WORLD.*/
-  id->job=JOB_INIT;
-  id->par=1;
-  id->sym=0;
-  id->comm_fortran=USE_COMM_WORLD;
-  dmumps_c(id);
-  /*Define the problem on the host*/
-  id->n = n;
-  id->nz = nnz;
-  id->irn = malloc(nnz*sizeof(int));
-  memcpy(id->irn, row, sizeof(int)*nnz);
-  id->jcn = malloc(nnz*sizeof(int));
-  memcpy(id->jcn, col, sizeof(int)*nnz);
-  /*No outputs*/
-  /*id->ICNTL(1)=-1;
-  id->ICNTL(2)=-1;
-  id->ICNTL(3)=-1;*/
-  /* output level */
-  id->ICNTL(4)=-1;
-
-  /* element format */
-  id->ICNTL(5)=0;
-  /* permutes the matrix to a zero-free diagonal */
-  id->ICNTL(6)=2;
-  /* reordering metis */
-  id->ICNTL(7)=5;
-  /* transpose matrix (local matrices c->f) */
-  id->ICNTL(9)=1;
-  /* no scaling */
-  id->ICNTL(8)=0;
-  /* When significant extra fill-in is caused by numerical pivoting, increasing ICNTL(14) may help */
-  /* 50 makes all the validation functional */
-  id->ICNTL(14)=50;
-  /* num of openmp threads */
-  id->ICNTL(16)=0;
-  return id;
-}
-
-void mumps_solve(DMUMPS_STRUC_C *id, double *a, double *rhs) {
-  int init = id->a == NULL;
-  id->a = a;
-  id->rhs = rhs;
-  id->job = (init ? 6 : 5);
-  dmumps_c(id);
-}
-
 DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, double *a_elt, double *rhs) {
   DMUMPS_STRUC_C *id = malloc(sizeof(DMUMPS_STRUC_C));
   int myid, ierr;
@@ -107,7 +56,7 @@ DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, dou
   /* 50 makes all the validation functional */
   id->ICNTL(14)=50;
   /* num of openmp threads */
-  id->ICNTL(16)=4;
+  id->ICNTL(16)=0;
   return id;
 }
 

python/_tools.py

@@ -12,13 +12,12 @@ signal.signal(signal.SIGINT, signal.SIG_DFL)
 
 timers = {}
 
-def _np2c(a,dtype=float) :
-    tmp = np.require(a,dtype,"C")
+def _np2c(a,dtype=float,order="C") :
+    tmp = np.require(a,dtype,order)
     r = ctypes.c_void_p(tmp.ctypes.data)
     r.tmp = tmp
     return r
 
-
 def timeit(func):
     def wrapper(*args, **kwargs):
         self = args[0]
@@ -59,7 +58,7 @@ except :
 
 def get_linear_system_package(choices=None):
     if choices is None:
-        choices = ["mumps","petsc","scipy"]
+        choices = ["mumps","petsc","petsc4py","scipy"]
     if isinstance(choices, str):
         choices = [choices]
     for choice in choices:

python/mumpslsys.py

@@ -25,7 +25,7 @@ import numpy as np
 import ctypes
 import ctypes.util
 import os
-from ._tools import timeit
+from ._tools import timeit, _np2c
 from .csr import CSR
 
 dir_path = os.path.dirname(os.path.realpath(__file__))
@@ -33,49 +33,6 @@ lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
 
 options_key = {"verbosity":4, "fill_in":14, "reorder_metis":7, "scaling":8, "openmp_threads" :16}
 
-def _np2c(a,dtype=float,order="C") :
-    tmp = np.require(a,dtype,order)
-    r = ctypes.c_void_p(tmp.ctypes.data)
-    r.tmp = tmp
-    return r
-
-class LinearSystem2 :
-    @timeit
-    def __init__(self, elements, n_fields, options=None, constraints=[]) :
-        self.localsize = elements.shape[1]*n_fields
-        self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
-        idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
-        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
-        self.csr = CSR(idx, self.idx,self.constraints)
-        self.val = np.zeros(self.csr.col.size)
-        row = np.repeat(np.arange(self.csr.size,dtype=np.int32),self.csr.row[1:]-self.csr.row[:-1])
-        lib2.mumps_new.restype = ctypes.c_void_p
-        self.mumps_ptr = ctypes.c_void_p(lib2.mumps_new(ctypes.c_int(self.csr.size),ctypes.c_int(self.csr.row[-1]),_np2c(row+1,np.int32),_np2c(self.csr.col+1,np.int32)))  
-        options = None if options=="" else options
-        if options is not None:
-            options_vec = []
-            for key, val in options.items():
-                if isinstance(key,str):
-                    key = options_key[key]
-                    options_vec.extend([key,val])
-            options_vec = np.array(options_vec, dtype=np.int32)
-            lib2.mumps_set_options(self.mumps_ptr, _np2c(options_vec, np.int32), ctypes.c_int(np.size(options_vec)//2))
-
-    @timeit
-    def local_to_global(self,localv,localm,u, constraints_value):
-        self.csr.assemble_mat(localm, constraints_value, self.val)
-        return self.csr.assemble_rhs(localv, u, constraints_value)
-
-    @timeit
-    def solve(self,rhs) :
-        r = np.copy(rhs)
-        lib2.mumps_solve(self.mumps_ptr,_np2c(self.val),_np2c(r))
-        return r[:self.csr.ndof]
-
-    def __delete__(self):
-        lib2.mumps_element_delete(self.mumps_ptr)
-
-
 class LinearSystem :
     @timeit
     def __init__(self, elements, n_fields, options=None, constraints=[]) :

python/petsc4pylsys.py

@@ -10,7 +10,7 @@ class LinearSystem:
 
     def __init__(self, elements, n_fields, options, constraints) :
         PETSc.Options().prefixPush("fluid_")
-        PETSc.Options().insertString("-pc_type lu -pc_factor_shift_type INBLOCKS -pc_factor_shift_amount 1e-16")
+        PETSc.Options().insertString("-pc_type ilu -pc_factor_levels 5 -pc_factor_shift_type INBLOCKS -pc_factor_shift_amount 1e-16")
         PETSc.Options().insertString(options)
         PETSc.Options().prefixPop()
         self.block = PETSc.Options().getBool("fluid_block_matrix",False)

python/petsclsys.py

@@ -20,7 +20,6 @@
 # see self.mat<http://www.gnu.org/licenses/>.
 import numpy as np
 import atexit
-import weakref
 
 from ._tools import timeit, _np2c
 from .csr import CSR
@@ -71,11 +70,11 @@ def get_option_bool(prefix, name, default):
 class LinearSystem:
 
     def __init__(self, elements, n_fields, options, constraints = []) :
-        insert_options("fluid", "-pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1e-16")
+        insert_options("fluid", "-pc_type ilu -pc_factor_levels 5 -pc_factor_shift_type INBLOCKS -pc_factor_shift_amount 1e-16")
         insert_options("fluid", options)
         self.block = get_option_bool("fluid","-block_matrix",0)
         self.ksp = PetscObj("KSPCreate", "KSPDestroy", COMM_WORLD)
-        lib.KSPSetOptionsPrefix(self.ksp,b"fluid_")
+        lib.KSPSetOptionsPrefix(self.ksp,b"fluid")
         lib.KSPSetFromOptions(self.ksp)
         if self.block:
             if len(constraints) != 0 :