Commit cc9c8895 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

fix reload

parent 6c70e9c6
Pipeline #3910 passed with stage
in 59 seconds
......@@ -1379,7 +1379,7 @@ void fluid_problem_after_import(FluidProblem *problem)
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload)
{
struct timespec tic, toc;
......@@ -1423,8 +1423,10 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_velocity[i*D+k] = velocity[i*D+k];
}
}
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
if (!reload){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
}
}
clock_get(&toc);
printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
......
......@@ -87,7 +87,7 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename);
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries, int n_fluids);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el);
......
......@@ -29,6 +29,7 @@
int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
{
const int n_fields = fluid_problem_n_fields(problem);
printf("n_fields = %d\n",n_fields);
MbXmlReader *r = mb_xml_reader_create(filename);
MbXmlElement fileelement,gridelement,pieceelement;
......@@ -102,7 +103,7 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
double *v = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
for (int i = 0; i < m->n_nodes; ++i)
for (int d = 0; d < DIMENSION; ++d)
problem->solution[n_fields*i+d] = v[i*DIMENSION+d];
problem->solution[n_fields*i+d] = v[i*3+d];
free(v);
}else if(!strcasecmp(a.name,"Concentration")) {
double *c = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
......
......@@ -31,7 +31,7 @@ import random
outputdir = "outputTest"
outputdir = "outputTestReload"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -57,14 +57,14 @@ rhop = 2500
#numerical parameters
lcmin = .1 # mesh size
dt = .00025 # time step
iReload = 3
iReload = 30
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,iReload)
outf = 10 # number of iterations between output files
outf = 1 # number of iterations between output files
ii = 0
t = 0
......@@ -75,7 +75,7 @@ def outerBndV(x) :
strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
# ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)
]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries,1)
fluid = fluid.fluid_problem("adapt/mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries,1)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
......@@ -83,12 +83,15 @@ fluid.set_weak_boundary("Injection","Inflow")
fluid.import_vtk(outputdir+"/fluid_%05d.vtu"%iReload)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
rel = 1
#set initial_condition
#fluid.solution()[:,3] = 1
ioutput = int(iReload) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),rel)
ii = iReload*outf+1
t = ii*dt
tic = time.clock()
......@@ -96,9 +99,11 @@ while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
if (ii%5==0 and ii != 0):
#fluid.adapt_mesh(1e-3,1e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),rel)
else :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),rel)
rel = 0
#Fluid solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
......
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