Merge remote-tracking branch 'origin/master' into margaux

AUTHORS.txt

@@ -3,5 +3,6 @@ MigFlow software has been developped by:
 
 Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
 Matthieu Constant <matthieu.constant@uclouvain.be>
+Nathan Coppin <nathan.coppin@uclouvain.be>
 Vincent Legat <vincent.legat@uclouvain.be>
 Frédéric Dubois <frederic.dubois@umontpellier.fr>

fluid/fluid_problem.c

@@ -86,6 +86,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
   const Mesh *mesh = problem->mesh;
   const double *porosity = problem->porosity;
   const double *solution = problem->solution;
+  int drag_in_stab = problem->drag_in_stab;
   int n_fields = fluid_problem_n_fields(problem);
   size_t local_size = N_SF*n_fields;
   double *s = malloc(sizeof(double)*n_fields);
@@ -106,8 +107,8 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
     double dxidx[D][D], dphi[N_SF][D];
     double detj = mesh_dxidx(mesh, iel, dxidx);
     grad_shape_functions(dxidx, dphi);
-    double *local_vector = &all_local_vector[local_size*iel];
-    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
+    double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
+    double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
     for (int i = 0; i < n_fields; ++i) {
       s[i] = 0;
       sold[i] = 0;
@@ -139,19 +140,23 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
     }
     double f[D],dfdu,dfddp;
     particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
+    if (!local_vector)
+      continue;
     double rho,nu;
     fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
-    double pspg = problem->taup[iel]/rho;
+    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
     for (int iD = 0; iD < D; ++iD) {
       for (int iphi = 0; iphi < N_SF; ++iphi){
         local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
         local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
         for (int jD = 0; jD < D; ++jD) {
-          double supg = sold[U+jD]*problem->taup[iel];
+          double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
           local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
         }
       }
     }
+    if (!local_matrix)
+      continue;
     for (int iD = 0; iD < D; ++iD){
       for (int iphi = 0; iphi < N_SF; ++iphi){
         for (int jphi = 0; jphi < N_SF; ++jphi){
@@ -160,7 +165,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
           LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
           LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
           for (int jD = 0; jD < D; ++jD) {
-            double supg = sold[U+jD]*problem->taup[iel];
+            double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
             LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
             LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
           }
@@ -173,7 +178,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
   free(sold);
 }
 
-static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data, double *f00)
+static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data, double *f00, double *f01)
 {
   const int vid = wbnd->vid;
   const int pid = wbnd->pid;
@@ -181,7 +186,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   const int aid = wbnd->aid;
   double p = s[P];
   const int n_fields = fluid_problem_n_fields(problem);
-  double u[D], du[D][D], dp[D], dun[D], dunt[D];
+  double u[D], c_du_o_c[D][D], dp[D];
   double s_c = 0;
   double unold = 0;
   double unext = 0;
@@ -189,18 +194,13 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   for (int iD = 0; iD < D; ++iD) {
     u[iD] = s[U+iD];
     dp[iD] = ds[P*D+iD];
-    dun[iD] = 0;
-    dunt[iD] = 0;
     unold += sold[U+iD]*n[iD];
     if(vid>=0)
     unext += data[vid+iD]*n[iD];
     dcn += dc[iD]*n[iD];
     s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
-    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
-    for (int jD = 0; jD < D; ++jD) {
-      dun[iD] += du[iD][jD]*n[jD];
-      dunt[iD] += du[jD][iD]*n[jD];
-    }
+    for (int jD = 0; jD < D; ++jD)
+      c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
   }
   double h = problem->element_size[eid];
   double sigma = (1+D)/(D*h)*mu*N_N;
@@ -217,11 +217,15 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     f0[P] = unext;
   }
   for (int id = 0; id < D; ++id) {
-    f0[U+id] = (pid<0?0:data[pid]-p)*n[id] + (vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id])); 
+    f0[U+id] = (pid<0?0:data[pid]-p)*n[id] + (vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]));
     f00[(U+id)*n_fields+P] += (pid<0?0:-n[id]);
     f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
     for (int jd = 0; jd <D; ++jd) {
-      f00[(U+id)*n_fields+U+jd] += (vid<0?0:n[id]*n[jd]*sigma);
+      f0[U+id] += mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
+      f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) - mu/c*dc[jd]*n[jd];
+      f00[(U+id)*n_fields+U+jd] -= mu/c*dc[id]*n[jd];
+      f01[(U+id)*n_fields*D+(U+id)*D+jd] += mu*n[jd];
+      f01[(U+id)*n_fields*D+(U+jd)*D+id] += mu*n[jd];
     }
   }
   double unbnd = unold;
@@ -280,7 +284,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
     a = fmin(1.,fmax(0.,a));
     fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
     double nu = mu/rho;
-    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)))*3;
+    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
     problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
   }
 }
@@ -349,17 +353,17 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
     gamma = gamma0*a + gamma1*(1-a);
   }
   else{
-    gamma = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
+    gamma = particle_drag_coeff(Due,mu,rho,vol,c);
   }
   double gammastar = gamma/(1+dt/mass*gamma);
   for (int i = 0; i < D; ++i) {
     double upstar = up[i]+dt/mass*(contact[i]+mass*G[i]-vol*dp[i]);
     double drag = gammastar*(upstar-u[i]/c);
     problem->particle_force[ip*D+i] = -drag+G[i]*mass-vol*dp[i];
-    f[U+i] = -drag;
+    f[U+i] = -drag-vol*dp[i];
   }
   *dfdu = gammastar/c;
-  *dfddp = gammastar*dt/mass*vol;
+  *dfddp = gammastar*dt/mass*vol-vol;
 }
 
 static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
@@ -381,7 +385,8 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
     }
     divu += du[iD][iD];
   }
-  f0[P] = (c-cold)/dt;
+  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
+  //f00[P*n_fields+P] = 1/dt*0.1;
   double G[D] = {0};
   G[1] = -problem->g;
   double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
@@ -397,9 +402,9 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
     // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
     //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
     for (int j = 0; j < D; ++j) {
-      f0[U+i] += rho/c*(uold[j]/c*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
+      f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
       f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
-      f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j]/c;
+      f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
     }
     for (int j = 0; j < D; ++j) {
       f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
@@ -445,6 +450,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
   double *dsold = malloc(sizeof(double)*n_fields*D);
   double *f0 = malloc(sizeof(double)*n_fields);
   double *f00 = malloc(sizeof(double)*n_fields*n_fields);
+  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
   double i_bnd = 0;
   double iv_bnd_up = 0;
   double iv_bnd_down = 0;
@@ -513,6 +519,9 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
         for (int i = 0; i < n_fields*n_fields; ++i) {
           f00[i] = 0;
         }
+        for (int i = 0; i < n_fields*n_fields*D; ++i) {
+          f01[i] = 0;
+        }
         double c = 0;
         double a = 0;
         for (int i = 0; i < DIMENSION; ++i) {
@@ -542,7 +551,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
             }
           }
         }
-        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00);
+        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
         for (int ifield = 0; ifield < n_fields; ++ifield) {
           for (int iphi = 0; iphi < D; ++iphi){
             local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
@@ -552,7 +561,15 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
           for (int ifield = 0; ifield < n_fields; ++ifield){
             for (int iphi = 0; iphi < N_LSF; ++iphi){
               for (int jphi = 0; jphi < N_LSF; ++jphi){
-                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*phi[jphi]*f00[ifield*n_fields+jfield];
+                double d = phi[jphi]*f00[ifield*n_fields+jfield];
+                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*d;
+              }
+              for (int jphi = 0; jphi < N_SF; ++jphi){
+                double d = 0;
+                for (int jD = 0; jD < D; ++jD) {
+                  d += dphi[jphi][jD]*f01[ifield*n_fields*D+jfield*D+jD];
+                }
+                LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += jw*phi[iphi]*d;
               }
             }
           }
@@ -1104,8 +1121,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
     for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
       problem->solution[j] -= dx[j];
     }
-    //break; /// equation is linear
+    node_force_volume(problem,problem->solution_explicit,dt,NULL,NULL);
+    break; /// equation is linear
   }
+  
   printf("total solve time : %g\n", totaltime);
   free(dx);
   free(rhs);
@@ -1155,7 +1174,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
   }
 }
 
-FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type)
+FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)
 {
   if (n_fluids != 1 && n_fluids != 2) {
     printf("error : only 1 or 2 fluids are supported\n");
@@ -1172,6 +1191,7 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
   problem->n_fluids = n_fluids;
   problem->sigma = sigma;
   problem->stab_param = 0;
+  problem->drag_in_stab = drag_in_stab;
   problem->reduced_gravity = 0;
   problem->strong_boundaries = NULL;
   problem->n_strong_boundaries = 0;

fluid/fluid_problem.h

@@ -97,6 +97,7 @@ struct FluidProblem {
   //parameters
   int reduced_gravity;
   double stab_param;
+  int drag_in_stab;
 };
 
 // complete force on the particle (including gravity)
@@ -105,7 +106,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
 void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload);
 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, const char *petsc_solver_type);
+FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab);
 void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
 
 void fluid_problem_adapt_mesh(FluidProblem *problem, double e_target, double lcmax, double lcmin, int n_el, double cmax, double cmin, double U_0, double P_0);

python/fluid.py

@@ -65,7 +65,7 @@ def _is_string(s) :
 class FluidProblem :
     """Create numerical representation of the fluid."""
 
-    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., petsc_solver_type="-pc_type lu"):
+    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., petsc_solver_type="-pc_type lu",drag_in_stab=0):
         """Build the fluid structure.
 
         Keyword arguments:
@@ -98,7 +98,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(c_double(g), n_fluids, np2c(mu), np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), petsc_solver_type.encode()))
+        self._fp = c_void_p(self._lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), petsc_solver_type.encode(),c_int(drag_in_stab)))
         if self._fp == None :
             raise NameError("Cannot create fluid problem.")
 

testcases/avalanch/avalanch1fluid/avalanch1fluidnofriction.py

@@ -23,7 +23,7 @@
 
 # TESTCASE DESCRIPTION:
 # Collapsing of a grain column in fluid. 
-# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/dep.py) to create a compact column as initial situation
+# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/depot.py) to create a compact column as initial situation
 # Positions, masses and radii of grains are setted in the depot.py file.
 from migflow import fluid
 from migflow import scontact
@@ -79,12 +79,12 @@ fluid.set_wall_boundary("Left")
 fluid.set_wall_boundary("Top",pressure=0)
 fluid.set_wall_boundary("Right")
 # Set locations of the grains in the mesh and compute the porosity in each computation cell.
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
 # Write output files for post-visualization.
 fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+
 
 # 
 # COMPUTATION LOOP
@@ -109,7 +109,7 @@ while t < tEnd :
         p.iterate(dt/nsub, forces)  
 
     t += dt
-    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
 
     # Output files writing.
     if ii%outf == 0 :

testcases/avalanch/avalanch2fluids/avalanch2fluids.py

@@ -126,9 +126,10 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
 s         =  fluid.solution()
 x         =  fluid.coordinates()
 #Set initial condition for the concentration
-for i in range(len(x[:,0])):
-    if (x[i,0])<lx:
-        s[i,3] = 1
+c = np.ndarray((fluid.n_nodes()))
+c[:] = 0
+c[x[:,0]<lx] = 1
+fluid.set_concentration_cg(c)
 fluid.export_vtk(outputdir,0,0)
 
 tic       =  time.time()
@@ -141,10 +142,6 @@ while t < tEnd :
     # Fluid solver
     fluid.implicit_euler(dt)
 
-    # Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
-    if (ii%15==0 and ii != 0):
-       fluid.adapt_mesh(1e-2,3e-3,10000)
-
     # Computation of the new velocities
     forces = fluid.compute_node_force(dt)
     vn     = p.velocity() + forces * dt / p.mass()

testcases/avalanch/avalanch2fluids/avalanch2fluidsnofriction.py

@@ -24,7 +24,7 @@
 # TESTCASE DESCRIPTION
 # Collapsing of a grain column in fluid. 
 # The grain column is initially mixed in sea water while the rest of the domain contain soft water. 
-# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/dep.py) to create a compact column as initial situation
+# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/depot.py) to create a compact column as initial situation
 # Positions, masses and radii of graisn are setted in the dep.py file.
 
 from migflow import fluid as fluid
@@ -50,8 +50,8 @@ rhom = 1050                                      # sea water density
 num = 1e-6                                       # sea water kinematic viscosity
 
 # Numerical parameters
-outf = 10                                        # number of iterations between output files
-dt = 1e-4                                        # time step
+outf = 100                                        # number of iterations between output files
+dt = 1e-5                                        # time step
 tEnd = 1                                         # final time
 
 # 
@@ -83,19 +83,19 @@ fluid.set_wall_boundary("Left")
 fluid.set_wall_boundary("Top",pressure=0)
 fluid.set_wall_boundary("Right")
 # Set location of the particles in the mesh and compute the porosity in each computation cell
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
 # Access to fluid fields and node coordiantes
 s = fluid.solution()
 x = fluid.coordinates()
 # Set initial condition for the concentration
 lx = 0.1
-for i in range(len(x[:,0])):
-    if (x[i,0])<lx:
-        s[i,3] = 1
+c = np.ndarray((fluid.n_nodes()))
+c[:] = 0
+c[x[:,0]<lx] = 1
+fluid.set_concentration_cg(c)
 fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
 
 #
 # COMPUTATION LOOP
@@ -104,10 +104,6 @@ while t < tEnd :
     # Fluid solver
     fluid.implicit_euler(dt)
 
-    # Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
-    if (ii%15==0 and ii != 0):
-       fluid.adapt_mesh(1e-2,3e-3,10000)
-
     # Computation of the new velocities
     forces = fluid.compute_node_force(dt)
     vn = p.velocity() + forces * dt / p.mass()
@@ -120,7 +116,7 @@ while t < tEnd :
         p.iterate(dt/nsub, forces)  
 
     t += dt
-    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
 
     # Output files writting
     if ii %outf == 0 :

testcases/avalanch/avalanchnofluid/avalanchnofluidfriction.py

@@ -23,7 +23,7 @@
 
 # TESTCASE DESCRIPTION:
 # Collapsing of a grain column in fluid. 
-# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/dep.py) to create a compact column as initial situation
+# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/depot.py) to create a compact column as initial situation
 # Positions, masses and radii of grains are setted in the depot.py file.
 from migflow import scontact
 from migflow import lmgc90Interface

testcases/avalanch/depot/depot.py

@@ -102,7 +102,7 @@ while t < tEnd :
     vn = p.velocity() + forces * dt / p.mass() 
     vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
     # Number of contact sub time step
-    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    nsub = max(3, int(np.ceil((vmax * dt * 4)/min(p.r()))))
     print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
     # NLGS iterations
     for i in range(nsub) :

testcases/couette-2d/mill.py

@@ -83,7 +83,7 @@ if not os.path.isdir(outputdir) :
 # Physical parameters
 g = 0                                          # gravity
 rho = 1.253e3                                  # fluid density
-rhop = 1000                                    # grains density
+rhop = 1500                                    # grains density
 nu = 1e-4                                      # kinematic viscosity
 
 # Numerical parameters
@@ -91,7 +91,7 @@ dt = 1e-2                                      # time step
 tEnd = 50                                      # final time
 
 # Geometry parameters
-outf = 5                                       # number of iterations between output files
+outf = 1                                       # number of iterations between output files
 rout = 0.0254                                  # outer radius
 rin = 0.0064                                   # inner radius
 r = 397e-6/2                                   # grains radius
@@ -106,10 +106,10 @@ if use_lmgc90 :
     lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
     p = lmgc90Interface.ParticleProblem(2)
 else :
-  p = scontact.ParticleProblem(2,True)
+  p = scontact.ParticleProblem(2,False)
   p.read_vtk(outputdir,0)
-  p.set_friction_coefficient(0.1,"Sand","Sand")#Particle-Particle
-  p.set_friction_coefficient(0.1,"Sand","Steel")#Particle-Wall
+  #p.set_friction_coefficient(0.1,"Sand","Sand")#Particle-Particle
+  #p.set_friction_coefficient(0.1,"Sand","Steel")#Particle-Wall
 
 # Initial time and iteration
 t = 0
@@ -121,14 +121,20 @@ ii = 0
 fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
 # Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
-fluid.set_wall_boundary("Inner",velocity=[0,0],pressure=0)
+fluid.set_wall_boundary("Inner",velocity=[0,0])
+#fluid.set_wall_boundary("BFix",velocity=[0,0],pressure=0)
 fluid.set_wall_boundary("Outer",velocity=[lambda x : -x[:, 1],lambda x : x[:, 0]])
+#fluid.set_strong_boundary("Inner",0,lambda x : -x[:, 1]*0.01);
+#fluid.set_strong_boundary("Inner",1,lambda x :  x[:, 0]*0.01);
+#fluid.set_strong_boundary("BFix",2,0);
+#fluid.set_strong_boundary("Outer",0,0.,0);
+#fluid.set_strong_boundary("Outer",1,0.,1);
 # Set location of the grains in the mesh and compute the porosity in each computation cell
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
 fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
 
 # 
 # COMPUTATION LOOP
@@ -154,7 +160,7 @@ while t < tEnd :
         tol = 1e-6                                # Numerical tolerance for grains intersection
         p.iterate(dt/nsub, forces, tol)
     # Localisation of the grains in the fluid
-    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces()*0)
     t += dt
     # Output files writting
     if ii %outf == 0 :

testcases/depot-2d/depot.py

@@ -34,7 +34,7 @@ import time
 import shutil
 import random
 
-def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
+def genInitialPosition(filename, r, H, ly, lx, rhop) :
     """Set all the particles centre positions and create the particles objects to add in the computing structure
     
     Keyword arguments:    
@@ -51,25 +51,11 @@ def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
     p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
 
     #Definition of the points where the particles are located
-    #x = np.arange(-lx/2+r, lx/2-r, 2*r)
-    #y = np.arange(0-r, -ly+r, -2*r)
-    #x, y = np.meshgrid(x, y)
-    #x = x.flat
-    #y = y.flat
-    
-    N = compacity*(lx*ly)/(np.pi*r**2)
-    e = 2*ly/(N)**(1./2.)
-    
-    print(N)
-
-    # Definition of the points where the particles are located
-    x = np.arange(lx/2-r, -lx/2+r, -e)
-    y = np.arange(-r, -ly+r, -e)
+    x = np.arange(-lx/2+r, lx/2-r, 2*r)
+    y = np.arange(H/2-r, H/2-ly+r, -2*r)
     x, y = np.meshgrid(x, y)
-    
     x = x.flat
     y = y.flat
-    
     # Add a grain at each centre position
     for i in range(len(x)) :
         p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
@@ -86,23 +72,22 @@ r = 1e-3                                        # particles radius
 rhop = 1500                                     # particles density
 rho = 1000                                      # fluid density
 nu = 1e-6                                       # kinematic viscosity
-compacity = 0.25 
 
 # Numerical parameters
 outf = 5                                        # number of iterations between output files
-dt = 5e-3                                       # time step
-tEnd = 10                                       # final time
+dt = 0.25e-2                                       # time step
+tEnd = 100                                      # final time
 
 # Geometrical parameters
-ly = 0.1                                        # particles area height
+ly = 5e-2                                       # particles area height
 lx = 4e-1                                       # particles area widht
-H = 0.6                                         # domain height
+H = 1.2                                         # domain height
 
 #
 # PARTICLE PROBLEM
 #
 # Initialise particles
-genInitialPosition(outputdir, r, H, ly, lx, rhop, compacity)
+genInitialPosition(outputdir, r, H, ly, lx, rhop)
 p = scontact.ParticleProblem(2)
 p.read_vtk(outputdir,0)
 
@@ -123,11 +108,11 @@ fluid.set_wall_boundary("Bottom")
 fluid.set_wall_boundary("Lateral")
 fluid.set_wall_boundary("Top",pressure=0)
 # Set location of the particles in the mesh and compute the porosity in each computation cell
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
 fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
 
 #
 # COMPUTATION LOOP
@@ -135,13 +120,7 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_
 while t < tEnd :
     # Fluid solver
     fluid.implicit_euler(dt)
-    
-    # Adaptation of the mesh.
-    # No adapt for this reference case !!! 
-    #if (ii%15==0 and ii != 0):   
-    #    fluid.adapt_mesh(500*r,8*r,6602*2)  # for the moment we cannot set the number of elements in the new mesh
-     
-        
+
     # Computation of the new velocities