Commit 7bdeb284 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

rewrite weak bnd

parent eae35a08
......@@ -44,8 +44,8 @@ SCONTACT_H=scontact/quadtree.h scontact/scontact.h src/mbxml.h src/yxml.h
libmbfluid2.so : ${FLUID_C} ${FLUID_H}
${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=2
libmbfluid3.so : ${FLUID_C} ${FLUID_H}
${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=3
#libmbfluid3.so : ${FLUID_C} ${FLUID_H}
# ${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=3
libscontact2.so : ${SCONTACT_C} ${SCONTACT_H}
${CC} ${SCONTACT_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=2
......
......@@ -37,12 +37,16 @@
#define newton_rtol 1e-5
#define D DIMENSION
#define deps 1e-8
#if D==2
#define N_SF 3
#define N_N 3
#define N_QP 3
#define N_LQP 2
static double LQP[] = {-5.773502691896257e-01, 5.773502691896257e-01};
static double LQW[] = {1,1};
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
static void shape_functions(double *uvw, double *sf)
......@@ -233,86 +237,31 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
#endif
}
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
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;
/*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];
const int i0 = bnd[1];
const int i1 = (bnd[1]+1) % 3;
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]};
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
for (int i = 0; i < D; ++i)
for (int j = 0; j < D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double dp[D]={0};
typedef void f_cb(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double dt);
const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {-dx[1]/l,dx[0]/l};
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;
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < D; ++j) {
dp[j] += dphi[i][j]*solution[el[i]*n_fields+3];
}
}
double p0 = solution[el[i0]*n_fields+P];
double p1 = solution[el[i1]*n_fields+P];
for (int k = 0; k < D; ++k) {
if (!forced) {
local_matrix[(i0*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/3*n[k];
local_matrix[(i0*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/3*n[k];
local_vector[i0+(U+k)*N_SF] += l*(p0/3+p1/6)*n[k];
local_vector[i1+(U+k)*N_SF] += l*(p0/6+p1/3)*n[k];
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
local_vector[i0+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
local_vector[i1+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
}
}
if (forced) {
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;
}
}
static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt)
{
int n_fluids = problem->n_fluids;
int P = n_fluids*(D+1);
const double epsilon = problem->epsilon;
double dp[D];
for (int id = 0; id < D; ++id) {
dp[id] = ds[P*D+id];
}
double p = s[P];
for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
int U = ifluid*(D+1);
int Q = ifluid*(D+1)+D;
double q = s[Q];
f[Q] = 0;
for (int id = 0; id < D; ++id) {
f[U+id] = q*p*n[id];
}
}
f[P] = 0;
for (int id = 0; id < D; ++id) {
f[P] += dp[id]*n[id]*epsilon;
}
}
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, double c, double dt) {
......@@ -322,8 +271,9 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double epsilonp = problem->epsilon;
double epsilonq = problem->epsilon/10;
f0[P] = 1-c;
f1[P*D+0] = 0;//epsilonp*dp[0];
f1[P*D+1] = 0;//epsilonp*dp[1];
for (int id = 0; id < D; ++id) {
f1[P*D+id] = epsilonp*dp[id];
}
for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
int Q = ifluid*(D+1)+D;
int U = ifluid*(D+1);
......@@ -348,59 +298,146 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
utau[i] += u[j]*tau[i][j];
}
}
f0[Q] = divu+(q-qold)/dt;
f0[Q] = (q-qold)/dt;
f0[P] -= q;
for (int i = 0; i < D; ++i) {
f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) - p*dq[i] + (i==(D-1) ? -rho*problem->g : 0);
for (int j = 0; j < D; ++j) {
f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
}
f1[Q*D+i] = epsilonq*dq[i]+epsilonp*q*dp[i];
f1[Q*D+i] = epsilonq*dq[i]-u[i];
}
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
{
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
double deps = 1e-8;
int n_fluids = problem->n_fluids;
int n_fields = fluid_problem_n_fields(problem);
const double *solution = problem->solution;
const double *compacity = problem->compacity;
const double epsilon = problem->epsilon;
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
for (size_t i = 0; i < N_E*local_size; ++i)
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
all_local_matrix[i] = 0;
const int n_fields = fluid_problem_n_fields(problem);
const size_t local_size = N_SF*n_fields;
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 *g0 = malloc(sizeof(double)*n_fields);
double *h0 = malloc(sizeof(double)*n_fields*D);
f_cb *f;
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
forced_row[i] = -1;
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys)
printf("Boundary tag \"%s\" not found.", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
if (strcmp(mesh->phys_name[iphys],"Bottom") == 0)
f = f_wall_pressure;
else if (strcmp(mesh->phys_name[iphys],"Right") == 0)
f = f_wall_pressure;
else if (strcmp(mesh->phys_name[iphys],"Left") == 0)
f = f_wall_pressure;
else if (strcmp(mesh->phys_name[iphys],"Top") == 0)
f = f_wall_pressure;
else
continue;
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];
const int i0 = bnd[1];
const int i1 = (bnd[1]+1) % 3;
const unsigned int *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]};
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
for (int i = 0; i < D; ++i)
for (int j = 0; j < D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {-dx[1]/l,dx[0]/l};
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double xi[2];
double a = (1+LQP[iqp])/2;
switch(i0) {
case 0 : xi[0] = a; xi[1] = 0; break;
case 1 : xi[0] = 1-a; xi[1] = a; break;
case 2 : xi[0] = 0; xi[1] = 1-a; break;
default:;
}
double dp[D]={0};
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double phi[N_SF];
shape_functions(xi,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];
s[j] += phi[i]*dof;
for (int k = 0; k < D; ++k)
ds[j*D+k] += dphi[i][k]*dof;
}
}
const double jw = LQW[iqp]*l/2;
f(problem,n,f0,s,ds,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 jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
f(problem,n,g0,s,ds,c,dt);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
f(problem,n,h0+jd*n_fields,s,ds,c,dt);
ds[jfield*D+jd] = store;
}
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 jd = 0; jd < D; ++jd) {
double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
m += jw*phi[iphi]*dphi[jphi][jd]*df01;
}
int IDX = (iphi*N2)*n_fields + jphi*n_fields;
LOCAL_MATRIX2(ifield,jfield) += m;
}
}
}
}
}
}
}
}
//compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *compacity = problem->compacity;
const double epsilon = problem->epsilon;
const double *solution = problem->solution;
int n_fluids = problem->n_fluids;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
......@@ -420,7 +457,6 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
int P = (D+1)*n_fluids;
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
......@@ -491,8 +527,62 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
}
free(s);
free(ds);
free(sold);
free(f0);
free(f1);
free(g0);
free(g1);
free(h0);
free(h1);
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
int n_fluids = problem->n_fluids;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
for (size_t i = 0; i < N_E*local_size; ++i)
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
all_local_matrix[i] = 0;
char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
forced_row[i] = -1;
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys)
printf("Boundary tag \"%s\" not found.", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
}
}
//compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
for (int iel=0; iel < mesh->n_elements; ++iel) {
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
for (int inode = 0; inode < N_SF; ++inode) {
for(int ifield = 0; ifield < n_fields; ++ifield) {
const unsigned int *el = &mesh->elements[iel*N_N];
char forced_field = forced_row[el[inode]*n_fields+ifield];
if (forced_field == -1) continue;
int i = inode*n_fields + ifield;
......@@ -508,19 +598,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
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(forced_row);
free(all_local_matrix);
free(all_local_vector);
free(forced_row);
}
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_fields, int nBnd, StrongBoundary *bnds)
......
......@@ -59,9 +59,9 @@ outf = 1 # number of iterations between ou
#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,1,0.),("Bottom",1,1,0.),("Right",1,1,0.),("Left",1,1,0.),
("Top",0,0,0.),("Bottom",0,0,0.),("Left",0,0,0),("Right",0,0,0),
("Top",2,3,0)]#,("Bottom",2,3,-rho*g)]#,("Left",2,3,lambda x : -rho*g*(1-x[:,1])),("Right",2,3,lambda x : -rho*g*(1-x[:,1]))]
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),
("Left",0,0,0),("Right",0,0,0),
("Top",2,3,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1)
ii = 0
t = 0
......@@ -70,7 +70,7 @@ t = 0
x = fluid.coordinates()
s = fluid.solution()
s[:,2] = 1.
s[:,3] = -rho*g*(1-x[:,1])
#s[:,3] = -rho*g*(1-x[:,1])
fluid.export_vtk(outputdir,0,0)
......@@ -96,4 +96,4 @@ print('Error', (vel.sum())**.5)
if (vel.sum())**.5<1.5e-3:
exit(1)
else:
exit(0)
\ No newline at end of file
exit(0)
......@@ -48,7 +48,7 @@ class Poiseuille(unittest.TestCase) :
#physical parameters
g = -9.81 # gravity
g = 0 # gravity
rho = 1000 # fluid density
nu0 = 1e-3 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
......
......@@ -47,7 +47,7 @@ class Poiseuille(unittest.TestCase) :
#physical parameters
g = -9.81 # gravity
g = 0 # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
......
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