Improve penalization with exact formulation

hysop/backend/host/python/operator/penalization.py

+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()

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"""

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)