petsc4py linear system

fluid/fluid_problem.c

@@ -1248,7 +1248,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
   }
 }
 
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, const char *petsc_solver_type, int drag_in_stab, double drag_coeff_factor, int temporal, int advection)
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, int temporal, int advection)
 {
   if (n_fluids != 1 && n_fluids != 2) {
     printf("error : only 1 or 2 fluids are supported\n");

fluid/fluid_problem.h

@@ -107,7 +107,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
 void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
 void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
 void fluid_problem_free(FluidProblem *problem);
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, const char *petsc_solver_type, int drag_in_stab, double drag_coeff_factor,int temporal, int advection);
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor,int temporal, int advection);
 void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
 void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume);
 int *fluid_problem_particle_element_id(FluidProblem *problem);

hxtlsys/hxt_linear_system_petsc.c

@@ -401,6 +401,7 @@ HXTStatus hxtLinearSystemPETScMapFromVec(HXTLinearSystemPETSc *system, Vec vec,
   HXT_PETSC_CHECK(VecRestoreArrayRead(vec, &vv));
   return HXT_STATUS_OK;
 }
+
 HXTStatus hxtLinearSystemPETScLocalToGlobal(HXTLinearSystemPETSc *sys, double *all_local_vector, double *all_local_matrix, double *rhs) {
   int local_size = sys->nNodesByElement*sys->nFields;
   for (int iel=0; iel <sys->nElements; ++iel) {

python/CMakeLists.txt

@@ -28,6 +28,8 @@ set(SRC
   scontact.py
   time_integration.py
   hxtlsys.py
+  petsclsys.py
+  scipylsys.py
   VTK.py
 )
 

python/fluid.py

@@ -34,7 +34,11 @@ import signal
 import os
 import sys
 from . import VTK
-from .hxtlsys import LinearSystem
+try :
+    from .petsclsys import LinearSystem
+except :
+    print("PETSc4py not available, falling back to scipy linear system")
+    from .scipylsys import LinearSystem
 
 dir_path = os.path.dirname(os.path.realpath(__file__))
 lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
@@ -69,10 +73,11 @@ def _np2c(a,dtype=float) :
     r.tmp = tmp
     return r
 
+    
 class FluidProblem :
     """Create numerical representation of the fluid."""
 
-    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="-pc_type lu", drag_in_stab=1, drag_coefficient_factor=1, temporal=True, advection=True):
+    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, temporal=True, advection=True):
         """Build the fluid structure.
 
         Keyword arguments:
@@ -93,8 +98,10 @@ class FluidProblem :
         ValueError -- if the dimension differs from 2 or 3
         NameError -- if C builder cannot be found
         """
+        self.solver_options = petsc_solver_type
         self.strong_cb_ref = []
         self.weak_cb_ref = []
+        self.sys = None
         if dim == 2 :
             self._lib = lib2
         elif dim == 3 :
@@ -105,7 +112,7 @@ class FluidProblem :
         n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
         self._n_fluids = n_fluids
         self._dim = dim
-        self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), petsc_solver_type.encode(),c_int(drag_in_stab),c_double(drag_coefficient_factor),c_int(temporal), c_int(advection)))
+        self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_int(temporal), c_int(advection)))
         if self._fp == None :
             raise NameError("Cannot create fluid problem.")
 
@@ -121,6 +128,7 @@ class FluidProblem :
         mesh_file_name -- name of the mesh.msh file containing information about the domain
         """
         self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
+        self.sys = None
 
     def set_wall_boundary(self, tag, pressure=None, velocity=None) :
         """Set the weak boundaries (=normal fluxes) for the fluid problem.
@@ -223,6 +231,7 @@ class FluidProblem :
         cmin -- optional argument to weigh minimum gradient used in the adaptation formula
         """
         self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
+        self.sys = None
 
     def _mesh_boundaries(self) :
         n = self._lib.fluid_problem_mesh_n_boundaries(self._fp)
@@ -295,6 +304,7 @@ class FluidProblem :
         self.porosity()[:] = data["porosity"]
         self.old_porosity()[:] = data["old_porosity"]
         self._lib.fluid_problem_after_import(self._fp)
+        self.sys = None
 
     def compute_node_force(self, dt) :
         """Compute the forces to apply on each grain resulting from the fluid motion.
@@ -334,20 +344,24 @@ class FluidProblem :
         rhs = np.zeros_like(solution_old)
         self._lib.fluid_problem_reset_boundary_force(self._fp)
         self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
-        sys = LinearSystem(self.elements(),self.n_fields())
+        if self.sys is None :
+            self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options
+                    )
+        sys = self.sys
         localv = np.ndarray([sys.localsize*self.n_elements()])
         localm = np.ndarray([sys.localsize**2*self.n_elements()])
         self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
         sys.local_to_global(localv,localm,rhs)
         if check_residual_norm > 0 :
-            norm0 = sys.rhs_norm(rhs)
+            norm0 = np.linalg.norm(rhs)
         self.solution()[:] -= sys.solve(rhs)
         if check_residual_norm > 0 :
             sys.zero_matrix()
             rhs[:] = 0
             self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
             sys.local_to_global(localv,localm,rhs)
-            norm = sys.rhs_norm(rhs)
+            norm = np.linalg.norm(rhs)
+            print("norm",norm)
             if norm > check_residual_norm :
                 raise ValueError("wrong derivative or linear system precision")
         self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)

python/hxtlsys.py

+from ctypes import *
+import numpy as np
+import os
+
+_dir_path = os.path.dirname(os.path.realpath(__file__))
+_lib = np.ctypeslib.load_library("libhxtlsys",_dir_path)
+
+def _np2c(a,dtype=float) :
+    tmp = np.require(a,dtype,"C")
+    r = c_void_p(tmp.ctypes.data)
+    r.tmp = tmp
+    return r
+
+class LinearSystem :
+    _initialized = False
+    def __init__(self, elements, n_fields) :
+        if not LinearSystem._initialized :
+             LinearSystem._initialized = True
+             _lib.hxtInitializeLinearSystems(c_int(0),None)
+        _lib.hxtLinearSystemCreatePETSc.restype = c_void_p
+        _lib.hxtPETScInsertOptions(b"-pc_type lu",b"fluid")
+        self._sysp = c_void_p()
+        _lib.hxtLinearSystemCreatePETSc(byref(self._sysp),c_int(elements.shape[0]),c_int(elements.shape[1]),c_int(n_fields),_np2c(elements,np.int32),b"fluid")
+        self.localsize = n_fields*elements.shape[1]
+        self._elements = elements
+        self.zero_matrix()
+
+    def zero_matrix(self) :
+        _lib.hxtLinearSystemZeroMatrix(self._sysp)
+
+    def local_to_global(self,localv,localm,rhs):
+        _lib.hxtLinearSystemLocalToGlobal(self._sysp,_np2c(localv),_np2c(localm),_np2c(rhs))
+
+    def rhs_norm(self,rhs) :
+        norm = c_double()
+        _lib.hxtLinearSystemGetRhsNorm(self._sysp,_np2c(rhs),byref(norm));
+        return norm.value
+
+    def solve(self,rhs) :
+        sol = np.ndarray(rhs.shape)
+        _lib.hxtLinearSystemSolve(self._sysp,_np2c(rhs),_np2c(sol))
+        return sol
+    
+    def __del__(self):
+        _lib.hxtLinearSystemDelete(byref(self._sysp))

python/petsclsys.py

+import petsc4py2
+petsc4py.init()
+from petsc4py import PETSc
+import numpy as np
+
+class LinearSystem :
+    def __init__(self, elements, n_fields, options) :
+        nnodes = np.max(elements)+1
+        nn = elements.shape[1]
+        pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
+        pairs['i0'] = elements[:,np.repeat(range(nn),nn-1)]
+        idx = np.hstack([[i for i in range(nn) if i!=j] for j in range(nn)]) 
+        pairs['i1'] = elements[:,idx]
+        pairs = np.unique(pairs)
+        nnz = (np.bincount(pairs["i0"])+1).astype(np.int32)
+        PETSc.Options().prefixPush("fluid_")
+        PETSc.Options().insertString(options)
+        PETSc.Options().prefixPop()
+        self.mat = PETSc.Mat().createBAIJ(nnodes*n_fields,n_fields,nnz)
+        self.ksp = PETSc.KSP().create()
+        self.ksp.setOptionsPrefix(b"fluid_")
+        self.ksp.setFromOptions()
+        self.mat.zeroEntries()
+        ###
+        self.localsize = elements.shape[1]*n_fields
+        self.elements = elements
+        self.idx = (self.elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
+
+    def local_to_global(self,localv,localm,rhs):
+        self.mat.zeroEntries()
+        np.add.at(rhs.reshape([-1]),self.idx,localv)
+        for e,m in zip(self.elements,localm.reshape([-1,self.localsize**2])) :
+            self.mat.setValuesBlocked(e,e,m,PETSc.InsertMode.ADD)
+        self.mat.assemble()
+
+    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.setOperators(self.mat)
+        self.ksp.solve(prhs,px)
+        return x

python/scipylsys.py

+import numpy as np
+import scipy.sparse
+import scipy.sparse.linalg
+
+class LinearSystem :
+    _initialized = False
+    def __init__(self, elements, n_fields,options) :
+        localsize = n_fields*elements.shape[1]
+        self.el = elements
+        idx = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,localsize])
+        self.rows = np.repeat(idx[:,:,None],localsize,axis=2)
+        self.rows = self.rows.reshape([-1])
+        self.cols = np.repeat(idx[:,None,:],localsize,axis=1)
+        self.cols = self.cols.reshape([-1])
+        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
+        self.localsize = localsize
+
+    def local_to_global(self,localv,localm,rhs):
+        np.add.at(rhs.reshape([-1]),self.idx,localv)
+        coo = scipy.sparse.coo_matrix((localm.reshape([-1]),(self.rows,self.cols)))
+        self.matrix = coo.tocsc()
+
+    def solve(self,rhs) :
+        return scipy.sparse.linalg.splu(self.matrix).solve(rhs.reshape([-1])).reshape(rhs.shape)
+

validation/darcy/darcy.py

@@ -105,7 +105,7 @@ class Darcy(unittest.TestCase) :
             p = scontact.ParticleProblem(2)
             p.read_vtk(outputdir,0)
             ii = 0
-            fluid = mbfluid.FluidProblem(2,g,mu,rho,drag_in_stab=True)
+            fluid = mbfluid.FluidProblem(2,g,mu,rho,drag_in_stab=True,petsc_solver_type="-pc_type ilu")
             fluid.load_msh("mesh.msh")
             V = 0.00001
             fluid.solution()[:,0] = V