cdP

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];
@@ -354,10 +354,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;//-vol*dp[i];
   }
   *dfdu = gammastar/c;
-  *dfddp = gammastar*dt/mass*vol-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) {
@@ -387,11 +387,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)