triangle surface

fluid/fluid_problem.c

@@ -418,8 +418,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
 
   for (int i = 0; i < D; ++i) {
     double dudt = (u[i]-uold[i])/dt;
-    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uudrho[i]/rho + (i==1 ? -(reduced_gravity ? 0 : 1)*c*problem->g : 0.);
-    f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? 0 : 1)*c*rho*problem->g : 0.);
+    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uudrho[i]/rho + (i==1 ? -(reduced_gravity ? (rho-problem->rho[0])/rho : 1)*c*problem->g : 0.);
+    f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? (rho-problem->rho[0])/rho : 1)*c*rho*problem->g : 0.);
     for (int j = 0; j < D; ++j) {
       f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + (stab_param==0 ? 1 : 0)*R*u[j]*rho*taup;
     }
@@ -760,7 +760,7 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
     invDD(dxdxi, dxidx);
   }
 
-  double sigma = 22.27*1e-9;
+  double sigma = 22.27*1e-3;
   int n_fields = fluid_problem_n_fields(problem);
   double *a_cg = malloc(sizeof(double)*mesh->n_nodes);
   for (int i = 0; i< mesh->n_nodes; ++i) {
@@ -791,7 +791,12 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
     for (int i=0; i<D; ++i)
       for (int j=0; j<D; ++j)
         dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
-    invDD(dxdxi, dxidx);
+    double det = invDD(dxdxi, dxidx);
+#if DIMENSION == 2
+    double vol = det/2;
+#else
+    double vol = det/6;
+#endif
     double dphi[N_N][D];
     grad_shape_functions(dxidx, dphi);
     double a = 0;
@@ -801,8 +806,10 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
       a += a_cg[el[i]];
     }
     a /= N_N;
-    //double fa = a>0.5 ? -1 : 1;
-    double fa = -1;
+    double c = 0;
+    for (int i = 0; i < N_N; ++i){
+      c += problem->porosity[el[i]]/N_N;
+    }
     for (int k=0; k<D; ++k){
       for (int i = 0; i < N_SF; ++i) {
         da[k] += a_cg[el[i]]*dphi[i][k];
@@ -814,9 +821,9 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
     for (int iphi = 0; iphi < N_N; ++iphi) {
       for (int id = 0; id < D; ++id) {
         for (int jd = 0; jd < D; ++jd) {
-          local_vector[(U+id)*N_N+iphi] += da[id]*da[jd]*dphi[iphi][jd]/nda*sigma*fa;
+          local_vector[(U+id)*N_N+iphi] -= c*vol*da[id]*da[jd]*dphi[iphi][jd]/nda*sigma;
         }
-        local_vector[(U+id)*N_N+iphi] -= nda*dphi[iphi][id]*sigma*fa;
+        local_vector[(U+id)*N_N+iphi] += c*vol*nda*dphi[iphi][id]*sigma;
       }
     }
   }