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

weak bnd 2 fluids

parent 31612c99
Pipeline #4767 passed with stage
in 30 seconds
......@@ -33,6 +33,8 @@
#define P D
#define U 0
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double rho, const double mu, const double dt, int eid, const double *v);
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
#if D==2
......@@ -313,7 +315,8 @@ static void f_wall(FluidProblem *problem,const double *n, double *f,const double
for (int i = 0; i < D; ++i) {
un += u[i]*n[i];
}
f[P] = 0;
double h = problem->element_size[eid];
f[P] = (v==NULL ? 0 : -(s[P]-v[0])/h*problem->taup[eid]);
for (int id = 0; id < D; ++id) {
f[U+id] = (v==NULL?p:v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
......@@ -324,10 +327,10 @@ static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
const Mesh *mesh = problem->mesh;
problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
problem->taua = realloc(problem->taua,sizeof(double)*mesh->n_elements);
problem->extraa = realloc(problem->extraa,sizeof(double)*mesh->n_elements);
double maxtaup = 0;
double mintaup = DBL_MAX;
......@@ -363,6 +366,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
#else
const double h = 2*cbrt(detDD(dxdxi)/(8*M_PI));
#endif
problem->element_size[iel] = h;
double rho, mu;
if(problem->n_fluids == 1){
rho = problem->rho[0];
......@@ -377,7 +381,6 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
problem->extraa[iel] = fmax(0.,fmax(-amin+0.01,amax-1+0.01))*problem->coeffStab;
}
}
......@@ -386,7 +389,6 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
double taua = problem->taua[iel];
double extraa = problem->extraa[iel];
double u[D], uold[D], du[D][D], dp[D];
for (int iD = 0; iD < D; ++iD) {
u[iD] = sol[U+iD];
......@@ -475,6 +477,15 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
wbnd->apply(bnd->n_nodes, x, v);
free(x);
}
f_cb *bndcb;
switch(wbnd->type) {
case BND_WALL : bndcb = f_wall; break;
case BND_PRESSURE : bndcb = f_pressure; break;
case BND_VELOCITY : bndcb = f_velocity; break;
default :
printf("Unknown boundary type : %i\n", wbnd->type);
exit(1);
};
double *vbnd = malloc(sizeof(double)*wbnd->n_value);
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
......@@ -509,25 +520,27 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
#endif
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double xi[D];
{
#if DIMENSION == 2
double a = LQP[iqp];
switch(icl) {
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 a = LQP[iqp];
switch(icl) {
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:;
}
#else
double a = LQP[iqp][0];
double b = LQP[iqp][1];
switch(icl) {
case 0 : xi[0] = a; xi[1] = b; xi[2] = 0; break;
case 1 : xi[0] = 0; xi[1] = a; xi[2] = b; break;
case 2 : xi[0] = b; xi[1] = 0; xi[2] = a; break;
case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break;
default:;
}
double a = LQP[iqp][0];
double b = LQP[iqp][1];
switch(icl) {
case 0 : xi[0] = a; xi[1] = b; xi[2] = 0; break;
case 1 : xi[0] = 0; xi[1] = a; xi[2] = b; break;
case 2 : xi[0] = b; xi[1] = 0; xi[2] = a; break;
case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break;
default:;
}
#endif
}
double dp[D]={0};
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
......@@ -545,6 +558,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
vbnd[i] = 0;
}
double dc[D] = {0};
double a = 0;
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
if (problem->n_fluids == 2) {
......@@ -571,7 +585,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
mu = problem->mu[0]*a + problem->mu[1]*(1-a);
}
const double jw = LQW[iqp]*detbnd;
wbnd->cb(problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0 ? NULL :vbnd);
bndcb(problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0 ? NULL :vbnd);
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;
......@@ -580,12 +594,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
wbnd->cb(problem,n,g0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0? NULL :vbnd);
bndcb(problem,n,g0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0? NULL :vbnd);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
wbnd->cb(problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value ==0? NULL : vbnd);
bndcb(problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value ==0? NULL : vbnd);
ds[jfield*D+jd] = store;
}
int N2 = n_fields*N_SF;
......@@ -978,7 +992,7 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
problem->taup = NULL;
problem->tauc = NULL;
problem->taua = NULL;
problem->extraa = NULL;
problem->element_size = NULL;
problem->porosity = NULL;
problem->oldporosity = NULL;
problem->solution = NULL;
......@@ -1018,7 +1032,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->taup);
free(problem->tauc);
free(problem->taua);
free(problem->extraa);
free(problem->element_size);
for (int i = 0; i < problem->n_strong_boundaries; ++i)
free(problem->strong_boundaries[i].tag);
for (int i = 0; i < problem->n_weak_boundaries; ++i)
......@@ -1472,7 +1486,7 @@ void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryC
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
p->weak_boundaries[p->n_weak_boundaries-1].tag = strdup(tag);
p->weak_boundaries[p->n_weak_boundaries-1].n_value = (pressurecb ? 1 : 0);
p->weak_boundaries[p->n_weak_boundaries-1].cb = f_wall;
p->weak_boundaries[p->n_weak_boundaries-1].type = BND_WALL;
p->weak_boundaries[p->n_weak_boundaries-1].apply = pressurecb;
}
......@@ -1488,20 +1502,20 @@ void fluid_problem_set_open_boundary(FluidProblem *p, const char *tag, const cha
p->n_weak_boundaries++;
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
p->weak_boundaries[p->n_weak_boundaries-1].tag = strdup(tag);
f_cb *f = NULL;
BoundaryType btype = 0;
if (strcasecmp(bndtype,"velocity") == 0){
f = f_velocity;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = D+p->n_fluids-1;
btype = BND_VELOCITY;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = D+p->n_fluids-1;
}
else if (strcasecmp(bndtype,"pressure") == 0){
f = f_pressure;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = p->n_fluids;
btype = BND_PRESSURE;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = p->n_fluids;
}
else {
printf("Unkown weak boundary type \"%s\".\n", bndtype);
exit(1);
}
p->weak_boundaries[p->n_weak_boundaries-1].cb = f;
p->weak_boundaries[p->n_weak_boundaries-1].type = btype;
p->weak_boundaries[p->n_weak_boundaries-1].apply = apply;
}
......
......@@ -35,6 +35,9 @@
#endif
typedef void BoundaryCallback(int n, const double *x, double *bnd);
typedef enum {BND_WALL=0,BND_VELOCITY,BND_PRESSURE} BoundaryType;
typedef struct {
char *tag;
int row;
......@@ -43,10 +46,9 @@ typedef struct {
} StrongBoundary;
typedef struct FluidProblem FluidProblem;
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double rho, const double mu, const double dt, int eid, const double *v);
typedef struct {
char *tag;
f_cb *cb;
BoundaryType type;
int n_value;
BoundaryCallback *apply;
} WeakBoundary;
......@@ -67,6 +69,7 @@ struct FluidProblem {
double *solution;
double *solution_explicit;
double *node_volume;
double *element_size;
int n_strong_boundaries;
StrongBoundary *strong_boundaries;
int n_weak_boundaries;
......@@ -84,7 +87,6 @@ struct FluidProblem {
double *taup;
double *tauc;
double *taua;
double *extraa;
int *particle_element_id;
};
......
......@@ -2,6 +2,8 @@
#include "fluid_problem_fem.h"
#include <math.h>
#include <float.h>
#include <string.h>
void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
const Mesh *mesh = problem->mesh;
const int n_fields = fluid_problem_n_fields(problem);
......@@ -102,6 +104,148 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
}
}
/*boundary_term*/
{
double u[D];
double *f0 = malloc(sizeof(double)*n_fields);
int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
#if DIMENSION == 2
int elbnd[3][2] = {{0,1},{1,2},{2,0}};
#else
int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
for (int i = 0; i < mesh->n_boundaries; ++i) {
if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
bndid = i;
}
}
if (bndid == -1) {
printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
}
MeshBoundary *bnd = &problem->boundaries[bndid];
double *v = malloc(wbnd->n_value*bnd->n_nodes*sizeof(double));
for (int i = 0; i < mesh->n_nodes; ++i) {
bnd_node_local_id[i] = -1;
}
if (wbnd->n_value != 0) {
double *x = (double*)malloc(bnd->n_nodes*sizeof(double)*D);
for (int i = 0; i < bnd->n_nodes; ++i){
int j = bnd->nodes[i];
bnd_node_local_id[j] = i;
for (int k = 0; k < D; ++k)
x[i*D+k] = mesh->x[j*3+k];
}
wbnd->apply(bnd->n_nodes, x, v);
free(x);
}
double *vbnd = malloc(sizeof(double)*wbnd->n_value);
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
const int iel = interface[0];
const int icl = interface[1];
const int *cl = elbnd[interface[1]];
int nodes[D];
double *x[D];
const unsigned int *el = &mesh->elements[iel*N_N];
for (int i = 0; i < D; ++i){
nodes[i] = el[cl[i]];
x[i] = &mesh->x[nodes[i]*3];
}
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);
#if DIMENSION == 2
const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {dx[1]/detbnd,-dx[0]/detbnd};
#else
const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
const double nn[3] = {
dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
const double n[3] = {-nn[0]/detbnd,-nn[1]/detbnd,-nn[2]/detbnd};
#endif
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double xi[D];
{
#if DIMENSION == 2
double a = LQP[iqp];
switch(icl) {
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:;
}
#else
double a = LQP[iqp][0];
double b = LQP[iqp][1];
switch(icl) {
case 0 : xi[0] = a; xi[1] = b; xi[2] = 0; break;
case 1 : xi[0] = 0; xi[1] = a; xi[2] = b; break;
case 2 : xi[0] = b; xi[1] = 0; xi[2] = a; break;
case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break;
default:;
}
#endif
}
double dp[D]={0};
double phi[N_SF];
shape_functions(xi,phi);
for (int i = 0; i < D; ++i) {
u[i] = 0;
}
for (int i = 0; i < wbnd->n_value; ++i){
vbnd[i] = 0;
}
double dc[D] = {0};
double a = 0;
for (int i = 0; i < N_SF; ++i) {
a += problem->concentration[iel*N_N+i]*phi[i];
int node_local_id = bnd_node_local_id[el[i]];
if (node_local_id != -1)
for (int j = 0; j < wbnd->n_value; j++){
vbnd[j] += v[node_local_id*wbnd->n_value+j]*phi[i];
}
for (int j = 0; j < D; ++j) {
u[j] += phi[i]*problem->solution[el[i]*n_fields+j];
}
}
double un = 0;
for (int i = 0; i< D; ++i) {
un += n[i]* (wbnd->type == BND_VELOCITY ? vbnd[i] : u[i]);
}
const double jw = LQW[iqp]*detbnd;
double f0;
switch(wbnd->type) {
case BND_WALL : f0 = 0; break;
case BND_PRESSURE : {
f0 = -un * (un > 0 ? a : 1);
break;
}
case BND_VELOCITY : {
f0 = -un * (un > 0 ? a : 1);
break;
}
}
for (int iphi = 0; iphi < N_SF; ++iphi){
rhs[iel*N_N+iphi] += phi[iphi]*f0*jw*dt;
}
}
}
free(vbnd);
free(v);
}
free(bnd_node_local_id);
free(f0);
}
/*resolution*/
#if DIMENSION == 2
......
......@@ -111,11 +111,11 @@ tEnd1 = 0.2
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],petsc_solver_type="-pc_type ilu -ksp_max_it 30")
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],petsc_solver_type="-pc_type ilu -ksp_max_it 50")
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Top")
fluid.set_wall_boundary("Top", pressure=0)
#fluid.set_wall_boundary("Top")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
......@@ -148,6 +148,7 @@ if init:
t += dt1
t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......@@ -157,7 +158,7 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#
while t < tEnd :
# Fluid solver
fluid.implicit_euler(dt)
fluid.implicit_euler(dt,newton_max_it=100)
# Adaptation of the mesh.
if (ii%15==0 and ii != 0):
......
......@@ -75,15 +75,13 @@ def outerBndV(x) :
print(.4*max(np.sin(t*np.pi*2./1),0))
return 0.4*max(np.sin(t*np.pi*2./1),0)
strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
#fluid.set_strong_boundary("Top",2,0)
#fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Top",pressure=0)
fluid.set_open_boundary("Bottom",velocity=[0,outerBndV])
ii = 0
t = 0
......@@ -111,4 +109,4 @@ while t < tEnd :
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -43,7 +43,6 @@ def genInitialPosition(filename, r, N, rhop) :
break
p.add_particle((xl, yl), r, r**2 * np.pi * rhop);
i += 1
print(len(p.mass()))
p.write_vtk(filename,0,0)
outputdir = "output"
......@@ -61,25 +60,19 @@ N = 10000
rhop = 1500
# numerical parameters
tEnd = 10 # final time
dt = .01 # time step
outf = 1 # number of iterations between output files
dt = .001 # time step
outf = 10 # number of iterations between output files
def outerBndV(x) :
print(.4*max(np.sin(t*np.pi*2./1),0))
if t > 0.5 :
return 0
return 0.4*max(np.sin(t*np.pi*2./1),0)
a = np.where(np.abs(x[:,0]-0.111) < 0.01, 0.4*max(np.sin(t*np.pi*2./1),0),0)
return a
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],petsc_solver_type="-pc_type ilu -ksp_rtol 1e-2 -ksp_max_it 100")
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_strong_boundary("Injection",3,1)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Top",pressure=0)
fluid.set_open_boundary("Bottom",velocity=[0,outerBndV])
ii = 0
t = 0
......@@ -95,4 +88,4 @@ while t < tEnd :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
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