Commit 0d71c513 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

weak wall pressure

parent 1a4e4463
Pipeline #4593 failed with stage
in 17 seconds
......@@ -633,15 +633,15 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
n_value = bnd->n_value;
if (n_value > 0){
v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
printf("tag : %s apply = %p %i\n", bnd->tag, bnd->apply,bnd->n_value);
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
}
break;
}
}
if (f == NULL) {
f = f_wall;
//printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
//exit(1);
printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
exit(1);
}
for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
......@@ -1565,6 +1565,54 @@ void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, i
problem->strong_boundaries[i].row = row;
}
void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
printf("Weak boundary is already defined for tag \"%s\".\n", tag);
exit(1);
}
}
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);
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].apply = pressurecb;
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
printf("Weak boundary is already defined for tag \"%s\".\n", tag);
exit(1);
}
}
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;
if (strcasecmp(bndtype,"wall") == 0){
f = f_wall;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = 0;
}
else if (strcasecmp(bndtype,"velocity") == 0){
f = f_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;
}
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].apply = apply;
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
......
......@@ -106,6 +106,7 @@ int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb);
int fluid_problem_private_n_physical(const FluidProblem *p);
void fluid_problem_private_physical(const FluidProblem *p, int i, char **phys_name, int *phys_dim, int *phys_id, int *phys_n_nodes, int **phys_nodes, int *phys_n_boundaries, int **phys_boundaries);
......
......@@ -32,6 +32,7 @@ from ctypes import *
import numpy
import signal
import os
import sys
from . import VTK
dir_path = os.path.dirname(os.path.realpath(__file__))
......@@ -41,6 +42,25 @@ lib3 = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
signal.signal(signal.SIGINT, signal.SIG_DFL)
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
class _Bnd :
def __init__(self, b, dim) :
self._b = b
self._dim = dim
def apply(self, n, xp, vp) :
nv = len(self._b)
x = numpy.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,-1])
v = numpy.frombuffer(cast(vp, POINTER(int(n)*nv*c_double)).contents).reshape([n,nv])
for i in range(nv):
if callable(self._b[i]) :
v[:,i] = self._b[i](x)
else :
v[:,i] = self._b[i]
def _is_string(s) :
if sys.version_info >= (3, 0):
return isinstance(s, str)
else :
return isinstance(s, basestring)
class FluidProblem :
"""Create numerical representation of the fluid."""
......@@ -98,15 +118,18 @@ class FluidProblem :
"""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, pressure, velocity
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.
tag -- physical tag (or list of tags), set in the mesh.geo file, of the wall boundaries
pressure -- the pressure value if imposed (callback or number)
"""
if not _is_string(tag) :
for t in tag :
self.set_wall_boundary(t,pressure)
return
pressurecb = None
if pressure :
raise NameError("Wall pressure not yet implemented")
self._lib.fluid_problem_set_weak_boundary(self._fp,tag.encode(),b"wall",None)
pressurecb = BNDCB(_Bnd(pressure,self._dim).apply)
self.weak_cb_ref.append(pressurecb)
self._lib.fluid_problem_set_wall_boundary(self._fp,tag.encode(), pressurecb)
def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1):
"""Set the weak boundaries (=normal fluxes) for the fluid problem.
......@@ -118,19 +141,6 @@ class FluidProblem :
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.
"""
class Bnd :
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)*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]
if porosity != 1 :
raise NameError("Only porosity=1 is implemented")
if (velocity is None and pressure is None) or (velocity is not None and pressure is not None) :
......@@ -143,7 +153,7 @@ class FluidProblem :
n_value = self._n_fluids # p, a
cb_or_value = [pressure,concentration] if self._n_fluids == 2 else [pressure]
bnd_type = "pressure"
bndcb = BNDCB(Bnd(cb_or_value, self._dim, n_value).apply)
bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
self.weak_cb_ref.append(bndcb)
self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bnd_type.encode(),bndcb)
......@@ -158,18 +168,7 @@ class FluidProblem :
"""
if row is None :
row = field
class Bnd :
def __init__(self, b, dim) :
self._b = b
self._dim = dim
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)
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)
......
......@@ -114,13 +114,12 @@ tEnd1 = 0.2
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_wall_boundary("Top",pressure=0)
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
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......@@ -186,4 +185,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.time() - tic))
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
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