small extraa and cropped a, aold

src/fluid_problem.c

@@ -181,6 +181,7 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
   }
   else{
     a = sol[A];
+    a = fmin(1.,fmax(0.,a));
     rho = problem->rho[0]*a + problem->rho[1]*(1-a);
     mu = problem->mu[0]*a + problem->mu[1]*(1-a);
   }
@@ -345,6 +346,7 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
   }
   else{
     double a = s[A];
+    a = fmin(1.,fmax(0.,a));
     rho = problem->rho[0]*a+problem->rho[1]*(1-a);
     mu = problem->mu[0]*a+problem->mu[1]*(1-a);
     f[A] = un*(un > 0 ? a : 0.);
@@ -385,6 +387,7 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
   }
   else{
     double a = s[A];
+    a = fmin(1.,fmax(0.,a));
     rho = problem->rho[0]*a+problem->rho[1]*(1-a);
     mu = problem->mu[0]*a+problem->mu[1]*(1-a);
     f[A] = un*(un > 0 ? a : 1.);
@@ -464,6 +467,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
       mu = problem->mu[0];
     }else{
       a /= N_N;
+      a = fmin(1.,fmax(0.,a));
       rho = problem->rho[0]*a + problem->rho[1]*(1-a);
       mu = problem->mu[0]*a + problem->mu[1]*(1-a);
     }
@@ -471,7 +475,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
     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(2/dt)+pow2(2*normu/h));
-    problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1));
+    problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*0.1;
     maxtaup = fmax(maxtaup, problem->taup[iel]);
     mintaup = fmin(mintaup, problem->taup[iel]);
   }
@@ -495,7 +499,9 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
   }
   else{
     a = sol[A];
+    a = fmin(1.,fmax(0.,a));
     aold = solold[A];
+    aold = fmin(1.,fmax(0.,aold));
     Ra = (a*c-aold*cold)/dt;
     rho = problem->rho[0]*a + problem->rho[1]*(1-a);
     rhoold = problem->rho[0]*aold + problem->rho[1]*(1-aold);