Commit 35b0def2 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

force implicit 2d

parent dbe00633
ALL: libmbfluid3.so libscontact2.so libscontact3.so
ALL: libmbfluid3.so libmbfluid2.so libscontact2.so libscontact3.so
CFLAGS=-Wno-unused-function -O3 -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt -Isrc
LDFLAGS=-shared -lm
......
......@@ -39,8 +39,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, lcmax, lcmin, n_el) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
def adapt_mesh(self, lcmax, lcmin, n_el, alpha, notComputeEpsilon) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(alpha), c_int(notComputeEpsilon))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
......
......@@ -77,6 +77,9 @@ class ParticleProblem :
def write_vtk(self, odir, i, t) :
lib.particleProblemWriteVtk(self._p, odir.encode(), i)
def read_vtk(self,dirname,i) :
lib.particleProblemReadVtk(self._p, dirname.encode(),i)
def add_boundary_disk(self, x0, r, tag) :
lib.particleProblemAddBoundaryDisk.restype = c_size_t
return lib.particleProblemAddBoundaryDisk(self._p, coord_type(*x0), c_double(r), tag.encode())
......
......@@ -6,7 +6,6 @@
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"
#define DIMENSION 3
#define n_fields (DIMENSION+1)
......@@ -109,12 +108,12 @@ static double particle_drag_coeff(double u[2], double mu, double rho, double vol
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u;
double f = pow(c,-1.8);
return = f*Cd_u*rho/2*d;
return f*Cd_u*rho/2*d;
}
#else
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
{
double d = 2*pow(3./4.*vol/M_PI,1./3.);
double d = 2*cbrt(3./4.*vol/M_PI);
double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
//Reynolds/|u_p-u_f|
double Re_O_u = d*c/mu*rho;
......@@ -138,8 +137,6 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
const double *solution = problem->solution;
for (int i = 0; i < problem->n_particles; ++i) {
int iel = problem->particle_element_id[i];
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
if(iel < 0){
continue;
}
......@@ -170,6 +167,8 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/c;
due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
}
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
double rhop = mass/vol;
double drho = rhop -problem->rho;
double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
......@@ -184,7 +183,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int d = 0; d < DIMENSION; ++d)
local_vector[iphi+N_SF*d] += fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*d] += fcoeff*gamma*(du[d]/*-dt/mass*vol*gradp[d]*/)*phi[iphi];
local_vector[iphi+N_SF*1] += fcoeff*gamma*(-dt*g)*phi[iphi];
}
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -490,6 +489,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
#if DIMENSION == 2
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#else
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
......@@ -772,6 +772,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
#if DIMENSION == 2
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
}
#else
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
......@@ -779,8 +780,8 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
#endif
if(notComputeEpsilon){
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->epsilon[i] = alpha;
}
problem->epsilon[i] = alpha;
}
}
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
......@@ -813,7 +814,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
}
problem->particle_position = malloc(sizeof(double)*DIMENSION*n);
problem->particle_mass = malloc(sizeof(double)*n);
problem->particle_force = malloc(sizeof(double)*n*3);
problem->particle_force = malloc(sizeof(double)*n*DIMENSION);
problem->particle_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
}
......
......@@ -21,7 +21,7 @@ def genInitialPosition(filename, r, rout, rhop) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.45
p.write(filename)
p.write_vtk(filename,0,0)
outputdir = "output"
......@@ -46,18 +46,18 @@ tEnd = 100
#numerical parameters
lcmin = 0.0001 # approx r*100 but should match the mesh size
dt = 5e-4
dt = 5e-4*10
alpha = 2.5e-2
epsilon = alpha*lcmin**2 /nu
print('epsilon',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)
genInitialPosition(outputdir, r, 1e-3, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
......@@ -73,24 +73,28 @@ fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%20==0 and ii != 0):
fluid.adapt_mesh(0.01,200,1e-4)
if (ii%5==0 and ii != 0):
#fluid.adapt_mesh(0.01,200,1e-4)
fluid.adapt_mesh(1e-3,1e-4,10000, epsilon, 1)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
else :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
print(forces)
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)
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+"/part-%05d" % ioutput)
p.write_boundary_vtk(outputdir, ioutput, t)
#p.write(outputdir+"/part-%05d" % ioutput)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -19,7 +19,7 @@ def genInitialPosition(filename, r, ly, rhop) :
for i in range(len(x)) :
rhop = random.choice([1100,1350,1600])
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.write(filename)
p.write_vtk(filename,0,0)
outputdir = "outputdifferentdensities"
......@@ -45,22 +45,23 @@ tEnd = 100
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-4
dt = 1e-4/5
alpha = 2.5e-6
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")
#p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, ly, rhop)
genInitialPosition(outputdir, r, ly, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 200
outf = 5
outf1 = 100000
ii = 0
......@@ -84,7 +85,9 @@ fluid.export(outputdir,0,0)
tic = time.clock()
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
......@@ -93,7 +96,6 @@ while t < tEnd :
for i in range(nsub) :
p.iterate(dt/nsub, forces)
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
......
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