Commit b853fd25 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

adapt

parent 334274bb
......@@ -193,12 +193,16 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
}
double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
double rhop = mass/vol;
double drho = rhop-rho;
double g = problem->g;
double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
double gamma = gamma0*a + gamma1*(1-a);
double gamma;
if(problem->tf){
double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
gamma = gamma0*a + gamma1*(1-a);
}
else{
gamma = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
}
double fcoeff = mass/(mass+dt*gamma);
for (int j = 0; j < D; ++j){
......@@ -351,7 +355,7 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a+problem->rho[1]*(1-a);
mu = problem->mu[0]*a+problem->mu[1]*(1-a);
f[A] = un*(un > 0 ? a : 0.);
f[A] = un*(un > 0 ? a : 1.);
}
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
......@@ -477,7 +481,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
problem->taup[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
problem->taua[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h));
problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*0.5;
problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*1e-1;
maxtaup = fmax(maxtaup, problem->taup[iel]);
mintaup = fmin(mintaup, problem->taup[iel]);
}
......@@ -1105,20 +1109,20 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
num_lc[i] = 0.;
}
double gradmax = 0.;
double gradPmax = 0.;
double gradmin = 1e8;
double gradPmin = 1e8;
double gradM = 0.;
double gradPM = 0.;
double gradAmax = 0.;
double gradAmin = 1e8;
double gradAM = 0.;
double volT = 0.;
int n_fields = fluid_problem_n_fields(problem);
double *porosity = problem->porosity;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double C[N_N],Un[D][N_N], Pn[N_N];
double C[N_N],Un[D][N_N],An[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
Pn[i] = solution[el[i]*n_fields+D];
if(problem->tf) An[i] = solution[el[i]*n_fields+D+1];
for (int j=0; j<D; ++j) {
Un[j][i] = solution[el[i]*n_fields+j];
}
......@@ -1149,32 +1153,38 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
volT += vol;
gradmax = fmax(gradmax,ngrad2);
gradmin = fmin(gradmin,ngrad2);
ngrad2 = 0;
for (int j=0; j<D; ++j){
double gradP = 0;
for (int k=0; k<N_SF; ++k){
gradP += dphi[k][j]*Pn[k];
if(problem->tf){
double ngradA = 0;
for (int j=0; j<D; ++j){
double gradA = 0;
for (int k=0; k<N_SF; ++k){
gradA += dphi[k][j]*An[k];
}
ngradA += gradA*gradA;
}
ngrad2 += gradP*gradP;
gradAM += sqrt(ngradA)*vol;
volT += vol;
gradAmax = fmax(gradAmax,ngradA);
gradAmin = fmin(gradAmin,ngradA);
}
gradPM += sqrt(ngrad2)*vol;
gradPmax = fmax(gradPmax,ngrad2);
gradPmin = fmin(gradPmin,ngrad2);
}
gradPM /= volT;
gradM /= volT;
gradPmin = sqrt(gradPmin)/gradPM;
gradPmax = sqrt(gradPmax)/gradPM;
gradmin = sqrt(gradmin)/gradM;
gradmax = sqrt(gradmax)/gradM;
printf("gradPmax=%e\n",gradPmax);
printf("gradPmin=%e\n",gradPmin);
printf("gradmax=%e\n",gradmax);
printf("gradmin=%e\n",gradmin);
gradAM /= volT;
gradAmin = sqrt(gradAmin)/gradAM;
gradAmax = sqrt(gradAmax)/gradAM;
printf("gradAmax=%e\n",gradAmax);
printf("gradAmin=%e\n",gradAmin);
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double C[N_N],Un[D][N_N], Pn[N_N];
double C[N_N],Un[D][N_N],An[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
Pn[i] = solution[el[i]*n_fields+D];
for (int j=0; j<D; ++j) {
Un[j][i] = solution[el[i]*n_fields+j];
}
......@@ -1196,31 +1206,25 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
ngrad2 += gradU*gradU;
}
}
/* double ngrad = pow(ngrad2, 0.25); */
/* for (int j=0; j<D; ++j){*/
/* double gradP = 0;*/
/* for (int k=0; k<N_SF; ++k){*/
/* gradP += dphi[k][j]*P[k];*/
/* }*/
/* ngrad2 += gradP*gradP;*/
/* }*/
double ngrad = sqrt(ngrad2)/gradM;
double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
ngrad2 = 0.;
for (int j=0; j<D; ++j){
double gradP = 0;
for (int k=0; k<N_SF; ++k){
gradP += dphi[k][j]*Pn[k];
if(problem->tf){
double ngradA = 0;
for (int j=0; j<D; ++j){
double gradA = 0;
for (int k=0; k<N_SF; ++k){
gradA += dphi[k][j]*An[k];
}
ngradA += gradA*gradA;
}
ngrad2 += gradP*gradP;
double ngrad = sqrt(ngradA/gradAM);
double lc2 = lcmin /fmin(1,fmax(ngrad/(gradAmax),gradAmin/(gradAmax)));
lc = fmin(lc,lc2);
}
ngrad = sqrt(ngrad2)/gradPM;
lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
for (int j = 0; j < N_N; ++j){
new_mesh_size[el[j]] += lc;
new_mesh_size[el[j]] += fmax(fmin(lc,lcmax),lcmin);
num_lc[el[j]] += 1.;
}
}
......@@ -1277,7 +1281,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
double *new_oldporosity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
clock_get(&tic);
int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
......@@ -1308,10 +1311,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
for (int k = 0; k < N_SF; ++k)
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
}
new_oldporosity[i] = 0;
for (int k=0; k < N_SF; ++k){
new_oldporosity[i] += problem->oldporosity[el[k]]*phi[k];
}
}
free(new_eid);
free(new_xi);
......@@ -1329,10 +1328,11 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(problem->porosity);
problem->porosity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->oldporosity);
problem->oldporosity = new_oldporosity;
problem->oldporosity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
{
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
mesh_free(problem->mesh);
problem->mesh = new_mesh;
......
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