Commit 3ba83798 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

début début début

parent 1b10a024
......@@ -30,7 +30,7 @@
#include "mesh_find.h"
#include "mbxml.h"
#define n_fields (D+2)
#define newton_max_it 10
#define newton_atol 5e-7
......@@ -290,6 +290,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
const Mesh *mesh = problem->mesh;
//a generaliser
double deps = 1e-8;
int n_fluids = problem->n_fluids;
int n_fields = problem->n_fields;
const double *solution = problem->solution;
const double *porosity = problem->porosity;
......@@ -317,204 +320,204 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
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 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) {
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];
}
}
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
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];
uold[j] += solution_old[el[i]*n_fields+j]*phi[i];
for (int nf = 0; nf < D; ++nf){
int P = (D+1)*n_fluids+1;
int Q = (D+1)*(nf+1);
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < D; ++j) {
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];
}
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];
}
const double jw = QW[iqp]*detj;
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[D]={0}, uold[D]={0}, p=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+nf*(D+1)]*phi[i];
uold[j] += solution_old[el[i]*n_fields+j+nf*(D+1)]*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];
}
const double jw = QW[iqp]*detj;
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[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_dq[D][D][D];
double fp0_du[D][D],fp1_dp[D][D],fp0_q,fq0_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];
double fp0_du[D][D],fp1_dp[D][D],fp0_q,fq0_q;
//discrete derivatives
//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);
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[l][j];
du[l][j] += deps;
f_p_0(&fp0e,du,q,qold,dt);
//discrete derivatives
//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);
du[l][j] = buf;
fp0_du[l][j] = (fp0e-fp0)/deps;
u[l] = buf;
for (int k = 0; k < D; ++k){
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;
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[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);
du[l][j] = buf;
fp0_du[l][j] = (fp0e-fp0)/deps;
for (int k = 0; k < D; ++k){
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;
}
}
}
}
}
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;
}
}
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);
buf = p;
p += deps;
f_u_1(fu1e,u,du,q,dq,p,mu);
dq[j] = buf;
p = 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;
fu1_p[k][i] = (fu1e[k][i]-fu1[k][i])/deps;
}
}
}
//filling
for (int iphi = 0; iphi < N_SF; ++iphi){
double phii = phi[iphi];
const double *dphii = dphi[iphi];
for (int k = 0; k < D; ++k){
double f = fu0[k]*phii;
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){
f += fu1[k][i]*dphii[i];
fp1_dp[i][j] = (fp1e[i]-fp1[i])/deps;
}
local_vector[iphi+N_SF*k] += jw*f;
}
double ff = fp0*phii;
for (int i = 0; i < D; ++i){
ff += fp1[i]*dphii[i];
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;
}
}
}
local_vector[iphi+N_SF*P] += jw*ff;
local_vector[iphi+N_SF*Q] += jw*fq0*phii;
for (int jphi = 0; jphi < N_SF; ++jphi){
double phij = phi[jphi];
const double *dphij = dphi[jphi];
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
//filling
for (int iphi = 0; iphi < N_SF; ++iphi){
double phii = phi[iphi];
const double *dphii = dphi[iphi];
for (int k = 0; k < D; ++k){
int U = k;
double mu_p = 0;
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;
double f = fu0[k]*phii;
for (int i = 0; i < D; ++i){
f += fu1[k][i]*dphii[i];
}
local_vector[iphi+N_SF*k] += jw*f;
}
double ff = fp0*phii;
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;
for (int jphi = 0; jphi < N_SF; ++jphi){
double phij = phi[jphi];
const double *dphij = dphi[jphi];
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
for (int k = 0; k < D; ++k){
int U = k;
double mu_p = 0;
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 += 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;
}
for (int i = 0; i < D; ++i){
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];
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_u += fu0_du[k][l][j]*phii*dphij[j];
mu_q += fu0_dq[k][j]*phii*dphij[j];
}
LOCAL_MATRIX(U,V) += jw*mu_u;
LOCAL_MATRIX(U,P) += jw*mu_p;
LOCAL_MATRIX(U,Q) += jw*mu_q;
}
for (int i = 0; i < D; ++i){
for (int l = 0; l < D; ++l){
int V = l;
double mp_u = 0;
for (int j = 0; j < D; ++j){
mu_q += fu1_dq[k][i][j]*dphii[i]*dphij[j];
mp_u += fp0_du[l][j]*dphij[j]*phii;
}
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(P,V) += jw*mp_u;
}
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];
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;
}
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;
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(all_local_matrix);
free(all_local_vector);
......@@ -533,6 +536,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
hxtLinearSystemSetRhsEntry(lsys, rhs,mesh->phys_nodes[iphys][i], bnd->field, 0.);
}
}
}
......@@ -660,7 +664,7 @@ static void fluid_problem_compute_porosity(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, int noGrains)
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int noGrains, int n_fluids)
{
static int initialized = 0;
if (!initialized) {
......@@ -675,9 +679,13 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
if (!mesh)
return NULL;
problem->mesh_tree = mesh_tree_create(mesh);
problem->rho = rho;
problem->mu = mu;
problem->n_fluids = n_fluids;
problem->n_fields = n_fluids*(D+1)+1;
for (int i = 0; i < n_fluids; ++i){
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
}
problem->g = g;
problem->mesh = mesh;
problem->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