fix csr

python/_tools.py

@@ -34,7 +34,6 @@ def timeprint(timers) :
         print(timers)
 
 atexit.register(timeprint, timers)
-from . import petsc4pylsys
 
 try :
     from . import mumpslsys
@@ -57,8 +56,6 @@ try :
 except :
     have_scipy = False
 
-print(have_petsc4py)
-
 def get_linear_system_package(choices=None):
     if choices is None:
         choices = ["mumps","petsc","petsc4py","scipy"]

python/csr.py

@@ -2,15 +2,14 @@ import numpy as np
 import time
 class CSR :
 
-    def __init__(self, idx,constraints):
+    def __init__(self, idx, rhsidx,constraints):
         shift = 2**32
         num = np.max(idx)
         self.ndof = num+1
-        idx = idx.astype('uint64')
-        pairs = (idx[:,None,:]*shift+idx[:,:,None]).flatten()
+        idx = rhsidx.reshape([idx.shape[0],-1]).astype('uint64')
+        pairs = (idx[:,:,None]*shift+idx[:,None,:]).flatten()
         allpairs = [pairs]
         self.constraints = constraints
-        self.idx = idx.reshape([-1]).astype('int64')
         for c in constraints :
             num += 1
             pairs = np.ndarray([c.size*2+1],np.uint64)
@@ -30,9 +29,10 @@ class CSR :
                               dtype=np.int32)])
         self.col = (pairs%shift).astype(np.int32)
         self.size = self.row.size-1
+        self.rhsidx = rhsidx
 
     def assemble_rhs(self, localv, u, constraints_value) :
-        rhs = np.bincount(self.idx,localv,self.size)
+        rhs = np.bincount(self.rhsidx,localv,self.size)
         for i, (c, cv) in enumerate(zip(self.constraints, constraints_value)):
             rhs[i+self.ndof] = cv[1] + np.sum(u[c]*cv[0])
         return rhs

python/petsc4pylsys.py

@@ -5,10 +5,10 @@ from petsc4py import PETSc
 import numpy as np
 from ._tools import timeit
 from .csr import CSR
-try :
-    from .dof_numbering import MeshCl, LagrangeSpace
-except :
-    print("dof_numbering not available !!!")
+# try :
+    # from .dof_numbering import MeshCl, LagrangeSpace
+# except :
+    # print("dof_numbering not available !!!")
 class LinearSystemAIJ:
 
     def __init__(self, idx, n_fields, options, constraints) :
@@ -17,6 +17,10 @@ class LinearSystemAIJ:
         self.idx = idx
         self.localsize = idx.shape[1]
         self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
+        elements = np.asarray(idx[:,:self.localsize//n_fields]/n_fields, np.int32)
+        idx_matrix = ((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.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
         # self.space_v = LagrangeSpace(MeshCl(elements), order=1).element_nodes
         # self.space_p = LagrangeSpace(MeshCl(elements), order=1).element_nodes
@@ -30,6 +34,7 @@ class LinearSystemAIJ:
         # pmap = dimension*nmax_v + self.space_p
         # self.idx = np.concatenate([self.idx, pmap], axis=1)
         
+        
         PETSc.Options().prefixPush("fluid_")
         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)
@@ -38,7 +43,7 @@ class LinearSystemAIJ:
         self.ksp = PETSc.KSP().create()
         self.ksp.setOptionsPrefix(b"fluid_")
         self.ksp.setFromOptions()
-        self.csr = CSR(self.idx,self.constraints)
+        self.csr = CSR(idx_matrix,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)
 
@@ -47,7 +52,6 @@ class LinearSystemAIJ:
     def local_to_global(self, localv, localm, u, constraints_value):
         self.csr.assemble_mat(localm, constraints_value, self.val)
         self.mat.assemble()
-        self.mat.view()
         return self.csr.assemble_rhs(localv, u, constraints_value)
 
     @timeit

validation/mbtests.py

@@ -5,7 +5,7 @@ all_tests = [
         # "inclined2d.inclined2d",
         # "inclined3d.inclined3d",
         "poiseuille.poiseuille",
-        # "cavity.cavity",
+        "cavity.cavity",
         # "oneGrain.oneGrain",
         # "darcy.darcy",
         # "periodic2d.periodic2d",

validation/poiseuille/poiseuille.py

@@ -27,6 +27,7 @@ import os
 import subprocess
 import time
 import shutil
+import random
 import unittest
 
 class Poiseuille(unittest.TestCase) :
@@ -39,7 +40,7 @@ class Poiseuille(unittest.TestCase) :
         if not os.path.isdir(outputdir) :
             os.makedirs(outputdir)
 
-        subprocess.call(["gmsh", "-2", "mesh_naive.geo","-clscale","2"])
+        subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","2"])
         t = 0
         ii = 0
 
@@ -49,31 +50,28 @@ class Poiseuille(unittest.TestCase) :
         rho = 1000                                      # fluid density
         nu = 1e-3                                   # kinematic viscosity
         mu = nu*rho                                     # dynamic viscosity
-        tEnd = 1000                                     # final time
+        tEnd = 100000                                     # final time
 
         #numerical parameters
         lcmin = .1                                  # mesh size
-        dt = 1                                       # time step
+        dt = 10                                       # time step
         alpha = 1e-4                                    # stabilization coefficient
         epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
         print('epsilon',epsilon)
 
-        shutil.copy("mesh.msh", outputdir +"/mesh_naive.msh")
+        shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
         outf = 1                                       # number of iterations between output files
 
         #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
         fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,solver="petsc4py")
-        fluid.load_msh("mesh_naive.msh")
-
-        # fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
-        # fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
-        fluid.set_open_boundary("Left",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
+        fluid.load_msh("mesh.msh")
+        fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
+        fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
         fluid.set_wall_boundary("Bottom",velocity=[0,0])
         fluid.set_wall_boundary("Top",velocity=[0,0])
-        fluid.set_open_boundary("Right",pressure=0)
-        # fluid.set_open_boundary("RightUp",pressure=0)
-        # fluid.set_open_boundary("RightDown",pressure=0)
+        fluid.set_open_boundary("RightUp",pressure=0)
+        fluid.set_open_boundary("RightDown",pressure=0)
 
 
         ii = 0
@@ -89,7 +87,7 @@ class Poiseuille(unittest.TestCase) :
         tic = time.perf_counter()
         while ii < 100 : 
             #Fluid solver
-            time_integration.iterate(fluid,None,dt)
+            time_integration.iterate(fluid,None,dt,check_residual_norm=1e-5)
             t += dt
             #Output files writting
             if ii %outf == 0 :
@@ -97,7 +95,6 @@ class Poiseuille(unittest.TestCase) :
                 fluid.write_vtk(outputdir, ioutput, t)
             ii += 1
             print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.perf_counter() - tic))
-            exit(0)
             
         s = fluid.solution()
         x = fluid.coordinates()