Commit 77aeb12d authored by Matthieu Constant's avatar Matthieu Constant
Browse files

New drag scheme

parent d0edf6d7
all :
touch *.c
python setup.py build_ext --inplace
python2 setup.py build_ext --inplace
petsc :
touch *.c
python setup_petsc.py build_ext --inplace
python2 setup_petsc.py build_ext --inplace
clean :
python setup.py clean
......@@ -17,8 +17,10 @@ cdef extern from "fluid_problem.h":
double rho
int n_particles
double *particle_force
double *node_force
pass
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int it);
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_boundary, StrongBoundary *boundaries);
......@@ -75,13 +77,15 @@ cdef class fluid_problem :
fluid_problem_adapt_mesh(self._fp, gradmin, gradmax, lcmin)
def export(self, output_dir, double t, int it) :
fluid_problem_export(self._fp, output_dir.encode(), t, it)
def implicit_euler(self, double dt) :
fluid_problem_implicit_euler(self._fp, dt)
def compute_node_force(self, double dt) :
fluid_problem_compute_node_particle_force(self._fp, dt)
forces = numpy.ndarray([self._fp.n_particles,2],"d")
for i in range(self._fp.n_particles) :
forces[i,0] = self._fp.particle_force[i*2+0]
forces[i,1] = self._fp.particle_force[i*2+1]
return forces
def implicit_euler(self, double dt) :
fluid_problem_implicit_euler(self._fp, dt)
def set_particles(self, mass, volume, position, velocity):
cdef double[:,::1] _mass = numpy.ascontiguousarray(mass)
......
......@@ -42,6 +42,8 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
return -1;
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, 0, 1))
return -1;
if(mesh_write_msh_vector(problem->mesh, f, "force", t, iter, problem->node_force, 2, 0, 1))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, 2))
return -1;
fclose(f);
......@@ -61,9 +63,13 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
drag[1] = GoU*u[1]/(1+dt/mass*GoU);
}
static void fluid_problem_add_particle_force(FluidProblem *problem, const double *solution, double dt) {
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
for (int i = 0; i < mesh->n_nodes*2; ++i) {
problem->node_force[i] = 0.0;
}
for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 2*i;
double phi[3] = {1-xi[0]-xi[1],xi[0],xi[1]};
......@@ -108,42 +114,17 @@ static void fluid_problem_add_particle_force(FluidProblem *problem, const double
double d = 2*sqrt(vol/M_PI);
double drag[2];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass);
double ddragdu[2],ddragdv[2];
double eps = 1e-8;
double duu[2] = {du[0]+eps, du[1]};
particle_drag(duu,problem->mu,problem->rho,d,c,ddragdu, dt, mass);
double duv[2] = {du[0], du[1]+eps};
particle_drag(duv,problem->mu,problem->rho,d,c,ddragdv, dt, mass);
ddragdu[0] = -(ddragdu[0]-drag[0])/(eps*c);
ddragdu[1] = -(ddragdu[1]-drag[1])/(eps*c);
ddragdv[0] = -(ddragdv[0]-drag[0])/(eps*c);
ddragdv[1] = -(ddragdv[1]-drag[1])/(eps*c);
problem->particle_force[i*2+0] = -drag[0]-vol*gradp[0];
problem->particle_force[i*2+1] = -drag[1]-9.81*(mass-problem->rho*vol)-vol*gradp[1];
double local_vector[3*n_fields];
double local_matrix[9*n_fields*n_fields];
for (int iphi = 0; iphi < 3; ++iphi) {
local_vector[iphi+0] = -drag[0]*phi[iphi];
local_vector[iphi+3] = -drag[1]*phi[iphi];
local_vector[iphi+6] = 0;
for (int jphi = 0; jphi < 3; ++jphi) {
local_matrix[(iphi+0)*9+jphi+0] = -phi[iphi]*ddragdu[0]*phi[jphi];
local_matrix[(iphi+0)*9+jphi+3] = -phi[iphi]*ddragdv[0]*phi[jphi];
local_matrix[(iphi+0)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][0];
local_matrix[(iphi+3)*9+jphi+0] = -phi[iphi]*ddragdu[1]*phi[jphi];
local_matrix[(iphi+3)*9+jphi+3] = -phi[iphi]*ddragdv[1]*phi[jphi];
local_matrix[(iphi+3)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][1];
local_matrix[(iphi+6)*9+jphi+0] = 0;
local_matrix[(iphi+6)*9+jphi+3] = 0;
local_matrix[(iphi+6)*9+jphi+6] = 0;
}
problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi]/problem->node_volume[tri[iphi]];
problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi]/problem->node_volume[tri[iphi]];
}
linear_system_add_to_matrix(problem->linear_system,iel,iel,local_matrix);
linear_system_add_to_rhs(problem->linear_system,iel,local_vector);
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt)
{
LinearSystem *lsys = problem->linear_system;
......@@ -151,6 +132,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
const double *solution = problem->solution;
const double *porosity = problem->porosity;
const double *old_porosity = problem->old_porosity;
const double *node_force = problem->node_force;
double *local_vector = malloc(sizeof(double)*n_fields*3);
double *local_matrix = malloc(sizeof(double)*n_fields*n_fields*9);
const double mu = problem->mu;
......@@ -184,6 +166,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
double DCDT[] = {(C[0]-OLDC[0])/dt,(C[1]-OLDC[1])/dt,(C[2]-OLDC[2])/dt};
double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
double FU[] = {node_force[tri[0]*2+0], node_force[tri[1]*2+0], node_force[tri[2]*2+0]};
double FV[] = {node_force[tri[0]*2+1], node_force[tri[1]*2+1], node_force[tri[2]*2+1]};
double U_old[] = {solution_old[tri[0]*3+0],solution_old[tri[1]*3+0],solution_old[tri[2]*3+0]};
double V_old[] = {solution_old[tri[0]*3+1],solution_old[tri[1]*3+1],solution_old[tri[2]*3+1]};
double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
......@@ -208,6 +192,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
const double phi[] = {1-xi-eta, xi, eta};
double u = phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2];
double v = phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2];
double fu = phi[0]*FU[0]+phi[1]*FU[1]+phi[2]*FU[2];
double fv = phi[0]*FV[0]+phi[1]*FV[1]+phi[2]*FV[2];
double u_old = phi[0]*U_old[0]+phi[1]*U_old[1]+phi[2]*U_old[2];
double v_old = phi[0]*V_old[0]+phi[1]*V_old[1]+phi[2]*V_old[2];
double p = phi[0]*P[0]+phi[1]*P[1]+phi[2]*P[2];
......@@ -222,13 +208,14 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
local_vector[iphi+0] += jw*(
2*mu*(dphii[0]*tau[0][0]+dphii[1]*(tau[0][1]+tau[1][0])*0.5)
+rho*phii*((u-u_old)/dt+(u*tau[0][0]+v*tau[0][1])/c)
-p*dphii[0]
-p*dphii[0]-fu*phii
);
local_vector[iphi+3] += jw*(
2*mu*(dphii[0]*(tau[0][1]+tau[1][0])*0.5+dphii[1]*tau[1][1])
+rho*phii*((v-v_old)/dt+(u*tau[1][0]+v*tau[1][1])/c)
-p*dphii[1]
-p*dphii[1]-fv*phii
);
local_vector[iphi+6] += jw*(
epsilon*(dphii[0]*dp[0]+dphii[1]*dp[1])
+phii*((du[0]+dv[1])+dcdt)
......@@ -258,7 +245,6 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
}
free(local_vector);
free(local_matrix);
fluid_problem_add_particle_force(problem, solution, dt);
}
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
......@@ -313,6 +299,30 @@ void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
free(solution_old);
}
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(problem->node_volume);
Mesh *mesh = problem->mesh;
problem->node_volume = malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->node_volume[i] = 0;
}
for (int iel = 0; iel < mesh->n_triangles; ++iel) {
const int *tri = &mesh->triangles[iel*3];
const double x[3][2] = {
{mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
{mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
{mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
const double dxdxi[2][2] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1]}};
const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
problem->node_volume[tri[0]] += detj/6;
problem->node_volume[tri[1]] += detj/6;
problem->node_volume[tri[2]] += detj/6;
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
FluidProblem *problem = malloc(sizeof(FluidProblem));
......@@ -331,6 +341,10 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
problem->strong_boundaries[i].field = strong_boundaries[i].field;
}
problem->porosity = malloc(mesh->n_nodes*sizeof(double));
problem->node_force = malloc(mesh->n_nodes*2*sizeof(double));
for (int i = 0; i < mesh->n_nodes*2; ++i) {
problem->node_force[i] = 0.0;
}
problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
// begin to remove
......@@ -341,6 +355,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
problem->solution[i*3+1] = 0;
problem->solution[i*3+2] = 0.;
}
problem->node_volume = malloc(0);
fluid_problem_compute_node_volume(problem);
// end to remove
problem->linear_system = linear_system_new(mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
problem->n_particles = 0;
......@@ -357,6 +373,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->node_force);
free(problem->porosity);
free(problem->old_porosity);
linear_system_free(problem->linear_system);
......@@ -376,25 +393,9 @@ void fluid_problem_free(FluidProblem *problem)
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
double *vertex_volume = malloc(sizeof(double)*mesh->n_nodes);
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 0;
vertex_volume[i] = 0;
}
for (int iel = 0; iel < mesh->n_triangles; ++iel) {
const int *tri = &mesh->triangles[iel*3];
const double x[3][2] = {
{mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
{mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
{mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
const double dxdxi[2][2] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1]}};
const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
vertex_volume[tri[0]] += detj/6;
vertex_volume[tri[1]] += detj/6;
vertex_volume[tri[2]] += detj/6;
}
for (int i = 0; i < problem->n_particles; ++i) {
const int *tri = &mesh->triangles[problem->particle_element_id[i]*3];
......@@ -404,9 +405,8 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
problem->porosity[tri[2]] += v*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = (1-problem->porosity[i]/vertex_volume[i]);
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
}
free(vertex_volume);
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
......@@ -418,8 +418,10 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
new_mesh_size[i] = DBL_MAX;
gradmax = sqrt(gradmax);
gradmin = sqrt(gradmin);
double *porosity = problem->porosity;
for (int iel = 0; iel < mesh->n_triangles; ++iel) {
int *tri = problem->mesh->triangles+iel*3;
double C[] = {porosity[tri[0]], porosity[tri[1]], porosity[tri[2]]};
double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
const double x[3][2] = {
......@@ -440,12 +442,12 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
{dxidx[1][0],dxidx[1][1]}
};
double gradU[2] = {
dphi[0][0]*U[0]+dphi[1][0]*U[1]+dphi[2][0]*U[2],
dphi[0][1]*U[0]+dphi[1][1]*U[1]+dphi[2][1]*U[2]
dphi[0][0]*U[0]/C[0]+dphi[1][0]*U[1]/C[1]+dphi[2][0]*U[2]/C[2],
dphi[0][1]*U[0]/C[0]+dphi[1][1]*U[1]/C[1]+dphi[2][1]*U[2]/C[2]
};
double gradV[2] = {
dphi[0][0]*V[0]+dphi[1][0]*V[1]+dphi[2][0]*V[2],
dphi[0][1]*V[0]+dphi[1][1]*V[1]+dphi[2][1]*V[2]
dphi[0][0]*V[0]/C[0]+dphi[1][0]*V[1]/C[1]+dphi[2][0]*V[2]/C[2],
dphi[0][1]*V[0]/C[0]+dphi[1][1]*V[1]/C[1]+dphi[2][1]*V[2]/C[2]
};
double ngrad = pow(gradU[0]*gradU[0]+gradU[1]*gradU[1] + gradV[0]*gradV[0] + gradV[1]*gradV[1], 0.25);
double lc = lcmin /fmin(1,fmax(ngrad/gradmax,gradmin/gradmax));
......@@ -499,12 +501,15 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->old_porosity);
problem->old_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->node_force);
problem->node_force = malloc(sizeof(double)*2*new_mesh->n_nodes);
mesh_free(problem->mesh);
problem->mesh = new_mesh;
mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_porosity(problem);
linear_system_free(problem->linear_system);
problem->linear_system = linear_system_new(problem->mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
problem->linear_system = linear_system_new(problem->mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
fluid_problem_compute_node_volume(problem);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
......
......@@ -12,6 +12,7 @@ typedef struct {
double *porosity;
double *old_porosity;
double *solution;
double *node_volume;
int n_strong_boundaries;
StrongBoundary *strong_boundaries;
LinearSystem *linear_system;
......@@ -23,10 +24,12 @@ typedef struct {
double *particle_velocity;
double *particle_mass;
double *particle_force; // complete force on the particle (including gravity)
double *node_force;
int *particle_element_id;
} FluidProblem;
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter);
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
......
#!/usr/bin/env python
import fluid as Fluid
from pyfefp import scontact2Interface
import scontact2
import numpy as np
import os
import time
import shutil
def genInitialPosition(filename, r, rout, rhop) :
p = scontact2.ParticleProblem()
scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
x = np.arange(rout, -rout, -2*r)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.13
p.write(filename)
outputdir = "outputtest41"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
r=25e-6/(2)
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = [0, -9.81]
rho = 1200
rhop = 2400
nu = 1e-4
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 2
#numerical parameters
h = 0.001 # approx r*100 but should match the mesh size
dt = 1e-4
c = h/dt * 7/50/2#50
#dt= 1e-6#0.5*h/c
epsilon = 1e-6#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
print("c = %g, eps = %g" % (c, epsilon))
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, 1e-3, rhop)
p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 100
outf1 = 100000
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
fluid = Fluid.fluid_problem("mesh.msh",nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
forces = g*p.mass()
while t < tEnd :
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
if ii%10==0 and ii != 0:
fluid.adapt_mesh(0.01,10,5e-4)
forces = fluid.compute_node_force(dt)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir, ioutput, t)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
Supports Markdown
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