Commit 9865d7dd authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

magic_pressure

parent 414d0f75
Pipeline #10140 failed with stages
in 62 minutes and 15 seconds
......@@ -386,7 +386,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
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];
double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D], dpold[D];
if(problem->stab_param != 0)
taup = tauc = 0;
double divu = 0;
......@@ -398,6 +398,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
nold += uold[iD]*uold[iD];
umesh[iD] = mesh_velocity[iD];
dp[iD] = dsol[P*D+iD];
dpold[iD] = dsolold[P*D+iD];
for (int jD = 0; jD < D; ++jD) {
du[iD][jD] = dsol[(U+iD)*D+jD];
duold[iD][jD] = dsolold[(U+iD)*D+jD];
......@@ -418,7 +419,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double drag = problem->volume_drag + problem->quadratic_drag*nold;
for (int i = 0; i < D; ++i) {
f0[U+i] =
+ c*(dp[i]- g[i]*rhoreduced- bf[i])
dp[i]+ c*(- g[i]*rhoreduced- bf[i])-(1-c)*dpold[i]
+ drag*u[i];
//+ 5.3e5*u[i];
f00[(U+i)*n_fields+U+i] = drag;//5.3e5;
......@@ -426,7 +427,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f0[U+i] += rho*(u[i]-uold[i])/dt;
f00[(U+i)*n_fields+U+i] += rho/dt;
}
f01[(U+i)*n_fields*D+P*D+i] = c;
f01[(U+i)*n_fields*D+P*D+i] = 1;
if (problem->advection) {
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
......
......@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0,solver="petsc",solver_options="-pc_type lu",usolid=True)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1,solver="petsc",solver_options="-pc_type lu",usolid=True)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
......@@ -118,7 +118,7 @@ tic = time.time()
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), use_predictor_corrector=False)
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), use_predictor_corrector=True)
t += dt
# Output files writting
......
L = 0.2;
H = 0.3;
L = 0.1;
H = 0.1;
y = 0;
lc = 0.01;
......
Supports Markdown
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