Skip to content
GitLab
Commit 74300d54 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

provisoire

parent 02115f3b
Hide whitespace changes
Inline Side-by-side
src/fluid_problem.c
View file @ 74300d54
......@@ -154,7 +154,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
double *porosity = problem->porosity;
double *compacity = problem->compacity;
Mesh *mesh = problem->mesh;
int n_fields = problem->nFields;
for (int i = 0; i < problem->n_particles; ++i) {
......@@ -171,7 +171,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
shape_functions(xi,phi);
uint32_t *el=problem->mesh->elements+iel*N_SF;
double *se = problem->solution_explicit;
double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
for (int k=0; k<DIMENSION; ++k)
for (int j=0; j<DIMENSION; ++j)
......@@ -180,12 +180,16 @@ 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){
double c=0,vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
double beta = 0;
for (int iphi=0; iphi < N_SF; ++iphi) {
c += porosity[el[iphi]]*phi[iphi];
c += compacity[el[iphi]]*phi[iphi];
beta += si[el[iphi]*n_fields+DIMENSION+ff*(DIMENSION+1)];
for (int j=0; j < DIMENSION; ++j) {
vf[j] += si[el[iphi]*n_fields+j]*phi[iphi];
vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+DIMENSION];
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];
gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+n_fields-1];
}
}
double du[DIMENSION], due[DIMENSION];
......@@ -197,19 +201,19 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double rhop = mass/vol;
double drho = rhop -problem->rho;
double g = problem->g*drho/rhop;
double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
double gamma = particle_drag_coeff(due,problem->mu[ff],problem->rho[ff],vol,c);
double fcoeff = 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+iel*n_fields*N_SF;
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
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)];
for (int iphi = 0; iphi < N_SF; ++iphi) {
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*1] -= fcoeff*gamma*dt*g*phi[iphi];
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];
}
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -219,11 +223,12 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double f = fcoeff*phi[iphi];
for (int d = 0; d < DIMENSION; ++d){
LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
LOCAL_MATRIX(d,d) -= -beta*f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) -= -beta*f*gamma*dt/mass*vol*dphi[jphi][d];
}
}
}
}
}
}
......@@ -235,11 +240,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
int n_fields = problem->nFields;
const double *compacity = problem->compacity;
const double *old_compacity = problem->old_compacity;
const double mu = problem->mu;
const double rho = problem->rho;
const double *mu = problem->mu;
const double *rho = problem->rho;
const double epsilon = problem->epsilon;
int n_fluids = problem->nFluids;
size_t local_size = N_SF*n_fields*problem->nFluids;
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));
......@@ -248,25 +253,38 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
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];
for (int ff = 0; ff < n_fluids; ++ff)
{
int P = n_fields;
int Q = DIMENSION+ff*(DIMENSION+1);
for (int i = 0; i < DIMENSION; ++i)
int P = n_fields;
double *local_vector = &all_local_vector[N_SF*n_fields*n_fluids*iel];
double c=0,p=0, dp[DIMENSION]={0};
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);
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)];
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0}, dcdq[DIMENSION]={0};
for (int iqp = 0; iqp< N_QP; ++iqp) {
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 *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel];
for (int ff = 0; ff < n_fluids; ++ff)
{
int Q = DIMENSION+ff*(DIMENSION+1);
double dq[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dqdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
dc[j] += dphi[i][j]*solution[el[i]*n_fields+P];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+Q];
dq[j] += dphi[i][j]*solution[el[i]*n_fields+Q];
dcdq[j] += dphi[i][j];
for (int k = 0; k < DIMENSION; ++k)
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k+ff*(DIMENSION+1)];
......@@ -275,19 +293,16 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0,q=0;
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, dqdt=0,q=0;
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];
}
p += solution[el[i]*n_fields+P]*phi[i];
q += solution[el[i]*n_fields+Q]*phi[i];
c += compacity[el[i]]*phi[i];
dcdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])/dt*phi[i];
dqdt += (q-solution_old[el[i]*n_fields+Q]*phi[i])/dt;
//dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
}
/////////////
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi) {
double phii = phi[iphi];
......@@ -307,17 +322,17 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int k = 0; k < DIMENSION; ++k) {
//utau += u[k]*dphii[k]*u[j]/c;
utau += u[k]*tau[j][k];
dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
dphisigma += 2*mu[ff]*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
}
local_vector[iphi+N_SF*j] += jw*(
dphisigma
//+rho*phii*dudt[j]-rho*utau
+rho*phii*(dudt[j]+u[j]/q*divu+utau/q)
+rho[ff]*phii*(dudt[j]+u[j]/q*divu+utau/q)
-p*dphii[j]
);
}
local_vector[iphi+N_SF*P] += jw*(epsilon*dphidp+phii*(divu+dcdt));
local_vector[iphi+N_SF*Q] += jw*phii*(c-q);
local_vector[iphi+N_SF*P] += jw*phii*(q);
local_vector[iphi+N_SF*Q] += jw*phii*(-dqdt+divu)+jw*epsilon*dphidp;
for (int jphi = 0; jphi < N_SF; ++jphi) {
double phij = phi[jphi];
const double *dphij = dphi[jphi];
......@@ -348,16 +363,16 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
int V = k;
utau += u[k]*tau[j][k];
dutaudq += u[k]*u[j]*(dc[k]*phij/(q*q)-dcdq[k]/q);
dphisigmadq += 2*mu*dphii[k]*0.5*(u[k]*dc[j]*phij/(q*q)+u[j]*dc[k]*phij/(q*q)-u[k]*dcdq[j]/q-u[j]*dcdq[k]/q);
dphisigmadq += 2*mu[ff]*dphii[k]*0.5*(u[k]*dc[j]*phij/(q*q)+u[j]*dc[k]*phij/(q*q)-u[k]*dcdq[j]/q-u[j]*dcdq[k]/q);
//LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
LOCAL_MATRIX(U,V) +=
jw*rho*phii*phij*tau[j][k]/q
+jw*2*mu*dphii[k]*dtau[j]*0.5
+jw*rho*phii*u[j]*dphij[k]/q;
jw*rho[ff]*phii*phij*tau[j][k]/q
+jw*2*mu[ff]*dphii[k]*dtau[j]*0.5
+jw*rho[ff]*phii*u[j]*dphij[k]/q;
}
//LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt+phij/q*divu+udtau/q);
LOCAL_MATRIX(U,Q) += jw*rho*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
LOCAL_MATRIX(U,U) += jw*mu[ff]*2*0.5*dphiidtau + jw*rho[ff]*phii*(phij/dt+phij/q*divu+udtau/q);
LOCAL_MATRIX(U,Q) += jw*rho[ff]*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
}
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(Q,Q) += -jw*phii*phij;
......
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