Commit 697a2850 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

false n_fluid

parent 46fe7aa1
...@@ -186,11 +186,11 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double ...@@ -186,11 +186,11 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double beta = 0; double beta = 0;
for (int iphi=0; iphi < N_SF; ++iphi) { for (int iphi=0; iphi < N_SF; ++iphi) {
c += compacity[el[iphi]]*phi[iphi]; c += compacity[el[iphi]]*phi[iphi];
beta += si[el[iphi]*n_fields+DIMENSION+ff*(DIMENSION+1)]; beta += si[el[iphi]*n_fields+DIMENSION+1+ff*(DIMENSION+1)];
for (int j=0; j < DIMENSION; ++j) { for (int j=0; j < DIMENSION; ++j) {
vf[j] += si[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi]; vf[j] += si[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi];
vfe[j] += se[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi]; vfe[j] += se[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi];
gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+n_fields-1]; gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+n_fields-2];
} }
} }
double du[DIMENSION], due[DIMENSION]; double du[DIMENSION], due[DIMENSION];
...@@ -210,12 +210,12 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double ...@@ -210,12 +210,12 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass; problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
double *local_vector = &all_local_vector[N_SF*n_fields*problem->nFluids*iel+ff*N_SF*(DIMENSION+1)]; 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*problem->nFluids*problem->nFluids*iel+ff*N_SF*N_E*(DIMENSION+1)*(DIMENSION+1)]; double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
for (int iphi = 0; iphi < N_SF; ++iphi) { for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int d = 0; d < DIMENSION; ++d) 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+ff*(DIMENSION+1)] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi]; local_vector[iphi+N_SF*1+ff*(DIMENSION+1)] -= fcoeff*gamma*dt*g*phi[iphi];
} }
for (int iphi = 0; iphi < N_SF; ++iphi) { for (int iphi = 0; iphi < N_SF; ++iphi) {
...@@ -225,8 +225,8 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double ...@@ -225,8 +225,8 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b] #define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double f = fcoeff*phi[iphi]; double f = fcoeff*phi[iphi];
for (int d = 0; d < DIMENSION; ++d){ for (int d = 0; d < DIMENSION; ++d){
LOCAL_MATRIX(d,d) -= -beta*f/(1-c)*phi[jphi]*gamma; LOCAL_MATRIX(d+ff*(DIMENSION+1),d+ff*(DIMENSION+1)) -= -beta*f/(1-c)*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) -= -beta*f*gamma*dt/mass*vol*dphi[jphi][d]; LOCAL_MATRIX(d+ff*(DIMENSION+1),DIMENSION) -= -beta*f*gamma*dt/mass*vol*dphi[jphi][d];
} }
} }
} }
...@@ -259,8 +259,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co ...@@ -259,8 +259,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iel=0; iel < mesh->n_elements; ++iel) { for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N]; const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION]; double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
int P = n_fields-1; int P = DIMENSION;
double *local_vector = &all_local_vector[N_SF*n_fields*n_fluids*iel]; double *local_vector = &all_local_vector[N_SF*n_fields*iel];
double c=0,p=0, dp[DIMENSION]={0}; double c=0,p=0, dp[DIMENSION]={0};
for (int i = 0; i < DIMENSION; ++i) for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j) for (int j = 0; j < DIMENSION; ++j)
...@@ -279,12 +279,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co ...@@ -279,12 +279,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int j = 0; j < DIMENSION; ++j) for (int j = 0; j < DIMENSION; ++j)
dp[j] += dphi[i][j]*solution[el[i]*n_fields+P]; dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
} }
}printf("c=%g,p=%g\n",c,p); }//printf("c=%g,p=%g\n",c,p);
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel]; double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
for (int ff = 0; ff < n_fluids; ++ff) for (int ff = 0; ff < n_fluids; ++ff)
{ {
int Q = DIMENSION+ff*(DIMENSION+1); int Q = n_fields-1;//DIMENSION+ff*(DIMENSION+1);
double dq[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dqdq[DIMENSION]={0}; double dq[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dqdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) { for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) { for (int j = 0; j < DIMENSION; ++j) {
...@@ -460,7 +460,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt) ...@@ -460,7 +460,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
rhs[i] = 0; rhs[i] = 0;
fluid_problem_assemble_system(problem, rhs, solution_old, dt); fluid_problem_assemble_system(problem, rhs, solution_old, dt);
//exit(0); //exit(0);
if (i==1) break; // if (i==1) break;
hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm); hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
printf("iter %i norm = %g\n", i, norm); printf("iter %i norm = %g\n", i, norm);
if (norm < newton_atol) if (norm < newton_atol)
...@@ -476,10 +476,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt) ...@@ -476,10 +476,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
for (int j=0; j<mesh->n_nodes*n_fields; ++j) { for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j]; problem->solution[j] -= dx[j];
} }
for (int j=0; j<mesh->n_nodes*n_fields; ++j){ /* for (int j=0; j<mesh->n_nodes*n_fields; ++j){
printf("sol=%g\n",dx[j]); printf("sol=%g\n",dx[j]);
} }
exit(0); exit(0);*/
} }
printf("total solve time : %g\n", totaltime); printf("total solve time : %g\n", totaltime);
free(dx); free(dx);
...@@ -538,7 +538,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu ...@@ -538,7 +538,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
hxtInitializeLinearSystems(0, NULL); hxtInitializeLinearSystems(0, NULL);
initialized = 1; initialized = 1;
#ifdef HXT_HAVE_PETSC #ifdef HXT_HAVE_PETSC
hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid"); hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
#endif #endif
} }
FluidProblem *problem = malloc(sizeof(FluidProblem)); FluidProblem *problem = malloc(sizeof(FluidProblem));
...@@ -573,7 +573,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu ...@@ -573,7 +573,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
for (int i = 0; i < problem->mesh->n_nodes; ++i) for (int i = 0; i < problem->mesh->n_nodes; ++i)
{ {
problem->compacity[i] = 0.; problem->compacity[i] = 0.;
problem->compacity[i] = 0.;
} }
// begin to remove // begin to remove
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){ for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
...@@ -930,7 +929,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou ...@@ -930,7 +929,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
} }
if (init) { if (init) {
//initial fluid_0+grains //initial fluid_0+grains
int Q = DIMENSION+0*(DIMENSION+1); int Q = DIMENSION+1+0*(DIMENSION+1);
for (int i = 0; i < problem->mesh->n_nodes; ++i){ for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->solution[i*n_fields+Q] = 1-problem->compacity[i]; problem->solution[i*n_fields+Q] = 1-problem->compacity[i];
// printf("compacity=%g\n",problem->compacity[i]); // printf("compacity=%g\n",problem->compacity[i]);
......
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