Commit b309e23e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

particle have mass

parent e2e53b36
This diff is collapsed.
......@@ -10,7 +10,7 @@ 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);
double particleProblemMaxDt(const ParticleProblem *p, double alert);
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r);
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);
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag);
#if DIMENSION==3
......@@ -18,6 +18,7 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIME
#endif
double *particleProblemDisk(ParticleProblem *p);
double *particleProblemSegment(ParticleProblem *p);
double *particleProblemParticle(ParticleProblem *p);
double *particleProblemVelocity(ParticleProblem *p);
double *particleProblemPosition(ParticleProblem *p);
double *particleProblemExternalForces(ParticleProblem *p);
......
......@@ -73,7 +73,7 @@ struct ParticleProblem{};
self->_this = particleProblemNew();
return self;
}
void addParticle(const double x[DIMENSION], double r) { particleProblemAddParticle($self->_this, x, r); }
void addParticle(const double x[DIMENSION], double r, double m) { particleProblemAddParticle($self->_this, x, r, m); }
size_t addBoundaryDisk(const double x0[DIMENSION], double r, int tag) {return particleProblemAddBoundaryDisk($self->_this, x0, r, tag);}
size_t addBoundarySegment(const double x0[DIMENSION], const double x1[DIMENSION], int tag) {return particleProblemAddBoundarySegment($self->_this, x0, x1, tag);}
#if DIMENSION==3
......@@ -90,6 +90,10 @@ struct ParticleProblem{};
PyObject *triangleTag() {return tagVector2PyArray(particleProblemTriangleTag($self->_this));}
#endif
PyObject *velocity() {return coordVector2PyArray(particleProblemVelocity($self->_this), particleProblemNParticle($self->_this), DIMENSION);}
PyObject *particles() {
double *d = particleProblemParticle($self->_this);
return coordVector2PyArray(d, vectorSize(d) / 2, 2);
}
PyObject *disks() {
double *d = particleProblemDisk($self->_this);
return coordVector2PyArray(d, vectorSize(d) / (2 * DIMENSION + 1), 2 * DIMENSION + 1);
......
......@@ -136,8 +136,8 @@ class ObjectCollection {
if (type == "P") {
readCoord(hf, coord, 1);
readCoord(hf, vel, 1);
double r;
hf >> r;
double r, m;
hf >> r >> m;
radius.push_back(r);
fixed.push_back(false);
}
......
#!/usr/bin/env python
import scontact
import scontact2
import numpy as np
import mesh
import os
......@@ -11,7 +11,7 @@ rin = 0.0064
r = 397e-6/2
compacity3d = 0.18
def genInitialPosition(p, N, r, rout, rin) :
def genInitialPosition(p, N, r, rout, rin, rhop) :
x = np.arange(rout, -rout, 2 * -r)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
......@@ -19,7 +19,7 @@ def genInitialPosition(p, N, r, rout, rin) :
x = x[keep][:N]
y = y[keep][:N]
for i in range(N) :
p.addParticle((x[i], y[i]), r);
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
addv = set()
def loadMeshBoundary(p, filename) :
......@@ -32,21 +32,21 @@ def loadMeshBoundary(p, filename) :
if v[3] in addv :
continue
addv.add(v[3])
p.addBoundaryDisk((v[0], v[1]), 0.000)
p.addBoundaryDisk((v[0], v[1]), 0.000, 0)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0], v0[1]), (v1[0], v1[1]))
p.addBoundarySegment((v0[0], v0[1]), (v1[0], v1[1]), 0)
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
surfacedense = N * 2 * np.sqrt(3) * r**2
surfacebulk = 0.34 * np.pi * (rout**2 - rin**2)
p = scontact.ParticleProblem2()
p = scontact2.ParticleProblem()
if meshBoundary :
loadMeshBoundary(p, "disk.msh")
else :
p.addBoundaryDisk((0, 0), rin)
p.addBoundaryDisk((0, 0), -rout)
genInitialPosition(p, N, r, rout, rin)
p.addBoundaryDisk((0, 0), rin, 0)
p.addBoundaryDisk((0, 0), -rout, 0)
genInitialPosition(p, N, r, rout, rin, rhop=1.18e3)
scale = (surfacebulk/surfacedense)**0.5
......@@ -63,8 +63,9 @@ os.makedirs(outputdir, exist_ok=True)
while t < tEnd :
p.externalForces()[:, :] = - 3 * p.velocity()
p.externalForces()[:, 1] += -g
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = min(0.001, p.maxdt(r)/2)
p.iterate(r, dt, tol=0.0005)
p.iterate(r, dt, tol=1e-5)
t += dt
i += 1
print("%.2g %.2g" % (t, tEnd))
......
import scontact
import scontact2
import numpy as np
import os
from dgpy import *
......@@ -28,23 +28,21 @@ class part2mesh :
vertexMomentum = self.particleOnMesh.particle2mesh(volume * velocity)/self.particleMesh.vertexVolume()
particleMesh2DGDof(self.particleMesh, vertexMomentum, self.momentum)
particles = scontact.ParticleProblem2()
particles = scontact2.ParticleProblem()
particles.load("deposited")
particles.velocity()[:, :] = 0.
N = particles.position().shape[0]
model = GModel()
model.load("disk.msh")
groups = dgGroupCollection(model, 2, -2)
r = 397e-6/2
#r = 0.000132333
rhop = 1.18e3
r = np.mean(particles.particles()[:, [0]])
volume = particles.particles()[:, [0]]**2 * np.pi
rhop = np.mean(particles.particles()[:, [1]] / volume)
rho = 1.253e3
mu = 0.588
volume = np.ndarray((N, 1), 'd')
volume[:, :] = np.pi * r * r
temporal = True
ns = False
ns = True
dt = 0.5/2
endt = 800
substeps=2
......@@ -99,7 +97,7 @@ for i in range (100000) :
solve.solve(sol)
gradp = t.particleOnMesh.function2particle(groups, sol.getFunctionGradient())[:, 6:8]
for j in range(substeps):
particles.externalForces()[:N,:] = -gradp/rhop
particles.externalForces()[:,:] = -gradp * volume
particles.iterate(r, dt/substeps)
t.particlesToMesh(particles.position(), volume, particles.velocity())
gradp = t.particleOnMesh.function2particle(groups, sol.getFunctionGradient())[:, 6:8]
......
......@@ -3,14 +3,14 @@ import scontact3
import numpy as np
import os
def genInitialPosition(p, N, r, rout, rin) :
def genInitialPosition(p, N, r, rout, rin, rhop) :
n = 0
for y in np.arange(rout, -rout, 2 * -r) :
for x in np.arange(rout, -rout, 2 * -r) :
for z in np.arange(-r, -rin+r, 2* -r) :
R2 = x**2 + y**2
if R2 < (rout - r)**2 and R2 > (rin + r) **2 :
p.addParticle((x, y, z), r);
p.addParticle((x, y, z), r, 4./3*r**3 * np.pi * rhop);
n +=1
if n == N:
return
......@@ -24,22 +24,23 @@ def loadmesh(p, filename):
for ent in m.entities :
if ent.dimension != 2 :
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], v[1], v[2]), 0.000)
p.addBoundaryDisk((v[0], v[1], v[2]), 0.000, tag)
for i in range(3) :
v0 = el.vertices[i]
v1 = el.vertices[(i + 1) % 3]
if set((v0[3], v1[3])) in adde :
continue
p.addBoundarySegment((v0[0], v0[1], v0[2]), (v1[0], v1[1], v1[2]))
p.addBoundarySegment((v0[0], v0[1], v0[2]), (v1[0], v1[1], v1[2]), tag)
v0 = el.vertices[0]
v1 = el.vertices[1]
v2 = el.vertices[2]
p.addBoundaryTriangle((v0[0], v0[1], v0[2]), (v1[0], v1[1], v1[2]), (v2[0], v2[1], v2[2]))
p.addBoundaryTriangle((v0[0], v0[1], v0[2]), (v1[0], v1[1], v1[2]), (v2[0], v2[1], v2[2]), tag)
outputdir ="output_depot"
if not os.path.isdir(outputdir) :
......@@ -49,7 +50,7 @@ rout = 0.0254
rin = 0.0064
p = scontact3.ParticleProblem()
loadmesh(p, "cyl.msh")
genInitialPosition(p, N = 8000, r = r, rout= rout, rin = rin)
genInitialPosition(p, N = 8000, r = r, rout= rout, rin = rin, rhop=1.18e3)
#p.useJacobi(True, relax = 0.8)
......@@ -65,6 +66,7 @@ p.write("%s/part-%05i" % (outputdir, 0))
while t < tEnd :
p.externalForces()[:, :] = - p.velocity()
p.externalForces()[:, 1] += -g
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = min(maxdt, p.maxdt(r) * 0.7)
mindt = min(dt, mindt)
print("%i T %.2g/%.2g, dt %.2g mindt %.2g" % (iter, t, tEnd, dt, mindt))
......
......@@ -34,12 +34,12 @@ N = particles.position().shape[0]
model = GModel()
model.load("cyl.msh")
groups = dgGroupCollection(model, 3, -2)
r = 397e-6
rhop = 1.18e3
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
volume = np.ndarray((N, 1), 'd')
volume[:, :] = 4./3 * np.pi * r * r * r
temporal = False
ns = False
......@@ -108,7 +108,7 @@ for i in range (100000) :
else:
solve.solve(sol)
gradp = t.particleOnMesh.function2particle(groups, sol.getFunctionGradient())[:, 9:12]
particles.externalForces()[:N,:] = -gradp/rhop
particles.externalForces()[:,:] = -gradp * volume
particles.iterate(r, dt/substeps)
if i % 10 == 0 :
print("output %i" % i)
......
......@@ -28,22 +28,23 @@ def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
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, N, r, L, H) :
x = np.arange(L - r, -L + r, 2.1 * -r)
y = np.arange(-H+r, H, 2.1 * r);
def genInitialPosition(p, N, rmin, rmax, L, H, rhop) :
x = np.arange(L - rmax, -L + rmax, 2 * -rmax)
y = -np.arange(-H+rmin, H, 2* rmax);
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
R = - np.random.random([N]) *0.2*r + r
print(x.shape)
R = rmin + (rmax - rmin ) * np.random.random([N])
x = x.flat[:N]
y = y.flat[:N]
for i in range(N) :
p.addParticle((x[i], y[i]), R[i])
p.addParticle((x[i], y[i]), R[i], R[i]**2 * np.pi * rhop)
return p
r = 0.015
r = 0.013
dr = 0.003
meshBoundary = True
N = 8000
genInitialPosition(p, N = N, r = r, L = 2, H = 1.5)
genInitialPosition(p, N = N, rmin = r -dr, rmax = r + dr , L = 2, H = 1.5, rhop=2)
p.useJacobi(False)
loadMeshBoundary(p, "mesh.msh", ["Box"])
......@@ -59,10 +60,9 @@ p.velocity()[:, 1] = 0;
while t < tEnd :
p.externalForces()[:, :] = - p.velocity()
p.externalForces()[:, 1] += -g
if (t> 2 and t < 8) :
p.externalForces()[:, :] += (np.random.random([N, 2]) - 0.5)
dt = min(0.001, p.maxdt(r)/2)
p.iterate(r, dt, 1e-5)
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = min(0.001, p.maxdt(r +dr)/2)
p.iterate(r + dr, dt, 1e-5)
t += dt
iter += 1
if iter %50 == 0 :
......
......@@ -32,15 +32,15 @@ def writep(p, odir, i) :
tmpoutputname = "%s/part-%05i_tmp" % (odir, i)
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
p.writeVtu("%s/part-%05i.vtu" % (odir, i))
tEnd = 18
tCycle = 3
r = 0.015
g = [0, -1]
r = np.mean(p.particles()[:, 0])
volume = np.pi * p.particles()[:, [0]]**2
rhop = np.mean(p.particles()[:, 1] / volume[:, 0])
dt = 0.1 * r * abs(tCycle)
outf = 10
rhop = 1
g = [0, -0.5]
mv = meshVelocity()
hydro = runf.hydro("mesh.msh", mv.cb, dt, r, rhop=rhop, g=g)
mv.V = 0.5/dt
......@@ -52,9 +52,6 @@ outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
writep(p, outputdir, 0)
N = p.position().shape[0]
volume = np.ndarray((N, 1), 'd')
volume[:, :] = np.pi * r * r
pileseg = np.where(p.segmentTag() == 2)[0]
piledisk = np.where(p.diskTag() == 2)[0]
......@@ -66,8 +63,7 @@ while t < tEnd :
hydro.moveMesh()
hydro.solve(volume, p.position(), p.velocity(), p.externalForces())
mv.update(dt)
p.externalForces()[:, :] += rhop * np.array([g]) * volume
p.externalForces()[:, :] /= rhop * volume
p.externalForces()[:, :] += p.particles()[:, [1]] * np.array([g])
p.iterate(r, dt, 1e-5)
p.disks()[piledisk, 1] += mv.V * dt
p.segments()[pileseg, 1] += mv.V * dt
......
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