Commit 50620881 authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

fix boundary

parent 2be74e44
Pipeline #5437 passed with stage
in 27 seconds
......@@ -216,11 +216,11 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
f00[(U+id)*n_fields+P] += (pid<0?0:-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*D+(U+id)*D+jd] += mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] += mu*n[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*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
}
}
double unbnd = unold;
......
......@@ -68,7 +68,7 @@ if not os.path.isdir(outputdir) :
# Physical parameters
g = -9.81 # gravity
r = 2.e-3 # particles radius
r = 3e-3 # particles radius
rhop = 2000 # particles density
rho = 1000 # fluid density
nu = 1.e-6 # kinematic viscosity
......@@ -119,7 +119,7 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_f
#
while t < tEnd :
# Fluid solver
fluid.implicit_euler(dt, newton_atol=5e-6, newton_rtol=1e-1, newton_max_it=35, reduced_gravity=0, stab_param=0.)
fluid.implicit_euler(dt, newton_atol=5e-6, newton_rtol=0.1, newton_max_it=100, reduced_gravity=0, stab_param=0.)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
......
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