Commit e2e53b36 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

clean hydro-quebec + tag for boundary objects

parent 159d021a
......@@ -11,23 +11,24 @@
\usepackage[round]{natbib}
\usetikzlibrary{arrows}
\newcommand{\re}{\mathit{R\kern-.15em e}}
\DeclareUnicodeCharacter{183}{{\,\cdot\,}}
\begin{document}
\bibliographystyle{plainnat}
\noindent
NB : all those papers deal with gaz-particle mixture in general and fluidized bed in particular. I'm not sure this can be applied to subsurface water.
\section*{Particle-fluid interaction forces}
While other forces can be seen in the literature \citep{zhu_discrete_2007}, the principal forces present in the gas-particle interactions are the drag force $F_d$ and the pressure gradient force $F_p$.
While other forces can be seen in the literature \citep{zhu_discrete_2007}, the principal forces present in the gas-particle interactions are the drag force $G$ and the pressure gradient force $F$.
\subsection*{Pressure gradient force}
The pressure gradient force can be easily computed by
\begin{equation}
F_p = -V_p∇p
F = -V_p∇p
\end{equation}
where $p$ is the fluid pressure field and $V_p$ the volume of the particle \citep{anderson_fluid_1967}.
\subsection*{Drag force}
There exists many parameterization of the drag force, \citet{li_gas-particle_2003} did a systematic study that indicates that they all possess similar predictive capability.
\citet{di_felice_fluid-particle_2012} uses the following simple formulation:
\begin{equation} \label{eq:fdp}
F_d = f(ε)C_d|u_p - u_f|(u_p - u_f)\frac{πd²__f}{8}
G = f(ε)C_d|u_p - u_f|(u_p - u_f)\frac{πd²__f}{8}
\end{equation}
where $C_d$ is the drag coefficient for an isolated particle and $f(ε)$ take into account the effect of the surrounding particles.
$C_d$ is a function of the Reynolds number of the particle $\re_p =d__f(u_p-u_f)εμ_f^{-1}$, one can use for example \citep{dallavalle_micromeritics:_1943} :
......@@ -46,9 +47,20 @@ dependency on the Reynolds number: for example, that proposed by \citet {di_feli
\section*{Discretization}
\begin{enumerate}
\item $ε$ is computed as the projection of the particles volume on the mesh nodes
\item $ε$ and $u_f$ are interpolated at the particles position
\item $F_d$ and $F_p$ are evaluated for each particle
\item $F_d$ and $F_p$ are projected on the mesh nodes for the fluid equation
\item[] This is exactly like a quadrature rule where the particles position would be the quadrature points and the particle volume the corresponding weights.
\item $ε$ and $u_f$ are interpolated at the particle positions
\item $F_d$ and $F_p$ are assembled on the residual.
\begin{align*}
F_i &= ∑_pF_p(P_p, (∇P)_p)φ_i(x_p)\\
\frac{∂F_i}{∂P_j}
&=∑_p\left(\frac{∂F_p}{∂P_p}\frac{∂P_p}{∂P_j}φ_i(x_p) + \frac{∂F_p}{∂(∇P)_p}·\frac{∂(∇P)_p}{∂P_j}\right_i(x_p)\\
&=∑_p\left(\frac{∂F_p}{∂P_p}φ_j(x_p) + \frac{∂F_p}{∂(∇P)_p}{·}∇φ_j(x_p)\right_i(x_p)
\end{align*}
\end{enumerate}
\section{Comparison with Darcy}
In our current implementation (Brinkman law), we have
\[
F = -ρ\frac{με²}{K₁}(u_p - u_f) + \dots\text{, with }K₁ = \frac{k₁ε³}{(1-ε)²}
\]
\bibliography{zotero.bib}
\end{document}
......@@ -19,11 +19,13 @@ struct _ParticleProblem{
Contact *contacts;
double *position, *velocity, *externalForces;
Disk *disks;
int *diskTag, *segmentTag;
Segment *segments;
Contact *contactParticleDisks, *contactParticleSegments;
#if DIMENSION == 3
Triangle *triangles;
Contact *contactParticleTriangles;
int *triangleTag;
#endif
};
......@@ -94,8 +96,9 @@ static void diskBoundingBox(const Disk *d, double *pmin, double *pmax) {
}
}
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r) {
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r, int tag) {
Disk *d = vectorPush(&p->disks);
*vectorPush(&p->diskTag) = tag;
d->r = r;
for (int i = 0; i < DIMENSION; ++i) {
d->x[i] = x[i];
......@@ -175,8 +178,9 @@ struct _Triangle {
double p[3][DIMENSION];
};
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3]) {
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag) {
Triangle *t = vectorPush(&p->triangles);
*vectorPush(&p->triangleTag) = tag;
for (int i = 0; i < DIMENSION; ++i) {
t->p[0][i] = x0[i];
t->p[1][i] = x1[i];
......@@ -518,8 +522,9 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
p->position[i] += p->velocity[i] * dt;
}
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION]) {
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag) {
Segment *s = vectorPush(&p->segments);
*vectorPush(&p->segmentTag) = tag;
for (int i = 0; i < DIMENSION; ++i) {
s->p[0][i] = x0[i];
s->p[1][i] = x1[i];
......@@ -551,19 +556,19 @@ void particleProblemWrite(const ParticleProblem *p, const char *filename) {
coordWrite(output, p->triangles[i].p[0]);
coordWrite(output, p->triangles[i].p[1]);
coordWrite(output, p->triangles[i].p[2]);
fprintf(output, "\n");
fprintf(output, " %i\n", p->triangleTag[i]);
}
#endif
for (size_t i = 0; i < vectorSize(p->segments); ++i){
fprintf(output, "S");
coordWrite(output, p->segments[i].p[0]);
coordWrite(output, p->segments[i].p[1]);
fprintf(output, "\n");
fprintf(output, " %i\n", p->segmentTag[i]);
}
for (size_t i = 0; i < vectorSize(p->disks); ++i) {
fprintf(output, "D");
coordWrite(output, p->disks[i].x);
fprintf(output, " %.16g\n", p->disks[i].r);
fprintf(output, " %.16g %i\n", p->disks[i].r, p->diskTag[i]);
}
for(size_t i = 0; i< vectorSize(p->particles); ++i) {
fprintf(output, "P");
......@@ -595,6 +600,7 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
{
FILE *input = fopen(filename, "r");
char type[128];
int tag;
while (fscanf(input, "%127s", type) == 1) {
if (!strcmp(type, "DIMENSION")) {
int dim;
......@@ -618,14 +624,15 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
double point[DIMENSION];
double r;
coordRead(input, point);
fscanf(input, "%le", &r);
particleProblemAddBoundaryDisk(p, point, r);
fscanf(input, "%le %d", &r, &tag);
particleProblemAddBoundaryDisk(p, point, r, tag);
}
else if (!strcmp(type, "S")) {
double x0[DIMENSION], x1[DIMENSION];
coordRead(input, x0);
coordRead(input, x1);
particleProblemAddBoundarySegment(p, x0, x1);
fscanf(input, "%d", &tag);
particleProblemAddBoundarySegment(p, x0, x1, tag);
}
#if DIMENSION == 3
else if (!strcmp(type, "T")) {
......@@ -633,7 +640,8 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
coordRead(input, x0);
coordRead(input, x1);
coordRead(input, x2);
particleProblemAddBoundaryTriangle(p, x0, x1, x2);
fscanf(input, "%d", &tag);
particleProblemAddBoundaryTriangle(p, x0, x1, x2, tag);
}
#endif
else {
......@@ -644,6 +652,11 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
}
int *particleProblemDiskTag(ParticleProblem *p) {return p->diskTag;}
int *particleProblemSegmentTag(ParticleProblem *p) {return p->segmentTag;};
#if DIMENSION == 3
int *particleProblemTriangleTag(ParticleProblem *p) {return p->triangleTag;};
#endif
double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
double *particleProblemDisk(ParticleProblem *p) {return (double*)&p->disks[0];}
double *particleProblemSegment(ParticleProblem *p) {return (double*)&p->segments[0];}
......@@ -653,7 +666,9 @@ void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax) {p->
void particleProblemDelete(ParticleProblem *p) {
vectorFree(p->particles);
vectorFree(p->disks);
vectorFree(p->diskTag);
vectorFree(p->segments);
vectorFree(p->segmentTag);
vectorFree(p->contactParticleDisks);
vectorFree(p->contactParticleSegments);
vectorFree(p->contacts);
......@@ -662,6 +677,7 @@ void particleProblemDelete(ParticleProblem *p) {
vectorFree(p->externalForces);
#if DIMENSION == 3
vectorFree(p->triangles);
vectorFree(p->triangleTag);
vectorFree(p->contactParticleTriangles);
#endif
free(p);
......@@ -676,12 +692,15 @@ ParticleProblem *particleProblemNew() {
p->contactParticleDisks = NULL;
p->contactParticleSegments = NULL;
p->disks = NULL;
p->diskTag = NULL;
p->segments = NULL;
p->segmentTag = NULL;
p->position = NULL;
p->velocity = NULL;
p->externalForces = NULL;
#if DIMENSION == 3
p->triangles = NULL;
p->triangleTag = NULL;
p->contactParticleTriangles = NULL;
#endif
return p;
......
......@@ -11,10 +11,10 @@ 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);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r);
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION]);
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
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION]);
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], int tag);
#endif
double *particleProblemDisk(ParticleProblem *p);
double *particleProblemSegment(ParticleProblem *p);
......@@ -25,4 +25,10 @@ unsigned long int particleProblemNParticle(ParticleProblem *p);
void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax);
void particleProblemWriteVtu(ParticleProblem *p, const char *filename);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
int *particleProblemSegmentTag(ParticleProblem *p);
#if DIMENSION == 3
int *particleProblemTriangleTag(ParticleProblem *p);
double *particleProblemTriangle(ParticleProblem *p);
#endif
#endif
......@@ -30,6 +30,13 @@
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
return array;
}
PyObject *tagVector2PyArray(int *fm)
{
npy_intp dims[1] = {vectorSize(fm)};
PyObject *array = PyArray_New(&PyArray_Type, 1, dims, NPY_INT, NULL, (void*)fm, 0, NPY_ARRAY_C_CONTIGUOUS|NPY_ARRAY_WRITEABLE, NULL);
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
return array;
}
int pySequenceToCoord(PyObject *o, double x[DIMENSION]) {
if (!PySequence_Check(o))
return 0;
......@@ -67,16 +74,21 @@ struct ParticleProblem{};
return self;
}
void addParticle(const double x[DIMENSION], double r) { particleProblemAddParticle($self->_this, x, r); }
size_t addBoundaryDisk(const double x0[DIMENSION], double r) {return particleProblemAddBoundaryDisk($self->_this, x0, r);}
size_t addBoundarySegment(const double x0[DIMENSION], const double x1[DIMENSION]) {return particleProblemAddBoundarySegment($self->_this, x0, x1);}
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
void addBoundaryTriangle(const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION]) {particleProblemAddBoundaryTriangle($self->_this, x0, x1, x2);}
void addBoundaryTriangle(const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], int tag) {particleProblemAddBoundaryTriangle($self->_this, x0, x1, x2, tag);}
#endif
void write(const char *filename) { particleProblemWrite($self->_this, filename); }
void iterate(double alert, double dt, double tol = 0.005) { particleProblemIterate($self->_this, alert, dt, tol); }
double maxdt(double alert) { return particleProblemMaxDt($self->_this, alert); }
void load(const char *filename) { particleProblemLoad($self->_this, filename); }
PyObject *diskTag() {return tagVector2PyArray(particleProblemDiskTag($self->_this));}
PyObject *segmentTag() {return tagVector2PyArray(particleProblemSegmentTag($self->_this));}
#if DIMENSION == 3
PyObject *triangleTag() {return tagVector2PyArray(particleProblemTriangleTag($self->_this));}
#endif
PyObject *velocity() {return coordVector2PyArray(particleProblemVelocity($self->_this), particleProblemNParticle($self->_this), DIMENSION);}
PyObject *disks() {
double *d = particleProblemDisk($self->_this);
......
......@@ -123,6 +123,7 @@ class ObjectCollection {
std::string type;
std::vector<float> vel;
dim = -1;
int tag;
while (hf>>type) {
if (type == "DIMENSION") {
hf >> dim;
......@@ -143,13 +144,17 @@ class ObjectCollection {
else if (type == "D") {
readCoord(hf, coord, 1);
double r;
hf >> r;
hf >> r >> tag;
radius.push_back(r);
fixed.push_back(true);
}
else if (type == "S") readCoord(hf, segcoord, 2);
else if (type == "S"){
readCoord(hf, segcoord, 2);
hf >> tag;
}
else if (type == "T"){
readCoord(hf, tricoord, 3);
hf >> tag;
}
else {
std::cout <<"Unknown type \""<< type << "\"in file \""<<filename<<"\"\n";
......
......@@ -6,25 +6,27 @@ import random
import os, time, sys
import shutil, random
outputdir ="output"
os.makedirs(outputdir, exist_ok=True)
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
addv = set()
p = scontact2.ParticleProblem()
def loadMeshBoundary(p, filename, tags) :
def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
m = mesh.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], v[1]), 0.000)
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], v0[1]), (v1[0], v1[1]))
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)
......@@ -70,4 +72,5 @@ while t < tEnd :
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
print("cpu time : %g" % (time.clock() - tic))
loadMeshBoundary(p, "mesh.msh", ["Pile"], [0, 0.5])
p.write("deposited")
#!/usr/bin/env python
import scontact2
import numpy as np
import mesh
import random
import os, time, sys
import shutil, random
import os, time, shutil
import runf
outputdir ="output_run"
os.makedirs(outputdir, exist_ok=True)
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
p = scontact2.ParticleProblem()
p.load("deposited")
def meshVelocity(cm, v, x) :
global V
l = 0.1
L = 2
h = 0.5
H = 1.5
fx = np.minimum(1 - (np.abs(x[:, 0]) - l) / (L - l), 1.)
HX = Y * fx
HH = np.where(x[:, 1] > HX, H - HX , H + HX)
fy = np.minimum(1 - (np.abs(x[:, 1] - HX) - h) / (HH - h), 1.)
v[:, 0] = 0
v[:, 1] = fx * fy * V
v[:, 2] = 0
def moveMesh(m, dt):
global Y
global V
for i, v in enumerate(m.vertices) :
x = np.array([[v[0], v[1], v[2]]])
vv = np.array([[0., 0., 0.]])
meshVelocity(None, vv, x)
m.vertices[i] = [x[0, 0] + vv[0, 0] * dt, x[0, 1] + vv[0, 1] * dt, v[2], v[3]]
Y += V * dt
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");
def loadMeshBoundary(p, filename, tags) :
m = mesh.mesh(filename)
addv = set()
pnum = [m.getPhysicalNumber(1, i) for i in tags]
piledisks = []
pilesegs = []
for ent in m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in addv :
continue
addv.add(v[3])
piledisks.append(p.addBoundaryDisk((v[0], v[1]), 0.))
v0 = el.vertices[1]
v1 = el.vertices[0]
pilesegs.append(p.addBoundarySegment((v0[0], v0[1]), (v1[0], v1[1])))
return piledisks, pilesegs
class meshVelocity :
def __init__(self) :
self.Y = 0
self.V = 0
def cb(self, cm, v, x) :
l = 0.1
L = 2
h = 0.5
H = 1.5
fx = np.minimum(1 - (np.abs(x[:, 0]) - l) / (L - l), 1.)
HX = self.Y * fx
HH = np.where(x[:, 1] > HX, H - HX , H + HX)
fy = np.minimum(1 - (np.abs(x[:, 1] - HX) - h) / (HH - h), 1.)
v[:, 0] = 0
v[:, 1] = fx * fy * self.V
def update(self, dt) :
self.Y += self.V * dt
def writep(p, odir, i) :
outputname = "%s/part-%05i" % (odir, i)
......@@ -72,51 +34,50 @@ def writep(p, odir, i) :
shutil.move(tmpoutputname, outputname)
p.writeVtu("%s/part-%05i.vtu" % (odir, i))
r = 0.015
#p.useJacobi(True)
Y = 0
m = mesh.mesh("mesh.msh")
global V
V = 0.5
moveMesh(m, 1)
m.write("meshI.msh")
piledisks, pilesegs = loadMeshBoundary(p, "meshI.msh", ["Pile"])
g = 0.5
tEnd = 18
tCycle = 3
r = 0.015
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
hydro.moveMesh()
mv.update(dt)
t = 0
ii = 0
tic = time.clock()
outputname = "%s/part-%05i" % (outputdir, 0)
p.velocity()[:, :] = 0;
p.write(outputname)
writep(p, outputdir, 0)
p.velocity()[:, :] = 0;
ysix = np.ix_(pilesegs, (1, 3))
dt = 0.1 * r * abs(tCycle)
outf = 1
hydro = runf.hydro("meshI.msh", meshVelocity, dt, r, rhop=2)
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]
tic = time.clock()
while t < tEnd :
f = (-1 if ((t % tCycle) < (tCycle / 2)) else 1) * 1/ (tCycle/2)
f = 0
V = f
p.disks()[piledisks, 3] = f
p.segments()[pilesegs, 5] = f
p.externalForces()[:, :] = 0
hydro.solve(p.position(), np.zeros_like(p.velocity()), p.externalForces())
p.externalForces()[:, 1] += -10
mv.V = (-1 if ((t % tCycle) < (tCycle / 2)) else 1) * 1/ (tCycle/2)
p.disks()[piledisk, 3] = mv.V
p.segments()[pileseg, 5] = mv.V
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.iterate(r, dt, 1e-5)
p.disks()[piledisks, 1] += f * dt
p.segments()[ysix] += f * dt
p.disks()[piledisk, 1] += mv.V * dt
p.segments()[pileseg, 1] += mv.V * dt
p.segments()[pileseg, 3] += mv.V * dt
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf)
writep(p, outputdir, ioutput)
hydro.writeMsh(outputdir, ioutput)
savepartposvec("%s/forces_%i.pos" % (outputdir, ioutput), "forces", p.position(), p.externalForces())
ii += 1
Y += V * dt
print("%.2g/%.2g" % (t, tEnd))
print("cpu time : %g" % (time.clock() - tic))
......@@ -5,7 +5,7 @@ import numpy as np
class hydro :
def bndMesh(self, cm) :
val = np.ndarray((cm.nEvaluationPoint(), 6))
val[:, 0:3] = cm.XXget(self.vmeshf)
val[:, 0:2] = cm.XXget(self.vmeshf)
val[:, 3:5] = 1
val[:, 5] = -1
return val
......@@ -22,51 +22,49 @@ class hydro :
val[:, 2] = 0.
return val
def __init__(self, fmesh, meshVelocityCB, dt, r, rhop = 1, rho = 1, mu = 1) :
def __init__(self, fmesh, meshVelocityCB, dt, r, rhop = 1, rho = 1, mu = 1, g = [0, -1]) :
self.model = GModel()
self.model.load(fmesh)
self.groups = dgGroupCollection(self.model, 2, -2)
self.vmeshf = functionNumpy(3, meshVelocityCB, [function.getCoordinates()])
self.updateMeshf = functorNumpy(lambda cm : function.getCoordinates().get(cm) + self.vmeshf.get(cm) * self.dt)
self.vmeshf = functionNumpy(2, meshVelocityCB, [function.getCoordinates()])
self.updateMeshf = functorNumpy(lambda cm : function.getCoordinates().get(cm)[:, 0:2] + self.vmeshf.get(cm) * self.dt)
self.porosity = dgDofContainer(self.groups, 1)
self.momentum = dgDofContainer(self.groups, 2)
self.r = r
self.rhop = rhop
self.rho = rho
self.dt = dt
self.g = np.array([g])
self.lsys = linearSystemPETScDouble()
self.lsys.setParameter("petscOptions", "-pc_type lu -ksp_atol 1e-16 -ksp_rtol 1e-8")
self.dof = dgDofManager.newCG(self.groups, 3, self.lsys)
self.law = brinkman.brinkman(
law = brinkman.brinkman(
mu = functionConstant(mu/rho),
k = functionConstant((self.r*2)**2/(150 * mu/rho)*1e6),
k = functionConstant((r*2)**2/(150 * mu/rho)),
porosity = self.porosity.getFunction(),
particleMomentum = self.momentum.getFunction(),
stab = dt/self.rhop,
stab = dt/rhop,
ns = False)
self.law.setP1P1Mini(self.dof)
self.law.bndBox = functorNumpy(self.bndMesh)
self.law.addStrongBoundaryCondition(1, "Pile", self.law.bndBox)
self.law.addStrongBoundaryCondition(1, "Box", self.law.bndBox)
self.law.ptFix = functionConstant([0, 0, 0, -1, -1, 1])
self.law.addStrongBoundaryCondition(0, "PtFix", self.law.ptFix)
self.law.setBoundary0Flux("Box")
self.law.setBoundary0Flux("Pile")
law.setP1P1Mini(self.dof)
law.bndBox = functorNumpy(self.bndMesh)
law.addStrongBoundaryCondition(1, "Pile", law.bndBox)
law.addStrongBoundaryCondition(1, "Box", law.bndBox)
law.ptFix = functionConstant([0, 0, 0, -1, -1, 1])
law.addStrongBoundaryCondition(0, "PtFix", law.ptFix)
law.setBoundary0Flux("Box")
law.setBoundary0Flux("Pile")
self.totalVelocityF = functorNumpy(self.totalVelocity)
self.solidVelocityF = functorNumpy(self.solidVelocity)
#self.solver = dgSteady(self.law, self.dof)
self.solver = dgDIRK(self.law, self.dof, 1)
#self.solver = dgSteady(law, self.dof)
self.solver = dgDIRK(law, self.dof, 1)
self.sol = dgDofContainer(self.groups, 3)
self.law = law
self.sol.setAll(0)
def solve(self, position, velocity, forces) :
def moveMesh(self) :
self.groups.updateCoordinates(self.updateMeshf)
self.groups.mesh().updateGModelPosition()
self.particleMesh = particleMesh(self.groups.getModel(), ["Domain"]);
N = position.shape[0]
volume = np.ndarray((N, 1), 'd')
volume[:, :] = np.pi * self.r * self.r
def solve(self, volume, position, velocity, forces) :
self.particleMesh = particleMesh(self.groups.getModel(), ["Domain"]);
self.particleOnMesh = particleList(self.particleMesh, position)
vertexCompacity = self.particleOnMesh.particle2mesh(volume) / self.particleMesh.vertexVolume()
porosity2d = 1 - vertexCompacity
......@@ -77,7 +75,8 @@ class hydro :
self.solver.iterate(self.sol, self.dt, 0)
gradp = self.particleOnMesh.function2particle(self.groups, self.sol.getFunctionGradient())[:, 6:8]
forces[:,:] = -gradp/self.rhop
#add hydrostatic pressure
forces[:,:] = -(gradp + self.g * self.rho) * volume
def writeMsh(self, odir, oiter) :
self.groups.exportFunctionMsh(self.law.pressure, "%s/f_%05i" % (odir, oiter), 0., oiter, "pressure", self.sol)
......@@ -85,5 +84,6 @@ class hydro :
self.groups.exportFunctionMsh(self.porosity.getFunction(), "%s/f_%05i" % (odir, oiter), 0., oiter, "porosity", self.sol)
self.groups.exportFunctionMsh(function.getCoordinates(), "%s/f_%05i" % (odir, oiter), 0., oiter, "xmesh")
self.groups.exportFunctionMsh(self.totalVelocityF, "%s/f_%05i" % (odir, oiter), 0., oiter, "totalVelocity")
self.groups.exportFunctionMsh(self.solidVelocityF, "%s/f_%05i" % (odir, oiter), 0., oiter, "solidVelocity")
self.groups.exportFunctionMsh(self.totalVelocityF, "%s/f_%05i" % (odir, oiter), 0., oiter, "totalVelocity")
self.groups.exportFunctionMsh(self.vmeshf, "%s/f_%05i" % (odir, oiter), 0., oiter, "meshVelocity")
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