Commit dd9b6cae authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

fix sub dt in boundary forces

parent 445d837c
Pipeline #7296 passed with stage
in 2 minutes and 43 seconds
......@@ -195,7 +195,7 @@ class ParticleProblem :
def contact_forces(self):
return (self._get_matrix("ContactForces",self._dim))
def get_boundary_forces(self,dt,tag="default") :
def get_boundary_impulsion(self,tag="default") :
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
def compute_fn_ft(contact_name,objects,tag) :
......@@ -214,7 +214,7 @@ class ParticleProblem :
f += f2
if self._dim == 3:
f+= compute_fn_ft("particle_triangle",self.triangles(),tag)
return f/dt
return f
def disks(self) :
return self._get_array("Disk",self._disktype)
......
import numpy as np
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0) :
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None) :
"""Contact solver for the grains
Keyword arguments:
......@@ -26,16 +26,21 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0) :
if iteration == 1:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
if after_sub_iter :
after_sub_iter(dt)
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)
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,iteration=iteration+1,after_sub_iter=after_sub_iter)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
else :
if after_sub_iter :
after_sub_iter(dt)
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5) :
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) :
"""Predictor-corrector scheme to solve fluid and grains.
Keyword arguments:
......@@ -83,11 +88,11 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.velocity()[:,:] = vp0
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol)
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), 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) :
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
......@@ -115,13 +120,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)
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter)
if fluid is not None and particles is not None and particles.r() is not None:
if fixed_grains:
fluid.implicit_euler(dt)
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces)
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter)
def particle_changed(fluid,p) :
co = fluid.old_porosity().copy()
......
......@@ -65,11 +65,11 @@ def genInitialPosition(filename, rmin, rmax, H, lx, ly, rhop) :
p.add_particle((x[i], y[i]), r[i], r[i]**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
def total_boundary_forces(p,dt):
fn, ft = p.get_boundary_forces(dt,"left")
fn2, ft2 = p.get_boundary_forces(dt,"right")
fn3, ft3 = p.get_boundary_forces(dt,"bottom")
fn4, ft4 = p.get_boundary_forces(dt,"top")
def total_boundary_impulsion(p):
fn, ft = p.get_boundary_impulsion("left")
fn2, ft2 = p.get_boundary_impulsion("right")
fn3, ft3 = p.get_boundary_impulsion("bottom")
fn4, ft4 = p.get_boundary_impulsion("top")
print(fn,ft)
print(fn2,ft2)
print(fn3,ft3)
......@@ -118,14 +118,18 @@ class Weight(unittest.TestCase) :
#
forces = g*p.mass()
while t < tEnd :
p.iterate(dt, forces, tol=1e-8)
bnd_forces = np.zeros((2,))
def accumulate(bnd_forces) :
bnd_forces += total_boundary_impulsion(p)/dt
time_integration.iterate(None,p,dt,min_nsub=3,external_particles_forces=forces,
after_sub_iter = lambda subdt : accumulate(bnd_forces))
#p.iterate(dt, forces, tol=1e-8)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
error = np.sqrt(sum(sum(forces) - total_boundary_forces(p,dt))**2)
error = np.sqrt(sum(sum(forces) - bnd_forces)**2)
error /= np.sqrt(sum(sum(forces)**2))
print(error*100,"%",p.n_particles())
self.assertLess(error,.5/100., "error is too large in weight")
Supports Markdown
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