Commit 056803bb authored by Matthieu Constant's avatar Matthieu Constant
Browse files

ajout derniere equation

parent 89d19005
......@@ -254,6 +254,11 @@ static void f_p_0(double *f, const double du[D][D], const double q, const double
*f = divu+(q-qold)/dt;
}
static void f_q_0(double *f, const double c, const double q)
{
*f = c-q;
}
//flux
static void f_u_1(double f[D][D], double u[D], const double du[D][D], const double q, const double dq[D], const double p, const double mu) {
......@@ -310,17 +315,18 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
grad_shape_functions(dxidx, dphi);
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*iel];
double dc[D]={0}, du[D][D]={{0}}, dp[D]={0};
double dq[D]={0}, du[D][D]={{0}}, dp[D]={0};
int P = D;
int Q = D+1;
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < D; ++j) {
dc[j] += dphi[i][j]*porosity[el[i]];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+D];
dq[j] += dphi[i][j]*solution[el[i]*n_fields+Q];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
for (int k = 0; k < D; ++k)
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
}
}
int P = D;
int Q = D+1;
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
......@@ -329,32 +335,33 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < D; ++j){
u[j] += solution[el[i]*n_fields+j]*phi[i];
uold[j] += solution_old[el[i]*n_fields+j])*phi[i];
uold[j] += solution_old[el[i]*n_fields+j]*phi[i];
}
p += solution[el[i]*n_fields+P]*phi[i];
q += solution[el[i]*n_fields+Q]*phi[i];
qold += solution_old[el[i]*n_fields+Q]*phi[i];
c += porosity[el[i]]*phi[i];
dcdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])/dt*phi[i];
//dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
}
const double jw = QW[iqp]*detj;
double fu0[D],fu0e[D],fu1[D][D],fu1e[D][D],*fp0,*fp0e,fp1[D],fp1e[D];
double fu0[D],fu0e[D],fu1[D][D],fu1e[D][D],fp0,fp0e,fp1[D],fp1e[D],fq0,fq0e;
f_u_0(fu0,u,du,q,dq,uold,dt,rho);
f_u_1(fu1,u,du,q,dq,p,mu);
f_p_0(&fp0,du,q,qold,dt);
f_p_1(fp1,dp,epsilon);
f_q_0(&fq0,c,q);
double fu0_u[D][D],fu1_u[D][D][D],fu0_du[D][D][D],fu1_du[D][D][D][D],fu1_p[D][D],fu0_q[D],fu1_q[D][D],fu0_dq[D][D],fu1[D][D][D];
//fu0_u[f][u],fu1_u[f][u][x],fu0_du[f][u][sfx],fu1_du[f][u][x][sfx],fu1_p[f][sfx],fu0_q[f],fu1_q[f][sfx],fu0_dq[f][x],fu1[f][x][sfx]
double fp0_du[D][D],fp1_dp[D][D],*fp0_q;
double fu0_u[D][D],fu1_u[D][D][D],fu0_du[D][D][D],fu1_du[D][D][D][D],fu1_p[D][D],fu0_q[D],fu1_q[D][D],fu0_dq[D][D],fu1_dq[D][D][D];
//fu0_u[f][u],fu1_u[f][u][x],fu0_du[f][u][sfx],fu1_du[f][u][x][sfx],fu1_p[f][sfx],fu0_q[f],fu1_q[f][sfx],fu0_dq[f][x],fu1_dq[f][x][sfx]
double fp0_du[D][D],fp1_dp[D][D],fp0_q,fq0_q;
//fp0_du[u][x],fp1_dp[x][sfx]
//discrete derivatives
for (int i = 0; i < D; ++i){
double buf = u[i];
u[i] += deps;
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
u[i] = buf;
buf = dp[i];
......@@ -384,7 +391,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
f_p_0(&fp0e,du,q,qold,dt);
fp0_du[i][j] = (*fp0e-*fp0)/deps;
fp0_du[i][j] = (fp0e-fp0)/deps;
du[i][j] = buf;
for (int k = 0; k < D; ++k){
fu0_du[k][i][j] = (fu0e[k]-fu0[k])/deps;
......@@ -409,9 +416,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
f_p_0(&fp0e,du,q+deps,qold,dt);
*fp0_q = (*fp0e- *fp0)/deps;
f_q_0(&fq0e,c,q+deps);
fp0_q = (fp0e- fp0)/deps;
fq0_q = (fq0e- fq0)/deps;
//remplissage
//filling
for (int iphi = 0; iphi < N_SF; ++iphi){
double phii = phi[iphi];
const double *dphii = dphi[iphi];
......@@ -422,11 +431,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
local_vector[iphi+N_SF*j] += jw*f;
}
double f = *fp0*phii;
double f = fp0*phii;
for (int j = 0; j < D; ++j){
f += fp1[j]*dphii[j];
}
local_vector[iphi+N_SF*P] += jw*f;
local_vector[iphi+N_SF*Q] += jw*fq0*phii;
for (int jphi = 0; jphi < N_SF; ++jphi){
......@@ -437,34 +447,36 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double mp_p = 0;
double mp_q = *fp0_q*phij*phii;
double mp_q = fp0_q*phij*phii;
double mq_q = fq0_q*phij*phii;
for (int j = 0; j < D; ++j){
int U = j;
double mp_u = 0;
double m_p = 0;
double m_q = fu0_q[j]*phii*phij;
double mu_p = 0;
double mu_q = fu0_q[j]*phii*phij;
for (int k = 0; k < D; ++k){
int V = k;
m_p += fu1_p[j][k]*dphii[k]*phij;
mu_p += fu1_p[j][k]*dphii[k]*phij;
mp_u += fp0_du[j][k]*phii*dphij[k];
m_q += fu0_dq[j][k]*phii*dphij[k]+fu1_q[j][k]*phij*dphii[k];
mu_q += fu0_dq[j][k]*phii*dphij[k]+fu1_q[j][k]*phij*dphii[k];
mp_p += fp1_dp[k][j]*dphii[j]*dphij[k];
double m = fu0_u[j][k]*phii*phij;
double mu_u = fu0_u[j][k]*phii*phij;
for (int i = 0; i < D; ++i){
m += fu0_du[j][k][i]*dphij[i]*phii+fu1_u[j][k][i]*phij*dphii[i];
m_q += fu1_dq[j][k][i]*dphij[k]*dphii[i];
mu_u += fu0_du[j][k][i]*dphij[i]*phii+fu1_u[j][k][i]*phij*dphii[i];
mu_q += fu1_dq[j][k][i]*dphij[k]*dphii[i];
for (int l = 0; l < D; ++l){
m += fu1_du[j][k][i][l]*dphij[i]*dphii[l];
mu_u+= fu1_du[j][k][i][l]*dphij[i]*dphii[l];
}
}
LOCAL_MATRIX(U,V) = jw*m;
LOCAL_MATRIX(U,V) = jw*mu_u;
}
LOCAL_MATRIX(P,U) = jw*mp;
LOCAL_MATRIX(U,P) = jw*m_p;
LOCAL_MATRIX(U,Q) = jw*m_q;
LOCAL_MATRIX(P,U) = jw*mp_u;
LOCAL_MATRIX(U,P) = jw*mu_p;
LOCAL_MATRIX(U,Q) = jw*mu_q;
}
LOCAL_MATRIX(P,P) = jw*mp_p;
LOCAL_MATRIX(P,Q) = jw*mp_q;
LOCAL_MATRIX(Q,Q) = jw*mq_q;
}
}
}
......
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