Skip to content

GitLab

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

broken

parent b717baec
Hide whitespace changes
Inline Side-by-side
src/fluid_problem.c
View file @ dcffb101
......@@ -123,7 +123,7 @@ static double invDD(const double m[3][3], double inv[3][3]) {
#endif
static int problem_n_fields(const FluidProblem *p) {return p->n_fluids*(D+1)+1;}
int fluid_problem_n_fields(const FluidProblem *p) {return p->n_fluids*(D+1)+1;}
#if D==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
......@@ -237,6 +237,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
static void f_u_0(
double f[D], const double u[D], const double du[D][D], const double q, const double dq[D],
const double *uold, const double dt, const double rho) {
double oq = 1./q;//1./fmax(q,1e-3);
double divu = 0;
for (int j = 0; j < D; ++j) {
divu += du[j][j];
......@@ -244,9 +245,9 @@ static void f_u_0(
for (int j = 0; j < D; ++j) {
double utau = 0;
for (int k = 0; k < D; ++k) {
utau += u[k]*(du[j][k]-u[j]*dq[k]/q);
utau += u[k]*(du[j][k]-u[j]*dq[k]*oq);
}
f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q);
f[j] = rho*((u[j]-uold[j])/dt + u[j]*oq*divu + utau*oq);
}
}
......@@ -262,22 +263,23 @@ static void f_q_0(double *f, const double du[D][D], const double q, const double
//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) {
double tau[D][D];
double oq = 1./q;//1./fmax(q,1e-3);
for (int i = 0; i < D; ++i) {
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dq[j]/q;
tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
}
}
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 ? -q*p : 0);
}
}
}
static void f_q_1(double f[D], const double dp[D], const double epsilon)
static void f_q_1(double f[D], const double dp[D], const double epsilon, const double q)
{
for (int j = 0; j < D; ++j) {
f[j] = epsilon*dp[j];
f[j] = epsilon*dp[j]*q;
}
}
......@@ -285,7 +287,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
{
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const int n_fields = problem_n_fields(problem);
const int n_fields = fluid_problem_n_fields(problem);
const size_t local_size = N_SF*n_fields;
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
int forced = 0;
......@@ -296,7 +298,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
const int iel = bnd[0];
const int i0 = bnd[1];
const int i1 = (bnd[1]+1) % 3;
const int *el = &mesh->elements[iel*N_N];
const uint32_t *el = &mesh->elements[iel*N_N];
const int nodes[2] = {el[i0],el[i1]};
const double *x[2] = {&mesh->x[nodes[0]*3], &mesh->x[nodes[1]*3]};
const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
......@@ -305,7 +307,6 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double *local_vector = &all_local_vector[local_size*iel];
const int U = 0;
const int Q = 2;
const int P = 3;
double p0 = solution[el[i0]*n_fields+P];
double p1 = solution[el[i1]*n_fields+P];
......@@ -320,13 +321,21 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}*/
if (forced) {
double a = -1;
local_vector[i0+Q*N_SF] += l*(p0/3+p1/6)*a;
local_vector[i1+Q*N_SF] += l*(p0/6+p1/3)*a;
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += l/3*a;
local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += l/6*a;
local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += l/6*a;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += l/3*a;
for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
const int Q = ifluid*(D+1)+D;
double q0 = solution[el[i0]*n_fields+Q];
double q1 = solution[el[i1]*n_fields+Q];
double a = -1;
local_vector[i0+Q*N_SF] += l*(p0/3+p1/6)*a*q0;
local_vector[i1+Q*N_SF] += l*(p0/6+p1/3)*a*q1;
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += l/3*a*q0;
local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += l/6*a*q0;
local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += l/6*a*q1;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += l/3*a*q1;
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+Q)] += l*(p0/3+p1/6)*a;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+Q)] += l*(p0/3+p1/6)*a;
}
}
}
}
......@@ -339,7 +348,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double deps = 1e-8;
int n_fluids = problem->n_fluids;
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
const double *solution = problem->solution;
const double *compacity = problem->compacity;
......@@ -424,11 +433,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
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);
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;
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){
......@@ -471,7 +480,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int j = 0; j < D; ++j){
buf = dp[j];
dp[j] += deps;
f_q_1(fq1e,dp,epsilon);
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;
......@@ -482,8 +491,12 @@ 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_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){
......@@ -572,6 +585,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
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];
}
LOCAL_MATRIX(Q,Q) += jw*mq_q;
LOCAL_MATRIX(P,Q) += -jw*phii*phij;
}
......@@ -593,6 +609,17 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
/*for (int i = 0; i < 7*3; ++i) {
printf("%8.1e | ", local_vector[i]);
for (int j = 0; j < 7*3; ++j) {
double m = local_matrix[i*21+j];
if (fabs(m) < 1e-8) printf(" ");
else printf("%8.1e ",m);
}
printf("\n");
}
printf("\n");
printf("\n");*/
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
......@@ -633,7 +660,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
static double totaltime = 0;
struct timespec tic, toc;
const Mesh *mesh = problem->mesh;
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
apply_boundary_conditions(mesh, problem->solution, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
......@@ -742,7 +769,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->mesh_tree = mesh_tree_create(mesh);
problem->n_fluids = n_fluids;
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
problem->mu = (double*)malloc(sizeof(double)*n_fluids);
problem->rho = (double*)malloc(sizeof(double)*n_fluids);
for (int i = 0; i < n_fluids; ++i){
......@@ -773,8 +800,10 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->solution[i*n_fields+D] = 1.;
for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
for (int i = 0; i < mesh->n_nodes; ++i){
problem->solution[i*n_fields+D+ifluid*(D+1)] = 0.5;
}
}
problem->node_volume = NULL;
fluid_problem_compute_node_volume(problem);
......@@ -836,7 +865,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double gradM = 0.;
double gradPM = 0.;
double volT = 0.;
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
double *compacity = problem->compacity;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
......@@ -991,7 +1020,6 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
system("gmsh -" xstr(D) " adapt/mesh.geo");
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
{
struct timespec tic, toc;
......@@ -1000,7 +1028,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
clock_get(&toc);
printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
clock_get(&tic);
......@@ -1069,7 +1097,7 @@ void fluid_problem_after_import(FluidProblem *problem)
{
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
free(problem->old_compacity);
problem->old_compacity = (double*)malloc(sizeof(double)*problem->mesh->n_nodes);
for (int i = 0; i < problem->n_particles; ++i)
......@@ -1090,7 +1118,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
{
struct timespec tic, toc;
int n_fields = problem_n_fields(problem);
int n_fields = fluid_problem_n_fields(problem);
if(problem->n_particles != n) {
problem->n_particles = n;
......
src/fluid_problem.h
View file @ dcffb101
......@@ -80,4 +80,5 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem);
int fluid_problem_n_fields(const FluidProblem *p);
#endif
src/fluid_problem_io.c
View file @ dcffb101
......@@ -25,29 +25,6 @@
#include "mbxml.h"
#include <string.h>
#define n_fields (DIMENSION+1)
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.msh",output_dir, iter);
FILE *f = fopen(file_name, "w");
if (!f){
printf("Cannot open file \"%s\" for writing.\n", file_name);
return -1;
}
if(mesh_write_msh(problem->mesh, f))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "compacity", t, iter, problem->compacity, 1, 0))
return -1;
int comp[]={0,1,2};
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, DIMENSION))
return -1;
fclose(f);
return 0;
}
int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
{
......@@ -177,6 +154,7 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
const int n_fields = fluid_problem_n_fields(problem);
Mesh *mesh = problem->mesh;
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.vtu",output_dir, iter);
......@@ -220,7 +198,7 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
#endif
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Cells>\n");
fprintf(f,"<PointData Scalars=\"Porosity Pressure\" Vectors=\"Velocity\">\n");
fprintf(f,"<PointData Scalars=\"Compacity Pressure\" Vectors=\"Velocity\">\n");
fprintf(f,"<DataArray Name=\"Compacity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->compacity[i]);
......@@ -228,17 +206,24 @@ 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*(DIMENSION+2)+DIMENSION]);
fprintf(f, "%g ", problem->solution[i*n_fields+problem->n_fluids*(DIMENSION+1)+1]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int j = 0; j < DIMENSION; ++j)
fprintf(f, "%g ", problem->solution[i*(DIMENSION+2)+j]);
if(DIMENSION==2)
fprintf(f, "0 ");
for (int ifluid = 0; ifluid < problem->n_fluids; ++ifluid) {
fprintf(f,"<DataArray Name=\"Concentration%i\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n",ifluid);
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*n_fields+DIMENSION+ifluid*(DIMENSION+1)]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Velocity%i\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n",ifluid);
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int j = 0; j < DIMENSION; ++j)
fprintf(f, "%g ", problem->solution[i*n_fields+j+ifluid*(DIMENSION+1)]);
if(DIMENSION==2)
fprintf(f, "0 ");
}
fprintf(f,"\n</DataArray>\n");
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</PointData>\n");
fprintf(f,"</Piece>\n");
fprintf(f,"<FieldData>\n");
......
testcases/poiseuille2fluids/poiseuille.py
View file @ dcffb101
......@@ -42,15 +42,15 @@ ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1 # kinematic viscosity
mu = nu*rho # dynamic viscosity
nu0 = 1e-3 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 5 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
epsilon = alpha*lcmin**2 /min(nu0,nu1) # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
......@@ -59,8 +59,8 @@ outf = 10 # number of iterations between o
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("LeftUp",1,0),("RightUp",1,0),("LeftDown",1,0),("RightDown",1,0),("Top",0,0),("Bottom",0,0),("Right",6,0),("LeftUp",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1])),("LeftDown",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1])),("LeftUp",2,1),("RightUp",2,1),("Top",2,1),("LeftDown",5,1),("RightDown",5,1),("Bottom",5,1),("LeftUp",5,0),("RightUp",5,0),("Top",5,0),("LeftDown",2,0),("RightDown",2,0),("Bottom",2,0)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries,1)
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),("LeftUp",1,1,0),("RightUp",1,1,0),("LeftDown",1,1,0),("RightDown",1,1,0),("Top",0,0,0),("Bottom",0,0,0),("RightUp",6,6,0),("RightDown",6,6,0),("LeftUp",0,0,lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),("LeftDown",0,0,lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1])),("LeftUp",2,2,1),("RightUp",2,2,1),("Top",2,2,0.5),("LeftDown",5,5,0.5),("RightDown",5,5,0.5),("Bottom",5,5,0.5),("LeftUp",5,5,0.5),("RightUp",5,5,0.5),("Top",5,5,0.5),("LeftDown",2,2,0.5),("RightDown",2,2,0.5),("Bottom",2,2,0.5)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,1,2)
ii = 0
t = 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