small oscillations

src/fluid_problem.c

@@ -233,56 +233,6 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
 #endif
 }
 
-// source term
-static void f_u_0(
-    double f[D], const double u[D], const double du[D][D], const double q, const double dq[D],
-    const double *uold, const double dt, const double rho) {
-  double oq = 1./q;//1./fmax(q,1e-3);
-  double divu = 0;
-  for (int j = 0; j < D; ++j) {
-    divu += du[j][j];
-  }
-  for (int j = 0; j < D; ++j) {
-    double utau = 0;
-    for (int k = 0; k < D; ++k) {
-      utau += u[k]*(du[j][k]-u[j]*dq[k]*oq);
-    }
-    f[j] = rho*((u[j]-uold[j])/dt + u[j]*oq*divu + utau*oq);
-  }
-}
-
-static void f_q_0(double *f, const double du[D][D], const double q, const double qold, const double dt)
-{
-  double divu = 0;
-  for (int j = 0; j < D; ++j) {
-    divu += du[j][j];
-  }
-  *f = divu+(q-qold)/dt;
-}
-
-//flux
-static void f_u_1(double f[D][D], double u[D], const double du[D][D], const double q, const double dq[D], const double p, const double mu) {
-  double tau[D][D];
-  double oq = 1./q;//1./fmax(q,1e-3);
-  for (int i = 0; i < D; ++i) {
-    for (int j = 0; j < D; ++j) {
-      tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
-    }
-  }
-  for (int j = 0; j < D; ++j) {
-    for (int k = 0; k < D; ++k) {
-      f[j][k] = mu*(tau[k][j]+tau[j][k]) + (k==j ? -q*p : 0);
-    }
-  }
-}
-
-static void f_q_1(double f[D], const double dp[D], const double epsilon, const double q)
-{
-  for (int j = 0; j < D; ++j) {
-    f[j] = epsilon*dp[j]*q;
-  }
-}
-
 static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
 {
   const Mesh *mesh = problem->mesh;
@@ -342,15 +292,11 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
 }
 
 static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, double c, double dt) {
-  double fu0[D];
-  double fu1[D][D];
-  double fq0;
-  double fq1[D];
   int P = problem->n_fluids*(D+1);
   double p = sol[P];
   double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
   double epsilon = problem->epsilon;
-  f0[P] = (1-c);
+  f0[P] = 1-c;
   f1[P*D+0] = 0;
   f1[P*D+1] = 0;
   for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
@@ -364,20 +310,28 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
     double u[D] = {sol[U],sol[U+1]};
     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]}};
-    f_u_0(fu0,u,du,q,dq,uold,dt,rho);
-    f_u_1(fu1,u,du,q,dq,p,mu);
-    f_q_0(&fq0,du,q,qold,dt);
-    f_q_1(fq1,dp,epsilon,q);
-    f0[U] = fu0[0];
-    f0[U+1] = fu0[1];
-    f0[Q] = fq0;
+
+    double oq = 1./q;//1./fmax(q,1e-3);
+    double divu = 0;
+    double tau[D][D];
+    double utau[D]={0.};
+    for (int i = 0; i < D; ++i) {
+      divu += du[i][i];
+      for (int j = 0; j < D; ++j) {
+        tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
+        utau[i] += u[j]*tau[i][j];
+      }
+    }
+    f0[Q] = divu+(q-qold)/dt;
     f0[P] -= q;
-    f1[U*D+0] = fu1[0][0];
-    f1[U*D+1] = fu1[0][1];
-    f1[U*D+D+0] = fu1[1][0];
-    f1[U*D+D+1] = fu1[1][1];
-    f1[Q*D+0] = fq1[0];
-    f1[Q*D+1] = fq1[1];
+    for (int i = 0; i < D; ++i) {
+      f0[Q] += epsilon*dp[i]*dq[i];
+      f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq*divu + utau[i]*oq) - 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] = epsilon*dp[i]*q;
+    }
   }
 }
 
@@ -716,9 +670,11 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
   for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
     problem->solution[i] = 0.;
   }
-  for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
-    for (int i = 0; i < mesh->n_nodes; ++i){
-      problem->solution[i*n_fields+D+ifluid*(D+1)] = 0.5;
+  for (int i = 0; i < mesh->n_nodes; ++i){
+    for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
+      double *x = mesh->x + i*3;
+      double a = (x[0]+5)/10;
+      problem->solution[i*n_fields+D+ifluid*(D+1)] = ifluid == 0 ? 0.8-0.6*a:0.2+0.6*a;
     }
   }
   problem->node_volume = NULL;

testcases/poiseuille2fluids/poiseuille.py

@@ -48,14 +48,14 @@ tEnd = 100000                                     # final time
 
 #numerical parameters
 lcmin = .1                                  # mesh size
-dt = 1                                       # time step
+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 = 10                                       # number of iterations between output files
+outf = 1                                       # 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)
@@ -65,27 +65,28 @@ p = 6
 strong_boundaries = [
         ("Bottom",u0,u0,0),
         ("Bottom",v0,v0,0.),
-        ("Bottom",q0,q0,0.5),
         ("Bottom",u1,u1,0),
         ("Bottom",v1,v1,0.),
-        #("Bottom",q1,q1,0.5),
+        #("Bottom",q0,q0,0.6),
         ("Top",u0,u0,0),
         ("Top",v0,v0,0.),
-        ("Top",q0,q0,0.5),
         ("Top",u1,u1,0),
+        #("Top",q0,q0,0.6),
         ("Top",v1,v1,0.),
-        #("Top",q1,q1,0.5),
-        ("Left",u0,u0,lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
-        ("Left",v0,v0,0),
-        ("Left",q0,q0,0.5),
-        ("Left",u1,u1,lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
-        ("Left",v1,v1,0),
-        #("Left",q1,q1,0.5),
-        ("Right",v0,v0,0),
-        ("Right",v1,v1,0),
-        #("Right",q0,q0,0.5),
-        #("Right",q1,q1,0.5),
-        #("Right",p,p,0)
+        ("LeftUp",u0,u0,lambda x : 0.8/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftUp",v0,v0,0),
+        ("LeftUp",q0,q0,0.8),
+        ("LeftUp",u1,u1,lambda x : 0.2/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftUp",v1,v1,0),
+        ("LeftDown",u0,u0,lambda x : 0.8/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftDown",v0,v0,0),
+        ("LeftDown",q0,q0,0.8),
+        ("LeftDown",u1,u1,lambda x : 0.2/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
+        ("LeftDown",v1,v1,0),
+        ("RightDown",v0,v0,0),
+        ("RightDown",v1,v1,0),
+        ("RightUp",v0,v0,0),
+        ("RightUp",v1,v1,0),
         ]
 fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,1,2)
 ii = 0