Commit 4edbd33d authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

imex

parent fca8368d
......@@ -25,8 +25,8 @@ class Brinkman :
_p_ = self.D
_U_ = range(self.D)
val = np.zeros_like(sol)
#self.g[None, None, :] * self.ca * self.rhof
val[..., _U_] = \
self.g[None, None, :] * self.ca * self.rhof\
- solgrad[..., _p_, :] \
+ self.mu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
val[..., _p_] = - solgrad[..., _u_, 0] - solgrad[..., _v_, 1] - (0 if self.D == 2 else solgrad[..., _w_, 2])
......@@ -36,12 +36,10 @@ class Brinkman :
def psiTerm1(self, sol, solgrad):
_p_ = self.D
_U_ = range(self.D)
val = np.ndarray(sol.shape + (self.D,))
val[..., _U_, :] = -solgrad[..., _U_, :] * self.mu
val[..., _p_, :] = 0#self.vs# + sol[..., _U_]
if self.stab != 0 :
pforces = solgrad[..., _p_, :]#self.particleForces(sol, solgrad)[..., _U_]
val[..., _p_, :] -= pforces * self.stab * (1-self.ca)
if self.ns :
val[..., _u_, :] += - sol[..., [_u_]] * sol[..., _U_] / self.ca * self.rhof
val[..., _v_, :] += - sol[..., [_v_]] * sol[..., _U_] / self.ca * self.rhof
......@@ -50,7 +48,7 @@ class Brinkman :
return val
# forces by the particle ON the fluid
def particleForces(self, sol, solgrad) :
def particleForces(self, sol, solgrad, fonpart = None) :
_p_ = self.D
_U_ = range(self.D)
val = np.zeros_like(sol)
......@@ -64,16 +62,22 @@ class Brinkman :
Cd_u = (0.63 * normU + 4.8/(Re_O_u**0.5))**2
f = self.ca**(-1.8)
F = solgrad[..., _p_, :]
G = f * Cd_u * self.rhof * U / 2
val[..., _U_] = F + G
GoU = f * Cd_u * self.rhof / 2
#mass = self.pf.volume * 3/2 * self.rhop
aext = self.g * (self.rhop - self.rhof) /self.rhop
UImp = U + aext * self.stab
if fonpart is not None:
fonpart[...] = (-F - GoU*UImp) / (1 + self.stab/self.pf.mass * GoU) + aext * self.pf.mass
val[..., _U_] = (F + GoU* UImp) / (1 + self.stab/self.pf.mass * GoU)
val[..., _p_] = 0
return val
def __init__(self, mu, pf, g, ns = False, rho = 1, stab = 0, D = 2) :
def __init__(self, mu, pf, g, rho, rhop, ns = False, stab = 0, D = 2) :
self.D = D
self.ns = ns
self.rhof = rho
self.rhop = rhop
self.stab = stab
self.mu = mu
self.pf = pf
......@@ -87,8 +91,7 @@ class Brinkman :
self.vs = self.pf.velocity
sol = solution.at_points(self.pf.eid, self.pf.uvw)
gradsol = solution.gradient_at_points(self.pf.eid, self.pf.uvw, meshJ)
pTerm = self.particleForces(sol, gradsol) * self.pf.volume
self.pf.forces[:, :] = -(pTerm[:, _U_])
pTerm = self.particleForces(sol, gradsol, self.pf.forces) * self.pf.volume
for i in range(dofm.n_fields()) :
np.add.at(res, dofm._fields[i][0][:, self.pf.eid], (dofm.get_field_basis(i).f(self.pf.uvw)*pTerm[:, [i]]).T)
d00 = derivate(lambda s : self.particleForces(s, gradsol), sol, 1) * self.pf.volume[..., None]
......
from .Brinkman import Brinkman
from .Steady import Steady
from .dofManager import DofManager
from .meshJacobian import MeshJacobian
from .legendre import *
from .particleFluid import ParticleFluid
from . import dofManager
import scontact2
import scontact3
import numpy as np
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.add_field(jac.dofManager.get_field_basis(0))
dofmanager2.complete()
self.momentum = dofmanager2.new_vector()
def set_particles(self, position, volume, velocity, forces, mass) :
self.position = position
self.volume = volume
self.velocity = velocity
self.forces = forces
self.mass = mass
els = self.jac.x
if self.jac.dim == 2 :
self.eid, self.uvw = scontact2.particleToMesh(position, els)
else :
self.eid, self.uvw = scontact3.particleToMesh(position, els)
self.psi = self.jac.dofManager.get_field_basis(0).f(self.uvw)
dof = self.jac.dofManager
vertexVolume = dof.new_vector()
if self.jac.dim == 2 :
np.add.at(vertexVolume, dof._fields[0][0], np.outer([1./3, 1./3, 1./3], self.jac.J/2))
else :
np.add.at(vertexVolume, dof._fields[0][0], np.outer([1./4, 1./4, 1./4, 1./4], self.jac.J/4))
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
self.momentum[:] = 0
for i in range(self.jac.dim) :
np.add.at(self.momentum, self.momentum.dof_manager()._fields[i][0][:, self.eid] , (self.psi * velocity[:, [i]] * volume/vollocal.T).T)
......@@ -656,13 +656,18 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) {
return alert / q;
}
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
_particleProblemGenContacts(p, alert);
_particleProblemSolveContacts(p, dt, tol, maxit);
}
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
for (size_t j = 0; j < vectorSize(p->particles); ++j)
for (size_t i = 0; i < DIMENSION; ++i)
p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
_particleProblemGenContacts(p, alert);
_particleProblemSolveContacts(p, dt, tol, maxit);
particleProblemSolve(p, alert, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i)
p->position[i] += p->velocity[i] * dt;
}
......
......@@ -9,6 +9,7 @@ void particleProblemDelete(ParticleProblem *p);
void particleProblemLoad(ParticleProblem *p, const char *filename);
void particleProblemWrite(const ParticleProblem *p, const char *filename);
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
double particleProblemMaxDt(const ParticleProblem *p, double alert);
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, int tag);
......
......@@ -102,6 +102,7 @@ struct ParticleProblem{};
#endif
void write(const char *filename) { particleProblemWrite($self->_this, filename); }
void iterate(double alert, double dt, double tol = 0.005, int maxit = -1) { particleProblemIterate($self->_this, alert, dt, tol, maxit); }
void solve(double alert, double dt, double tol = 0.005, int maxit = -1) { particleProblemSolve($self->_this, alert, dt, tol, maxit); }
double maxdt(double alert) { return particleProblemMaxDt($self->_this, alert); }
void load(const char *filename) { particleProblemLoad($self->_this, filename); }
......
......@@ -15,9 +15,10 @@ class Hydro :
mu = mu,
pf = self._pf,
ns = False,
g = np.array([0, 0]),
rho = rho
)
g = self.g, #np.array([0, 0]),
rho = rho,
rhop = rhop,
stab = dt)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
d.add_field(pyfefp.TriangleP1Mini)
......@@ -30,16 +31,21 @@ class Hydro :
d.complete()
self.rhop = rhop
self._dof = d
self.dt = dt
self.law = law
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac, ilu=False)
def set_particles(self, pos, vol, vel, forces) :
self._pf.set_particles(pos, vol, vel, forces)
self._solverEx = pyfefp.Steady(law, d, self._jac, ilu=False)
def solve(self) :
tic = time.clock()
v0 = self._pf.velocity.copy()
self.law.stab = self.dt *0.5
self._sol[:] = self._solver.solve()
self._pf.forces[:] += self._pf.volume * (self.rhop - self.rho) * self.g
self._pf.velocity += self._pf.forces * self.dt*0.5 / self._pf.mass
self.law.stab = 0
self._sol[:] = self._solverEx.solve()
self._pf.velocity[:] = v0
print("cpu : ", time.clock() - tic)
def write_solution(self, odir, oiter, time, filetype="msh") :
......
from pyfefp.mesh import Mesh
from hydro import Hydro
import scontact2
import numpy as np
import shutil
import time
import os
import sys
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
addv = set()
m = Mesh(filename)
pnum = [m.getPhysicalNumber(1, i) for i in tags]
for ent in m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
tag = 0 if not ent.physicals else ent.physicals[0]
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in addv :
continue
addv.add(v[3])
p.addBoundaryDisk((v[0] + shift[0], v[1] + shift[1]), 0.000, tag)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), tag)
def genInitialPosition(p, r, lx, ly, rhop) :
x = np.arange(lx-r, -lx, 2 * -r)
y = np.arange(ly-r, -ly, 2 * -r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
for i in range(len(x)) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
dt = 0.01
tEnd = 1000
outputdir = "output/fine_water"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
rhop=1100
genInitialPosition(p, 0.01, 2, 1, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 4
g = [0, -9.81]
r = np.mean(p.particles()[:, 0])
volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
#rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
print("RHOP = %g" % rhop)
outf = 5
m = Mesh("mesh.msh")
mu = 1e-3
rho = 1000
hydro = Hydro(m, dt, r, rhop=rhop, rho = rho, g=g, mu = mu)
t = 0
ii = 0
outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
p.r = r
p.dt = dt
writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
while t < tEnd :
hydro.solve()
p.iterate(r, dt, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
writep(p, outputdir, ioutput)
hydro.write_solution(outputdir, ioutput, t, filetype="msh")
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -32,40 +32,49 @@ def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), tag)
def genInitialPosition(p, r, lx, ly, rhop) :
x = np.arange(lx-r, -lx, 2 * -r)
y = np.arange(ly-r, -ly, 2 * -r)
x = np.arange(lx-r, -lx + r*0.999, 2 * -r)
y = np.arange(ly-r, -ly + r*0.999, 2 * -r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
for i in range(len(x)) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
dt = 0.005 #fine : 0.01
tEnd = 1500
dt = 0.0025
tEnd = 15000
outputdir = "output/fine3"
outputdir = "output/fine5-imex-dt-00025-outf5"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
rhop=1100
genInitialPosition(p, 0.01, 2, 1, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 4
d = 0.01
mu = 0.001
rho = 1000
ii = 0
if ii == 0 :
genInitialPosition(p, d, 2, 1, rhop)
loadMeshBoundary(p, "mesh.msh", ["Top", "Box"])
p.position()[:, 1] += 4
p.velocity()[:, :] = 0;
else :
p.load(outputdir+"/part-%05d"%ii)
mass = p.particles()[:, [1]]
g = [0, -9.81]
r = np.mean(p.particles()[:, 0])
volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
print("r = %g, m = %g\n" % (r, p.particles()[0, 1]))
print("dtmax = %g Re_p/u = %g" % (p.particles()[0, 1] / (150 * mu * d), d * rho * 0.3/mu))
#rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
print("RHOP = %g" % rhop)
outf = 20
outf = 5
m = Mesh("mesh.msh")
mu = 10
rho = 1000
hydro = Hydro(m, dt, r, rhop=rhop, rho = rho, g=g, mu = mu)
t = 0
ii = 0
outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
p.r = r
p.dt = dt
......@@ -73,14 +82,20 @@ writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
while t < tEnd :
hydro.solve()
p.iterate(r, dt, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
nsub = 20
vn = p.velocity() + p.externalForces() * dt / mass
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 2)/r)))
print("NSUB", nsub, "VMAX * dt", vmax * dt, "r", r)
for i in range(nsub) :
p.iterate(r, dt/nsub, 1e-5)
hydro._pf.set_particles(p.position(), volume, p.velocity(), p.externalForces(), mass)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
This diff is collapsed.
......@@ -4,19 +4,22 @@ import pyfefp
class Hydro :
def __init__(self, mesh, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
def __init__(self, mesh, dt, r, rhop, rho, mu, g) :
self._jac = pyfefp.MeshJacobian(mesh, 2)
self._pf = pyfefp.ParticleFluid(self._jac)
self._dt = dt
self.g = np.array(g)
self.rho = rho
self.rhop = rhop
#todo add vmesh in law
law = pyfefp.Brinkman(
self.law = pyfefp.Brinkman(
mu = mu,
pf = self._pf,
ns = False,
g = np.array([0, 0]),
rho = rho
g = self.g,
rho = rho,
rhop = rhop,
stab = 0
)
d = pyfefp.DofManager(mesh)
d.add_field(pyfefp.TriangleP1Mini)
......@@ -28,18 +31,25 @@ class Hydro :
d.add_boundary_condition("Top", 1, 0.)
d.add_boundary_condition("Top", 2, 0.)
d.complete()
self.rhop = rhop
self._dof = d
self._sol = d.new_vector()
self._solver = pyfefp.Steady(law, d, self._jac)
self._solver = pyfefp.Steady(self.law, d, self._jac)
self._solverex = pyfefp.Steady(self.law, d, self._jac)
self.dt = dt
def set_particles(self, pos, vol, vel, forces) :
self._pf.set_particles(pos, vol, vel, forces)
def set_particles(self, pos, vol, vel, forces, mass) :
self._pf.set_particles(pos, vol, vel, forces, mass)
def solve(self) :
tic = time.clock()
self._sol[:] = self._solver.solve()
self._pf.forces[:] += self._pf.volume * (self.rhop - self.rho) * self.g
v0 = self._pf.velocity.copy()
#self.law.stab = self.dt
#self._solver.solve()
#self._pf.velocity[:] += self._pf.forces * self.dt/self._pf.mass
self.law.stab = 0
self._sol[:] = self._solverex.solve()
self._pf.velocity[:] = v0
#self._pf.forces[:] += self._pf.volume * self.rhop * self.g
print("cpu : ", time.clock() - tic)
def write_solution(self, odir, oiter, time, filetype="msh") :
......
......@@ -44,7 +44,7 @@ def genInitialPosition(p, r, rout, rhop) :
dt = 0.005
tEnd = 10
outputdir = "output_run"
outputdir = "output_run_%g" % dt
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
......@@ -74,14 +74,19 @@ writep(p, outputdir, 0)
pileseg = np.where(p.segmentTag() == 3)[0]
piledisk = np.where(p.diskTag() == 3)[0]
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces(), p.particles()[:, [1]])
hydro.write_solution(outputdir, 0, 0, filetype="msh")
tic = time.clock()
v0 = 0
while t < tEnd :
hydro.solve()
p.iterate(r, dt, 1e-5)
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces())
p.velocity()[:] += p.externalForces() * dt/ p.particles()[:, [1]]
nsub = 1
for i in range(nsub) :
p.solve(dt/nsub, 1e-5)
p.position()[:] += p.velocity() * dt/nsub
hydro.set_particles(p.position(), volume, p.velocity(), p.externalForces(), p.particles()[:, [1]])
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
This diff is collapsed.
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