mostly ok

src/fluid_problem.c

@@ -298,8 +298,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
   double epsilonp = problem->epsilon;
   double epsilonq = problem->epsilon/10;
   f0[P] = 1-c;
-  f1[P*D+0] = 0;
-  f1[P*D+1] = 0;
+  f1[P*D+0] = epsilonp*dp[0];
+  f1[P*D+1] = epsilonp*dp[1];
   for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
     int Q = ifluid*(D+1)+D;
     int U = ifluid*(D+1);
@@ -312,7 +312,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
     double uold[D] = {solold[U],solold[U+1]};
     double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
 
-    double oq = 1./fmax(q,0.1);
+    double oq2 = 1./fmax(q,1e-5);
+    double oq = 1./fmax(q,1e-2);
     double divu = 0;
     double tau[D][D];
     double utau[D]={0.};
@@ -326,11 +327,11 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
     f0[Q] = divu+(q-qold)/dt;
     f0[P] -= q;
     for (int i = 0; i < D; ++i) {
-      f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq*divu + utau[i]*oq) - p*dq[i];
+      f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) - p*dq[i];
       for (int j = 0; j < D; ++j) {
         f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
       }
-      f1[Q*D+i] = epsilonp*dp[i]*fmax(q,0.1)+epsilonq*dq[i];
+      f1[Q*D+i] = epsilonq*dq[i];
     }
   }
 }

testcases/poiseuille2fluids/poiseuille.py

@@ -52,7 +52,7 @@ dt = 10                                       # time step
 alpha = 1e-4                                    # stabilization coefficient
 epsilon = alpha*lcmin**2 /min(nu0,nu1)                    # stabilization parametre
 print('epsilon',epsilon)
-6
+
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
 outf = 1                                       # number of iterations between output files
@@ -94,8 +94,8 @@ t = 0
 #set initial_condition
 x = fluid.coordinates()
 s = fluid.solution()
-alpha = (x[:,0]+5)/10
-s[:,2] = 1-alpha*0.8
+a = (x[:,0]+5)/10
+s[:,2] = 1-a*0.8
 s[:,5] = 1-s[:,2]
 
 fluid.export_vtk(outputdir,0,0)