Commit 44f6ef0e authored by Matthieu Constant's avatar Matthieu Constant
Browse files

bug bug bug

parent 89112986
Pipeline #6604 failed with stage
in 1 minute and 3 seconds
......@@ -337,6 +337,8 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
static int crado = 0;
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip, double *upred, double *vold) {
const double *contact = &problem->particle_contact[ip*D];
double u[D], uold[D], dp[D];
......@@ -350,7 +352,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
double Due[D];
const double *up = &problem->particle_velocity[ip*D];
for (int j = 0; j < D; ++j) {
Due[j] = problem->particle_velocity[ip*D+j]-upred[j]/c;
Due[j] = vold[j]-upred[j]/c;
}
double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
......@@ -369,16 +371,17 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
else{
gamma = particle_drag_coeff(Due,mu,rho,vol,c);
}
double aaa = 0;
gamma = gamma*problem->drag_coeff_factor;
double gammastar = gamma/(1+dt/mass*gamma);
double gammastar = gamma/(1+dt*(1+crado*aaa)/mass*gamma);
for (int i = 0; i < D; ++i) {
double upstar = vold[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
double upstar = vold[i]+dt*(1+crado*aaa)/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
double drag = gammastar*(upstar-u[i]/c);
problem->particle_force[ip*D+i] = -drag-vol*dp[i];
f[U+i] = -drag;//-vol*dp[i];
}
*dfdu = gammastar/c;
*dfddp = gammastar*dt/mass*vol;//-vol;
*dfddp = gammastar*dt*(1+crado*aaa)/mass*vol;//-vol;
}
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
......@@ -1255,6 +1258,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double *ufpre
node_force_volume(problem,problem->solution_explicit,ufpred,vpold,dt,NULL,NULL);
//break; /// equation is linear
}
crado = 1-crado;
printf("total solve time : %g\n", totaltime);
free(dx);
......
......@@ -41,28 +41,29 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.implicit_euler(dt,ufpred=sf0[:,:-1],vpold=vp0)
sf1 = np.copy(fluid.solution())
f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
_advance_particles(particles,f0,dt,min_nsub,contact_tol)
_advance_particles(particles,f0,dt/2,min_nsub,contact_tol)
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
particles.position()[:,:] = pp0
particles.contact_forces()[:,:] = cp0
#particles.position()[:,:] = pp0
#particles.contact_forces()[:,:] = cp0
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces(), reload=True)
#particles.position()[:,:] = pp0
# Corrector
#
# Compute the fluid solution from the previous time-step one using the predicted solid velocities in the fluid-grain interaction force
fluid.solution()[:,:] = sf0
fluid.implicit_euler(dt,ufpred=sf1[:,:-1],vpold=vp0)
fluid.implicit_euler(dt,ufpred=sf0[:,:-1],vpold=particles.velocity())
# Fluid solution is a weighted average of the predicted one and the corrected one.
# The alpha parametre is the weight
fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
#fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
f1 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
# For two fluids flows, advance the concentration field with the fluid velocity.
# The number of sub-time steps for the advection is automatically computed.
if fluid.n_fluids() == 2 :
fluid.advance_concentration(dt)
# Reset solid velocities
particles.velocity()[:,:] = vp0
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol)
_advance_particles(particles,f1,dt/2,min_nsub,contact_tol)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces())
......
......@@ -63,7 +63,7 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
outputdir = "output1"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment