new adapt + testcase

marblesbag/fluid.py

@@ -39,8 +39,8 @@ class fluid_problem :
         if(self._fp != None) :
             lib.fluid_problem_free(self._fp)
 
-    def adapt_mesh(self, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el) :
-        lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(gradPmin), c_double(gradPmax), c_double(lcmin), c_double(n_el))
+    def adapt_mesh(self, lcmax, lcmin, n_el) :
+        lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
 
     def export(self, output_dir, t, it) :
         lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))

marblesbag/fluid3.py

@@ -16,7 +16,7 @@ class StrongBoundary(Structure):
 
 class fluid_problem :
 
-    def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries):
+    def __init__(self, mesh_file_name, g, mu, rho, alpha, strong_boundaries, notComputeEpsilon):
         nsb = len(strong_boundaries)
         class Bnd :
             def __init__(self, b) :
@@ -32,7 +32,7 @@ class fluid_problem :
         asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).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), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb))
+        self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(alpha), nsb, asb, c_int(notComputeEpsilon)))
         if self._fp == None :
             raise NameError("cannot create fluid problem")
 
@@ -40,8 +40,8 @@ class fluid_problem :
         if(self._fp != None) :
             lib.fluid_problem_free(self._fp)
 
-    def adapt_mesh(self, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el) :
-        lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(gradPmin), c_double(gradPmax), c_double(lcmin), c_double(n_el))
+    def adapt_mesh(self, lcmax, lcmin, n_el, alpha, notComputeEpsilon) :
+        lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(alpha), c_int(notComputeEpsilon))
 
     def export(self, output_dir, t, it) :
         lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))

marblesbag/src/fluid_problem.c

@@ -574,7 +574,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
   const double *old_porosity = problem->old_porosity;
   const double mu = problem->mu;
   const double rho = problem->rho;
-  const double epsilon = problem->epsilon;
+  const double *epsilon = problem->epsilon;
   size_t local_size = N_SF*n_fields;
   size_t N_E = mesh->n_elements;
   double *all_local_vector = malloc(N_E*local_size*sizeof(double));
@@ -643,7 +643,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
           dphidp += dphii[k]*dp[k];
           divu += du[k][k];
         }
-        local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon*dphidp+phii*(divu+dcdt));
+        local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon[el[iphi]]*dphidp+phii*(divu+dcdt));
         for (int jphi = 0; jphi < N_SF; ++jphi) {
           double phij = phi[jphi];
           const double *dphij = dphi[jphi];
@@ -673,7 +673,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
             }
             LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
           }
-          LOCAL_MATRIX(P,P) +=  jw*epsilon*dphiidphij;
+          LOCAL_MATRIX(P,P) +=  jw*epsilon[el[iphi]]*dphiidphij;
         }
       }
     }
@@ -813,7 +813,7 @@ static void fluid_problem_compute_porosity(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)
+FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, int n_strong_boundaries, StrongBoundary *strong_boundaries, int notComputeEpsilon)
 {
   static int initialized = 0;
   if (!initialized) {
@@ -828,7 +828,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
   if (!mesh)
     return NULL;
   problem->rho = rho;
-  problem->epsilon = epsilon;
+
   problem->mu = mu;
   problem->g = g;
   problem->mesh = mesh;
@@ -840,6 +840,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
     problem->strong_boundaries[i].field = strong_boundaries[i].field;
   }
   problem->porosity = malloc(mesh->n_nodes*sizeof(double));
+  problem->epsilon = malloc(mesh->n_nodes*sizeof(double));
   problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
   problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
   // begin to remove
@@ -848,6 +849,21 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
   }
   problem->node_volume = malloc(0);
   fluid_problem_compute_node_volume(problem);
+#if DIMENSION == 2
+  for (int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
+#else
+  for (int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
+    printf("eps = %e\n",problem->epsilon[i]);
+  }
+#endif
+  if(notComputeEpsilon){
+    for (int i = 0; i < problem->mesh->n_nodes; ++i){
+      problem->epsilon[i] = alpha;
+      printf("eps = %e\n",problem->epsilon[i]);
+    }
+  }
   // end to remove
   problem->linear_system = NULL;
 #ifdef HXT_HAVE_PETSC
@@ -870,6 +886,7 @@ void fluid_problem_free(FluidProblem *problem)
 {
   free(problem->solution);
   free(problem->porosity);
+  free(problem->epsilon);
   free(problem->old_porosity);
   hxtLinearSystemDelete(&problem->linear_system);
   free(problem->particle_uvw);
@@ -885,7 +902,7 @@ void fluid_problem_free(FluidProblem *problem)
   mesh_free(problem->mesh);
 }
 
-void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
+void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
 {
   double ngradM = 0.;
   Mesh *mesh = problem->mesh;
@@ -896,9 +913,70 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
     new_mesh_size[i] = 0.;
     num_lc[i] = 0.;
   }
-  gradmax = sqrt(gradmax);
-  gradmin = sqrt(gradmin);
+  double gradmax = 0.;
+  double gradPmax = 0.;
+  double gradmin = 1e8;
+  double gradPmin = 1e8;
+  double gradM = 0.;
+  double gradPM = 0.;
+  double volT = 0.;
   double *porosity = problem->porosity;
+  for (int iel = 0; iel < mesh->n_elements; ++iel) {
+    const unsigned int *el = &mesh->elements[iel*N_N];
+    double C[N_N],U[DIMENSION][N_N], P[N_N];
+    for (int i = 0; i < N_N; ++i){
+      C[i] = porosity[el[i]];
+      P[i] = solution[el[i]*n_fields+DIMENSION];
+      for (int j=0; j<DIMENSION; ++j) {
+        U[j][i] = solution[el[i]*n_fields+j];
+      }
+    }
+    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
+    for (int i=0; i<DIMENSION; ++i)
+      for (int j=0; j<DIMENSION; ++j)
+        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
+#if DIMENSION == 2
+    const double vol = detDD(dxdxi)/2;
+#else        
+    const double vol = detDD(dxdxi)/6;
+#endif    
+    invDD(dxdxi, dxidx);
+    double dphi[N_SF][DIMENSION];
+    grad_shape_functions(dxidx, dphi);
+    double ngrad2 = 0;
+    for (int i=0; i<DIMENSION; ++i){
+      for (int j=0; j<DIMENSION; ++j){
+        double gradU = 0;
+        for (int k=0; k<N_SF; ++k){
+          gradU += dphi[k][j]*U[i][k]/C[k];
+        }
+        ngrad2 += gradU*gradU;
+      }
+    }
+    gradM += sqrt(ngrad2)*vol;
+    volT += vol;
+    gradmax = fmax(gradmax,ngrad2);
+    gradmin = fmin(gradmin,ngrad2);
+    ngrad2 = 0;
+    for (int j=0; j<DIMENSION; ++j){
+      double gradP = 0;
+      for (int k=0; k<N_SF; ++k){
+        gradP += dphi[k][j]*P[k];
+      }
+      ngrad2 += gradP*gradP;
+    }
+    gradPM += sqrt(ngrad2)*vol;
+    gradPmax = fmax(gradPmax,ngrad2);
+    gradPmin = fmin(gradPmin,ngrad2);
+  }
+  gradPM /= volT;
+  gradM /= volT;
+  gradPmin = 20*sqrt(gradPmin)/gradPM;
+  gradPmax = .1*sqrt(gradPmax)/gradPM;
+  gradmin = 20*sqrt(gradmin)/gradM;
+  gradmax = .1*sqrt(gradmax)/gradM;
+  printf("gradPmax=%e\n",gradPmax);
+  printf("gradPmin=%e\n",gradPmin);
   for (int iel = 0; iel < mesh->n_elements; ++iel) {
     const unsigned int *el = &mesh->elements[iel*N_N];
     double C[N_N],U[DIMENSION][N_N], P[N_N];
@@ -934,10 +1012,10 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
 /*      }*/
 /*      ngrad2 += gradP*gradP;*/
 /*    }*/
-    double ngrad = sqrt(ngrad2);
+    double ngrad = sqrt(ngrad2)/gradM;
     double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
     
-    ngrad = pow(ngrad2, 0.25); 
+    ngrad2 = 0.; 
     for (int j=0; j<DIMENSION; ++j){
       double gradP = 0;
       for (int k=0; k<N_SF; ++k){
@@ -945,8 +1023,8 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
       }
       ngrad2 += gradP*gradP;
     }
-    ngrad = sqrt(ngrad2);
-    lc = fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc);
+    ngrad = sqrt(ngrad2)/gradPM;
+    lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
     
     
     for (int j = 0; j < N_N; ++j){
@@ -972,6 +1050,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
     nOfEl += nOfElByNode;
   }
   double ratioEl = n_el/nOfEl;
+  printf("RatioEl = %e\n",ratioEl);
   for (int i = 0; i < mesh->n_nodes; ++i){
     new_mesh_size[i] = fmax(new_mesh_size[i]/cbrt(ratioEl),lcmin);
   }
@@ -996,20 +1075,33 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
 }
 
 
-void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
+void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double alpha, int notComputeEpsilon)
 {
-  fluid_problem_adapt_gen_mesh(problem, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el);
+  struct timespec tic, toc;
+  fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
+  clock_get(&tic);
   Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
+  clock_get(&toc);
+  printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
   double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
   double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
+  clock_get(&tic);
   int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
   for (int i = 0; i < new_mesh->n_nodes; ++i)
     new_eid[i] = -1;
+  clock_get(&toc);
+  printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
+  clock_get(&tic);
   double *newx = malloc(sizeof(double)*DIMENSION*new_mesh->n_nodes);
   for (int i = 0;i < new_mesh->n_nodes;++i) 
     for (int k = 0; k < DIMENSION; ++k)
       newx[i*DIMENSION+k] = new_mesh->x[i*3+k];
+  clock_get(&toc);
+  printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
+  clock_get(&tic);
   mesh_particles_to_mesh(problem->mesh, new_mesh->n_nodes, newx, new_eid, new_xi);
+  clock_get(&toc);
+  printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
   for (int i = 0; i < new_mesh->n_nodes; ++i) {
     unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
     if(new_eid[i] < 0) {
@@ -1029,6 +1121,8 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
   free(newx);
   free(problem->solution);
   problem->solution = new_solution;
+  free(problem->epsilon);
+  problem->epsilon = malloc(new_mesh->n_nodes*sizeof(double));
   free(problem->porosity);
   problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
   free(problem->old_porosity);
@@ -1039,6 +1133,19 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
     problem->particle_element_id[i] = -1;
   mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
   fluid_problem_compute_node_volume(problem);
+#if DIMENSION == 2
+  for (int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
+#else
+  for (int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
+  }
+#endif
+  if(notComputeEpsilon){
+    for (int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->epsilon[i] = alpha;
+  }
+  }
   fluid_problem_compute_porosity(problem);
   hxtLinearSystemDelete(&problem->linear_system);
   #ifdef HXT_HAVE_PETSC
@@ -1083,6 +1190,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
   mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw);
   clock_get(&toc);
   printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
+  clock_get(&tic);
   for (int i = 0; i < n; ++i) {
     problem->particle_mass[i] = mass[i];
     problem->particle_volume[i] = volume[i];
@@ -1091,10 +1199,15 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
       problem->particle_velocity[i*DIMENSION+k] = velocity[i*DIMENSION+k];
     }
   }
+  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_porosity[i] = problem->porosity[i];
   }
+  clock_get(&tic);
   fluid_problem_compute_porosity(problem);
+  clock_get(&toc);
+  printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
 }
 
 int *fluid_problem_particle_element_id(FluidProblem *problem)

marblesbag/src/fluid_problem.h

@@ -11,7 +11,7 @@ typedef struct {
 }StrongBoundary;
 
 typedef struct {
-  double epsilon;
+  double *epsilon;
   double rho;
   double mu;
   double alpha;
@@ -42,8 +42,8 @@ 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);
 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);
-void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el);
+FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, int n_strong_boundaries, StrongBoundary *strong_boundaries, int notComputeEpsilon);
+void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double alpha, int notComputeEpsilon);
 int *fluid_problem_particle_element_id(FluidProblem *problem);
 int fluid_problem_n_particles(FluidProblem *problem);
 #endif

testcases/drop-3d/metzger/drop.py

@@ -27,7 +27,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
     p.write(filename)
 
 
-outputdir = "MetzgerB1"
+outputdir = "MetzgerB"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
@@ -36,7 +36,7 @@ ii = 0
 
 r=154e-6
 R = 3.3e-3
-compacity = 0.2
+compacity = .20
 drho = 35.4
 p = scontact3.ParticleProblem()
 #R = np.random.uniform(45e-06, 90e-06, len(x))
@@ -48,7 +48,7 @@ nu = 1.17/1030.
 rho = 1030#rhop-drho/compacity
 V = 0.5 # todo : estimate V base on limit velocity
 print('V',V)
-tEnd = 300
+tEnd = 1000
 
 #numerical parameters
 lcmin = 0.0004 # approx r*100 but should match the mesh size
@@ -68,16 +68,16 @@ p = scontact3.ParticleProblem(outputdir+"/part-00000")
 
 print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
 print("RHOP = %g" % rhop)
-outf = 1
+outf = 10
 outf1 = 100000
 
 #strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
-strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
-#strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
-fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
+strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries, 1)
 
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
-fluid.export_vtk(outputdir,0,0)
+
 
 ii = 0
 jj = 0
@@ -87,34 +87,39 @@ print(p.velocity()[1,1])
 print(1-compacity)
 
 
-#for jj in range(3):
-#    dt1 = dt
-#    t = 0
-#    if jj!=0:
-#        fluid.adapt_mesh(0.01,50,1,50,4e-4,100000)
-#    while t<tEnd1:
-#        forces = fluid.compute_node_force(dt1)
-#        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
-#        niter = fluid.implicit_euler(dt1)
-#        t += dt1
-#        if niter < 5 and dt1 < 20:
-#            dt1 *= 3
-#        print(dt1)
-
+for jj in range(4):
+    dt1 = dt
+    t = 0
+    if jj!=0:
+        fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1)
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    
+    while t < tEnd1 : 
+        niter = fluid.implicit_euler(dt1)
+        forces = fluid.compute_node_force(dt1)
+        vn = p.velocity() + forces * dt1 / p.mass()
+        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+        if niter < 5 and dt1 < 20:
+            dt1 *= 3
+        elif niter > 5 and dt1 > dt:
+            dt1 /= 3
+        print(dt1)
+        t += dt1
 #        if ii %outf == 0 :
-#            ioutput = int(ii/outf) + 1
+#                ioutput = int(ii/outf) + 1