Commit 1b10a024 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

equations with q

parent 010232f1
......@@ -241,7 +241,7 @@ static void f_u_0(
for (int k = 0; k < D; ++k) {
utau += u[k]*(du[j][k]-u[j]*dq[k]/q);
}
f[j] = rho*((u[j]-uold[j])/dt + 0*u[j]/q*divu + 0*utau/q);
f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q);
}
}
......@@ -251,7 +251,7 @@ static void f_p_0(double *f, const double du[D][D], const double q, const double
for (int j = 0; j < D; ++j) {
divu += du[j][j];
}
*f = divu+0*(q-qold)/dt;
*f = divu+(q-qold)/dt;
}
static void f_q_0(double *f, const double c, const double q)
......@@ -270,7 +270,7 @@ static void f_u_1(double f[D][D], double u[D], const double du[D][D], const doub
}
for (int j = 0; j < D; ++j) {
for (int k = 0; k < D; ++k) {
f[j][k] = mu*(tau[k][j]+tau[j][k]) + (k==j ? -p : 0);
f[j][k] = mu*(tau[k][j]+tau[j][k]) + (k==j ? -p : 0);
}
}
}
......@@ -289,7 +289,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
//a generaliser
double deps = 1e-12;
double deps = 1e-8;
const double *solution = problem->solution;
const double *porosity = problem->porosity;
......@@ -353,95 +353,106 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
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_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,p,mu);
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];
dq[i] += deps;
//k:equation; i:dphii; l:variable; j:dphij
double buf;
for (int l = 0; l < D; ++l){
buf = u[l];
u[l] += deps;
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;
u[l] = buf;
for (int k = 0; k < D; ++k){
fu0_u[k][l] = (fu0e[k]-fu0[k])/deps;
for (int i = 0; i < D; ++i){
fu1_u[k][i][l] = (fu1e[k][i]-fu1[k][i])/deps;
}
}
for (int j = 0; j < D; ++j){
buf = du[i][j];
du[i][j] += deps;
buf = du[l][j];
du[l][j] += deps;
f_p_0(&fp0e,du,q,qold,dt);
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;
du[l][j] = buf;
fp0_du[l][j] = (fp0e-fp0)/deps;
for (int k = 0; k < D; ++k){
fu0_du[k][i][j] = (fu0e[k]-fu0[k])/deps;
for (int l = 0; l < D; ++l){
fu1_du[k][i][j][l] = (fu1e[k][l]-fu1[k][l])/deps;
fu0_du[k][l][j] = (fu0e[k]-fu0[k])/deps;
for (int i = 0; i< D; ++i){
fu1_du[k][i][l][j] = (fu1e[k][i]-fu1[k][i])/deps;
}
}
}
}
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;
buf = p;
p += deps;
f_u_1(fu1e,u,du,q,dq,p,mu);
p = buf;
for (int k = 0; k < D; ++k){
for (int i = 0; i < D; ++i){
fu1_p[k][i] = (fu1e[k][i]-fu1[k][i])/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;
for (int j = 0; j < D; ++j){
buf = dp[j];
dp[j] += deps;
f_p_1(fp1e,dp,epsilon);
dp[j] = buf;
for (int i = 0; i < D; ++i){
fp1_dp[i][j] = (fp1e[i]-fp1[i])/deps;
}
}
buf = q;
q += deps;
f_q_0(&fq0e,c,q);
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);
q = buf;
fq0_q = (fq0e-fq0)/deps;
fp0_q = (fp0e-fp0)/deps;
for (int k = 0; k < D;++k){
fu0_q[k] = (fu0e[k]-fu0[k])/deps;
for (int i = 0; i < D; ++i){
fu1_q[k][i] = (fu1e[k][i]-fu1[k][i])/deps;
}
}
for (int j = 0; j < D; ++j){
buf = dq[j];
dq[j] += deps;
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
dq[j] = buf;
for (int k = 0; k < D; ++k){
fu0_dq[k][j] = (fu0e[k]-fu0[k])/deps;
for (int i = 0; i < D; ++i){
fu1_dq[k][i][j] = (fu1e[k][i]-fu1[k][j])/deps;
}
}
}
f_p_0(&fp0e,du,q+deps,qold,dt);
f_q_0(&fq0e,c,q+deps);
fp0_q = (fp0e- fp0)/deps;
fq0_q = (fq0e- fq0)/deps;
//filling
for (int iphi = 0; iphi < N_SF; ++iphi){
double phii = phi[iphi];
const double *dphii = dphi[iphi];
for (int j = 0; j < D; ++j){
double f = fu0[j]*phii;
for (int k = 0; k < D; ++k){
f += fu1[j][k]*dphii[k];
for (int k = 0; k < D; ++k){
double f = fu0[k]*phii;
for (int i = 0; i < D; ++i){
f += fu1[k][i]*dphii[i];
}
local_vector[iphi+N_SF*j] += jw*f;
local_vector[iphi+N_SF*k] += jw*f;
}
double ff = fp0*phii;
for (int j = 0; j < D; ++j){
ff += fp1[j]*dphii[j];
for (int i = 0; i < D; ++i){
ff += fp1[i]*dphii[i];
}
local_vector[iphi+N_SF*P] += jw*ff;
local_vector[iphi+N_SF*Q] += jw*fq0*phii;
//printf("iphi=%d,iel=%d: local_vectorU=%g,local_vectorV=%g,local_vectorP=%g,local_vectorQ=%g\n",iphi,iel,local_vector[iphi+N_SF*0],local_vector[iphi+N_SF*1],local_vector[iphi+N_SF*P],local_vector[iphi+N_SF*Q]);
for (int jphi = 0; jphi < N_SF; ++jphi){
double phij = phi[jphi];
const double *dphij = dphi[jphi];
......@@ -449,36 +460,55 @@ 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;
double mq_q = fq0_q*phij*phii;
for (int j = 0; j < D; ++j){
int U = j;
double mp_u = 0;
for (int k = 0; k < D; ++k){
int U = k;
double mu_p = 0;
double mu_q = fu0_q[j]*phii*phij;
for (int k = 0; k < D; ++k){
int V = k;
mu_p += fu1_p[j][k]*dphii[k]*phij;
mp_u += fp0_du[j][k]*phii*dphij[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 mu_u = fu0_u[j][k]*phii*phij;
double mu_q = fu0_q[k]*phii*phij;
for (int l = 0; l < D; ++l){
int V = l;
double mu_u = fu0_u[k][l]*phii*phij;
for (int i = 0; i < D; ++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){
mu_u += fu1_du[j][k][i][l]*dphij[i]*dphii[l];
mu_u += fu1_u[k][i][l]*dphii[i]*phij;
for (int j = 0; j < D; ++j){
mu_u += fu1_du[k][i][l][j]*dphii[i]*dphij[j];
}
}
for (int j = 0; j < D; ++j){
mu_u += fu0_du[k][l][j]*phii*dphij[j];
}
LOCAL_MATRIX(U,V) += jw*mu_u;
}
LOCAL_MATRIX(P,U) += jw*mp_u;
for (int i = 0; i < D; ++i){
for (int j = 0; j < D; ++j){
mu_q += fu1_dq[k][i][j]*dphii[i]*dphij[j];
}
mu_p += fu1_p[k][i]*phij*dphii[i];
mu_q += fu1_q[k][i]*phij*dphii[i];
}
for (int j = 0; j < D; ++j){
mu_q += fu0_dq[k][j]*phii*dphij[j];
}
LOCAL_MATRIX(U,P) += jw*mu_p;
LOCAL_MATRIX(U,Q) += jw*mu_q;
}
for (int l = 0; l < D; ++l){
int V = l;
double mp_u = 0;
for (int j = 0; j < D; ++j){
mp_u += fp0_du[l][j]*dphij[j]*phii;
}
LOCAL_MATRIX(P,V) += jw*mp_u;
}
double mp_p = 0;
for (int i = 0; i < D; ++i){
for (int j = 0; j < D; ++j){
mp_p += fp1_dp[i][j]*dphii[i]*dphij[j];
}
}
LOCAL_MATRIX(P,P) += jw*mp_p;
double mp_q = fp0_q*phii*phij;
LOCAL_MATRIX(P,Q) += jw*mp_q;
double mq_q = fq0_q*phij*phii;
LOCAL_MATRIX(Q,Q) += jw*mq_q;
}
}
......
L = 3;
L = 5;
H = 1;
y = 0;
lc = .1;
lc = .05;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, -H, 0, lc};
Point(3) = {L, -H, 0, lc};
Point(2) = {-L, 0, 0, lc};
Point(3) = {L, 0, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
......
......@@ -44,11 +44,11 @@ g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 1000 # final time
tEnd = 100000 # final time
#numerical parameters
lcmin = 0.1 # mesh size
dt = 5e-2 # time step
lcmin = .1 # mesh size
dt = 5 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......
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