Commit f0e50faa authored by Michel Henry's avatar Michel Henry
Browse files

fix state

parent b60343a1
......@@ -61,7 +61,7 @@ H = 0.6
#numerical parameters
lcmin = .1 # mesh size
dt = 2.5e-3 # time step
dt = 1e-3 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......@@ -70,7 +70,7 @@ shutil.copy("meshDepot.msh", outputdir +"/mesh.msh")
p = scontact.ParticleProblem(2)
p.load_msh_boundaries("meshDepot.msh", ["Bottom", "Top", "Right", "Left"])
for x in np.arange(2*r,L-2*r,2.1*r):
for x in np.arange(2*r+0.1,L-2*r-0.1,2.1*r):
for y in np.arange(H-0.2,H-r,2.1*r):
R = r*(1-np.random.random()/4)
p.add_particle((x+R,y),R,R**2*np.pi*rhop)
......@@ -91,16 +91,16 @@ fluid.set_wall_boundary("Top",velocity=[0,0])
# fluid.set_strong_boundary("Bottom",0,0)
# if strong boundary on periodic line, it should be forced on both sides
fluid.set_strong_boundary("Right",2,0)
fluid.set_strong_boundary("Left",2,0)
# fluid.set_strong_boundary("Right",2,0)
# fluid.set_strong_boundary("Left",2,0)
ii = 0
t = 0
#set initial_condition
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.write_vtk(outputdir,0,0)
G = np.zeros((p.n_particles(),2))
G[:,1] = p.mass()[:,0]*g[1]
......@@ -109,22 +109,19 @@ ii = 0
tic = time.time()
while t < tEnd :
#Fluid solver
oldstate = np.copy(p.state())
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
state = p.state()
state.x[:,0] = np.fmod(state.x[:,0],L)
ind = np.where(state.x[:,0] <= 0)
state.x[ind,0] += L
p.set_state(state)
# state = p.state()
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = np.where(p.position()[:,0] <= 0)
p.position()[ind,0] += L
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.set_particles(p.mass(), p.volume(), oldstate,p.contact_forces(),init=False)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=False)
#Output files writting
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
ii += 1
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