Skip to content

GitLab

Commit c7b86890 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

newton tolerance

parent f8aad78e
Pipeline #4535 failed with stage
in 20 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ c7b86890
......@@ -164,7 +164,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
static void node_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip) {
static void particle_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip, int in_derivative) {
double p = sol[P];
double u[D], uold[D], dp[D];
......@@ -193,8 +193,7 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
}
double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
double rhop = mass/vol;
double g = problem->g*(rhop-problem->rho)/rhop;
double g = problem->g;
double gamma;
if(problem->n_fluids == 2){
......@@ -207,10 +206,12 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
}
double fcoeff = mass/(mass+dt*gamma);
for (int j = 0; j < D; ++j){
problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
if(!in_derivative)
problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
}
problem->particle_force[ip*D+1] += fcoeff*g*mass;
if(!in_derivative)
problem->particle_force[ip*D+1] += fcoeff*g*mass;
f0[U+1] += -fcoeff*gamma*dt*g;
}
......@@ -280,7 +281,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
node_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip);
particle_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield];
......@@ -289,12 +290,12 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
node_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip);
particle_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
node_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip);
particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1);
ds[jfield*D+jd] = store;
}
int N2 = n_fields*N_SF;
......@@ -420,7 +421,7 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
f[P] = 0;
for (int id = 0; id < D; ++id) {
f[U+id] = 0*p*n[id];
f[U+id] = p*n[id];
}
if(problem->n_fluids == 2){
......@@ -546,7 +547,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
}
}
f0[P] = (c-cold)/dt+divu;
f0[P] = (c-cold)/dt;
if(problem->n_fluids == 2){
f0[A] = (a*c-aold*cold)/dt;
......@@ -555,12 +556,12 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
for (int i = 0; i < D; ++i) {
double dudt = (u[i]-uold[i])/dt;
double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -c*problem->g : 0.);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -0*c*rho*problem->g : 0.);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -c*rho*problem->g : 0.);
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup*0;
}
f1[(U+i)*D+i] += 0*divu*tauc*rho;
f1[P*D+i] = -u[i]*0 + 1e-6*dp[i];
f1[P*D+i] = -u[i] + 1e-6*(dp[i]+(i==1 ? -rho*problem->g : 0.));
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
......@@ -1138,7 +1139,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double C[N_N],Un[D][N_N],An[N_N],Pn[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
Pn[i] = solution[el[i]*n_fields+D]-0*problem->g *problem->rho[0]*mesh->x[el[i]*3+1];
Pn[i] = solution[el[i]*n_fields+D]-problem->g*problem->rho[0]*mesh->x[el[i]*3+1];
if(problem->n_fluids ==2) An[i] = solution[el[i]*n_fields+D+1];
for (int j=0; j<D; ++j) {
Un[j][i] = solution[el[i]*n_fields+j];
......@@ -1212,7 +1213,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double C[N_N],Un[D][N_N],An[N_N],Pn[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
Pn[i] = solution[el[i]*n_fields+D]-0*problem->g *problem->rho[0]*mesh->x[el[i]*3+1];
Pn[i] = solution[el[i]*n_fields+D]-problem->g *problem->rho[0]*mesh->x[el[i]*3+1];
if(problem->n_fluids ==2) An[i] = solution[el[i]*n_fields+D+1];
for (int j=0; j<D; ++j) {
Un[j][i] = solution[el[i]*n_fields+j];
......
testcases/drop-3d/InteractingDrops/Vert/2VertDrops.py
View file @ c7b86890
......@@ -46,6 +46,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
rhop -- particles density
compacity -- initial compacity in the cloud
"""
random.seed(0)
# Particles structure builder
p = scontact.ParticleProblem(3)
# Load mesh.msh file specifying physical boundaries names
......@@ -81,7 +82,7 @@ if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Enable computation of the initial velocity
init = 1
init = False
# Physical parameters
g = -9.81 # gravity
......@@ -112,7 +113,7 @@ p.read_vtk(outputdir,0)
t = 0
ii = 0
jj = 0
tEnd1 = 0.2
tEnd1 = 100
#
# FLUID PROBLEM
......@@ -133,16 +134,16 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
# The initial mesh is obtained by adapting for times the mesh on the converged fields
if init:
for jj in range(4):
dt1 = dt/100
dt1 = dt
t = 0
if jj!=0:
# Adaptation of the mesh.
fluid.adapt_mesh(8e-3,8e-4,250000,cmax=.1,cmin=20)
fluid.adapt_mesh(8e-3,8e-4,50000,cmax=.1,cmin=20)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd1 :
# Fluid solver
niter = fluid.implicit_euler(dt1)
niter = fluid.implicit_euler(dt1,newton_atol=1e-8)
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
......@@ -156,20 +157,23 @@ if init:
t += dt1
t = 0
fluid.export_vtk(outputdir,0,0)
tic = time.time()
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#fluid.solution()[:,3] = (fluid.coordinates()[:,1]-0.25)*g*rho
ioutput = 0
fluid.export_vtk(outputdir,ioutput,0)
#
# COMPUTATION LOOP
#
while t < tEnd :
# Fluid solver
fluid.implicit_euler(dt)
fluid.implicit_euler(dt,newton_atol=1e-8,newton_max_it=1000)
# Adaptation of the mesh. Args are maximal mesh radius, minimal mesh radius and number of elements. Optional arguments can be given to weigh the max and min gradient used in the adapt formula
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(8e-3,8e-4,250000,cmin=20,cmax=0.1)
fluid.adapt_mesh(8e-3,8e-4,150000,cmin=20,cmax=0.1)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
......@@ -187,8 +191,8 @@ while t < tEnd :
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
ioutput += 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
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