wip

fluid/fluid_problem_fem.h

@@ -31,6 +31,7 @@ void fe_fields_eval_grad(const FEFields *fields, const double *xi, const double
 void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]);
 void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, const double f1[][D], const double *sf, const double dsf[][D], double jw, double *local_vector);
 static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *f00, const double f01[][D], const double f10[][D], const double f11[][D][D], const double *sf, const double dsf[][D], double jw, double *local_matrix) {
+  printf("ndof_local = %6d\n", fields->local_size); // dof_local
   int jc = 0;
   int n_fields = fields->n;
   for (int jfield = 0; jfield < fields->n; ++jfield) {
@@ -39,6 +40,7 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
       for (int ifield = 0; ifield < fields->n; ++ifield){
         for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
           double m = 0;
+          // int index = ifield*fields->n+jfield;
           if (f00) {
             m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
           }
@@ -59,13 +61,23 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
               m += jw*sf[ic]*dsf[jc][jd]*f01[ifield*fields->n+jfield][jd];
             }
           }
-          local_matrix[(ifield*n_fields+iphi)*n_fields*n_fields+jfield*n_fields+jphi] += m; // to change
+          // local_matrix[] += m; // elt*dof_local 
+          local_matrix[(ifield*n_fields+iphi)*n_fields*n_fields+jfield*n_fields+jphi] += m;
           ic ++;
         }
       }
       jc ++;
     }
   }
+  // for(int jfield=0; jfields < fields->n; ++jfield){
+    // for(int jphi=0; jphi < fields->element[jfield]->nlocal; ++jphi){
+      // for(int ifields=0; ifields < fields->n; ++ifields){
+        // for(int iphi=0; iphi < fields->element[ifields]->nlocal; ++iphi){
+          // local_matrix[]
+        // }
+      // }
+    // }
+  // }
 }
 void fe_fields_free(FEFields *fields);
 

python/_tools.py

@@ -57,6 +57,8 @@ 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,7 +2,7 @@ import numpy as np
 import time
 class CSR :
 
-    def __init__(self, idx, rhsidx,constraints):
+    def __init__(self, idx,constraints):
         shift = 2**32
         num = np.max(idx)
         self.ndof = num+1
@@ -10,7 +10,7 @@ class CSR :
         pairs = (idx[:,:,None]*shift+idx[:,None,:]).flatten()
         allpairs = [pairs]
         self.constraints = constraints
-        self.idx = idx.reshape([-1])
+        self.idx = idx.reshape([-1]).astype('int64')
         for c in constraints :
             num += 1
             pairs = np.ndarray([c.size*2+1],np.uint64)
@@ -32,12 +32,30 @@ class CSR :
         self.size = self.row.size-1
 
     def assemble_rhs(self, localv, u, constraints_value) :
+        localv_print = localv.reshape([4,-1])
+        ii = 0
+        for elt in range(4):
+            print("-------- elt : ------------ ", elt)
+            for i in range(9):
+                print("\t:", "%3d"%self.idx[ii], "->%2.4f" % localv_print[elt,i], end="")
+                ii += 1
+            print("")
+
         rhs = np.bincount(self.idx,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
 
     def assemble_mat(self, localm, constraints_value, v) :
+        localm_print = localm.reshape([4,9,-1])
+        ii = 0
+        for elt in range(4):
+            print("-------- elt : ------------ ", elt)
+            for i in range(9):
+                for j in range(9):
+                    print("\t:", "%4d"%self.map[0][ii], "->%2.5f" % localm_print[elt,i,j], end="")
+                    ii += 1
+                print("")
         v[:] = np.bincount(self.map[0],localm,self.col.size)
         for cmap, cv in zip (self.map[1:], constraints_value) :
             v += np.bincount(cmap, np.concatenate((cv[0],cv[0],[0])), self.col.size)

python/petsc4pylsys.py

@@ -122,3 +122,5 @@ class LinearSystemBAIJ :
             return x
         else:
             return x[:self.csr.ndof]
+
+LinearSystem = LinearSystemAIJ
\ No newline at end of file

validation/poiseuille/poiseuille.py

@@ -39,7 +39,7 @@ class Poiseuille(unittest.TestCase) :
         if not os.path.isdir(outputdir) :
             os.makedirs(outputdir)
 
-        subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","2"])
+        subprocess.call(["gmsh", "-2", "mesh_naive.geo","-clscale","2"])
         t = 0
         ii = 0
 
@@ -58,20 +58,22 @@ class Poiseuille(unittest.TestCase) :
         epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
         print('epsilon',epsilon)
 
-        shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+        shutil.copy("mesh.msh", outputdir +"/mesh_naive.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)
-        fluid.load_msh("mesh.msh")
+        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("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.set_wall_boundary("Bottom",velocity=[0,0])
         fluid.set_wall_boundary("Top",velocity=[0,0])
-        fluid.set_open_boundary("RightUp",pressure=0)
-        fluid.set_open_boundary("RightDown",pressure=0)
+        fluid.set_open_boundary("Right",pressure=0)
+        # fluid.set_open_boundary("RightUp",pressure=0)
+        # fluid.set_open_boundary("RightDown",pressure=0)
 
 
         ii = 0