Commit 89efe011 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

ajout cas test automatique

parent a74c070d
Pipeline #7156 passed with stage
in 2 minutes and 29 seconds
......@@ -194,6 +194,20 @@ class ParticleProblem :
def contact_forces(self):
return (self._get_matrix("ContactForces",self._dim))
def get_boundary_forces(self,dt,tag="default") :
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
odk,vdk,basisdk =self.get_contacts("particle_disk")
osg,vsg,basissg =self.get_contacts("particle_segment")
otr,vtr,basistr =self.get_contacts("particle_triangle")
o,v,basis = np.concatenate((odk,osg,otr)),np.concatenate((vdk,vsg,vtr)),np.concatenate((basisdk,basissg,basistr))
m = self.mass()[o[o[:,0]==tag,1]]
fn = np.sum(m*(basis[o[:,0]==tag,0,:]*v[o[:,0]==tag,:][:,[0]])/dt,axis=0)
ft = np.sum(m*(basis[o[:,0]==tag,1,:]*v[o[:,0]==tag,:][:,[1]])/dt,axis=0)
if self._dim == 3:
ft += np.sum(m*(basis[o[:,0]==tag,1,:]*v[o[:,0]==tag,:][:,[2]])/dt,axis=0)
return fn, ft
def set_boundary_velocity(self, tag, v) :
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
......
import numpy as np
def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iteration=0) :
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0) :
"""Contact solver for the grains
Keyword arguments:
......@@ -9,7 +9,6 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
dt -- time step
min_nsub -- minimal nsub for the particles iterations
contact_tol -- tolerance for the contact solver
momentum -- array to be filled with the forces on the boundaries due to the contacts
iteration -- internal argument -- used to keep track of the number of time step splitting
"""
# Compute free velocities of the grains
......@@ -32,16 +31,11 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
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,momentum=momentum,iteration=iteration+1)
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,iteration=iteration+1)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
if momentum is not None:
#print(particles.momentum()*nsub/dt)
momentum += particles.momentum()
if momentum is not None and iteration == 0:
momentum /= dt
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,momentum=None) :
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5) :
"""Predictor-corrector scheme to solve fluid and grains.
Keyword arguments:
......@@ -52,7 +46,6 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
contact_tol -- tolerance for the contact solver
external_particles_forces -- external forces applied on the particles
alpha -- parametre of the predictor-corrector scheme [alpha*f(n)+(1-alpha)*f(n+1)]
momentum -- array to be filled with the forces on the boundaries due to the contacts
"""
if external_particles_forces is None:
external_particles_forces = np.zeros_like(particles.velocity())
......@@ -90,11 +83,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,momentum=momentum)
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol)
# 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, momentum=None) :
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
......@@ -105,7 +98,6 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
contact_tol -- tolerance for the contact solver
external_particles_forces -- vector of external forces applied on the particles
fixed_grains -- boolean variable specifying if the grains are fixed in the fluid
momentum -- array to be filled with the forces on the boundaries due to the contacts
Raises:
ValueError -- fluid and particles cannot be both None
......@@ -123,13 +115,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,momentum=momentum)
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol)
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,momentum=momentum)
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces)
def particle_changed(fluid,p) :
co = fluid.old_porosity().copy()
......
import unittest
import sys
all_tests = [
all_tests = ["weight.weight",
"poiseuille.poiseuille",
"cavity.cavity",
"oneGrain.oneGrain",
......
......@@ -20,20 +20,17 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Bidimensional particles sedimentation in fluid
from migflow import fluid
from migflow import fluid as mbfluid
from migflow import scontact
from migflow import time_integration
import numpy as np
import os
import subprocess
import time
import shutil
import random
import unittest
def genInitialPosition(filename, rmin, rmax, H, lx, ly, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
......@@ -68,64 +65,64 @@ 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)
# Define output directory
outputdir = "depot"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = 9.81 # gravity
r = 5e-3 # particles radius
rhop = 1500 # particles density
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 0.2e-3 # time step
tEnd = 100 # final time
# Geometrical parameters
ly = 1
lx = 0.25
H = 0.4
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r*0.7,r, H, lx, ly, rhop)
p = scontact.ParticleProblem(2,True,True)
p.set_friction_coefficient(0.2)
p.read_vtk(outputdir,0)
p.velocity()[:,0] = np.random.uniform(0,1e-3,p.velocity().shape[0])
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
# Initial time and iteration
t = 0
ii = 0
#
# COMPUTATION LOOP
#
forces = np.zeros_like(p.position())
forces[:,[1]] = -g*p.mass()
def get_momentum_boundary() :
o,v,basis = p.get_contacts("particle_segment")
m = p.mass()[o[:,1]]
fn = np.sum(m*(basis[:,0,:]*v[:,[0]])/0.2e-3)
ft = np.sum(m*(basis[:,1,:]*v[:,[1]])/0.2e-3)
return fn, ft
while t < tEnd :
vn2 = np.copy(p.velocity())
p.iterate(dt, forces, tol=1e-8)
print(np.sum(g*p.mass()))
t += dt
fn,ft = get_momentum_boundary()
print(f"momentum : {fn+ft}")
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
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")
return fn + fn2 + fn3 + fn4 + ft + ft2 + ft3 + ft4
class Weight(unittest.TestCase) :
def runTest(self) :
dir_path = os.path.dirname(os.path.realpath(__file__))
os.chdir(dir_path)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
outf = 100
outputdir = "depot"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81])
r = 6e-3 # particles radius
rhop = 1500 # particles density
# Numerical parameters
dt = 2e-3 # time step
tEnd = 1 # final time
# Geometrical parameters
ly = 0.5
lx = 0.125
H = 0.15
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r*0.7,r, H, lx, ly, rhop)
p = scontact.ParticleProblem(2,True,True)
p.set_friction_coefficient(0)
p.read_vtk(outputdir,0)
p.velocity()[:,0] = np.random.uniform(0,1e-3,p.velocity().shape[0])
# Initial time and iteration
t = 0
ii = 0
#
# COMPUTATION LOOP
#
forces = g*p.mass()
while t < tEnd :
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)**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