mumps elements with constraints

fluid/mumps_solver.c

 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include "mpi.h"
 #include "dmumps_c.h"
 #include "string.h"
@@ -14,6 +15,57 @@ 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)=4;
+  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;
@@ -58,6 +110,7 @@ DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, dou
   id->ICNTL(16)=4;
   return id;
 }
+
 void mumps_element_solve(DMUMPS_STRUC_C *id, double *a_elt, double *rhs) {
   int init = id->a_elt == NULL;
   id->a_elt = a_elt;

python/CMakeLists.txt

@@ -24,9 +24,11 @@ set(SRC
   fluid.py
   _tools.py
   lmgc90Interface.py
+  csr.py
   scontact.py
   time_integration.py
   petsclsys.py
+  petsc4pylsys.py
   scipylsys.py
   VTK.py
 )

python/_tools.py

@@ -37,6 +37,11 @@ try :
     have_petsc = True
 except :
     have_petsc = False
+try :
+    from . import petsc4pylsys
+    have_petsc4py = True
+except :
+    have_petsc4py = False
 try :
     from . import scipylsys
     have_scipy = True
@@ -58,5 +63,8 @@ def get_linear_system_package(choices=None):
         if have_scipy and choice=="scipy":
             print("using scipy linear system")
             return scipylsys.LinearSystem
+        if have_petsc4py and choice=="petsc4py":
+            print("using petsc4py linear system")
+            return petsc4pylsys.LinearSystem
     raise ValueError("linear system not available")
 

python/mumpslsys.py

@@ -25,6 +25,8 @@ import numpy as np
 import ctypes
 import ctypes.util
 import os
+from ._tools import timeit
+from .csr import CSR
 
 dir_path = os.path.dirname(os.path.realpath(__file__))
 lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
@@ -37,14 +39,61 @@ def _np2c(a,dtype=float,order="C") :
     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 :
-    def __init__(self, elements, n_fields, options=None) :
+    @timeit
+    def __init__(self, elements, n_fields, options=None, constraints=[]) :
         n_nodes = np.max(elements)+1
-        self.globalsize = n_nodes*n_fields
+        self.globalsize = n_nodes*n_fields+len(constraints)
         self.localsize = elements.shape[1]*n_fields
         self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
         idxf = ((elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1])+1).astype(np.int32)
         eltptr = np.arange(1,elements.size*n_fields+2,self.localsize)
+        self.ndof = n_nodes*n_fields
+        num = self.ndof
+        self.n_fields = n_fields
+        self.constraints = constraints
+        for constraint in constraints:
+            idxc = np.c_[np.full(constraint.shape[0], num+1), constraint[:,0]*n_fields+constraint[:,1]+1].flatten()
+            idxf = np.concatenate((idxf,idxc))
+            eltptr = np.concatenate((eltptr, eltptr[-1]+np.arange(2,constraint.shape[0]*2+2,2)))
+            num += 1
         lib2.mumps_element_new.restype = ctypes.c_void_p
         self.mumps_ptr = ctypes.c_void_p(lib2.mumps_element_new(ctypes.c_int(self.globalsize),ctypes.c_int(eltptr.size-1),_np2c(eltptr,np.int32),_np2c(idxf,np.int32)))  
         options = None if options=="" else options
@@ -54,19 +103,27 @@ class LinearSystem :
                 if isinstance(key,str):
                     key = options_key[key]
                     options_vec.extend([key,val])
-                    print(options_vec)
             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))
 
-
-    def local_to_global(self,localv,localm,rhs):
-        rhs[:] = np.bincount(self.idx,localv,self.globalsize)
+    @timeit
+    def local_to_global(self,localv,localm,u, constraints_values):
+        rhs = np.bincount(self.idx,localv,self.globalsize)
+        for i, (c, cv) in enumerate(zip(self.constraints,constraints_values)):
+            rhs[i+self.ndof] = cv[1] + np.sum(u[c[:,0]*self.n_fields+c[:,1]]*cv[0])
         self.v = localm.reshape([-1])
+        for constraint, cv  in constraints_values:
+            localc = np.zeros((constraint.shape[0],4))
+            localc[:,1] = constraint
+            localc[:,2] = constraint
+            self.v = np.concatenate([self.v,localc.flatten()])
+        return rhs
 
+    @timeit
     def solve(self,rhs) :
         r = np.copy(rhs)
         lib2.mumps_element_solve(self.mumps_ptr,_np2c(self.v),_np2c(r))
-        return r
+        return r[:self.ndof]
 
     def __delete__(self):
         lib2.mumps_element_delete(self.mumps_ptr)

python/petsc4pylsys.py

@@ -3,74 +3,41 @@ import sys
 petsc4py.init(sys.argv)
 from petsc4py import PETSc
 import numpy as np
-import atexit
 from ._tools import timeit
-
-def gen_csr(idx) :
-    pairs = np.ndarray([idx.shape[0],idx.shape[1],idx.shape[1]],dtype=([('i0',np.int32),('i1',np.int32)]))
-    pairs['i0'][:,:,:] = idx[:,:,None]
-    pairs['i1'][:,:,:] = idx[:,None,:]
-    pairs, csrmap = np.unique(pairs,return_inverse=True)
-    csr_rowidx = np.hstack([np.array([0],dtype=np.int32),
-                            np.cumsum(np.bincount(pairs["i0"]),
-                           dtype=np.int32)])
-    csr_col = pairs['i1'].reshape(-1)
-    return csr_rowidx, csr_col, csrmap
+from .csr import CSR
 
 class LinearSystemAIJ:
 
-    def __init__(self, elements, n_fields, options, constraints = []) :
+    def __init__(self, elements, n_fields, options, constraints) :
         idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
         self.localsize = elements.shape[1]*n_fields
         self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
-        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1,9]).astype(np.int32)
-        self.idx2 = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,9]).astype(np.int32)
-        options="-ksp_monitor -pc_type lu -info"
+        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
         PETSc.Options().prefixPush("fluid_")
+        PETSc.Options().insertString("-pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1e-16")
         PETSc.Options().insertString(options)
         PETSc.Options().prefixPop()
         self.ksp = PETSc.KSP().create()
         self.ksp.setOptionsPrefix(b"fluid_")
         self.ksp.setFromOptions()
-        self.csr_rowidx, self.csr_col, self.csrmap = gen_csr(idx)
-        self.val = np.zeros(self.csr_col.size)
-        self.mat = PETSc.Mat()
-        self.mat.createAIJWithArrays(self.csr_rowidx.shape[0]-1,(self.csr_rowidx,self.csr_col,self.val),bsize=3)
-
-    #@timeit
-    #def local_to_global(self,localv,localm, u, constraints_value = []):
-    #    msize = self.csr_rowidx.shape[0]-1
-    #    rhs = np.bincount(self.idx,localv,msize)
-    #    self.val[:] = np.bincount(self.csrmap,localm,self.csr_col.size)
-    #    return rhs
+        self.csr = CSR(idx, self.idx,self.constraints)
+        self.val = np.zeros(self.csr.col.size)
+        self.mat = PETSc.Mat().createAIJWithArrays(self.csr.size,(self.csr.row,self.csr.col,self.val),bsize=1)
 
     @timeit
-    def local_to_global2(self,localv,localm, u, constraints_value = []):
-        self.mat.zeroEntries()
-        for e,m in zip(self.idx2,localm.reshape([-1,self.localsize**2])) :
-            self.mat.setValues(e,e,m,PETSc.InsertMode.ADD)
+    def local_to_global(self, localv, localm, u, constraints_value):
+        self.csr.assemble_mat(localm, constraints_value, self.val)
         self.mat.assemble()
-
-    def local_to_global(self,localv,localm, u, constraints_value = []):
-        print(self.idx.shape)
-        msize = self.csr_rowidx.shape[0]-1
-        rhs = np.zeros((msize,))
-        #rhs = np.zeros((self.ndof + len(constraints_value,)))
-        np.add.at(rhs.reshape([-1]),self.idx.reshape([-1]),localv)
-        self.local_to_global2(localv,localm,u,constraints_value)
-        return rhs
+        return self.csr.assemble_rhs(localv, u, constraints_value)
 
     @timeit
     def solve(self,rhs) :
         x = np.ndarray(rhs.shape)
         prhs = PETSc.Vec().createWithArray(rhs.reshape([-1]))
         px = PETSc.Vec().createWithArray(x.reshape([-1]))
-        #self.ksp.getPC().setOperators(None)
-        #self.ksp.getPC().setOperators(self.mat)
-        #self.ksp.setOperators(None)
         self.ksp.setOperators(self.mat)
         self.ksp.solve(prhs,px)
-        return x
+        return x[:self.csr.ndof]
 
 
 class LinearSystemBAIJ :

python/petsclsys.py

@@ -23,15 +23,13 @@ import atexit
 import weakref
 
 from ._tools import timeit
+from .csr import CSR
 
 import ctypes
-import ctypes.util
 import os
 
-print(os.environ["LD_LIBRARY_PATH"])
 petscpath = os.environ["PETSC_DIR"]+"/"+os.environ["PETSC_ARCH"]+"/lib"
 os.environ["LD_LIBRARY_PATH"] += ":"+petscpath
-print(os.environ["LD_LIBRARY_PATH"])
 lib = np.ctypeslib.load_library("libpetsc",petscpath)
 
 lib.PetscInitialize(None, None, None)
@@ -70,9 +68,10 @@ class KSP(PetscObj) :
 
     def __init__(self, options) :
         super().__init__("KSPCreate","KSPDestroy",COMM_WORLD)
+
         lib.PetscOptionsPrefixPush(None, b"fluid_")
+        lib.PetscOptionsInsertString(None, b"-pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1e-16")
         lib.PetscOptionsInsertString(None, options.encode())
-        lib.PetscOptionsInsertString(None, b"-pc_type lu -ksp_monitor")
         lib.PetscOptionsPrefixPop(None)
         lib.KSPSetOptionsPrefix(self,b"fluid_")
         lib.KSPSetFromOptions(self)
@@ -104,25 +103,6 @@ class MatSeqBAIJ(PetscObj) :
         lib.MatAssemblyEnd(self,0)
 
 
-class MatSeqAIJ(PetscObj) :
-
-    def __init__(self, csr) :
-        super().__init__("MatCreate","MatDestroy",COMM_WORLD)
-        lib.MatSetType(self, b"seqaij")
-        lib.MatSetSizes(self,
-                ctypes.c_int(csr.size),ctypes.c_int(csr.size),
-                ctypes.c_int(csr.size),ctypes.c_int(csr.size))
-        lib.MatSeqAIJSetPreallocationCSR(self, _np2c(csr.row,np.int32), _np2c(csr.col,np.int32), None)
-
-    def set(self, v) :
-        aptr = ctypes.POINTER(ctypes.c_double)()
-        lib.MatSeqAIJGetArray(self, ctypes.byref(aptr))
-        np.ctypeslib.as_array(aptr,shape=(v.size,))[:] = v
-        lib.MatSeqAIJRestoreArray(self, ctypes.byref(aptr))
-        lib.MatAssemblyBegin(self,0)
-        lib.MatAssemblyEnd(self,0)
-
-
 class Vec(PetscObj) :
 
     def __init__(self, v) :
@@ -130,53 +110,6 @@ class Vec(PetscObj) :
                 ctypes.c_int(v.shape[-1]), ctypes.c_int(v.size),ctypes.c_void_p(v.ctypes.data))
 
 
-class CSR :
-
-    def __init__(self, idx, rhsidx,constraints):
-        pairs = np.ndarray([idx.shape[0],idx.shape[1],idx.shape[1]],dtype=([('i0',np.int32),('i1',np.int32)]))
-        pairs['i0'][:,:,:] = idx[:,:,None]
-        pairs['i1'][:,:,:] = idx[:,None,:]
-        pairs = pairs.reshape([-1])
-        allpairs = [pairs.reshape([-1])]
-        num = np.max(idx)
-        self.constraints = constraints
-        for c in constraints :
-            num += 1
-            pairs = np.ndarray([c.size*2],dtype=([('i0',np.int32),('i1',np.int32)]))
-            pairs['i0'][:c.size] = c
-            pairs['i1'][:c.size] = num
-            pairs['i0'][c.size:] = num
-            pairs['i1'][c.size:] = c
-            allpairs.append(pairs)
-        pairs = np.vstack(allpairs)
-        pairs, pmap = np.unique(pairs,return_inverse=True)
-        self.map = []
-        count = 0
-        for p in allpairs :
-            self.map.append(pmap[count:count+p.size])
-            count += p.size
-        self.row = np.hstack([np.array([0],dtype=np.int32),
-                              np.cumsum(np.bincount(pairs["i0"]),
-                              dtype=np.int32)])
-        self.col = pairs['i1']
-        self.size = self.row.size-1
-        self.rhsidx = rhsidx
-
-    def assemble_rhs(self, localv, u, constraints_value) :
-        rhs = np.bincount(self.rhsidx,localv,self.size)
-        for i, (c, cv) in enumerate(zip(self.constraints, constraints_value)):
-            localm = np.hstack((localm, cv[0], cv[0]))
-            rhs[i+self.ndof] = cv[1] + np.sum(u[c]*cv[0])
-        return rhs
-
-    def assemble_mat(self, localm, constraints_value) :
-        m = np.bincount(self.map[0],localm,self.col.size)
-        for cmap, cv in zip (self.map[1:], constraints_value) :
-            m += np.bincount(np.hstack(cmap, [cv[0],cv[0]], self.col.size))
-        return m
-
-
-
 class LinearSystemAIJ:
 
     def __init__(self, elements, n_fields, options, constraints = []) :
@@ -185,17 +118,22 @@ class LinearSystemAIJ:
         self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
         self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
         self.ksp = KSP(options)
-        self.csr = CSR(idx, self.idx,[])
-        self.mat = MatSeqAIJ(self.csr)
+        self.csr = CSR(idx, self.idx,self.constraints)
+        self.val = np.zeros(self.csr.col.size)
+        self.mat = PetscObj("MatCreateSeqAIJWithArrays", "MatDestroy", COMM_WORLD,
+                ctypes.c_int(self.csr.size), ctypes.c_int(self.csr.size), 
+                _np2c(self.csr.row,np.int32), _np2c(self.csr.col,np.int32), _np2c(self.val))
 
     @timeit
     def local_to_global(self,localv,localm, u, constraints_value = []):
-        self.mat.set(self.csr.assemble_mat(localm, constraints_value))
+        self.csr.assemble_mat(localm, constraints_value, self.val)
+        lib.MatAssemblyBegin(self.mat,0)
+        lib.MatAssemblyEnd(self.mat,0)
         return self.csr.assemble_rhs(localv, u , constraints_value)
 
     @timeit
     def solve(self,rhs) :
-        return self.ksp.solve(self.mat, rhs.reshape([-1])).reshape(rhs.shape)
+        return self.ksp.solve(self.mat, rhs.reshape([-1])).reshape(rhs.shape)[:self.csr.ndof]
 
 
 class LinearSystemBAIJ :
@@ -206,7 +144,7 @@ class LinearSystemBAIJ :
         nnodes = np.max(elements)+1
         nn = elements.shape[1]
         self.ksp = KSP(options)
-        self.csr = CSR(elements,[])
+        self.csr = CSR(elements,elements,[])
         self.mat = MatSeqBAIJ(self.csr.size*n_fields, n_fields, self.csr.row, self.csr.col)
         self.localsize = elements.shape[1]*n_fields
         self.elements = elements
@@ -215,6 +153,7 @@ class LinearSystemBAIJ :
         self.n_fields = n_fields
         csrmap = ((self.csr.map[0]*9)[:,None,None]+np.arange(0,9,3)[None,None,:]+np.arange(0,3)[None,:,None])
         self.csrmap = np.copy(np.swapaxes(csrmap.reshape([elements.shape[0],nn,nn, n_fields,n_fields]),2,3).reshape([-1]))
+        self.val = np.zeros(self.csr.col.size)
 
     @timeit
     def local_to_global(self,localv,localm, u, constraints_value = []):

python/scipylsys.py

@@ -18,13 +18,15 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program (see COPYING and COPYING.LESSER files).  If not,
 # see <http://www.gnu.org/licenses/>.
+
 import numpy as np
 import scipy.sparse
 import scipy.sparse.linalg