oscillate

python/marblesbag/fluid.py

@@ -36,7 +36,7 @@ class StrongBoundary(Structure):
 
 class fluid_problem :
 
-    def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries, noGrains, n_fluids):
+    def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries, n_fluids=1):
         nsb = len(strong_boundaries)
         class Bnd :
             def __init__(self, b) :
@@ -56,7 +56,7 @@ class fluid_problem :
         asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2],BNDCB(Bnd(i[3]).apply)) for i in strong_boundaries])
         self.asb = asb
         lib.fluid_problem_new.restype = c_void_p
-        self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), np2c(mu), np2c(rho), c_double(epsilon), nsb, asb, c_int(noGrains), c_int(n_fluids)))
+        self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), np2c(mu), np2c(rho), c_double(epsilon), nsb, asb, c_int(n_fluids)))
         if self._fp == None :
             raise NameError("cannot create fluid problem")
 
@@ -93,6 +93,24 @@ class fluid_problem :
         self.n_particles = mass.shape[0]
         lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(init))
 
+    def _get_matrix(self, f_name, nrow, ncol) :
+        f = getattr(lib,"fluid_problem_"+f_name)
+        f.restype = POINTER(c_double)
+        return numpy.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
+
+    def solution(self) :
+        return self._get_matrix("solution", self.n_nodes(), self.n_fields())
+
+    def coordinates(self) :
+        return self._get_matrix("coordinates",self.n_nodes(), 3)
+    
+    def n_fields(self):
+        lib.fluid_problem_n_fields.restype = c_int
+        return lib.fluid_problem_n_fields(self._fp);
+
+    def n_nodes(self):
+        lib.fluid_problem_n_nodes.restype = c_int
+        return lib.fluid_problem_n_nodes(self._fp);
 
     def particle_element_id(self) :
         f = getattr(lib,"fluid_problem_particle_element_id")

src/fluid_problem.c

@@ -295,7 +295,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
   int P = problem->n_fluids*(D+1);
   double p = sol[P];
   double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
-  double epsilon = problem->epsilon;
+  double epsilonp = problem->epsilon;
+  double epsilonq = problem->epsilon/10;
   f0[P] = 1-c;
   f1[P*D+0] = 0;
   f1[P*D+1] = 0;
@@ -311,7 +312,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
     double uold[D] = {solold[U],solold[U+1]};
     double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
 
-    double oq = 1./q;//1./fmax(q,1e-3);
+    double oq = 1./fmax(q,0.1);
     double divu = 0;
     double tau[D][D];
     double utau[D]={0.};
@@ -325,12 +326,11 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
     f0[Q] = divu+(q-qold)/dt;
     f0[P] -= q;
     for (int i = 0; i < D; ++i) {
-      f0[Q] += epsilon*dp[i]*dq[i];
       f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq*divu + utau[i]*oq) - p*dq[i];
       for (int j = 0; j < D; ++j) {
         f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
       }
-      f1[Q*D+i] = epsilon*dp[i]*q;
+      f1[Q*D+i] = epsilonp*dp[i]*fmax(q,0.1)+epsilonq*dq[i];
     }
   }
 }
@@ -347,7 +347,6 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
   
   const double *solution = problem->solution;
   const double *compacity = problem->compacity;
-  const double *old_compacity = problem->old_compacity;
   const double epsilon = problem->epsilon;
   size_t local_size = N_SF*n_fields;
   size_t N_E = mesh->n_elements;
@@ -622,14 +621,14 @@ static void fluid_problem_compute_compacity(FluidProblem *problem)
   }
 }
 
-FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int noGrains, int n_fluids)
+FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int n_fluids)
 {
   static int initialized = 0;
   if (!initialized) {
     hxtInitializeLinearSystems(0, NULL);
     initialized = 1;
 #ifdef HXT_HAVE_PETSC
-    hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30 -ksp_monitor", NULL);
+    hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
 #endif
   }
   FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
@@ -658,28 +657,15 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
     problem->strong_boundaries[i].row = strong_boundaries[i].row;
   }
   problem->compacity = (double*)malloc(mesh->n_nodes*sizeof(double));
-  problem->old_compacity = (double*)malloc(mesh->n_nodes*sizeof(double));
-  if(noGrains){
-    for(int i = 0; i < mesh->n_nodes; ++i){
-      problem->compacity[i] = 0.;
-      problem->old_compacity[i] = 0.;
-    }
+  for(int i = 0; i < mesh->n_nodes; ++i){
+    problem->compacity[i] = 0.;
   }
   problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
-  // begin to remove
   for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
     problem->solution[i] = 0.;
   }
-  for (int i = 0; i < mesh->n_nodes; ++i){
-    for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
-      double *x = mesh->x + i*3;
-      double a = (x[0]+5)/10;
-      problem->solution[i*n_fields+D+ifluid*(D+1)] = ifluid == 0 ? 0.8-0.6*a:0.2+0.6*a;
-    }
-  }
   problem->node_volume = NULL;
   fluid_problem_compute_node_volume(problem);
-  // end to remove
   problem->linear_system = NULL;
 #ifdef HXT_HAVE_PETSC
     hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
@@ -703,7 +689,6 @@ void fluid_problem_free(FluidProblem *problem)
   free(problem->rho);
   free(problem->solution);
   free(problem->compacity);
-  free(problem->old_compacity);
   hxtLinearSystemDelete(&problem->linear_system);
   free(problem->particle_uvw);
   free(problem->particle_element_id);
@@ -941,12 +926,9 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
   problem->solution = new_solution;
   free(problem->compacity);
   problem->compacity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
-  free(problem->old_compacity);
-  problem->old_compacity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
   for (int i = 0; i < new_mesh->n_nodes; ++i)
   {
     problem->compacity[i] = 0.;
-    problem->old_compacity[i] = 0.;
   }
   mesh_free(problem->mesh);
   problem->mesh = new_mesh;
@@ -970,8 +952,6 @@ void fluid_problem_after_import(FluidProblem *problem)
   mesh_tree_free(problem->mesh_tree);
   problem->mesh_tree = mesh_tree_create(problem->mesh);
   int n_fields = fluid_problem_n_fields(problem);
-  free(problem->old_compacity);
-  problem->old_compacity = (double*)malloc(sizeof(double)*problem->mesh->n_nodes);
   for (int i = 0; i < problem->n_particles; ++i)
     problem->particle_element_id[i] = -1;
   mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
@@ -1032,9 +1012,6 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
   }
   clock_get(&toc);
   printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
-  for (int i = 0; i < problem->mesh->n_nodes; ++i){
-    problem->old_compacity[i] = problem->compacity[i];
-  }
   clock_get(&tic);
   fluid_problem_compute_compacity(problem);
   clock_get(&toc);
@@ -1050,3 +1027,18 @@ int fluid_problem_n_particles(FluidProblem *problem)
 {
   return problem->n_particles;
 }
+
+double *fluid_problem_solution(const FluidProblem *p)
+{
+  return p->solution;
+}
+
+double *fluid_problem_coordinates(const FluidProblem *p)
+{
+  return p->mesh->x;
+}
+
+int fluid_problem_n_nodes(const FluidProblem *p)
+{
+  return p->mesh->n_nodes;
+}

src/fluid_problem.h

@@ -49,7 +49,6 @@ typedef struct {
   Mesh *mesh;
   MeshTree *mesh_tree;
   double *compacity;
-  double *old_compacity;
   double *solution;
   double *solution_explicit;
   double *node_volume;
@@ -75,10 +74,13 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
 int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
 void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init);
 void fluid_problem_free(FluidProblem *problem);
-FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int noGrains, int n_fluids);
+FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int n_fluids);
 void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el);
 int *fluid_problem_particle_element_id(FluidProblem *problem);
 int fluid_problem_n_particles(FluidProblem *problem);
 void fluid_problem_after_import(FluidProblem *problem);
 int fluid_problem_n_fields(const FluidProblem *p);
+double *fluid_problem_solution(const FluidProblem *p);
+double *fluid_problem_coordinates(const FluidProblem *p);
+int fluid_problem_n_nodes(const FluidProblem *p);
 #endif

testcases/poiseuille2fluids/poiseuille.py

@@ -67,31 +67,36 @@ strong_boundaries = [
         ("Bottom",v0,v0,0.),
         ("Bottom",u1,u1,0),
         ("Bottom",v1,v1,0.),
-        #("Bottom",q0,q0,0.6),
         ("Top",u0,u0,0),
         ("Top",v0,v0,0.),
         ("Top",u1,u1,0),
-        #("Top",q0,q0,0.6),
         ("Top",v1,v1,0.),
-        ("LeftUp",u0,u0,lambda x : 0.8/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftUp",u0,u0,lambda x : 1./(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
         ("LeftUp",v0,v0,0),
-        ("LeftUp",q0,q0,0.8),
-        ("LeftUp",u1,u1,lambda x : 0.2/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftUp",q0,q0,0.999),
+        ("LeftUp",u1,u1,lambda x : 0/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
         ("LeftUp",v1,v1,0),
-        ("LeftDown",u0,u0,lambda x : 0.8/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftDown",u0,u0,lambda x : 1./(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
         ("LeftDown",v0,v0,0),
-        ("LeftDown",q0,q0,0.8),
-        ("LeftDown",u1,u1,lambda x : 0.2/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftDown",q0,q0,0.999),
+        ("LeftDown",u1,u1,lambda x : 0/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
         ("LeftDown",v1,v1,0),
         ("RightDown",v0,v0,0),
         ("RightDown",v1,v1,0),
         ("RightUp",v0,v0,0),
         ("RightUp",v1,v1,0),
         ]
-fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,1,2)
+
+fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,2)
 ii = 0
 t = 0
 
+#set initial_condition
+x = fluid.coordinates()
+s = fluid.solution()
+alpha = (x[:,0]+5)/10
+s[:,2] = 1-alpha*0.8
+s[:,5] = 1-s[:,2]
 
 fluid.export_vtk(outputdir,0,0)