Commit 015397a0 authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

marche pas

parent 1d5e3609
Pipeline #5451 passed with stage
in 28 seconds
......@@ -361,7 +361,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*dfddp = gammastar*dt/mass*vol-vol;
}
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double *u_mesh,const double *u_mesh_old, const double *du_mesh_vec, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc,const double cold, const double *dcold,const double *u_mesh,const double *u_mesh_old, const double *du_mesh_vec, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
......@@ -385,7 +385,9 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
divu_mesh -= u_mesh[iD]*dc[iD];
}
//f0[P] = (c-cold)/dt ;// + (sol[P]-solold[P])/dt*0.1;
f0[P] = (c-cold)/dt + divu_mesh;
f0[P] = (c-cold)/dt+ divu_mesh;
//f0[P] = divu_mesh;
//f00[P*n_fields+P] = 1/dt*0.1;
double G[D] = {0};
G[1] = -problem->g;
......@@ -401,13 +403,17 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
// = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
double fo[D]={0},fo1[2*D]={0},foo[D]={0};
for (int j = 0; j < D; ++j) {
f0[U+i] += rho/c*((uold[j]-c*u_mesh[j])*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
//f0[U+i] += rho/c*(u[j]-c*u_mesh[j])*du[i][j];
//f0[U+i] += rho*((uold[j]/cold-u_mesh[j])*du[i][j] + u[i]*(duold[j][j]/cold-du_mesh[j][j]-uold[j]/(cold*cold)*dcold[j]));
fo[U+i] += rho*((uold[j]/cold)*du[i][j] + u[i]*(duold[j][j]/cold-uold[j]/(cold*cold)*dcold[j]));
//f00[(U+i)*n_fields+U+i] += rho*(duold[j][j]/cold-du_mesh[j][j]-uold[j]/(cold*cold)*dcold[j]);
f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
//f00[(U+i)*n_fields+U+i] += rho/c*du[i][j];
//f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(u[j]-c*u_mesh[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(uold[j]-c*u_mesh[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(u[j]-c*u_mesh[j]);
//f01[(U+i)*n_fields*D+(U+i)*D+j] += rho*(uold[j]/cold-u_mesh[j]);
fo1[i+j] += rho*(uold[j]/cold);
foo[U+i] += rho*(duold[j][j]/cold-uold[j]/(cold*cold)*dcold[j]);
}
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
......@@ -419,14 +425,17 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double supg = uold[j]*taup;
f1[(U+i)*D+j] += supg*f0[U+i];
//f1[(U+i)*D+j] += supg*fo[U+i];
f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
//f10[((U+i)*D+j)*n_fields+U+i] += supg*foo[i];
f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i];
for (int k = 0; k < D; ++k)
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
//f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*fo1[i+k];
}
double lsic = tauc*rho;
f1[(U+i)*D+i] += ((c-cold)/dt +divu_mesh +divu)*lsic;
f1[(U+i)*D+i] += ((c-cold)/dt+divu_mesh +divu)*lsic;
//f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
for (int j = 0; j < D; ++j) {
f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
......@@ -437,10 +446,13 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
// PSPG
double pspg = taup/rho;
f1[P*D+i] += pspg*f0[U+i] + problem->stab_param*dp[i];
//f1[P*D+i] += pspg*fo[U+i] + problem->stab_param*dp[i];
f11[(P*D+i)*(n_fields*D)+P*D+i] += pspg*f01[(U+i)*n_fields*D+P*D+i]+problem->stab_param;
f10[(P*D+i)*n_fields+U+i] += pspg*f00[(U+i)*n_fields+(U+i)];
//f10[(P*D+i)*n_fields+U+i] += pspg*foo[i];
for (int j = 0; j < D; ++j) {
f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*f01[(U+i)*n_fields*D+(U+i)*D+j];
//f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*fo1[i+j];
}
}
}
......@@ -780,7 +792,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
du_mesh_vec[i*D+j] = 0;
}
}
double dc[D] = {0}, da[D]={0};
double dc[D] = {0}, da[D]={0},dcold[D]={0};
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
......@@ -799,6 +811,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
dcold[k] += problem->oldporosity[el[i]]*dphi[i][k];
}
if (problem->n_fluids==2) {
for (int k=0; k<D; ++k){
......@@ -839,9 +852,10 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
u_mesh[j] += dofmesh*phi[i];
u_mesh_old[j] += dofmesh_old*phi[i];
}
/*
if (problem->volume_old[el[i]]>1.e-30){
double rapport = problem->volume_new[el[i]]/problem->volume_old[el[i]];
double rapport = problem->volume_old[el[i]]/problem->volume_new[el[i]];
//c += (1-rapport*(1-problem->porosity[el[i]]))*phi[i];
cold += rapport*problem->oldporosity[el[i]]*phi[i];
}
......@@ -851,7 +865,6 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}*/
cold += problem->oldporosity[el[i]]*phi[i];
//printf("\n %.5f \n\n\n",cold);
c += problem->porosity[el[i]]*phi[i];
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
......@@ -868,7 +881,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
const double jw = QW[iqp]*detj;
problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,cold,u_mesh,u_mesh_old,du_mesh_vec,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,cold,dcold,u_mesh,u_mesh_old,du_mesh_vec,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -1154,7 +1167,30 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
problem->node_volume[el[i]] += vol;
}
}
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 0;
}
for (int i = 0; i < problem->n_particles; ++i) {
if(problem->particle_element_id[i] == -1) continue;
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*D],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j)
problem->porosity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
problem->porosity[i] = 0.;
}
else{
problem->porosity[i] = 1-problem->porosity[i]/problem->node_volume[i];
}
}
}
int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it)
{
......@@ -1175,11 +1211,12 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
volume_old[i] = problem->node_volume[i];
//fluid_problem_compute_node_volume(problem);
for (int i=0; i<mesh->n_nodes; ++i)
volume_new[i] = problem->node_volume[i];
problem->volume_new = volume_new;
problem->volume_old = volume_old;
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
solution_old[i] = problem->solution[i];
......@@ -1187,8 +1224,13 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
problem->solution_explicit = solution_old;
double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0;
int i;
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
if (-problem->mesh->x[3*i]/10+0.55 > 1) problem->porosity[i] = 1;
else problem->porosity[i] = -problem->mesh->x[3*i]/10+0.55;
}
compute_stab_parameters(problem,dt);
int i;
for (i=0; i<newton_max_it; ++i) {
double norm = 0;
hxtLinearSystemZeroMatrix(problem->linear_system);
......@@ -1218,8 +1260,6 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
//break; /// equation is linear
}
printf("total solve time : %g\n", totaltime);
free(dx);
free(rhs);
......@@ -1230,30 +1270,6 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
return i;
}
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 0;
}
for (int i = 0; i < problem->n_particles; ++i) {
if(problem->particle_element_id[i] == -1) continue;
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*D],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j)
problem->porosity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
problem->porosity[i] = 0.;
}
else{
problem->porosity[i] = 1-problem->porosity[i]/problem->node_volume[i];
}
}
}
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)
{
......@@ -1598,8 +1614,10 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
for(int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
if (-problem->mesh->x[3*i]/10+0.55 > 1){problem->porosity[i] = 1;problem->oldporosity[i] = 1;}
else {problem->porosity[i] = -problem->mesh->x[3*i]/10+0.55;problem->oldporosity[i] = -problem->mesh->x[3*i]/10+0.55;}
//problem->porosity[i] = 1.;
//problem->oldporosity[i] = 1.;
problem->velocity_mesh[i*D] = 0.;
problem->velocity_mesh[i*D+1] = 0.;
}
......
......@@ -70,9 +70,9 @@ struct FluidProblem {
double *concentration;
double *solution;
double *solution_explicit;
double *volume_new;
double *volume_old;
double *velocity_mesh_old;
double *volume_old;
double *volume_new;
double *velocity_old;
double *node_volume;
double *element_size;
......
......@@ -68,7 +68,7 @@ if not os.path.isdir(outputdir) :
# Physical parameters
g = -9.81 # gravity
r = 3e-3 # particles radius
r = 1e-3 # particles radius
rhop = 2000 # particles density
rho = 1000 # fluid density
nu = 1.e-6 # kinematic viscosity
......
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