Commit ee884dd9 authored by mozul's avatar mozul

couette3d : lmgc90 use instead of scontact

parent f7f20d4a
......@@ -45,12 +45,12 @@ def loadmesh(p, filename):
outputdir ="output_depot"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
r = 397e-6
r = 397e-6 * 2
rout = 0.0254
rin = 0.0064
p = scontact3.ParticleProblem()
loadmesh(p, "cyl.msh")
genInitialPosition(p, N = 8000, r = r, rout= rout, rin = rin, rhop=1.18e3)
genInitialPosition(p, N = 1000, r = r, rout= rout, rin = rin, rhop=1.18e3)
#p.useJacobi(True, relax = 0.8)
......
import numpy as np
from pylmgc90 import chipy
class LmgcParticles(object):
"""
A class definition holding data of lmgc90
in a compatible form to use coupling with
gmsh
Attributs are :
- number of disks
- volume of diskx
- reference on diskx2rbdy2 map
The methods are :
- iterate() : do one time iteration
- position() : give back position of all diskx
- velocity() : give back position of all diskx
- externalForces() : update external forces on particles
- getMeanRadius() : return the mean radius of particles
- getVolume() : give back the volume of all diskx
- writeHeader() : write header file for plot
"""
def __init__(self,space_dim, dt, theta) :
"""
Initialize LMGC90
"""
### Set dimension in chipy for dummies ###
chipy.SetDimension(space_dim)
chipy.Initialize()
### definition des parametres du calcul ###
chipy.utilities_logMes('INIT TIME STEPPING')
chipy.TimeEvolution_SetTimeStep(dt)
chipy.Integrator_InitTheta(theta)
### lecture du modele ###
chipy.utilities_logMes('READ BEHAVIOURS')
chipy.ReadBehaviours()
chipy.utilities_logMes('READ BODIES')
chipy.ReadBodies()
chipy.utilities_logMes('LOAD BEHAVIOURS')
chipy.LoadBehaviours()
chipy.utilities_logMes('READ INI DOF')
chipy.ReadIniDof()
chipy.utilities_logMes('READ DRIVEN DOF')
chipy.ReadDrivenDof()
chipy.utilities_logMes('LOAD TACTORS')
chipy.LoadTactors()
chipy.utilities_logMes('READ INI Vloc Rloc')
chipy.ReadIniVlocRloc()
##### ecriture paranoiaque du modele ###
##utilities_logMes('WRITE BODIES')
##WriteBodies()
##utilities_logMes('WRITE BEHAVIOURS')
##WriteBehaviours()
##utilities_logMes('WRITE DRIVEN DOF')
##WriteDrivenDof()
chipy.OpenPostproFiles()
chipy.OpenDisplayFiles()
chipy.WriteDisplayFiles()
chipy.ComputeMass()
self._space_dim = space_dim
if space_dim == 2:
self._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._getCoor = chipy.DISKx_GetContactorCoor
self._getVect = chipy.RBDY2_GetBodyVector
self._setVect = chipy.RBDY2_PutBodyVector
nb_dof = 3
else :
self._nbDisk = chipy.SPHER_GetNbSPHER()
self._d2bMap = chipy.SPHER_GetPtrSPHER2RBDY3()
self._getCoor = chipy.SPHER_GetContactorCoor
self._getVect = chipy.RBDY3_GetBodyVector
self._setVect = chipy.RBDY3_PutBodyVector
nb_dof = 6
self._position = np.empty([self._nbDisk,nb_dof], 'd')
self._velocity = np.empty([self._nbDisk,nb_dof], 'd')
self._externalF= np.empty([self._nbDisk,nb_dof], 'd')
self._volume = np.empty([self._nbDisk,1], 'd')
if space_dim == 2:
for i in xrange(self._nbDisk):
self._volume[i] = np.pi * chipy.DISKx_GetContactorRadius(i+1)**2
else:
for i in xrange(self._nbDisk):
self._volume[i] = 4. * np.pi * chipy.SPHER_GetContactorRadius(i+1)**3 / 3.
def iterate(self,freq_detect,freq_write,freq_display,ref_radius,**solver_args):
"""
Do one step of a LMGC90 computation.
Last argument is a dictionnary holding nlgs_solver list of arguments.
"""
#
chipy.IncrementStep()
chipy.ComputeFext()
for i in xrange(self._nbDisk):
self._setVect('Fext_', int(self._d2bMap[i,0]), self._externalF[i,:])
chipy.ComputeBulk()
chipy.ComputeFreeVelocity()
chipy.SelectProxTactors(freq_detect)
chipy.RecupRloc()
chipy.ExSolver(**solver_args)
chipy.StockRloc()
chipy.ComputeDof()
if freq_write > 0:
chipy.WriteOutDof(freq_write)
chipy.WriteOutVlocRloc(freq_write)
else:
chipy.WriteLastDof()
chipy.WriteLastVlocRloc()
chipy.WritePostproFiles()
if freq_display > 0:
chipy.WriteDisplayFiles(freq_display,ref_radius)
chipy.UpdateStep()
def volume(self):
""" Return the volume of particles """
return self._volume
def position(self):
""" Get current position of contactors and return it """
for i in xrange(self._nbDisk):
self._position[i,:] = self._getCoor(i+1)
self._position[i,self._space_dim] = 0.
return self._position[:,:self._space_dim]
def velocity(self):
""" Get current velocity of body of contactor and return it
Beware : it should not work very well with clusters !
"""
for i in xrange(self._nbDisk):
self._velocity[i,:] = self._getVect('V____',int(self._d2bMap[i,0]))
self._velocity[i,self._space_dim:] = 0.
return self._velocity[:,:self._space_dim]
def externalForces(self):
""" Get an external forces array
"""
return self._externalF[:,:self._space_dim]
def __del__(self):
"""
"""
chipy.ClosePostproFiles()
chipy.CloseDisplayFiles()
chipy.Finalize()
import numpy as np
import os
from dgpy import *
from brinkman import brinkman
nb_steps = 100000
def savepartposvec(filename, viewname, position, field) :
with open(filename, "w") as output :
output.write("View \"%s\" {\n" %viewname)
for i in range(position.shape[0]) :
output.write("VP (%g, %g, 0) {%g, %g, 0};\n" % (
position[i, 0], position[i, 1], field[i, 0], field[i, 1]))
output.write("};\n");
class part2mesh :
def __init__(self, groups, tag) :
self.groups = groups
self.particleMesh = particleMesh(self.groups.getModel(), tag);
self.porosity = dgDofContainer(self.groups, 1)
self.momentum = dgDofContainer(self.groups, 3)
def particlesToMesh (self, position, volume, velocity) :
self.particleOnMesh = particleList(self.particleMesh, position)
vertexCompacity = self.particleOnMesh.particle2mesh(volume) / self.particleMesh.vertexVolume()
porosity3d = 1 - vertexCompacity
particleMesh2DGDof(self.particleMesh, porosity3d, self.porosity)
vertexMomentum = self.particleOnMesh.particle2mesh(volume * velocity)/self.particleMesh.vertexVolume()
particleMesh2DGDof(self.particleMesh, vertexMomentum, self.momentum)
rho_particle = 1.18e3
#from .geo file
rin = 0.0064;
rout = 0.0254;
fric = 0.
from pylmgc90 import pre_lmgc
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
f = open('deposited')
l = f.readline()
space_dim = int(l.split()[1])
x = []
v = []
r = []
for l in f.readlines():
if l[0] == 'P':
l = l.split()
x.append( map(float,l[1:space_dim+1]) )
v.append( map(float,l[space_dim+1:2*space_dim+1]) )
r.append( float(l[2*space_dim+1]) )
x = np.array(x)
r = np.array(r)
r_min = np.min(r)
avs = pre_lmgc.avatars()
mat = pre_lmgc.materials()
mod = pre_lmgc.models()
svs = pre_lmgc.see_tables()
tac = pre_lmgc.tact_behavs()
mater = pre_lmgc.material(name='STONE',type='RIGID',density=rho_particle)
model = pre_lmgc.model(name='rigid', type='MECAx', element='Rxx3D',dimension=space_dim)
mat.addMaterial(mater)
mod.addModel(model)
for i in xrange(r.size) :
P = pre_lmgc.rigidSphere( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
#P.imposeInitValue(component=[1,2,3],value=v[i])
avs.addAvatar(P)
up = pre_lmgc.rigidPlan(axe1=rout, axe2=rout, axe3=r_min, center=[0., 0., r_min], model=model, material=mater,color='OUTxx')
up.rotate(theta=np.pi, center=up.nodes[1].coor)
down= pre_lmgc.rigidPlan(axe1=rout, axe2=rout, axe3=r_min, center=[0., 0., -rin-r_min], model=model, material=mater,color='OUTxx')
cyl = pre_lmgc.rigidCylinder( rin, rin, [0.,0.,-rin/2.], model, mater, 'OUTxx' )
lyc = pre_lmgc.rigidCylinder( rout, rin, [0.,0.,-rin/2.], model, mater, 'OUTxx', is_Hollow=True )
up.imposeDrivenDof(component=range(1,7),dofty='vlocy')
down.imposeDrivenDof(component=range(1,7),dofty='vlocy')
lyc.imposeDrivenDof(component=range(1,7),dofty='vlocy')
cyl.imposeDrivenDof(component=range(1,6),dofty='vlocy')
cyl.imposeDrivenDof(component=6,dofty='vlocy',ct=1.)
avs.addAvatar(cyl)
avs.addAvatar(lyc)
avs.addAvatar(down)
avs.addAvatar(up)
clb_fric = pre_lmgc.tact_behav(name='iqsc0',type='IQS_CLB',fric=fric)
tac += clb_fric
IN = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="SPHER", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.max())
CYL = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="CYLND", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.min())
LYC = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="DNLYC", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.min())
OUT = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="PLANx", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.min())
svs += IN
svs += CYL
svs += LYC
svs += OUT
post = pre_lmgc.postpro_commands()
solv = pre_lmgc.postpro_command(type='SOLVER INFORMATIONS', step=1)
post.addCommand(solv)
# file writting
pre_lmgc.writeBodies(avs,chemin=datbox_path)
pre_lmgc.writeDrvDof(avs,chemin=datbox_path)
pre_lmgc.writeDofIni(avs,chemin=datbox_path)
pre_lmgc.writeModels(mod,chemin=datbox_path)
pre_lmgc.writeBulkBehav(mat,chemin=datbox_path,gravy=[0.,0.,0.])
pre_lmgc.writeTactBehav(tac,svs,chemin=datbox_path)
pre_lmgc.writeVlocRlocIni(chemin=datbox_path)
pre_lmgc.writePostpro(commands=post, parts=avs, path=datbox_path)
from lmgc_particles import *
t = 0.
dt = 0.1
endt = 800
theta = 0.5
freq_write = 10
freq_display = 10
ref_radius = 0.1e-1
freq_detect = 1
solver_params = { 'type' :'Stored_Delassus_Loops ',
'norm' : 'Quad ',
'conv' : 1e-5,
'relax' : 1.0,
'gsit1' : 20,
'gsit2' : 20000
}
#chipy.nlgs_SetWithQuickScramble()
chipy.checkDirectories()
particles = LmgcParticles(space_dim,dt,theta)
model = GModel()
model.load("cyl.msh")
groups = dgGroupCollection(model, 3, -2)
#r = np.mean(particles.particles()[:, 0])
#volume = particles.particles()[:, [0]]**3 * 4./3 * np.pi
#rhop = np.mean(particles.particles()[:, [1]] / volume)
rho = 1.253e3
mu = 0.588
temporal = False
ns = False
substeps=1
#substeps=2
outputdir ="loutput_" + ("ns" if ns else "stokes") + ("" if temporal else "_steady")
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
with open(outputdir + "/runinfo", "w") as log :
log.write("%g\n" % dt);
t = part2mesh(groups, ["Domain"])
law = brinkman(
mu = functionConstant(mu/rho),
k = functionConstant((np.mean(r)*2)**2/(150 * mu/rho)),
porosity = t.porosity.getFunction(),
particleMomentum = t.momentum.getFunction(),
stab = dt/rho_particle * 1,
ns = ns)
lsys = linearSystemPETScDouble()
lsys.setParameter("petscOptions", "-pc_type lu -ksp_atol 1e-16 -ksp_rtol 1e-8 -ksp_monitor")
dof = dgDofManager.newCG(groups, law.getNbFields(), lsys)
sol = dgDofContainer(groups, law.getNbFields())
law.setP1P1Mini(dof)
############# boundary conditions ################
law.bndOut = functionConstant([0, 0, 0, 0, 1, 1, 1, -1])
def bndInF(cmap, val, xyz):
val[:, 0] = -xyz[:, 1]
val[:, 1] = xyz[:, 0]
val[:, 2] = 0
val[:, 3] = 0
val[:, 4:7] = 1
val[:, 7] = -1
def lateralF(cmap, val, valin):
val[:, [0, 1, 3]] = valin[:, [0, 1, 3]]
val[:, 2] = -valin[:, 2]
law.bndIn = functionNumpy(8, bndInF, [function.getCoordinates()])
law.addStrongBoundaryCondition(2, "Inner", law.bndIn)
law.addStrongBoundaryCondition(2, "Outer", law.bndOut)
law.ptFix = functionConstant([0, 0, 0, 0, -1, -1, -1, 1])
law.addStrongBoundaryCondition(1, "PtFix", law.ptFix)
law.setBoundary0Flux("Inner")
law.setBoundary0Flux("Outer")
bndLatF = functionNumpy(4, lateralF, [function.getSolution()])
bndLat = law.newOutsideValueBoundary(bndLatF)
law.addBoundaryCondition(["Top", "Bottom"], bndLat)
################################################
if temporal :
solve = dgDIRK(law, dof, 1)
else :
solve = dgSteady(law, dof)
T = 0
for i in range(nb_steps) :
t.particlesToMesh(particles.position(), particles.volume(), particles.velocity())
if temporal :
solve.iterate(sol, dt, T)
else:
solve.solve(sol)
gradp = t.particleOnMesh.function2particle(groups, sol.getFunctionGradient())[:, 9:12]
particles.externalForces()[:,:] = -gradp * particles.volume()
particles.iterate(freq_detect,freq_write,freq_display,ref_radius,**solver_params)
if i % freq_display == 0 :
print("output %i" % i)
groups.exportFunctionVtk(law.pressure, "%s/couette" % (outputdir), i*dt, i, "pressure", sol)
groups.exportFunctionVtk(law.velocity, "%s/couette" % (outputdir), i*dt, i, "velocity", sol)
groups.exportFunctionVtk(t.porosity.getFunction(), "%s/couette" % (outputdir), i*dt, i, "porosity", sol)
groups.exportFunctionVtk(t.momentum.getFunction(), "%s/couette" % (outputdir), i*dt, i, "velocityB", sol)
if (i * dt > endt) :
break
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