Commit 9e584e2d authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

smooth velocity

parent 8dbdb427
......@@ -45,11 +45,11 @@ class fluid_problem :
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
def compute_node_force(self, dt) :
def compute_node_force(self, dt, lsmooth) :
def np2c(a) :
return c_void_p(a.ctypes.data)
forces = numpy.ndarray([self.n_particles,2],"d")
lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), np2c(forces))
lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_double(lsmooth), np2c(forces))
return forces
def implicit_euler(self, dt) :
......
......@@ -56,9 +56,87 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
drag[0] = GoU*(rhop*u[0]-gradp[0]*dt)/(rhop+dt*GoU*h/(mass/rhop));
}
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
static double mesh_get_min_edge_length(Mesh *mesh) {
double emin = DBL_MAX;
for (int i= 0; i < mesh->n_triangles; ++i) {
const int *tri = mesh->triangles+i*3;
for (int j = 0; j < 3; j++) {
double *x0 = mesh->x + tri[j]*3;
double *x1 = mesh->x + tri[(j+1)%3]*3;
double dx[] = {x0[0]-x1[0],x0[1]-x1[1],x0[2]-x1[2]};
double l2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
emin = fmin(l2,emin);
}
}
return sqrt(emin);
}
static void fluid_problem_smooth_velocity(FluidProblem *problem, double lsmooth, double *smooth_velocity) {
Mesh *mesh = problem->mesh;
double emin = mesh_get_min_edge_length(mesh);
int n = ceil(lsmooth*lsmooth/(emin*emin));
double alpha = lsmooth*lsmooth/n;
double *rhs = malloc(mesh->n_nodes*2*mesh->n_nodes);
double *mass = malloc(mesh->n_nodes*mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i){
smooth_velocity[i*2+0] = problem->solution[i*3+0];
smooth_velocity[i*2+1] = problem->solution[i*3+1];
}
for (int iter = 0; iter < n; ++iter) {
for (int i = 0; i < mesh->n_nodes; ++i){
rhs[i*2+0] = 0;
rhs[i*2+1] = 0;
mass[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];
const double dxidx[2][2] = {
{dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
{-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
};
const double dphi[][2] = {
{-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
{dxidx[0][0],dxidx[0][1]},
{dxidx[1][0],dxidx[1][1]}
};
double U[] = {smooth_velocity[tri[0]*2+0],smooth_velocity[tri[1]*2+0],smooth_velocity[tri[2]*2+0]};
double V[] = {smooth_velocity[tri[0]*2+1],smooth_velocity[tri[1]*2+1],smooth_velocity[tri[2]*2+1]};
double du[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]
};
double dv[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]
};
for (int i = 0; i < 3; ++i) {
rhs[tri[i]*2+0] -= (dphi[i][0]*du[0] + dphi[i][1]*du[1])*detj/2;
rhs[tri[i]*2+1] -= (dphi[i][0]*dv[0] + dphi[i][1]*dv[1])*detj/2;
mass[tri[i]] += detj/6;
}
}
for (int i = 0; i < mesh->n_nodes; ++i) {
smooth_velocity[i*2+0] += rhs[i*2+0]/mass[i]*alpha;
smooth_velocity[i*2+1] += rhs[i*2+1]/mass[i]*alpha;
}
}
free(mass);
free(rhs);
}
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double lsmooth, double *particle_force) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
double *smoothv = malloc(mesh->n_nodes*2*sizeof(double));
fluid_problem_smooth_velocity(problem, lsmooth, smoothv);
const double *solution = problem->solution;
for (int i = 0; i < mesh->n_nodes*2; ++i) {
problem->node_force[i] = 0.0;
......@@ -74,11 +152,10 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
if(iel < 0){
continue;
}
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]};
double U[] = {smoothv[tri[0]*2+0],smoothv[tri[1]*2+0],smoothv[tri[2]*2+0]};
double V[] = {smoothv[tri[0]*2+1],smoothv[tri[1]*2+1],smoothv[tri[2]*2+1]};
double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
double vf[2] = {
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2],
......@@ -131,6 +208,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
problem->node_force[i*2+0] /= problem->node_volume[i];
problem->node_force[i*2+1] /= problem->node_volume[i];
}
free(smoothv);
}
static void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt)
......
......@@ -38,7 +38,7 @@ typedef struct {
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter);
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double lsmooth, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
......
......@@ -74,7 +74,7 @@ forces = g*p.mass()
while t < tEnd :
# if (ii%50==0 and ii != 0):
# fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon)
forces = fluid.compute_node_force(dt)
forces = fluid.compute_node_force(dt,r*10)
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()))))
......
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