Skip to content

GitLab

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

debug

parent fa7ec368
Hide whitespace changes
Inline Side-by-side
src/fluid_problem.c
View file @ d788dfc7
......@@ -157,6 +157,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double *compacity = problem->compacity;
Mesh *mesh = problem->mesh;
int n_fields = problem->nFields;
int N_E = problem->mesh->n_elements;
for (int i = 0; i < problem->n_particles; ++i) {
int iel = problem->particle_element_id[i];
double mass = problem->particle_mass[i];
......@@ -180,7 +181,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double dphi[N_SF][DIMENSION];
grad_shape_functions(dxidx, dphi);
double *si = problem->solution;
for (int ff=0;ff<problem->n_fluids;++ff){
for (int ff=0;ff<problem->nFluids;++ff){
double c=0,vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
double beta = 0;
for (int iphi=0; iphi < N_SF; ++iphi) {
......@@ -199,21 +200,22 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
}
double vol = problem->particle_volume[i];
double rhop = mass/vol;
double drho = rhop -problem->rho;
double drho = rhop -problem->rho[ff];
double g = problem->g*drho/rhop;
double gamma = particle_drag_coeff(due,problem->mu[ff],problem->rho[ff],vol,c);
double fcoeff = mass/(mass+dt*gamma);
double fcoeff = beta*mass/(mass+dt*gamma);
for (int d = 0; d < DIMENSION; ++d)
problem->particle_force[i*DIMENSION+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
double *local_vector = &all_local_vector[N_SF*n_fields*n_fluids*iel+ff*N_SF*(DIMENSION+1)];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel+ff*N_SF*N_E*(DIMENSION+1)*(DIMENSION+1)];
////pas bon
double *local_vector = &all_local_vector[N_SF*n_fields*problem->nFluids*iel+ff*N_SF*(DIMENSION+1)];
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)];
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int d = 0; d < DIMENSION; ++d)
local_vector[iphi+N_SF*d] -= beta*fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*1] -= beta*fcoeff*gamma*dt*g*phi[iphi];
local_vector[iphi+N_SF*d] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi];
}
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -263,17 +265,19 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
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);
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
const double jwl = QW[iqp]*detj;
double phi[N_SF];
for (int i = 0; i < N_SF; ++i) {
phii = phi[i];
p += solution[el[i]*n_fields+P]*phi[i];
c += compacity[el[i]]*phi[i];
local_vector[iphi+N_SF*P] += jw*phii*(c-1);
for (int j = 0; j < DIMENSION; ++j)
dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
double phii = phi[i];
p += solution[el[i]*n_fields+P]*phi[i];
c += compacity[el[i]]*phi[i];
local_vector[i+N_SF*P] += jwl*phii*(c-1);
for (int j = 0; j < DIMENSION; ++j)
dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
}
}
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel];
......@@ -297,7 +301,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j){
u[j] += solution[el[i]*n_fields+j+ff*(DIMENSION+1)]*phi[i];
dudt[j] += (solution[el[i]*n_fields+j+ff*(DIMENSION+1]-solution_old[el[i]*n_fields+j+ff*(DIMENSION+1])/dt*phi[i];
dudt[j] += (solution[el[i]*n_fields+j+ff*(DIMENSION+1)]-solution_old[el[i]*n_fields+j+ff*(DIMENSION+1)])/dt*phi[i];
}
q += solution[el[i]*n_fields+Q]*phi[i];
dqdt += (q-solution_old[el[i]*n_fields+Q]*phi[i])/dt;
......@@ -310,7 +314,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double tau[DIMENSION][DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
tau[i][j] = du[i][j]-u[i]*dc[j]/q;
tau[i][j] = du[i][j]-u[i]*dq[j]/q;
double dphidp = 0, divu = 0;
for (int k = 0; k < DIMENSION; ++k){
dphidp += dphii[k]*dp[k];
......@@ -556,8 +560,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i = 0; i < problem->mesh->n_nodes; ++i)
{
problem->porosity[i] = 1.;
problem->old_porosity[i] = 1.;
problem->compacity[i] = 0.;
problem->compacity[i] = 0.;
}
// begin to remove
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
......@@ -621,7 +625,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double gradM = 0.;
double gradPM = 0.;
double volT = 0.;
int n_fields = nFields;
int n_fields = problem->nFields;
double *compacity = problem->compacity;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
......@@ -683,7 +687,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
const unsigned int *el = &mesh->elements[iel*N_N];
double C[N_N],U[DIMENSION][N_N], P[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
C[i] = compacity[el[i]];
P[i] = solution[el[i]*n_fields+DIMENSION];
for (int j=0; j<DIMENSION; ++j) {
U[j][i] = solution[el[i]*n_fields+j];
......@@ -854,8 +858,8 @@ void fluid_problem_after_import(FluidProblem *problem)
{
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
free(problem->old_porosity);
problem->old_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
free(problem->old_compacity);
problem->old_compacity = malloc(sizeof(double)*problem->mesh->n_nodes);
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
......@@ -863,7 +867,7 @@ void fluid_problem_after_import(FluidProblem *problem)
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,problem->nFields,problem->mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
#endif
......@@ -900,7 +904,6 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_element_id[i]=elid[i];
}
}
clock_get(&tic);
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
......@@ -918,7 +921,6 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
for (int i = 0; i < problem->mesh->n_nodes; ++i)
problem->solution[i*n_fields+Q] = 1-problem->compacity[i];
}
clock_get(&tic);
fluid_problem_compute_porosity(problem);
}
......
src/fluid_problem.h
View file @ d788dfc7
......@@ -77,7 +77,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_fluids, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
......
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