cavity ok but poiseuille error slightly larger

src/fluid_problem.c

@@ -242,37 +242,53 @@ typedef void f_cb(FluidProblem *problem, const double *n, double *f, const doubl
 static void f_inflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
   int n_fluids = problem->n_fluids;
   int P = n_fluids*(D+1);
+  const double epsilon = problem->epsilon;
   double p = s[P];
+  double dp[D];
+  for (int id = 0; id < D; ++id) {
+    dp[id] = ds[P*D+id];
+  }
   for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
     int U = ifluid*(D+1);
     int Q = ifluid*(D+1)+D;
     double q = s[Q];
     f[Q] = 0;
     for (int id = 0; id < D; ++id) {
-      f[Q] -= s[U+id]*n[id];
-      //f[Q] += dp[id]*n[id]*epsilon*q;
+      f[Q] += -s[U+id]*n[id];
+      f[Q] += dp[id]*n[id]*epsilon*q;
       f[U+id] = 0;//-p*q*n[id];
     }
   }
   f[P] = 0;
+  for (int id = 0; id < D; ++id) {
+    //f[P] += dp[id]*n[id]*epsilon;
+  }
 }
 
 static void f_outflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
   int n_fluids = problem->n_fluids;
   int P = n_fluids*(D+1);
+  const double epsilon = problem->epsilon;
   double p = s[P];
+  double dp[D];
+  for (int id = 0; id < D; ++id) {
+    dp[id] = ds[P*D+id];
+  }
   for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
     int U = ifluid*(D+1);
     int Q = ifluid*(D+1)+D;
     double q = s[Q];
     f[Q] = 0;
     for (int id = 0; id < D; ++id) {
-      f[Q] -= s[U+id]*n[id];
-      //f[Q] += dp[id]*n[id]*epsilon*q;
+      f[Q] += -s[U+id]*n[id];
+      f[Q] += dp[id]*n[id]*epsilon*q;
       f[U+id] = 0;//-p*q*n[id];
     }
   }
   f[P] = 0;
+  for (int id = 0; id < D; ++id) {
+    //f[P] += dp[id]*n[id]*epsilon;
+  }
 }
 
 static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt)
@@ -291,13 +307,13 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
     double q = s[Q];
     f[Q] = 0;
     for (int id = 0; id < D; ++id) {
-      //f[Q] += dp[id]*n[id]*epsilon*q;
+      f[Q] += dp[id]*n[id]*epsilon*q;
       f[U+id] = q*p*n[id];
     }
   }
   f[P] = 0;
   for (int id = 0; id < D; ++id) {
-    f[P] += dp[id]*n[id]*epsilon;
+    //f[P] += dp[id]*n[id]*epsilon;
   }
 }
 
@@ -309,7 +325,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
   double epsilonq = problem->epsilon/10;
   f0[P] = 1-c;
   for (int id = 0; id < D; ++id) {
-    f1[P*D+id] = epsilonp*dp[id];
+    f1[P*D+id] = 0;//epsilonp*dp[id];
   }
   for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
     int Q = ifluid*(D+1)+D;
@@ -343,7 +359,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
       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] = epsilonq*dq[i]-u[i];
+      f1[Q*D+i] = epsilonq*dq[i]-u[i]+q*epsilonp*dp[i];
     }
   }
 }

testcases/separation/cavity.py

@@ -40,7 +40,7 @@ ii = 0
 
 
 #physical parameters
-g =  -9.81                                      # gravity
+g =  0                                      # gravity
 rho = 1000                                      # fluid density
 nu = 1e-3                                   # kinematic viscosity
 mu = nu*rho                                     # dynamic viscosity
@@ -48,8 +48,8 @@ tEnd = 100000                                     # final time
 
 #numerical parameters
 lcmin = .1                                  # mesh size
-dt = 5                                       # time step
-alpha = 1e-7                                    # stabilization coefficient
+dt = 0.1                                       # time step
+alpha = 1e-6                                    # stabilization coefficient
 epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
 print('epsilon',epsilon)
 
@@ -61,13 +61,14 @@ outf = 1                                       # number of iterations between ou
 #Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
 strong_boundaries = [
         ("Top",0,0,1),
-        ("PtFix",2,3,0),
+        ("Top",1,1,0),
+        ("PtFix",3,3,0),
         ("Bottom",0,0,0),
         ("Bottom",1,1,0),
         ("Left",0,0,0),
         ("Left",1,1,0),
         ("Right",0,0,0),
-        ("Right",1,1,0)
+        ("Right",1,1,0),
         ]
 fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1)
 ii = 0

testcases/separation/mesh.geo

@@ -20,4 +20,4 @@ Physical Line("Right") = {3};
 Physical Line("Bottom") = {2};
 Physical Line("Top") = {4};
 Physical Surface("Domain") = {1};
-Physical Point("PtFix") = {1};
+Physical Point("PtFix") = {2};

validation/poiseuille/poiseuille.py

@@ -62,7 +62,7 @@ class Poiseuille(unittest.TestCase) :
 
         shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
-        outf = 100                                       # number of iterations between output files
+        outf = 10                                       # number of iterations between output files
 
         #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
         #Format: strong_boundaries = [(Boundary tag, Fluid field, Value)