Revert "continuous grad a in surface tension"

This reverts commit 50602a91.

fluid/fluid_problem.c

@@ -760,7 +760,7 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
     invDD(dxdxi, dxidx);
   }
 
-  double sigma = 22.27*1e-3*1000;
+  double sigma = 22.27*1e-9;
   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) {
@@ -785,23 +785,24 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
   for (int i = 0; i< mesh->n_nodes; ++i) {
     a_cg[i] /= problem->node_volume[i];
   }
-  double *da_cg = malloc(sizeof(double)*D*mesh->n_nodes);
   for (int iel = 0; iel < mesh->n_elements; ++iel) {
     const int *el = &mesh->elements[iel*N_N];
     double dxdxi[D][D],dxidx[D][D];
     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];
-    double det = invDD(dxdxi, dxidx);
-#if DIMENSION == 2
-    double vol = det/6;
-#else
-    double vol = det/24;
-#endif
+    invDD(dxdxi, dxidx);
     double dphi[N_N][D];
     grad_shape_functions(dxidx, dphi);
+    double a = 0;
     double da[D] = {0};
     double nda = 0;
+    for (int i = 0; i < N_SF; ++i) {
+      a += a_cg[el[i]];
+    }
+    a /= N_N;
+    //double fa = a>0.5 ? -1 : 1;
+    double fa = -1;
     for (int k=0; k<D; ++k){
       for (int i = 0; i < N_SF; ++i) {
         da[k] += a_cg[el[i]]*dphi[i][k];
@@ -809,53 +810,17 @@ static void fluid_problem_surface_tension(FluidProblem *problem, double *all_loc
       nda += da[k]*da[k];
     }
     nda = fmax(sqrt(nda),1e-8);
-    for (int iphi = 0; iphi < N_N; ++iphi) {
-      for (int id = 0; id< D; ++id) {
-        da_cg[el[iphi]*D+id] += da[id]*vol/N_N;
-      }
-    }
-  }
-  for (int iel = 0; iel < mesh->n_elements; ++iel) {
-    const int *el = &mesh->elements[iel*N_N];
-    double dxdxi[D][D],dxidx[D][D];
-    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];
-    double det = invDD(dxdxi, dxidx);
-#if DIMENSION == 2
-    double vol = det/6;
-#else
-    double vol = det/24;
-#endif
-    double dphi[N_N][D];
-    grad_shape_functions(dxidx, dphi);
     double *local_vector = &all_local_vector[iel*n_fields*N_N];
-    
-    for (int iqp = 0; iqp < N_QP; ++iqp) {
-      double phi[N_N];
-      shape_functions(QP[iqp],phi);
-      double jw = QW[iqp]*det;
-      double da[D] = {0};
-      double nda = 0;
-      for (int k=0; k<D; ++k){
-        for (int i = 0; i < N_SF; ++i) {
-          da[k] += da_cg[el[i]*D+k]*phi[i];
-        }
-        nda += da[k]*da[k];
-      }
-      nda = fmax(sqrt(nda),1e-8);
-      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*jw;
-          }
-          local_vector[(U+id)*N_N+iphi] += nda*dphi[iphi][id]*sigma*jw;
+    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] -= nda*dphi[iphi][id]*sigma*fa;
       }
     }
   }
   free(a_cg);
-  free(da_cg);
 }
 
 static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, int reduced_gravity, double stab_param)

testcases/goutte/liqDrop.py

@@ -56,7 +56,7 @@ print("rhom = %g, num = %g" % (rhom,num))
 
 #numerical parameters
 tEnd = 200                                     # final time
-dt = 1e-1                                       # time step
+dt = 1e-2                                       # time step
 outf = 5                                       # number of iterations between output files
 
 def outerBndV(x) :