Commit 55510a90 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

merge in progress

parent d6163da8
......@@ -1096,7 +1096,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
break; /// equation is linear
//break; /// equation is linear
}
printf("total solve time : %g\n", totaltime);
free(dx);
......@@ -1440,11 +1440,11 @@ void fluid_problem_adapt_gen_mesh_in_advance(FluidProblem *problem, Mesh *mesh,
//sprintf(porosity_filename, "/Users/margauxboxho/migflow/testcases/Boycott/Porosity/Porosity_%d.dat",global_count);
sprintf(porosity_filename, "/Users/margauxboxho/migflow/testcases/depot-2d/Porosity/Porosity_%d.dat",global_count);
FILE *file = fopen(porosity_filename, "w");
/*FILE *file = fopen(porosity_filename, "w");
if(file == NULL){
printf("Error!\n");
exit(1);
}
}*/
// These definition should be more general
double h_max = lcmax;
......@@ -1560,13 +1560,13 @@ void fluid_problem_adapt_gen_mesh_in_advance(FluidProblem *problem, Mesh *mesh,
bool3 = (h_target<h_k*alpha) && (h_target>h_k/alpha);
h_target = bool1*fmax(h_k/alpha, h_min) + bool2*fmin(h_k*alpha, h_max) + bool3*h_k;
if (global_count%10==0) {
/*if (global_count%10==0) {
fprintf(file, "%d, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.5f, %.10f, %.5f, %.5f\n", iel, mesh->x[el[0]*3+0], mesh->x[el[0]*3+1], mesh->x[el[1]*3+0], mesh->x[el[1]*3+1], mesh->x[el[2]*3+0], mesh->x[el[2]*3+1], Pn[0], Pn[1], Pn[2], dPndx, dPndy, C[0], C[1], C[2], dCdx, dCdy, Un[0][0], Un[0][1], Un[0][2], Un[1][0], Un[1][1], Un[1][2], dUdx, dUdy, dVdx, dVdy, gradient_x[iel*3+0], gradient_x[iel*3+1], gradient_x[iel*3+2], gradient_y[iel*3+0], gradient_y[iel*3+1], gradient_y[iel*3+2], eta, h_target, h_k);
}
}*/
//mesh->min_h_target[iel] = fmin(h_target, mesh->min_h_target[iel]);
mesh->min_h_target[iel] = fmin(h_target, 1.1*mesh->min_h_target[iel]);
}
fclose(file);
//fclose(file);
free(edges_vec);
free(gradient_x);
free(gradient_y);
......
......@@ -124,7 +124,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
......@@ -135,7 +135,7 @@ if init:
#if jj!=0:
# Adaptation of the mesh.
# fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
while t < tEnd1 :
# Fluid solver
......@@ -143,7 +143,7 @@ if init:
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
# Increase time step if convergence
if niter < 5 and dt1 < 20:
dt1 *= 3
......@@ -156,7 +156,7 @@ t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
#
# COMPUTATION LOOP
......@@ -181,7 +181,7 @@ while t < tEnd :
p.iterate(dt/nsub, forces)
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
# Output files writting
if ii %outf == 0 :
......
......@@ -122,7 +122,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
......@@ -134,7 +134,7 @@ if init:
# Adaptation of the mesh.
#fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.adapt_gen_mesh(8*r, 4*r, 6602*2, "adapt/lc.pos")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
while t < tEnd1 :
# Fluid solver
......@@ -142,7 +142,7 @@ if init:
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
# Increase time step if convergence
if niter < 5 and dt1 < 20:
dt1 *= 3
......@@ -154,7 +154,7 @@ if init:
t = 0
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
#
# COMPUTATION LOOP
......@@ -188,7 +188,7 @@ while t < tEnd :
# load the save past state
fluid.import_vtk("output/fluid_00000.vtu")
p.read_vtk(outputdir,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
# create the .msh on the adapt/lc.pos
fluid.adapt_load_mesh("adapt/mesh.msh")
#fluid.implicit_euler(dt)
......@@ -208,7 +208,7 @@ while t < tEnd :
p.iterate(dt/nsub, forces)
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
# Output files writting
if ii %outf == 0 :
......
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