flag stab

fluid/fluid_problem.c

@@ -174,80 +174,6 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
   }
 }
 
-static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double ds[][D], const double sold[D], const double dsold[][D], const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double f01[][D])
-{
-  const int vid = wbnd->vid;
-  const int pid = wbnd->pid;
-  const int cid = wbnd->cid;
-  const int aid = wbnd->aid;
-  double p = s[P];
-  const int n_fields = fluid_problem_n_fields(problem);
-  double u[D], dp[D], uext[D];
-  double s_c = 0;
-  double unold = 0;
-  double cext = cid<0?c:data[cid];
-  double unext = 0;
-  double unmesh = 0;
-  double dcn = 0;
-  for (int iD = 0; iD < D; ++iD) {
-    u[iD] = s[U+iD];
-    dp[iD] = ds[P][iD];
-    unold += sold[U+iD]*n[iD];
-    unmesh += mesh_velocity[iD]*n[iD];
-    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]*c);
-    unext += uext[iD]*n[iD];
-    dcn += dc[iD]*n[iD];
-    s_c += (u[iD]-uext[iD])*n[iD];
-  }
-  double h = problem->element_size[eid];
-  if (wbnd->type == BND_WALL && pid >= 0) {
-      f0[P] = -(p-data[pid])/h*problem->taup[eid];
-      f00[P*n_fields+P] += -1/h*problem->taup[eid];
-  }
-  if(wbnd->type == BND_OPEN) {
-    f0[P] = unext;
-    for (int i = 0; i < D; ++i) {
-      f00[P*n_fields+U+i] += vid<0?n[i]:0;
-    }
-  }
-  for (int id = 0; id < D; ++id) {
-    f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
-    f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
-  }
-  double c_du_o_c[D][D];
-  double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
-  for (int iD = 0; iD < D; ++iD) {
-    for (int jD = 0; jD < D; ++jD)
-      c_du_o_c[iD][jD] = ds[U+iD][jD]-u[iD]/c*dc[jD];
-  }
-  if(wbnd->compute_viscous_term == 1){
-    for (int id = 0; id < D; ++id) {
-      f0[U+id] += sigma*(u[id]-uext[id]+s_c*n[id]);
-      f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
-      for (int jd = 0; jd <D; ++jd) {
-        f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
-        f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
-        f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
-        f01[(U+id)*n_fields+U+id][jd] -= mu*n[jd];
-        f01[(U+id)*n_fields+U+jd][id] -= mu*n[jd];
-      }
-    }
-  }
-  if (problem->advection) {
-    double unbnd = unold;
-    if (wbnd->type == BND_WALL) unbnd = unold/2;
-    else if(vid>0) unbnd = (unold+unext)/2;
-    if (unbnd - unmesh < 0) {
-      for (int id = 0; id < D; ++id) {
-        f0[U+id] += ((unbnd - unmesh)*(vid<0?0:uext[id])-(unold - unmesh)*u[id])*rho/c;
-        f00[(U+id)*n_fields+U+id] -= (unold-unmesh)*rho/c;
-      }
-    }
-  }
-}
-
-
-#if 0
 static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double ds[][D], const double sold[D], const double dsold[][D], const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double f01[][D])
 {
   const int vid = wbnd->vid;
@@ -321,9 +247,9 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     }
   }
 }
-#endif
 
 static double pow2(double a) {return a*a;}
+
 void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
   const Mesh *mesh = problem->mesh;
   problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
@@ -455,108 +381,6 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
   *dfddp = gammastar*dt/mass*vol;//-vol;
 }
 
-static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double f1[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {
-  double p = sol[P];
-  double taup = problem->taup[iel];
-  double tauc = problem->tauc[iel];
-  const int n_fields = fluid_problem_n_fields(problem);
-  double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
-  for (int i = 0; i < n_fields; ++i) {
-    for (int j = 0; j < n_fields; ++j) {
-      int r = i*n_fields+j;
-      f00[r] = 0;
-      for (int id = 0; id < D; ++id) {
-        f01[r][id] = 0;
-        f10[r][id] = 0;
-        for (int jd = 0; jd < D; ++jd) {
-          f11[r][id][jd] = 0;
-        }
-      }
-    }
-  }
-  if(problem->stab_param != 0)
-    taup = tauc = 0;
-  double divu = 0;
-  double nold = 0;
-  for (int iD = 0; iD < D; ++iD) {
-    u[iD] = sol[U+iD];
-    uold[iD] = solold[U+iD];
-    nold += uold[iD]*uold[iD];
-    umesh[iD] = mesh_velocity[iD];
-    dp[iD] = dsol[P][iD];
-    for (int jD = 0; jD < D; ++jD) {
-      du[iD][jD] = dsol[U+iD][jD];
-      duold[iD][jD] = dsolold[U+iD][jD];
-    }
-    divu += du[iD][iD];
-  }
-  nold = sqrt(nold);
-  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
-  //f00[P*n_fields+P] = 1/dt*0.1;
-  double *g = problem->g;
-  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
-  double drag = problem->volume_drag + problem->quadratic_drag*nold;
-  // f0[P] = p;
-  // f00[P*n_fields+P] = 1;
-  for (int i = 0; i < D; ++i) {
-    f0[U+i] =
-      + c*dp[i]
-      - g[i]*rhoreduced*c
-      - bf[i]*c
-      + drag*u[i];
-      // + 5.3e5*u[i];
-    f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
-    f01[(U+i)*n_fields+P][i] += c;
-    if(problem->temporal) {
-      f0[U+i] += rho*(u[i]-uold[i])/dt;// - rho;
-      f00[(U+i)*n_fields+U+i] += rho/dt;
-    }
-    if (problem->advection) {
-      // advection :
-      // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
-      //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
-      for (int j = 0; j < D; ++j) {
-        f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
-        f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
-        f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
-        f01[(U+i)*n_fields+(U+i)][j] += rho/c*uold[j];
-        f01[(U+i)*n_fields+(U+i)][j] -= rho/c*umesh[j];
-      }
-    }
-    for (int j = 0; j < D; ++j) {
-      f1[U+i][j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
-      f10[(U+i)*n_fields+U+j][j] += -mu*dc[i]/c;
-      f10[(U+i)*n_fields+U+i][i] += -mu*dc[j]/c;
-      f11[(U+i)*n_fields+U+j][j][i] += mu;
-      f11[(U+i)*n_fields+U+i][j][j] += mu;
-
-      // SUPG
-      double supg = uold[j]*taup;
-      f1[U+i][j] += supg*f0[U+i];
-      f10[(U+i)*n_fields+U+i][j] += supg*f00[(U+i)*n_fields+U+i];
-      f11[(U+i)*n_fields+P][j][i] += supg*f01[(U+i)*n_fields+P][i];
-      for (int k = 0; k < D; ++k)
-        f11[(U+i)*n_fields+(U+i)][j][k] += supg*f01[(U+i)*n_fields+(U+i)][k];
-    }
-    double lsic = tauc*rho;
-    f1[U+i][i] += ((c-cold)/dt+divu)*lsic;
-    for (int j = 0; j < D; ++j) {
-      f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
-    }
-    f1[P][i] = -u[i];
-    f10[P*n_fields+U+i][i] += -1;
-    // PSPG
-    double pspg = taup/rho;
-    f1[P][i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
-    f11[P*n_fields+P][i][i] += pspg*(f01[(U+i)*n_fields+P][i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
-    f10[P*n_fields+U+i][i] += pspg*f00[(U+i)*n_fields+(U+i)];
-    for (int j = 0; j < D; ++j) {
-      f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
-    }
-  }
-}
-
-#if 0
 static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double f1[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {
   double p = sol[P];
   double taup = problem->stability?problem->taup[iel]:0;
@@ -663,7 +487,6 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
     }
   }
 }
-#endif
 
 void weak_boundary_values(const FEFields *fields, const Mesh *mesh, const MeshBoundary *bnd, const WeakBoundary *wbnd, double *data);
 static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)

testcases/depot-2d/depot.py

@@ -102,7 +102,7 @@ ii        =  0
 #
 # FLUID PROBLEM
 #
-fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)
 # Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
 fluid.set_wall_boundary("Bottom")
@@ -120,9 +120,6 @@ tic = time.time()
 #
 while t < tEnd :
     time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), check_residual_norm=1e-3)
-    # fluid.implicit_euler(dt, check_residual_norm=1e-3)
-    # time_integration._advance_particles(p, fluid, dt, 5, 1e-8)
-    # fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
 
     t += dt
     # Output files writting

validation/darcy/darcy.py

@@ -104,7 +104,7 @@ class Darcy(unittest.TestCase) :
             p = scontact.ParticleProblem(2)
             p.read_vtk(outputdir,0)
             ii = 0
-            fluid = mbfluid.FluidProblem(2,g,mu,rho,order=[1,1,1], solver="petsc4py", solver_options="-pc_type lu", stability=True, drag_in_stab=True)
+            fluid = mbfluid.FluidProblem(2,g,mu,rho,order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)#, stability=True, drag_in_stab=True)
             fluid.load_msh("mesh.msh")
             U = 0.00001                                             # averaged fluid velocity
             V = U/(1-compacity)                                     # fluid velocity

validation/poiseuille/poiseuille.py

@@ -64,7 +64,7 @@ class Poiseuille(unittest.TestCase) :
         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)
-        fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[1,1,1], solver="petsc4py", solver_options="-pc_type lu", stability=True)
+        fluid = mbfluid.FluidProblem(2,g,nu*rho,rho, order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)
         fluid.load_msh("mesh.msh")
         fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
         fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
@@ -95,7 +95,7 @@ class Poiseuille(unittest.TestCase) :
             print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.perf_counter() - tic))
             
         s = fluid.velocity()
-        x = fluid.coordinates(order=1)
+        x = fluid.coordinates(order=2)
         vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
         vS = np.sum((1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2)
         print('Error', (vel.sum())**.5)