correction sign boundary term of f(q)

src/fluid_problem.c

@@ -254,9 +254,9 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
     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[U+id] = 0;//-p*q*n[id];
+      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;
@@ -280,9 +280,9 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
     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[U+id] = 0;//-p*q*n[id];
+      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;
@@ -307,7 +307,7 @@ 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];
     }
   }
@@ -357,7 +357,6 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
       //f0[Q] += epsilonp*dq[i]*dp[i];
       f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) -p*dq[i]+ (i==(D-1) ? -q*rho*problem->g : 0);
       for (int j = 0; j < D; ++j) {
-        f0[U+i] += mu*(tau[i][j]+tau[j][i])*dq[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]+q*epsilonp*dp[i];
@@ -417,7 +416,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
       invDD(dxdxi, dxidx);
       grad_shape_functions(dxidx, dphi);
       const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
-      const double n[2] = {-dx[1]/l,dx[0]/l};
+      const double n[2] = {dx[1]/l,-dx[0]/l};
 
       for (int iqp = 0; iqp < N_LQP; ++iqp) {
         double xi[2];