Skip to content

GitLab

Commit c98a0a38 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

may not be broken anymore

parent 881f8934
Hide whitespace changes
Inline Side-by-side
src/fluid_problem.c
View file @ c98a0a38
......@@ -291,8 +291,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
const size_t local_size = N_SF*n_fields;
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
int forced = 0;
if (strcmp(mesh->phys_name[iphys],"Right") == 0)
forced = 1;
/*if (strcmp(mesh->phys_name[iphys],"Right") == 0)
forced = 1;*/
for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
const int iel = bnd[0];
......@@ -341,6 +341,47 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, double c, double dt) {
double fu0[D];
double fu1[D][D];
double fq0;
double fq1[D];
int P = problem->n_fluids*(D+1);
double p = sol[P];
double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
double epsilon = problem->epsilon;
f0[P] = (1-c);
f1[P*D+0] = 0;
f1[P*D+1] = 0;
for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
int Q = ifluid*(D+1)+D;
int U = ifluid*(D+1);
double rho = problem->rho[ifluid];
double mu = problem->mu[ifluid];
double q = sol[Q];
double qold = solold[Q];
double dq[D] = {dsol[Q*D+0],dsol[Q*D+1]};
double u[D] = {sol[U],sol[U+1]};
double uold[D] = {solold[U],solold[U+1]};
double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
f_u_0(fu0,u,du,q,dq,uold,dt,rho);
f_u_1(fu1,u,du,q,dq,p,mu);
f_q_0(&fq0,du,q,qold,dt);
f_q_1(fq1,dp,epsilon,q);
f0[U] = fu0[0];
f0[U+1] = fu0[1];
f0[Q] = fq0;
f0[P] -= q;
f1[U*D+0] = fu1[0][0];
f1[U*D+1] = fu1[0][1];
f1[U*D+D+0] = fu1[1][0];
f1[U*D+D+1] = fu1[1][1];
f1[Q*D+0] = fq1[0];
f1[Q*D+1] = fq1[1];
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
HXTLinearSystem *lsys = problem->linear_system;
......@@ -382,6 +423,15 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
//compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
double *f0 = malloc(sizeof(double)*n_fields);
double *f1 = malloc(sizeof(double)*n_fields*D);
double *g0 = malloc(sizeof(double)*n_fields);
double *g1 = malloc(sizeof(double)*n_fields*D);
double *h0 = malloc(sizeof(double)*n_fields*D);
double *h1 = malloc(sizeof(double)*n_fields*D*D);
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
......@@ -396,200 +446,68 @@ 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);
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
sold[i] = 0;
for (int j = 0; j < D; ++j) {
ds[i*D+j] = 0;
}
}
double c = 0;
for (int i = 0; i < N_SF; ++i) {
c += problem->compacity[el[i]]*phi[i];
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
double dofold = solution_old[el[i]*n_fields+j];
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
for (int k = 0; k < D; ++k)
ds[j*D+k] += dphi[i][k]*dof;
}
}
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi){
#define LOCAL_VECTOR(a) local_vector[iphi+N_SF*(a)]
double phii = phi[iphi];
LOCAL_VECTOR(P) += (1-c)*phii*jw;
}
for (int nf = 0; nf < problem->n_fluids; ++nf){
double dq[D]={0}, du[D][D]={{0}}, dp[D]={0};
double rho = problem->rho[nf];
double mu = problem->mu[nf];
int Q = (D+1)*nf+D;
int U = (D+1)*nf;
double u[D]={0}, uold[D]={0}, p=0, qold=0, q=0;
for (int i = 0; i< N_SF; ++i) {
const double *sol = &solution[el[i]*n_fields];
const double *sol_old = &solution_old[el[i]*n_fields];
for (int j = 0; j < D; ++j) {
dq[j] += dphi[i][j]*sol[Q];
dp[j] += dphi[i][j]*sol[P];
for (int k = 0; k < D; ++k)
du[j][k] += dphi[i][k]*sol[U+j];
u[j] += phi[i]*sol[U+j];
uold[j] += phi[i]*sol_old[U+j];
}
p += phi[i]*sol[P];
q += phi[i]*sol[Q];
qold += phi[i]*sol_old[Q];
}
double fu0[D],fu0e[D],fu1[D][D],fu1e[D][D],fq0,fq0e,fq1[D],fq1e[D];
f_u_0(fu0,u,du,q,dq,uold,dt,rho);
f_u_1(fu1,u,du,q,dq,p,mu);
f_q_0(&fq0,du,q,qold,dt);
f_q_1(fq1,dp,epsilon,q);
//discrete derivatives
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 fq0_du[D][D],fq1_dp[D][D],fq0_q,fq1_q[D];
//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_q_0(&fq0e,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;
fq0_du[l][j] = (fq0e-fq0)/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_q_1(fq1e,dp,epsilon,q);
dp[j] = buf;
for (int i = 0; i < D; ++i){
fq1_dp[i][j] = (fq1e[i]-fq1[i])/deps;
}
}
buf = q;
q += deps;
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
f_q_0(&fq0e,du,q,qold,dt);
f_q_1(fq1e,dp,epsilon,q);
q = buf;
fq0_q = (fq0e-fq0)/deps;
for (int i = 0; i < D; ++i){
fq1_q[i] = (fq1e[i] - fq1[i])/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;
fluid_problem_f(problem,f0,f1,s,ds,sold,c,dt);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
for (int id = 0; id < D; ++id) {
local_vector[iphi+N_SF*ifield] += dphi[iphi][id]*f1[ifield*D+id]*jw;
}
}
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;
}
}
}
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
fluid_problem_f(problem,g0,g1,s,ds,sold,c,dt);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dt);
ds[jfield*D+jd] = store;
}
//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 i = 0; i < D; ++i){
f += fu1[k][i]*dphii[i];
}
LOCAL_VECTOR(U+k) += jw*f;
}
double fq = fq0*phii;
for (int i = 0; i < D; ++i){
fq += fq1[i]*dphii[i];
}
LOCAL_VECTOR(Q) += jw*fq;
LOCAL_VECTOR(P) += -q*jw*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)*n_fields + jphi*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*(a)+(b)]
for (int k = 0; k < D; ++k){
double mu_p = 0;
double mu_q = fu0_q[k]*phii*phij;
for (int l = 0; l < D; ++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+k,U+l) += jw*mu_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];
int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
for (int ifield = 0; ifield < n_fields; ++ifield){
for (int iphi = 0; iphi < N_SF; ++iphi){
for (int jphi = 0; jphi < N_SF; ++jphi){
double df00 = (g0[ifield]-f0[ifield])/deps;
double m = jw*phi[iphi]*phi[jphi]*df00;
for (int id = 0; id < D; ++id) {
double df10 = (g1[ifield*D+id]-f1[ifield*D+id])/deps;
m += jw*dphi[iphi][id]*phi[jphi]*df10;
for (int jd = 0; jd < D; ++jd) {
double df11 = (h1[jd*n_fields*D+ifield*D+id]-f1[ifield*D+id])/deps;
m += jw*dphi[iphi][id]*dphi[jphi][jd]*df11;
}
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];
for (int jd = 0; jd < D; ++jd) {
double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
m += jw*phi[iphi]*dphi[jphi][jd]*df01;
}
LOCAL_MATRIX(U+k,P) += jw*mu_p;
LOCAL_MATRIX(U+k,Q) += jw*mu_q;
}
for (int l = 0; l < D; ++l){
double mq_u = 0;
for (int j = 0; j < D; ++j){
mq_u += fq0_du[l][j]*dphij[j]*phii;
}
LOCAL_MATRIX(Q,U+l) += jw*mq_u;
}
double mq_p = 0;
for (int i = 0; i < D; ++i){
for (int j = 0; j < D; ++j){
mq_p += fq1_dp[i][j]*dphii[i]*dphij[j];
}
}
LOCAL_MATRIX(Q,P) += jw*mq_p;
double mq_q = fq0_q*phii*phij;
for (int i = 0; i < D; ++i){
mq_q += fq1_q[i]*phij*dphii[i];
int IDX = (iphi*N2)*n_fields + jphi*n_fields;
LOCAL_MATRIX2(ifield,jfield) += m;
}
LOCAL_MATRIX(Q,Q) += jw*mq_q;
LOCAL_MATRIX(P,Q) += -jw*phii*phij;
}
}
}
......@@ -623,6 +541,15 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(s);
free(ds);
free(sold);
free(f0);
free(f1);
free(g0);
free(g1);
free(h0);
free(h1);
free(all_local_matrix);
free(all_local_vector);
free(forced_row);
......@@ -690,10 +617,6 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
break;
clock_get(&tic);
hxtLinearSystemSolve(problem->linear_system, rhs, dx);
for (int i = 0; i < mesh->n_nodes*n_fields; ++i)
printf("%g\n", dx[i]);
printf("\n");
exit(0);
clock_get(&toc);
//break;
totaltime += clock_diff(&tic, &toc);
......
src/fluid_problem_io.c
View file @ c98a0a38
......@@ -206,7 +206,7 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*n_fields+problem->n_fluids*(DIMENSION+1)+1]);
fprintf(f, "%g ", problem->solution[i*n_fields+problem->n_fluids*(DIMENSION+1)]);
}
fprintf(f,"\n</DataArray>\n");
for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
......
testcases/poiseuille2fluids/poiseuille.py
View file @ c98a0a38
......@@ -48,7 +48,7 @@ tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 5 # time step
dt = 1 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /min(nu0,nu1) # stabilization parametre
print('epsilon',epsilon)
......@@ -65,38 +65,27 @@ p = 6
strong_boundaries = [
("Bottom",u0,u0,0),
("Bottom",v0,v0,0.),
#("Bottom",q0,q0,0.5),
("Bottom",q0,q0,0.5),
("Bottom",u1,u1,0),
("Bottom",v1,v1,0.),
#("Bottom",q1,q1,0.5),
("Top",u0,u0,0),
("Top",v0,v0,0.),
#("Top",q0,q0,0.5),
("Top",q0,q0,0.5),
("Top",u1,u1,0),
("Top",v1,v1,0.),
#("Top",q1,q1,0.5),
("LeftDown",u0,u0,lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
("LeftDown",v0,v0,0),
("LeftDown",q0,q0,0.5),
("LeftDown",u1,u1,lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
("LeftDown",v1,v1,0),
("LeftDown",q1,q1,0.5),
("LeftUp",u0,u0,lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
("LeftUp",v0,v0,0),
("LeftUp",q0,q0,0.5),
("LeftUp",u1,u1,lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
("LeftUp",v1,v1,0),
("LeftUp",q1,q1,0.5),
("RightDown",v0,v0,0),
("RightDown",v1,v1,0),
#("RightDown",q0,q0,0.5),
#("RightDown",q1,q1,0.5),
("RightDown",p,p,0),
("RightUp",v0,v0,0),
("RightUp",v1,v1,0),
#("RightUp",q0,q0,0.5),
#("RightUp",q1,q1,0.5),
("RightUp",p,p,0)
("Left",u0,u0,lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
("Left",v0,v0,0),
("Left",q0,q0,0.5),
("Left",u1,u1,lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
("Left",v1,v1,0),
#("Left",q1,q1,0.5),
("Right",v0,v0,0),
("Right",v1,v1,0),
#("Right",q0,q0,0.5),
#("Right",q1,q1,0.5),
#("Right",p,p,0)
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,1,2)
ii = 0
......
Markdown is supported
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