Commit 2ebded12 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

true dt in prepro

parent af3ad5f7
......@@ -9,8 +9,8 @@
#define n_fields (DIMENSION+1)
#define newton_max_it 20
#define newton_atol 1e-12
#define newton_max_it 10
#define newton_atol 1e-6
#define newton_rtol 1e-5
#if DIMENSION == 2
......@@ -710,7 +710,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
hxtInitializeLinearSystems(0, NULL);
initialized = 1;
#ifdef HXT_HAVE_PETSC
hxtPETScInsertOptions("-pc_type ilu", "fluid");
hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid");
#endif
}
FluidProblem *problem = malloc(sizeof(FluidProblem));
......
......@@ -23,7 +23,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop);
p.position()[:, 1] += 0.52
print(len(p.position()[:, 1] ))
print("nb particles : ", len(p.position()[:, 1] ))
p.write(filename)
......@@ -68,12 +68,11 @@ p = scontact3.ParticleProblem(outputdir+"/part-00000")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 200
outf1 = 100000
outf = 100
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
#strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......@@ -88,10 +87,9 @@ for jj in range(5):
dt1 = dt
t = 0
if jj!=0:
print("ADAPT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n")
fluid.adapt_mesh(0.01,50,1,50,5e-4,50000)
while t<tEnd1:
forces = fluid.compute_node_force(dt)
forces = fluid.compute_node_force(dt1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
niter = fluid.implicit_euler(dt1)
t += dt1
......@@ -99,13 +97,14 @@ for jj in range(5):
dt1 *= 3
print(dt1)
# if ii %outf == 0 :
# ioutput = int(ii/outf) + 1
# p.write_vtk(outputdir, ioutput, t)
# fluid.export_vtk(outputdir, t, ioutput)
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
print("export %i" % ioutput)
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
ii = 0
#ii = 0
t = 0
tic = time.clock()
forces = g*p.mass()
......@@ -126,6 +125,7 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
print("export %i"% ioutput)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - 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