Commit 1298f8e1 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

new formula for drag

parent e2e2d4eb
ALL: libmbfluid2.so libmbfluid3.so libscontact2.so libscontact3.so
ALL: libmbfluid3.so libscontact2.so libscontact3.so
CFLAGS=-Wno-unused-function -O3 -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt
LDFLAGS=-shared -lm
......
......@@ -6,6 +6,7 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#define DIMENSION 3
#define n_fields (DIMENSION+1)
......@@ -127,8 +128,6 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
int comp[]={0,1,2};
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
return -1;
if(mesh_write_msh_vector(problem->mesh, f, "force", t, iter, problem->node_force, DIMENSION, comp))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, DIMENSION))
return -1;
fclose(f);
......@@ -299,9 +298,6 @@ 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*DIMENSION; ++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]};
......@@ -374,6 +370,10 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
#else
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
for (int i = 0; i < problem->n_particles*DIMENSION; ++i) {
particle_force[i] = problem->particle_force[i];
}
/*
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
......@@ -405,11 +405,6 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
phi[0]*W[0]+phi[1]*W[1]+phi[2]*W[2]+phi[3]*W[3]
};
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2]+phi[3]*C[3];
/*printf("c=%e\n",c);
printf("phi[0]=%e\n",phi[0]);
printf("phi[1]=%e\n",phi[1]);
printf("phi[2]=%e\n",phi[2]);
printf("phi[3]=%e\n",phi[3]);*/
const double x[4][3] = {
{mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
{mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
......@@ -461,7 +456,112 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
problem->node_force[i*3+0] /= problem->node_volume[i];
problem->node_force[i*3+1] /= problem->node_volume[i];
problem->node_force[i*3+2] /= problem->node_volume[i];
}*/
}
static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 3*i;
double phi[4] = {1-xi[0]-xi[1]-xi[2],xi[0],xi[1],xi[2]};
int iel = problem->particle_element_id[i];
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
if(iel < 0){
continue;
}
uint32_t *tet = problem->mesh->elements+iel*4;
double C[] = {porosity[tet[0]],porosity[tet[1]],porosity[tet[2]],porosity[tet[3]]};
double U[] = {solution[tet[0]*4+0],solution[tet[1]*4+0],solution[tet[2]*4+0],solution[tet[3]*4+0]};
double V[] = {solution[tet[0]*4+1],solution[tet[1]*4+1],solution[tet[2]*4+1],solution[tet[3]*4+1]};
double W[] = {solution[tet[0]*4+2],solution[tet[1]*4+2],solution[tet[2]*4+2],solution[tet[3]*4+2]};
double P[] = {solution[tet[0]*4+3],solution[tet[1]*4+3],solution[tet[2]*4+3],solution[tet[3]*4+3]};
double vf[3] = {
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2]+phi[3]*U[3],
phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]+phi[3]*V[3],
phi[0]*W[0]+phi[1]*W[1]+phi[2]*W[2]+phi[3]*W[3]
};
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2]+phi[3]*C[3];
const double x[4][3] = {
{mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
{mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
{mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
{mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]}};
const double dxdxi[3][3] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0],x[3][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1],x[3][1]-x[0][1]},
{x[1][2]-x[0][2],x[2][2]-x[0][2],x[3][2]-x[0][2]}};
const double detj = dxdxi[0][0]*(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])-dxdxi[1][0]*(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])+dxdxi[2][0]*(dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2]);
const double dxidx[3][3] = {
{(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])/detj,-(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])/detj, (dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2])/detj},
{-(dxdxi[1][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[1][2])/detj,(dxdxi[0][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[0][2])/detj,-(dxdxi[0][0]*dxdxi[1][2]-dxdxi[1][0]*dxdxi[0][2])/detj },
{(dxdxi[1][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[1][1])/detj,-(dxdxi[0][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[0][1])/detj, (dxdxi[0][0]*dxdxi[1][1]-dxdxi[1][0]*dxdxi[0][1])/detj}
};
const double dphi[][3] = {
{-dxidx[0][0]-dxidx[1][0]-dxidx[2][0],-dxidx[0][1]-dxidx[1][1]-dxidx[2][1],-dxidx[0][2]-dxidx[1][2]-dxidx[2][2]},
{dxidx[0][0],dxidx[0][1],dxidx[0][2]},
{dxidx[1][0],dxidx[1][1],dxidx[1][2]},
{dxidx[2][0],dxidx[2][1],dxidx[2][2]}
};
double gradp[3] = {
dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2]+dphi[3][0]*P[3],
dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]+dphi[3][1]*P[3],
dphi[0][2]*P[0]+dphi[1][2]*P[1]+dphi[2][2]*P[2]+dphi[3][2]*P[3]
};
double du[3] = {
problem->particle_velocity[i*3+0]-vf[0]/c,
problem->particle_velocity[i*3+1]-vf[1]/c,
problem->particle_velocity[i*3+2]-vf[2]/c
};
problem->particle_force[i*3+0] = -du[0]*mass/dt;
problem->particle_force[i*3+1] = -du[1]*mass/dt;
problem->particle_force[i*3+2] = -du[2]*mass/dt;
double rhop = mass/vol;
double drho = rhop -problem->rho;
double drag[3] = {
mass*du[0]/dt - gradp[0]*vol,
mass*du[1]/dt +problem->g*vol*drho - gradp[1]*vol,
mass*du[2]/dt - gradp[2]*vol
};
double *local_vector = all_local_vector+iel*n_fields*N_SF;
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
for (int iphi = 0; iphi < 4; ++iphi) {
local_vector[iphi+N_SF*0] -= drag[0]*phi[iphi];
local_vector[iphi+N_SF*1] -= drag[1]*phi[iphi];
local_vector[iphi+N_SF*2] -= drag[2]*phi[iphi];
}
for (int iphi = 0; iphi < 4; ++iphi) {
for (int jphi = 0; jphi < 4; ++jphi) {
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
LOCAL_MATRIX(0,0) += phi[jphi]*phi[iphi]/(c*dt)*mass;
LOCAL_MATRIX(1,1) += phi[jphi]*phi[iphi]/(c*dt)*mass;
LOCAL_MATRIX(2,2) += phi[jphi]*phi[iphi]/(c*dt)*mass;
LOCAL_MATRIX(0,3) += phi[iphi]*vol*dphi[jphi][0];
LOCAL_MATRIX(1,3) += phi[iphi]*vol*dphi[jphi][1];
LOCAL_MATRIX(2,3) += phi[iphi]*vol*dphi[jphi][2];
}
}
/*
for (int iphi = 0; iphi < 4; ++iphi) {
problem->node_force[tet[iphi]*3+0] += drag[0]*phi[iphi];
problem->node_force[tet[iphi]*3+1] += drag[1]*phi[iphi];
problem->node_force[tet[iphi]*3+2] += drag[2]*phi[iphi];
}*/
}
/*for (int i = 0; i < mesh->n_nodes; ++i) {
problem->node_force[i*3+0] /= problem->node_volume[i];
problem->node_force[i*3+1] /= problem->node_volume[i];
problem->node_force[i*3+2] /= problem->node_volume[i];
}*/
}
#endif
......@@ -472,10 +572,18 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
const double *solution = problem->solution;
const double *porosity = problem->porosity;
const double *old_porosity = problem->old_porosity;
const double *node_force = problem->node_force;
const double mu = problem->mu;
const double rho = problem->rho;
const double epsilon = problem->epsilon;
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
double *all_local_vector = malloc(N_E*local_size*sizeof(double));
double *all_local_matrix = malloc(N_E*local_size*local_size*sizeof(double));
for (size_t i = 0; i < N_E*local_size; ++i)
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
all_local_matrix[i] = 0;
compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
......@@ -484,8 +592,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double local_vector[n_fields*N_SF] = {0};
double local_matrix[n_fields*n_fields*N_SF*N_SF] = {0};
double *local_vector = &all_local_vector[N_SF*n_fields*iel];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
......@@ -503,7 +611,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int j = 0; j < DIMENSION; ++j){
u[j] += solution[el[i]*n_fields+j]*phi[i];
dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*phi[i];
fu[j] += node_force[el[i]*DIMENSION+j]*phi[i];
//fu[j] += node_force[el[i]*DIMENSION+j]*phi[i];
}
p += solution[el[i]*n_fields+DIMENSION]*phi[i];
c += porosity[el[i]]*phi[i];
......@@ -572,6 +680,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(all_local_matrix);
free(all_local_vector);
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
......@@ -730,10 +840,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->strong_boundaries[i].field = strong_boundaries[i].field;
}
problem->porosity = malloc(mesh->n_nodes*sizeof(double));
problem->node_force = malloc(mesh->n_nodes*DIMENSION*sizeof(double));
for (int i = 0; i < mesh->n_nodes*DIMENSION; ++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
......@@ -754,6 +860,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->particle_element_id = malloc(0);
problem->particle_position = malloc(0);
problem->particle_mass = malloc(0);
problem->particle_force = malloc(0);
problem->particle_volume = malloc(0);
problem->particle_velocity = malloc(0);
return problem;
......@@ -762,7 +869,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->node_force);
free(problem->porosity);
free(problem->old_porosity);
hxtLinearSystemDelete(&problem->linear_system);
......@@ -770,6 +876,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
for (int i = 0; i < problem->n_strong_boundaries; ++i)
......@@ -926,8 +1033,6 @@ 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)*DIMENSION*new_mesh->n_nodes);
mesh_free(problem->mesh);
problem->mesh = new_mesh;
for (int i = 0; i < problem->n_particles; ++i)
......@@ -955,6 +1060,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
problem->particle_uvw = malloc(sizeof(double)*DIMENSION*n);
......@@ -964,6 +1070,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_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
}
......
......@@ -31,7 +31,7 @@ typedef struct {
double *particle_volume;
double *particle_velocity;
double *particle_mass;
double *node_force;
double *particle_force;
int *particle_element_id;
} FluidProblem;
......
......@@ -52,7 +52,7 @@ tEnd = 300
#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-5
dt = 5e-5*1000
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
......@@ -111,10 +111,12 @@ ii = 0
t = 0
tic = time.clock()
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
fluid.adapt_mesh(0.01,50,1,50,4e-4,100000)
forces = fluid.compute_node_force(dt)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
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]))
......@@ -123,7 +125,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