Commit 320a778d authored by Michel Henry's avatar Michel Henry
Browse files

fix testcases, remove state

parent f0e50faa
Pipeline #9343 passed with stages
in 3 minutes and 52 seconds
......@@ -39,16 +39,13 @@ os.chdir(dir_path)
outputdir = "outputDepotNP"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "meshDepot.geo","-clscale","1"])
#subprocess.call(["gmsh", "-2", "meshDepot.geo","-clscale","1"])
t = 0
ii = 0
#physical parameters
# g = np.array([0.001,-0.001])
# g = np.array([0.1,-0.5]) # gravity
g = np.array([0.,-9.81]) # gravity
rho = 1000 # fluid density
rhop = 1500 # particle density
......@@ -61,7 +58,7 @@ H = 0.6
#numerical parameters
lcmin = .1 # mesh size
dt = 1e-3 # time step
dt = 2.5e-3 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......@@ -70,7 +67,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+0.1,L-2*r-0.1,2.1*r):
for x in np.arange(2*r,L-2*r,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)
......@@ -111,7 +108,6 @@ while t < tEnd :
#Fluid solver
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
# state = p.state()
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = np.where(p.position()[:,0] <= 0)
p.position()[ind,0] += L
......@@ -123,6 +119,7 @@ while t < tEnd :
ioutput = int(ii/outf) + 1
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))
......
......@@ -96,8 +96,8 @@ 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 = p.mass()*g
ii = 0
......@@ -105,22 +105,18 @@ tic = time.time()
while t < tEnd :
# while ii < 10:
#Fluid solver
oldstate = np.copy(p.state())
time_integration.iterate(fluid,p,dt,external_particles_forces=G)
t += dt
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)
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)
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = p.position()[:,0] <= 0
p.position()[ind,0] += L
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#Output files writting
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