Commit 28c86d8e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

pyfefp as a library

parent cf7f249c
import numpy as np
import quadrature
from . import quadrature
D = 2
_u_ = 0
......@@ -21,7 +21,7 @@ def derivate(f, x, nd = 1, eps = 1e-8) :
def dot(a, b) :
return a[..., 0] * b[..., 0] + a[..., 1] * b[..., 1] + (0 if D == 2 else a[..., 2] * b[..., 2])
class brinkman :
class Brinkman :
def psiTerm0(self, sol, solgrad):
val = np.zeros_like(sol)
......@@ -66,7 +66,6 @@ class brinkman :
self.rhof = 1
self.mu = mu
self.pf = pf
self.strongBC = {}
def particleInteractionTerm(self, meshJ, solution, res, jac) :
dofm = solution.dof_manager()
......
......@@ -3,7 +3,7 @@ import time
import scipy.sparse
import scipy.sparse.linalg
class steady :
class Steady :
def __init__(self, law, dof, jac) :
self.law = law
self.meshJacobian = jac
......@@ -43,7 +43,7 @@ class steady :
if newtonIter > 10:
break
tic = time.clock()
matrix = jac
matrix = jac.tocsc()
if not self.preconditioner :
P = scipy.sparse.linalg.spilu(matrix)
self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
......
import mesh
from . import mesh
import numpy as np
import scipy.sparse
class dofManager :
class DofManager :
def getN(self) :
return self._ndof
......@@ -28,7 +28,7 @@ class dofManager :
def n_fields(self) :
return len(self._fields)
def addField(self, basis) :
def add_field(self, basis) :
field = []
elist = []
ifield = len(self._fields)
......@@ -54,7 +54,7 @@ class dofManager :
self._fields.append((np.transpose(np.array(field, np.int)), elist, basis))
def setStrongBC(self, tag, field, val) :
def set_boundary_condition(self, tag, field, val) :
vdof = self._dof_by_dimension[0]
edof = self._dof_by_dimension[self._fields[field][2].dimension()]
for en in self._mesh.entities :
......@@ -69,7 +69,7 @@ class dofManager :
for dof in edof.get((field, e.tag), ()) :
self._fixedval[self._fixed[dof]] = val
def addStrongBC(self, tag, field, val) :
def add_boundary_condition(self, tag, field, val) :
vdof = self._dof_by_dimension[0]
edof = self._dof_by_dimension[self._fields[field][2].dimension()]
for en in self._mesh.entities :
......
import gmshType
import time
class mesh() :
class Mesh() :
class element():
def __init__(self, etype, vertices, partition, tag):
self.etype = etype
......@@ -13,11 +13,11 @@ class mesh() :
v = self.vertices
order = self.etype.order
if self.etype.baseType == gmshType.TYPE_TRI :
return [mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 3]] + v[3 + i * (order-1):3 + (i+1) * (order-1)], self.partition, -1) for i in range(3)]
return [Mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 3]] + v[3 + i * (order-1):3 + (i+1) * (order-1)], self.partition, -1) for i in range(3)]
elif self.etype.baseType == gmshType.TYPE_QUA :
return [mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 4]] + v[4 + i * (order-1):4 + (i+1) * (order-1)], self.partition, -1) for i in range(4)]
return [Mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 4]] + v[4 + i * (order-1):4 + (i+1) * (order-1)], self.partition, -1) for i in range(4)]
elif self.etype.baseType == gmshType.TYPE_LIN :
return [mesh.element(gmshType.MSH_PNT, [v[i]], self.partition, -1) for i in range(2)]
return [Mesh.element(gmshType.MSH_PNT, [v[i]], self.partition, -1) for i in range(2)]
class entity():
def __init__(self, tag, dim, physicals) :
......@@ -108,7 +108,7 @@ class mesh() :
for i in range(n) :
l = fin.readline().split()
j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
self.entities.append(mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
self.entities.append(Mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
entitymap[(dim, j)] = self.entities[-1]
fin.readline()
elif w[0] == "$Nodes" :
......@@ -136,10 +136,10 @@ class mesh() :
edim = gmshType.Type[t].baseType.dimension
entity = entitymap.get((edim, e), None)
if not self.useFormat3 and not entity:
entity = mesh.entity(e, edim, [p])
entity = Mesh.entity(e, edim, [p])
self.entities.append(entity)
entitymap[(edim, e)] = entity
entity.elements.append(mesh.element(gmshType.Type[t], vertices, partition, j))
entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
self.maxeid = max(self.maxeid, j)
l = fin.readline()
......@@ -162,11 +162,11 @@ class mesh() :
def newElement(self, t, partition, entity, vertices) :
self.maxeid += 1
entity.elements.append(mesh.element(t, vertices, partition, self.maxeid))
entity.elements.append(Mesh.element(t, vertices, partition, self.maxeid))
return entity.elements[-1]
def newEntity(self, entityDim, physicalTag) :
self.entities.append(mesh.entity(len(self.entities) + 1, entityDim, [self.getPhysicalNumber(entityDim, physicalTag)]))
self.entities.append(Mesh.entity(len(self.entities) + 1, entityDim, [self.getPhysicalNumber(entityDim, physicalTag)]))
return self.entities[-1]
def genFaces(self, dim, physicalTag) :
......
import dofManager
import integrationMatrices
from . import dofManager
from . import legendre
import numpy as np
import scipy.sparse
import legendre
D = 2
class meshJacobian :
class MeshJacobian :
def update(self, x) :
self.x = x
dim = self.dim
......@@ -75,9 +74,9 @@ class meshJacobian :
def __init__(self, mesh, dim) :
d = dofManager.dofManager(mesh)
d = dofManager.DofManager(mesh)
self.dofManager = d
d.addField(legendre.TriangleP1)
d.add_field(legendre.TriangleP1)
d.complete()
self.dim = dim
self.mesh = mesh
......
#!/usr/bin/env python
import scontact2
import numpy as np
import mesh
......@@ -40,8 +39,8 @@ def genInitialPosition(p, N, rmin, rmax, L, H, rhop) :
p.addParticle((x[i], y[i]), R[i], R[i]**2 * np.pi * rhop)
return p
r = 0.013*2
dr = 0#0.003*2
r = 0.025
dr = 0.002
N = 2000
genInitialPosition(p, N = N, rmin = r -dr, rmax = r + dr , L = 2, H = 1.5, rhop=2)
p.useJacobi(False)
......
import Brinkman
import Steady
import numpy as np
import dofManager
import meshJacobian
import scontact2
import time
import legendre
import scipy.sparse
import pyfefp
class particleFluid :
def __init__(self, jac) :
self.jac = jac
self.porosity = jac.dofManager.new_vector()
self.compacity = jac.dofManager.new_vector()
self.mesh = jac.dofManager._mesh
dofmanager2 = dofManager.dofManager(self.mesh)
for i in range(self.jac.dim) :
dofmanager2.addField(jac.dofManager.get_field_basis(0))
dofmanager2.complete()
self.momentum = dofmanager2.new_vector()
class Hydro :
def setParticles(self, position, volume, velocity, forces) :
self.position = position
self.volume = volume
self.velocity = velocity
self.forces = forces
els = self.jac.x
self.eid, self.uvw = scontact2.particleToMesh(position, els)
xi, eta = self.uvw[:, [0]], self.uvw[:, [1]]
self.psi = np.hstack([1 - xi - eta, xi, eta])
dof = self.jac.dofManager
vertexVolume = dof.new_vector()
np.add.at(vertexVolume, dof._fields[0][0], np.outer([1./3, 1./3, 1./3], self.jac.J/2))
vollocal = vertexVolume[dof._fields[0][0][:,self.eid]]
self.compacity[:] = 0
np.add.at(self.compacity, dof._fields[0][0][:, self.eid] , (self.psi * volume/vollocal.T).T)
self.porosity[:] = 1 - self.compacity
#todo 3D
self.momentum[:] = 0
np.add.at(self.momentum, self.momentum.dof_manager()._fields[0][0][:, self.eid] , (self.psi * velocity[:, [0]] * volume/vollocal.T).T)
np.add.at(self.momentum, self.momentum.dof_manager()._fields[0][0][:, self.eid] , (self.psi * velocity[:, [1]] * volume/vollocal.T).T)
class hydro :
def __init__(self, mesh, meshVelocityCB, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
self.mesh = mesh
self.jac = meshJacobian.meshJacobian(mesh, 2)
self.pf = particleFluid(self.jac)
self.meshVelocity = meshVelocityCB
self.rho = rho
self.dt = dt
self.g = np.array([g])
def __init__(self, mesh, mesh_velocity_cb, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
self._jac = pyfefp.MeshJacobian(mesh, 2)
self._pf = pyfefp.ParticleFluid(self._jac)
self._mesh_velocity = mesh_velocity_cb
self._dt = dt
#todo add vmesh in law
law = Brinkman.brinkman(
law = pyfefp.Brinkman(
mu = mu,
pf = self.pf,
pf = self._pf,
ns = False,
rhof = 1)
d = dofManager.dofManager(mesh)
d.addField(legendre.TriangleP1Mini)
d.addField(legendre.TriangleP1Mini)
d.addField(legendre.TriangleP1)
d.addStrongBC("Pile", 0, 0.)
d.addStrongBC("Pile", 1, 0.)
d.addStrongBC("Box", 0, 0.)
d.addStrongBC("Box", 1, 0.)
d.addStrongBC("Top", 0, 0.)
d.addStrongBC("Top", 1, 0.)
d.addStrongBC("Top", 2, 0.)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1)
d.add_boundary_condition("Pile", 0, 0.)
d.add_boundary_condition("Pile", 1, 0.)
d.add_boundary_condition("Box", 0, 0.)
d.add_boundary_condition("Box", 1, 0.)
d.add_boundary_condition("Top", 0, 0.)
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.dof = d
self.sol = np.zeros((d.getN()))
self.solver = Steady.steady(law, d, self.jac)
self.law = law
self._dof = d
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac)
def set_particles(self, pos, vol, vel, forces) :
self._pf.set_particles(pos, vol, vel, forces)
def setVelocity(self, v) :
self.law.strongBC["Pile"] = ((0, 0.), (1, v))
self.solver.dof.setStrongBC("Pile", 1, v)
def set_velocity(self, v) :
self._dof.set_boundary_condition("Pile", 1, v)
def moveMesh(self) :
dx = self.meshVelocity.cb(self.jac.x.reshape([-1, 2])) * self.dt
self.jac.update(self.jac.x + dx.reshape([-1, 3, 2]))
def move_mesh(self) :
dx = self._mesh_velocity.cb(self._jac.x.reshape([-1, 2])) * self._dt
self._jac.update(self._jac.x + dx.reshape([-1, 3, 2]))
def solve(self, volume, position, velocity, forces) :
def solve(self) :
tic = time.clock()
self.sol[:] = self.solver.solve()
self._sol[:] = self._solver.solve()
print("cpu : ", time.clock() - tic)
def writeSolution(self, odir, oiter, time, filetype="msh") :
def write_solution(self, odir, oiter, time, filetype="msh") :
print(oiter, time)
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self.jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self.jac.dofManager.getField(0, self.pf.porosity))
self.jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self.dof.getField(2, self.sol))
X = np.zeros((9, self.jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self.jac.x.reshape([-1, 6]).T
self.jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
X[0::3, :] = self.dof.getField(0, self.sol)[:-1,:]
X[1::3, :] = self.dof.getField(1, self.sol)[:-1,:]
self.jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
elif filetype == "vtk" :
print("vtk output not implemented")
#exportfun = self.groups.exportFunctionVtk
#basename = odir+"/f"
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(2, self._sol))
X = np.zeros((9, self._jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
X[0::3, :] = self._dof.getField(0, self._sol)[:-1,:]
X[1::3, :] = self._dof.getField(1, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
return
import numpy as np
class integrationMatrices :
def __init__(self, integ) :
N = integ.psi().shape[1]
D = integ.dPsi().shape[0]//integ.psi().shape[0]
self.psiPsiW = np.einsum("nmq->qnm", integ.psiPsiW().reshape(N, N, -1))
self.dPsiPsiW = np.einsum("nmxq->qxnm", integ.dPsiPsiW().reshape(N, N, D, -1))
self.psiDPsiW = np.einsum("nmqx->qxnm", integ.psiDPsiW().reshape(N, N, -1, D))
self.dPsiDPsiW = np.einsum("nmxyq->qxynm", integ.dPsiDPsiW().reshape(N, N, D, D, -1))
self.dPsi = integ.dPsi().reshape(-1, D, N)
self.dPsiW = integ.dPsiW().reshape(N, -1, D)
self.psiW = integ.psiW().T #Q N
L = 2;
l = 0.1;
H = 1.5;
h = 0.5;
y = 0;
lc = 0.1;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, -H, 0, lc};
Point(3) = {L, -H, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1:4};
Point(11) = {-l, y + h, 0, lc};
Point(12) = {l, y + h, 0, lc};
Point(13) = {l, y - h, 0, lc};
Point(14) = {-l, y - h, 0, lc};
Line(11) = {11, 12};
Line(12) = {12, 13};
Line(13) = {13, 14};
Line(14) = {14, 11};
Line Loop(2) = {11:14};
Plane Surface(1) = {1, 2};
Physical Line("Box") = {1,2,3};
Physical Line("Top") = {4};
Physical Line("Pile") = {11:14};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
import mesh
from pyfefp.mesh import Mesh
from hydro import Hydro
import scontact2
import numpy as np
import hydro as Hydro
import shutil, time, os
import shutil
import time
import os
reloadit = -1
class meshVelocity :
class MeshVelocity :
def __init__(self) :
self.Y = 0
self.V = 0
......@@ -38,8 +40,7 @@ if not os.path.isdir(outputdir) :
p = scontact2.ParticleProblem()
p.load("deposited" if reloadit == -1 else "output_run/part-%05d" % reloadit)
tEnd = 18
tEnd = 1.5
tEnd = 9
tCycle = 3
g = [0, -1]
r = np.mean(p.particles()[:, 0])
......@@ -47,13 +48,13 @@ volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
dt = 0.025 * r * abs(tCycle)
outf = 10
mv = meshVelocity()
m = mesh.mesh("mesh.msh")
hydro = Hydro.hydro(m, mv, dt, r, rhop=rhop, g=g, mu = 0.1)
mv = MeshVelocity()
m = Mesh("mesh.msh")
hydro = Hydro(m, mv, dt, r, rhop=rhop, g=g, mu = 0.1)
mv.V = 0.5/dt
if reloadit != -1 :
mv.V -= (reloadit * dt) / (tCycle/2) / dt
hydro.moveMesh()
hydro.move_mesh()
mv.update(dt)
t = 0
ii = 0
......@@ -66,20 +67,20 @@ writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
hydro.pf.setParticles(p.position(), volume, p.velocity(), p.externalForces())
hydro.writeSolution(outputdir, 0, 0, filetype="msh")
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
while t < tEnd :
mv.V = (-1 if ((t % tCycle) < (tCycle / 2)) else 1) * 1/ (tCycle/2)
hydro.setVelocity(mv.V)
hydro.set_velocity(mv.V)
p.disks()[piledisk, 3] = mv.V
p.segments()[pileseg, 5] = mv.V
hydro.solve(volume, p.position(), p.velocity(), p.externalForces())
hydro.solve()
mv.update(dt)
p.iterate(r, dt, 1e-5)
hydro.moveMesh()
hydro.pf.setParticles(p.position(), volume, p.velocity(), p.externalForces())
hydro.move_mesh()
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
p.disks()[piledisk, 1] += mv.V * dt
p.segments()[pileseg, 1] += mv.V * dt
p.segments()[pileseg, 3] += mv.V * dt
......@@ -87,10 +88,8 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
writep(p, outputdir, ioutput)
hydro.writeSolution(outputdir, ioutput, t, filetype="msh")
hydro.write_solution(outputdir, ioutput, t, filetype="msh")
ii += 1
if (ii == 10) :
break
print("%i : %.2g/%.2g" % (ii, t, tEnd))
......
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