Commit de313e53 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

updating splitting in particle solver

parent 24fb090b
Pipeline #8768 passed with stages
in 5 minutes and 26 seconds
......@@ -22,7 +22,7 @@ import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None,max_split=1,nsubt=[]) :
def _advance_particles(particles, f, dt, min_nsub,contact_tol,max_nsub=None,after_sub_iter=None,n_divide=None) :
"""Contact solver for the grains
Keyword arguments:
......@@ -31,9 +31,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
dt -- time step
min_nsub -- minimal nsub for the particles iterations
contact_tol -- tolerance for the contact solver
iteration -- internal argument -- used to keep track of the number of time step splitting
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
# Compute free velocities of the grains
v = particles.velocity()
......@@ -47,27 +45,32 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
else:
nsub = 1
vmax = 0
if max_nsub is None:
max_nsub = min_nsub*8
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()) if (particles.r() is not None and particles.r().size !=0) else 0 )
nsubt2 = nsubt + [nsub]
if n_divide is None:
n_divide = nsub
else:
n_divide *= nsub
#If time step was split too many times, ignore the convergence and move the grains
if iteration == max_split:
if min_nsub > max_nsub:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
if after_sub_iter :
after_sub_iter(nsubt2)
after_sub_iter(n_divide)
return
# For each sub-time step solve contacts (do NLGS iterations) and move the grains
for i in range(nsub) :
if not particles.iterate(dt/nsub, f, tol=contact_tol):
print("Splitting time-step to level %d..."%(iteration+1))
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,iteration=iteration+1,after_sub_iter=after_sub_iter,nsubt=nsubt2)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
print("Splitting time-step")
_advance_particles(particles,f,dt/nsub,min_nsub*2,contact_tol/2,max_nsub/nsub,after_sub_iter=after_sub_iter,n_divide=n_divide)
print("Convergence reached for subinvervals")
else :
if after_sub_iter :
after_sub_iter(nsubt2)
after_sub_iter(n_divide)
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_split=1) :
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_nsub=None) :
"""Predictor-corrector scheme to solve fluid and grains.
Keyword arguments:
......@@ -99,7 +102,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.implicit_euler(dt)
sf1 = np.copy(fluid.solution())
f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_nsub=max_nsub)
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
fluid.particle_velocity()[:] = particles.velocity()
......@@ -119,12 +122,12 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.restore_state()
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
# Give to the fluid the solid information
fluid.move_particles(particles.position(), particles.velocity(), particles.contact_forces())
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_nsub=None) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
......@@ -154,13 +157,13 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if fluid.n_fluids() == 2 :
fluid.advance_concentration(dt)
if particles is not None and fluid is None:
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
if fluid is not None and particles is not None and particles.r() is not None and particles.r().size !=0:
if fixed_grains:
fluid.implicit_euler(dt)
fluid.move_particles(particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
class FreeSurface():
......@@ -359,4 +362,4 @@ class FreeSurface():
def corr_un(self, vbox):
integral_un = self.compute_integral_un()
print(f"Integral un before correction: {integral_un - vbox[-1]}")
self.fluid.solution()[self.fs_nodes,1] += (vbox[-1]-integral_un)
\ No newline at end of file
self.fluid.solution()[self.fs_nodes,1] += (vbox[-1]-integral_un)
......@@ -42,7 +42,7 @@ dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
outf = 1000 #iterations between data frames
# Define output directory
outputdir = "notOkTrial/" + vitesse[index]
if not os.path.isdir(outputdir) :
......@@ -53,14 +53,14 @@ filename = outputdir + "/Drag.csv"
#
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/outputNotOk",10)
p.read_vtk("depot/outputFluid",10)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.61))
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/outputNotOk",10)
p2.read_vtk("depot/outputFluid",10)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel")
p2.load_msh_boundaries("meshP.msh", ["Inner","Left","Right"],material = "Steel")
p2.set_friction_coefficient(0.3,"Sand","Sand")
p2.set_friction_coefficient(0.3,"Sand","Steel")
p2.write_vtk(outputdir, 0, 0)
......@@ -82,6 +82,7 @@ p2.write_vtk(outputdir, 0, 0)
def accumulate(F,nsub):
F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/nsub
while t < tEnd :
print(vit)
p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
if t > 2:
p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
......@@ -93,7 +94,7 @@ while t < tEnd :
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.55))
F = np.zeros(2)
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub))
time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub),max_nsub=4)
with open(filename,"a") as file1:
#F += fluid.boundary_force()[0,:]
file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
......
......@@ -113,12 +113,13 @@ class Weight(unittest.TestCase) :
# COMPUTATION LOOP
#
forces = g*p.mass()
def accumulate(bnd_forces,nsub) :
bnd_forces += total_boundary_force(p)/nsub
print("coucou")
def accumulate(bnd_forces,n_divide) :
bnd_forces += total_boundary_force(p)/n_divide
while t < tEnd :
bnd_forces = np.zeros((2,))
time_integration.iterate(None,p,dt,min_nsub=1,external_particles_forces=forces,
after_sub_iter = lambda nsub : accumulate(bnd_forces,nsub))
after_sub_iter = lambda n_divide : accumulate(bnd_forces,n_divide))
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
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