diff --git a/hysop/backend/host/python/operator/penalization.py b/hysop/backend/host/python/operator/penalization.py
index dedb44151c37570533c68aead17714fa7f614af5..c11f92a9ce09412076fcc822c7f119255ce32b3c 100644
--- a/hysop/backend/host/python/operator/penalization.py
+++ b/hysop/backend/host/python/operator/penalization.py
@@ -1,8 +1,10 @@
+from hysop.constants import PenalizationFormulation
from hysop.backend.host.host_operator import HostOperator
from hysop.tools.types import check_instance, InstanceOf
from hysop.tools.decorators import debug
from hysop.fields.continuous_field import Field
from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.parameters.tensor_parameter import TensorParameter
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.methods import SpaceDiscretization
from hysop.core.memory.memory_request import MemoryRequest
@@ -28,6 +30,7 @@ class PythonPenalizeVorticity(HostOperator):
dm = super(PythonPenalizeVorticity, cls).default_method()
dm.update(cls.__default_method)
return dm
+
@classmethod
def available_methods(cls):
am = super(PythonPenalizeVorticity, cls).available_methods()
@@ -37,15 +40,17 @@ class PythonPenalizeVorticity(HostOperator):
@debug
def __init__(self, obstacles, variables,
velocity, vorticity,
- dt, coeff=None, **kwds):
+ dt, coeff=None, ubar=None, formulation=None, **kwds):
check_instance(velocity, Field)
check_instance(vorticity, Field)
check_instance(variables, dict, keys=Field,
values=CartesianTopologyDescriptors)
check_instance(dt, ScalarParameter)
- check_instance(coeff, (ScalarParameter, float, type(lambda x:x)))
+ check_instance(coeff, (ScalarParameter, float, type(lambda x: x)), allow_none=True)
+ check_instance(ubar, TensorParameter, allow_none=True)
+ check_instance(formulation, PenalizationFormulation, allow_none=True)
check_instance(obstacles, (tuple, dict), values=Field,
- keys=(ScalarParameter, float), check_kwds=False)
+ keys=(ScalarParameter, float, type(lambda x: x)), check_kwds=False)
input_fields = {velocity: variables[velocity],
vorticity: variables[vorticity]}
@@ -54,14 +59,22 @@ class PythonPenalizeVorticity(HostOperator):
if isinstance(coeff, ScalarParameter):
self.coeff = lambda o: coeff()*o
- if isinstance(coeff, type(lambda x:x)):
+ elif isinstance(coeff, type(lambda x: x)):
self.coeff = coeff
- else:
+ elif not coeff is None:
c = ScalarParameter("penal_coeff", initial_value=coeff)
self.coeff = lambda o: c()*o
if isinstance(obstacles, dict):
- obs = obstacles
+ obs = {}
+ for c, o in obstacles.iteritems():
+ if isinstance(c, ScalarParameter):
+ obs[lambda x: c()*x] = o
+ elif isinstance(c, type(lambda x: x)):
+ obs[c] = o
+ elif not c is None:
+ _ = ScalarParameter("penal_coeff", initial_value=c)
+ obs[lambda x: _()*x] = o
else:
obs = {}
for o in obstacles:
@@ -69,6 +82,16 @@ class PythonPenalizeVorticity(HostOperator):
for o in obs.values():
assert o.nb_components == 1
input_fields[o] = variables[o]
+ self._ubar = ubar
+ if ubar is None:
+ self._ubar = TensorParameter(
+ name="ubar", shape=(velocity.dim,),
+ quiet=True, dtype=velocity.dtype,
+ initial_value=(0.,)*velocity.dim)
+ if formulation is None or formulation is PenalizationFormulation.IMPLICIT:
+ self._compute_penalization = self._compute_penalization_implicit
+ else:
+ self._compute_penalization = self._compute_penalization_exact
self.velocity = velocity
self.vorticity = vorticity
@@ -94,7 +117,8 @@ class PythonPenalizeVorticity(HostOperator):
@debug
def get_field_requirements(self):
- requirements = super(PythonPenalizeVorticity, self).get_field_requirements()
+ requirements = super(PythonPenalizeVorticity,
+ self).get_field_requirements()
stencil = self.stencil[0]
G = max(max(a, b) for a, b in zip(stencil.L, stencil.R))
for is_input, (field, td, req) in requirements.iter_requirements():
@@ -122,16 +146,18 @@ class PythonPenalizeVorticity(HostOperator):
self.dobstacles = {}
for c, o in self.obstacles.iteritems():
o_df = self.input_discrete_fields[o]
- self.dobstacles[c] = o_df.data[0][o_df.local_slices(ghosts=(G, G, G))]
+ self.dobstacles[c] = o_df.data[0][o_df.local_slices(
+ ghosts=(G, G, G))]
for s, dx in zip(self.stencil, dw.space_step):
s.replace_symbols({s.dx: dx})
- self.ghost_exchanger = dw.build_ghost_exchanger()
+ self.w_ghost_exchanger = dw.build_ghost_exchanger()
+ self.v_ghost_exchanger = dv.build_ghost_exchanger()
@debug
def get_work_properties(self):
- requests = super(PythonPenalizeVorticity, self).get_work_properties()
+ requests = super(PythonPenalizeVorticity, self).get_work_properties()
buffers = MemoryRequest.empty_like(
a=self.V[0], nb_components=4, backend=self.dvelocity.backend)
requests.push_mem_request('Wtmp', buffers)
@@ -151,18 +177,29 @@ class PythonPenalizeVorticity(HostOperator):
def supports_mpi(cls):
return True
+ def _compute_penalization_implicit(self):
+ dt = self.dt()
+ self.tmp[...] = 0.
+ for c, o in self.dobstacles.iteritems():
+ self.tmp[...] += (-dt) * c(o) / (1.0 + dt * c(o))
+
+ def _compute_penalization_exact(self):
+ dt = self.dt()
+ self.tmp[...] = 1.
+ for c, o in self.dobstacles.iteritems():
+ self.tmp[...] *= np.exp(-dt*c(o))
+ self.tmp[...] -= 1.
+
@op_apply
def apply(self, **kwds):
super(PythonPenalizeVorticity, self).apply(**kwds)
- dt = self.dt()
V, W = self.V, self.W
+ self.v_ghost_exchanger()
- # Penalize velocity
- self.tmp[...] = 0.
- for c, o in self.dobstacles.iteritems():
- self.tmp[...] += - dt * c(o) / (1.0 + dt * c(o))
- for v, tmp in zip(V, self.tmp_v):
- tmp[...] = v * self.tmp
+ # # Penalize velocity
+ self._compute_penalization()
+ for ubar, v, tmp in zip(self._ubar(), V, self.tmp_v):
+ tmp[...] = (v - ubar) * self.tmp
# compute penalized vorticity:
# (dvz/dy - dvy/dz)
@@ -186,4 +223,4 @@ class PythonPenalizeVorticity(HostOperator):
self.tmp = self.stencil[2](a=self.tmp_v[0], out=self.tmp, axis=0)
W[1][...] += self.tmp
- self.ghost_exchanger()
+ self.w_ghost_exchanger()
diff --git a/hysop/constants.py.in b/hysop/constants.py.in
index cf9c1aeaf7b6d58b88cf5dfff284575386c3e5fb..12742767f0ac49b3c26b960e184f43b156577998 100644
--- a/hysop/constants.py.in
+++ b/hysop/constants.py.in
@@ -189,6 +189,11 @@ StretchingFormulation = EnumFactory.create('StretchingFormulation',
['GRAD_UW', 'GRAD_UW_T', 'MIXED_GRAD_UW', 'CONSERVATIVE'])
"""Stretching formulations"""
+PenalizationFormulation = EnumFactory.create(
+ 'PenalizationFormulation',
+ ['IMPLICIT', 'EXACT'])
+"""Penalization formulations"""
+
SpaceDiscretization = EnumFactory.create('SpaceDiscretization',
['FDC2', 'FDC4', 'FDC6', 'FDC8'])
"""Space discretization for stencil generation"""
diff --git a/hysop/operator/penalization.py b/hysop/operator/penalization.py
index f829168b17439c0a4e30df996524d28eb06405ea..94d750e7392f588948447ad1f053fdb54d3541c4 100644
--- a/hysop/operator/penalization.py
+++ b/hysop/operator/penalization.py
@@ -7,12 +7,13 @@
See details in :ref:`penalisation` section of HySoP user guide.
"""
-from hysop.constants import Implementation
+from hysop.constants import Implementation, PenalizationFormulation
from hysop.tools.types import check_instance, to_list
from hysop.tools.decorators import debug
from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.parameters.tensor_parameter import TensorParameter
from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
from hysop.backend.host.python.operator.penalization import PythonPenalizeVorticity
@@ -26,7 +27,7 @@ class PenalizeVorticity(ComputationalGraphNodeFrontend):
using penalization.
"""
__implementations = {
- Implementation.PYTHON: PythonPenalizeVorticity
+ Implementation.PYTHON: PythonPenalizeVorticity
}
@classmethod
@@ -40,7 +41,7 @@ class PenalizeVorticity(ComputationalGraphNodeFrontend):
@debug
def __init__(self, obstacles, variables,
velocity, vorticity,
- dt, coeff=None,
+ dt, coeff=None, ubar=None, formulation=None,
implementation=None, **kwds):
"""
Parameters
@@ -53,6 +54,10 @@ class PenalizeVorticity(ComputationalGraphNodeFrontend):
output vorticity
coeff : ScalarParameter, optional
penalization factor (\f\lambda\f) applied to all geometries.
+ ubar : TensorParameter, optional
+ Solid velocity (default to 0)
+ formulation: PenalizationFormulation
+ Solving penalization either with IMPLICIT scheme or EXACT solution
variables: dict
dictionary of fields as keys and topologies as values.
dt: ScalarParameter
@@ -82,24 +87,29 @@ class PenalizeVorticity(ComputationalGraphNodeFrontend):
Warning : coeff as a function is not yet implemented!!
"""
- obstacles = to_list(obstacles)
- assert len(set(obstacles)) == len(obstacles)
- obstacles = tuple(obstacles)
+ if not isinstance(obstacles, dict):
+ obstacles = to_list(obstacles)
+ assert len(set(obstacles)) == len(obstacles)
+ obstacles = tuple(obstacles)
check_instance(velocity, Field)
check_instance(vorticity, Field)
check_instance(variables, dict, keys=Field,
values=CartesianTopologyDescriptors)
check_instance(dt, ScalarParameter)
- check_instance(coeff, (ScalarParameter, float, type(lambda x:x)), allow_none=True)
+ check_instance(coeff, (ScalarParameter, float, type(lambda x: x)), allow_none=True)
+ check_instance(formulation, PenalizationFormulation, allow_none=True)
+ check_instance(ubar, TensorParameter, allow_none=True)
check_instance(obstacles, (tuple, dict), values=Field,
- keys=(ScalarParameter, float), check_kwds=False)
+ keys=(ScalarParameter, float, type(lambda x: x)), check_kwds=False)
super(PenalizeVorticity, self).__init__(
velocity=velocity,
vorticity=vorticity,
coeff=coeff,
+ ubar=ubar,
obstacles=obstacles,
dt=dt,
+ formulation=formulation,
variables=variables,
implementation=implementation,
**kwds)