c dP

fluid/fluid_problem.c

@@ -212,8 +212,8 @@ 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]));
-    f00[(U+id)*n_fields+P] += (pid<0?0:-n[id]);
+    f0[U+id] = c*(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] += c*(pid<0?0:-n[id]);
     f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
     for (int jd = 0; jd <D; ++jd) {
       f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
@@ -355,10 +355,10 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
     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-vol*dp[i];
+    f[U+i] = -drag;
   }
   *dfdu = gammastar/c;
-  *dfddp = gammastar*dt/mass*vol-vol;
+  *dfddp = gammastar*dt/mass*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 *dcold,const double *u_mesh,const double *u_mesh_old, const double *du_mesh_vec, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
@@ -395,11 +395,11 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
   for (int i = 0; i < D; ++i) {
     f0[U+i] = 
       rho*(u[i]-uold[i])/dt
-      + dp[i] 
+      + c*dp[i] 
       + G[i]*rhoreduced*c
       + problem->volume_drag*u[i]*mu;
     f00[(U+i)*n_fields+U+i] = rho/dt + problem->volume_drag*mu;
-    f01[(U+i)*n_fields*D+P*D+i] = 1;
+    f01[(U+i)*n_fields*D+P*D+i] = c;
     // advection :
     // 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)
@@ -410,7 +410,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
       fo[U+i] += rho*((uold[j]/cold)*du[i][j] + u[i]*(duold[j][j]/cold-uold[j]/(cold*cold)*dcold[j]));
       //f00[(U+i)*n_fields+U+i] += rho*(duold[j][j]/cold-du_mesh[j][j]-uold[j]/(cold*cold)*dcold[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*(u[j]-c*u_mesh[j]);
+      f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(uold[j]-c*u_mesh[j]);
       //f01[(U+i)*n_fields*D+(U+i)*D+j] += rho*(uold[j]/cold-u_mesh[j]);
       fo1[i+j] +=  rho*(uold[j]/cold);
       foo[U+i] += rho*(duold[j][j]/cold-uold[j]/(cold*cold)*dcold[j]);
@@ -1224,11 +1224,17 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
   problem->solution_explicit = solution_old;
   double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
   double norm0 = 0;
-  for (int i=0; i<problem->mesh->n_nodes; ++i){
+  /*for (int i=0; i<problem->mesh->n_nodes; ++i){
       problem->oldporosity[i] = problem->porosity[i];
       if (-problem->mesh->x[3*i]/10+0.55 > 1) problem->porosity[i] = 1;
       else problem->porosity[i] = -problem->mesh->x[3*i]/10+0.55;
-    }
+    }*/
+  for (int i=0; i<mesh->n_nodes; ++i){
+    double x = mesh->x[i*3+0];
+    problem->porosity[i] = (x+5.1)/10.1;
+    double oldx = x-dt*problem->velocity_mesh[i*2+0];
+    problem->oldporosity[i] = (oldx+5.1)/10.1;
+  }
   compute_stab_parameters(problem,dt);
   int i;
   for (i=0; i<newton_max_it; ++i) {
@@ -1614,10 +1620,10 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
   free(problem->oldporosity);
   problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
   for(int i = 0; i < mesh->n_nodes; ++i){
-    if (-problem->mesh->x[3*i]/10+0.55 > 1){problem->porosity[i] = 1;problem->oldporosity[i] = 1;}
-    else {problem->porosity[i] = -problem->mesh->x[3*i]/10+0.55;problem->oldporosity[i] = -problem->mesh->x[3*i]/10+0.55;}
-    //problem->porosity[i] = 1.;
-    //problem->oldporosity[i] = 1.;
+    //if (-problem->mesh->x[3*i]/10+0.55 > 1){problem->porosity[i] = 1;problem->oldporosity[i] = 1;}
+    //else {problem->porosity[i] = -problem->mesh->x[3*i]/10+0.55;problem->oldporosity[i] = -problem->mesh->x[3*i]/10+0.55;}
+    problem->porosity[i] = 1.;
+    problem->oldporosity[i] = 1.;
     problem->velocity_mesh[i*D] = 0.;
     problem->velocity_mesh[i*D+1] = 0.;
   }

testcases/hugo/adaptation_mailage/adaptation1.py

@@ -280,7 +280,7 @@ nb1 = int(nbtot/(2*(np.pi+2)))
 nb2 = nbtot -2*nb1
 print((nb1,nb2))
 
-liste_depl = geo_ligne(nb1,rayon,[1,0])+ + 1*geo_cercle(nb2,rayon) + geo_ligne(nb1,rayon,[-1,0]) 
+liste_depl = geo_ligne(nb1,rayon,[1,0]) + 1*geo_cercle(nb2,rayon) + geo_ligne(nb1,rayon,[-1,0]) 
 #liste_depl =geo_ligne(1,rayon,[1,0])+geo_cercle(nb2,rayon)
 #liste_depl = 1*(geo_ligne(nb2,rayon,[0,-5])+geo_ligne(200,rayon,[2,1]))
 #liste_depl = geo_ligne(nb1,rayon,[0,-1])
@@ -303,7 +303,8 @@ cl_vitesse = [0.,0.]
 depl[0] += liste_depl[1][0] 
 depl[1] += liste_depl[1][1]
 
-for k in range(2,n):
+fluid.velocity_mesh()[:] = 0
+for k in range(n):
     #norm_max = calcul_norm(para,liste_depl[k],pos_centre)
     for j in range(N):
         #[poubelle,fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)