Commit c4807f3c authored by Matthieu Constant's avatar Matthieu Constant
Browse files

generalized weak boundaries

parent 65c3f712
Pipeline #4569 passed with stage
in 1 minute and 52 seconds
......@@ -328,7 +328,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
free(h0);
}
static void f_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, double v){
static void f_pressure(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
double p = s[P];
double u[D], du[D][D], dp[D];
for (int iD = 0; iD < D; ++iD) {
......@@ -361,14 +361,14 @@ static void f_pressure(FluidProblem *problem,const double *n, double *f,const do
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a+problem->rho[1]*(1-a);
mu = problem->mu[0]*a+problem->mu[1]*(1-a);
f[A] = un*(un > 0 ? a : 1.);
f[A] = un*(un > 0 ? a : v[1]);
}
for (int id = 0; id < D; ++id) {
f[U+id] = v*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
f[U+id] = v[0]*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
}
}
static void f_velocity(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid,double v){
static void f_velocity(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
double p = s[P];
double u[D], du[D][D], dp[D];
for (int iD = 0; iD < D; ++iD) {
......@@ -380,12 +380,12 @@ static void f_velocity(FluidProblem *problem,const double *n, double *f,const do
double divu = 0;
double tau[D][D];
double utau[D]={0.};
double rho, mu;
double rho, mu, rho_ext;
double unint = 0;
double un = 0;
for (int i = 0; i < D; ++i) {
un += v*n[i];
un += v[i]*n[i];
unint += u[i]*n[i];
divu += du[i][i];
for (int j = 0; j < D; ++j) {
......@@ -397,21 +397,24 @@ static void f_velocity(FluidProblem *problem,const double *n, double *f,const do
if(problem->n_fluids == 1){
mu = problem->mu[0];
rho = problem->rho[0];
rho_ext = rho;
}
else{
double a = s[A];
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a+problem->rho[1]*(1-a);
rho_ext = problem->rho[0]*v[D]+problem->rho[1]*(1-v[D]);
mu = problem->mu[0]*a+problem->mu[1]*(1-a);
f[A] = un*(un > 0 ? a : 1.);
f[A] = un*(un > 0 ? a : v[D]);
//printf("%g\n",v[D]);
}
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*unint/c : rho*un*un*n[id]/c);
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*unint/c : rho_ext*un*un*n[id]/c);
}
}
static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc,const double dt,int eid,double v)
static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc, const double dt, int eid, const double *v)
{
double p = s[P];
double u[D], du[D][D], dp[D];
......@@ -613,7 +616,6 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
bnd_node_local_id[i] = -1;
if(mesh->phys_dimension[iphys] != D-1) continue;
double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
int j = mesh->phys_nodes[iphys][i];
bnd_node_local_id[j] = i;
......@@ -622,11 +624,17 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
f = NULL;
int n_value;
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *bnd = problem->weak_boundaries + ibnd;
if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
f = bnd->cb;
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
n_value = bnd->n_value;
if (n_value > 0){
v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
}
break;
}
}
......@@ -701,13 +709,18 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
double c = 0;
double vbnd = 0;
double *vbnd = malloc(sizeof(double)*n_value);
for (int i = 0; i < n_value; ++i){
vbnd[i] = 0;
}
double dc[D] = {0};
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
int node_local_id = bnd_node_local_id[el[i]];
if (node_local_id != -1)
vbnd += v[node_local_id]*phi[i];
for (int j = 0; j < n_value; j++){
vbnd[j] += v[node_local_id*n_value+j]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
s[j] += phi[i]*dof;
......@@ -753,6 +766,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
}
free(vbnd);
}
}
free(x);
......@@ -1558,12 +1572,18 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const cha
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;
if (strcasecmp(bndtype,"wall") == 0)
if (strcasecmp(bndtype,"wall") == 0){
f = f_wall;
else if (strcasecmp(bndtype,"velocity") == 0)
p->weak_boundaries[p->n_weak_boundaries-1].n_value = 0;
}
else if (strcasecmp(bndtype,"velocity") == 0){
f = f_velocity;
else if (strcasecmp(bndtype,"pressure") == 0)
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;
}
else {
printf("Unkown weak boundary type \"%s\".\n", bndtype);
exit(1);
......
......@@ -43,10 +43,11 @@ 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 dt, int eid, double v);
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 dt, int eid, const double *v);
typedef struct {
char *tag;
f_cb *cb;
int n_value;
StrongBoundaryCallback *apply;
} WeakBoundary;
......
......@@ -75,6 +75,7 @@ class FluidProblem :
raise ValueError("dim should be 2 or 3.")
self._lib.fluid_problem_new.restype = c_void_p
n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
self._dim = dim
self._fp = c_void_p(self._lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), c_double(coeff_stab), petsc_solver_type.encode()))
if self._fp == None :
......@@ -93,26 +94,41 @@ class FluidProblem :
"""
self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
def set_weak_boundary(self, tag, bnd_type, callback_or_value=0) :
def set_weak_boundary(self, tag, bnd_type, callback_or_value=[]) :
"""Set the weak boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
tag -- physical tag (set in the mesh.geo file) of the boundary on which the flux type is added
bnd_type -- type of the boundary flux: wall, inflow, outflow
bnd_type -- type of the boundary flux: wall, pressure, velocity
callback_or_value -- list of values (functions) to specify the weak boundary values
Raises:
ValueError -- if the number of values set in the list callback_or_value do not match what is expected for the boundary type : wall do not take any values; pressure takes a value for the pressure and concentration if any; velocity takes values for each velocity dimension and the concentration if any.
"""
if bnd_type.lower() == 'wall' :
self._n_value = 0
elif bnd_type.lower() == 'velocity' :
self._n_value = self._dim + self._n_fluids-1 # u, v, w, p, a
elif bnd_type.lower() == 'pressure' :
self._n_value = self._n_fluids # p, a
if len(callback_or_value) != self._n_value :
raise ValueError("Number of weak boundary values don't match the weak boundary type.")
class Bnd :
def __init__(self, b, dim) :
def __init__(self, b, dim, n_value) :
self._b = b
self._dim = dim
self._n_value = n_value
def apply(self, n, xp, vp) :
x = numpy.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,self._dim])
v = numpy.frombuffer(cast(vp, POINTER(int(n)*c_double)).contents)
if callable(self._b) :
v[:] = self._b(x)
else :
v[:] = self._b
bndcb = BNDCB(Bnd(callback_or_value,self._dim).apply)
v = numpy.frombuffer(cast(vp, POINTER(int(n)*self._n_value*c_double)).contents).reshape([n,self._n_value])
for i in range(self._n_value):
if callable(self._b[i]) :
v[:,i] = self._b[i](x)
else :
v[:,i] = self._b[i]
bndcb = BNDCB(Bnd(callback_or_value, self._dim, self._n_value).apply)
self.weak_cb_ref.append(bndcb)
self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bnd_type.encode(),bndcb)
......@@ -138,7 +154,7 @@ class FluidProblem :
v[:] = self._b(x)
else :
v[:] = self._b
bndcb = BNDCB(Bnd(callback_or_value,self._dim).apply)
bndcb = BNDCB(Bnd(callback_or_value, self._dim).apply)
self.strong_cb_ref.append(bndcb)
self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
......
......@@ -66,30 +66,28 @@ lx = .222
ly = .14*2
rhop = 2500
# numerical parameters
dt = .0001 # time step
dt = .0005 # time step
genInitialPosition(outputdir, N, r, lx, ly, rhop)
#genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction=0.3
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
# number of iterations between output files
outf = 25
outf = 1
ii = 0
t = 0
def outerBndV(x) :
return max(1*t,0.1)
return min((6*t**5-15*t**4+10*t**3),1)*0.1
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
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_strong_boundary("Injection",0,0)
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_weak_boundary("Bottom","wall")
fluid.set_weak_boundary("Lateral","wall")
fluid.set_weak_boundary("Top","pressure",[0, 1])
fluid.set_weak_boundary("Injection","velocity",[0, outerBndV, 1])
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......@@ -106,7 +104,7 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
#Fluid solver
fluid.implicit_euler(dt)
fluid.implicit_euler(dt, newton_max_it=20)
if (ii%65==0 and ii != 0):
fluid.adapt_mesh(8e-2,3e-3,10000)
forces = fluid.compute_node_force(dt)
......
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