fix adapt

src/fluid_problem.c

@@ -33,8 +33,8 @@
 
 
 #define newton_max_it 10
-#define newton_atol 5e-7
-#define newton_rtol 1e-5
+#define newton_atol 5e-6
+#define newton_rtol 1e-4
 #define D DIMENSION
 
 #define  deps 1e-8
@@ -478,10 +478,10 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
       mu = problem->mu[0]*a + problem->mu[1]*(1-a);
     }
     double nu = mu/rho;
-    problem->taup[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
+    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);
-    problem->taua[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h));
-    problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*0.03;
+    problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
+    problem->extraa[iel] = fmax(0.,fmax(-amin+0.01,amax-1+0.01))*0.01;
     maxtaup = fmax(maxtaup, problem->taup[iel]);
     mintaup = fmin(mintaup, problem->taup[iel]);
   }
@@ -556,6 +556,9 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
       f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
     }
   }
+  if(problem->tf){
+    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
+  }
 }
 
 static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
@@ -1170,12 +1173,12 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
     }
   }
   gradM /= volT;
-  gradmin = sqrt(gradmin)/gradM;
+  gradmin = sqrt(gradmin)*sqrt(gradmax)/gradM;
   gradmax = sqrt(gradmax)/gradM;
   printf("gradmax=%e\n",gradmax);
   printf("gradmin=%e\n",gradmin);
   gradAM /= volT;
-  gradAmin = sqrt(gradAmin)/gradAM;
+  gradAmin = sqrt(gradAmin)*sqrt(gradAmax)/gradAM;
   gradAmax = sqrt(gradAmax)/gradAM;
   printf("gradAmax=%e\n",gradAmax);
   printf("gradAmin=%e\n",gradAmin);
@@ -1185,6 +1188,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
     double C[N_N],Un[D][N_N],An[N_N];
     for (int i = 0; i < N_N; ++i){
       C[i] = porosity[el[i]];
+      An[i] = solution[el[i]*n_fields+D+1];
       for (int j=0; j<D; ++j) {
         Un[j][i] = solution[el[i]*n_fields+j];
       }
@@ -1210,17 +1214,22 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
     double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
     
     if(problem->tf){
-      double ngradA = 0;
+      double ngradA2 = 0;
+      //printf("2fluids\n");
       for (int j=0; j<D; ++j){
         double gradA = 0;
         for (int k=0; k<N_SF; ++k){
           gradA += dphi[k][j]*An[k];
         }
-        ngradA += gradA*gradA;
+        ngradA2 += gradA*gradA;
       }
-      double ngrad = sqrt(ngradA/gradAM);
+      double ngrad = sqrt(ngradA2)/gradAM;
       double lc2 = lcmin /fmin(1,fmax(ngrad/(gradAmax),gradAmin/(gradAmax)));
+      //if(ngradA2>1e-5)
+      //printf("lc2=%g,  gradAm=%g,  ngradA=%g,  gradnAmin=%g,  gradAmax=%g\n",lc2,gradAM,ngradA2,gradAmin,gradAmax);
+      //printf("before lc=%g\n",lc);
       lc = fmin(lc,lc2);
+      //printf("after lc=%g\n",lc);
     }
 
     for (int j = 0; j < N_N; ++j){