Commit 89d19005 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

derivees de u par rapport a p et q + derivee de p par rapport a u et p

parent e71261a1
......@@ -245,6 +245,15 @@ static void f_u_0(
}
}
static void f_p_0(double *f, const double du[D][D], const double q, const double qold, const double dt)
{
double divu = 0;
for (int j = 0; j < D; ++j) {
divu += du[j][j];
}
*f = divu+(q-qold)/dt;
}
//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) {
......@@ -261,6 +270,15 @@ static void f_u_1(double f[D][D], double u[D], const double du[D][D], const doub
}
}
static void f_p_1(double f[D], const double dp[D], const double epsilon)
{
for (int j = 0; j < D; ++j) {
f[j] = epsilon*dp[j];
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
HXTLinearSystem *lsys = problem->linear_system;
......@@ -306,44 +324,67 @@ 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[D]={0}, dudt[D]={0}, p=0, c=0, dcdt=0,q=0;
double u[D]={0}, uold[D]={0}, p=0, c=0, qold=0, q=0;
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];
dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*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],fu1[D][D];
double fu0[D],fu0e[D],fu1[D][D],fu1e[D][D],*fp0,*fp0e,fp1[D],fp1e[D];
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);
double fu0_u[D][D],fu1_u[D][D][D],fu0_du[D][D][D],fu1_du[D][D][D][D];
//fu0_u[f][u],fu1_u[f][u][x],fu0_du[f][u][sfx],fu1_du[f][u][x][sfx]
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;
//fp0_du[u][x],fp1_dp[x][sfx]
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);
u[i] = buf;
buf = dp[i];
dp[i] += deps;
f_p_1(fp1e,dp,epsilon);
dp[i] = buf;
for (int j = 0; j < D; ++j){
fu0_u[j][i] = (fu0e[j]-fu0[j])/deps;
fp1_dp[i][j] = (fp1e[j]-fp1[j])/deps;
for (int k = 0; k < D; ++k){
fu1_u[j][i][k] = (fu1e[j][k]-fu1[j][k])/deps;
}
}
buf = dq[i];
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
dq[i] = buf;
for (int j = 0; j < D; ++j){
fu0_dq[j][i] = (fu0e[j]-fu0[j])/deps;
for (int k = 0; k < D; ++k){
fu1_dq[j][i][k] = (fu1e[j][k]-fu1[j][k])/deps;
}
}
for (int j = 0; j < D; ++j){
buf = du[i][j];
du[i][j] += deps;
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;
du[i][j] = buf;
for (int k = 0; k < D; ++k){
fu0_du[k][i][j] = (fu0e[k]-fu0[k])/deps;
......@@ -353,6 +394,24 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
}
f_u_1(fu1e,u,du,q,dq,p+deps,mu);
for (int i = 0; i < D; ++i){
for (int j = 0; j < D; ++j){
fu1_p[i][j] = (fu1e[i][j]-fu1[i][j])/deps;
}
}
f_u_0(fu0e,u,du,q+deps,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q+deps,dq,p,mu);
for (int i = 0; i < D; ++i){
fu0_q[i] = (fu0e[i]-fu0[i])/deps;
for (int j = 0; j < D; ++j){
fu1_q[i][j] = (fu1e[i][j]-fu1[i][j])/deps;
}
}
f_p_0(&fp0e,du,q+deps,qold,dt);
*fp0_q = (*fp0e- *fp0)/deps;
//remplissage
for (int iphi = 0; iphi < N_SF; ++iphi){
double phii = phi[iphi];
const double *dphii = dphi[iphi];
......@@ -363,6 +422,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
local_vector[iphi+N_SF*j] += jw*f;
}
double f = *fp0*phii;
for (int j = 0; j < D; ++j){
f += fp1[j]*dphii[j];
}
local_vector[iphi+N_SF*P] += jw*f;
for (int jphi = 0; jphi < N_SF; ++jphi){
double phij = phi[jphi];
......@@ -371,20 +436,35 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
int N2 = n_fields*N_SF;
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;
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;
for (int k = 0; k < D; ++k){
int V = k;
m_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];
mp_p += fp1_dp[k][j]*dphii[j]*dphij[k];
double m = 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];
for (int l = 0; l < D; ++l){
m += fu1_du[j][k][i][l]*dphij[i]*dphii[l];
}
}
LOCAL_MATRIX(U,V) = jw*m;
}
LOCAL_MATRIX(P,U) = jw*mp;
LOCAL_MATRIX(U,P) = jw*m_p;
LOCAL_MATRIX(U,Q) = jw*m_q;
}
LOCAL_MATRIX(P,P) = jw*mp_p;
LOCAL_MATRIX(P,Q) = jw*mp_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