diff --git a/hysop/backend/host/python/operator/penalization.py b/hysop/backend/host/python/operator/penalization.py
index fddd21532dec3e3dd81a05ad11092098efe7e6ac..747f2118473e2544bfb4b0bbdb21868fe9c71fa6 100644
--- a/hysop/backend/host/python/operator/penalization.py
+++ b/hysop/backend/host/python/operator/penalization.py
@@ -13,38 +13,11 @@ from hysop.numerics.stencil.stencil_generator import StencilGenerator, CenteredS
import numpy as np
-class PythonPenalizeVorticity(HostOperator):
- """
-
- """
- __default_method = {
- SpaceDiscretization: SpaceDiscretization.FDC4
- }
-
- __available_methods = {
- SpaceDiscretization: (InstanceOf(SpaceDiscretization), InstanceOf(int))
- }
-
- @classmethod
- def default_method(cls):
- dm = super(PythonPenalizeVorticity, cls).default_method()
- dm.update(cls.__default_method)
- return dm
-
- @classmethod
- def available_methods(cls):
- am = super(PythonPenalizeVorticity, cls).available_methods()
- am.update(cls.__available_methods)
- return am
-
+class CommonPenalization(object):
@debug
- def __init__(self, obstacles, variables,
- velocity, vorticity,
+ def __init__(self, obstacles,
dt, coeff=None, ubar=None, formulation=None, **kwds):
- check_instance(velocity, Field)
- check_instance(vorticity, Field)
- check_instance(variables, dict, keys=Field,
- values=CartesianTopologyDescriptors)
+ super(CommonPenalization, self).__init__()
check_instance(dt, ScalarParameter)
check_instance(coeff, (ScalarParameter, float, type(lambda x: x)), allow_none=True)
check_instance(ubar, TensorParameter, allow_none=True)
@@ -52,11 +25,6 @@ class PythonPenalizeVorticity(HostOperator):
check_instance(obstacles, (tuple, dict), values=Field,
keys=(ScalarParameter, float, type(lambda x: x)), check_kwds=False)
- input_fields = {velocity: variables[velocity],
- vorticity: variables[vorticity]}
- output_fields = {vorticity: variables[vorticity]}
- input_params = {dt.name: dt}
-
if isinstance(coeff, ScalarParameter):
self.coeff = lambda o: coeff()*o
elif isinstance(coeff, type(lambda x: x)):
@@ -81,7 +49,6 @@ class PythonPenalizeVorticity(HostOperator):
obs[self.coeff] = o
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(
@@ -93,13 +60,69 @@ class PythonPenalizeVorticity(HostOperator):
else:
self._compute_penalization = self._compute_penalization_exact
+ self.obstacles = obs
+
+ 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.
+
+
+class PythonPenalizeVorticity(HostOperator, CommonPenalization):
+ """
+
+ """
+ __default_method = {
+ SpaceDiscretization: SpaceDiscretization.FDC4
+ }
+
+ __available_methods = {
+ SpaceDiscretization: (InstanceOf(SpaceDiscretization), InstanceOf(int))
+ }
+
+ @classmethod
+ def default_method(cls):
+ dm = super(PythonPenalizeVorticity, cls).default_method()
+ dm.update(cls.__default_method)
+ return dm
+
+ @classmethod
+ def available_methods(cls):
+ am = super(PythonPenalizeVorticity, cls).available_methods()
+ am.update(cls.__available_methods)
+ return am
+
+ @debug
+ def __init__(self, obstacles, variables,
+ velocity, vorticity,
+ dt, coeff=None, ubar=None, formulation=None, **kwds):
+ check_instance(velocity, Field)
+ check_instance(vorticity, Field)
+ check_instance(variables, dict, keys=Field,
+ values=CartesianTopologyDescriptors)
+
+ input_fields = {velocity: variables[velocity],
+ vorticity: variables[vorticity]}
+ output_fields = {vorticity: variables[vorticity]}
+ input_params = {dt.name: dt}
+ for o in (obstacles.values() if isinstance(obstacles, dict) else obstacles):
+ input_fields[o] = variables[o]
+
self.velocity = velocity
self.vorticity = vorticity
- self.obstacles = obs
self.dt = dt
super(PythonPenalizeVorticity, self).__init__(
input_fields=input_fields, output_fields=output_fields,
- input_params=input_params, **kwds)
+ input_params=input_params,
+ obstacles=obstacles, dt=dt, coeff=coeff, ubar=ubar, formulation=formulation, **kwds)
@debug
def handle_method(self, method):
@@ -185,19 +208,6 @@ 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)
@@ -205,7 +215,7 @@ class PythonPenalizeVorticity(HostOperator):
# Penalize velocity
self._compute_penalization()
- for ubar, v, tmp in zip(self._ubar(), self.V, self.tmp_v):
+ for ubar, v, tmp in zip(self._ubar()[::-1], self.V, self.tmp_v):
tmp[...] = (v - ubar) * self.tmp
self._compute_vorticity()
@@ -246,7 +256,7 @@ class PythonPenalizeVorticity(HostOperator):
self.W[0][...] += -self.tmp
-class PythonPenalizeVelocity(HostOperator):
+class PythonPenalizeVelocity(HostOperator, CommonPenalization):
"""
"""
@@ -277,108 +287,32 @@ class PythonPenalizeVelocity(HostOperator):
check_instance(velocity, 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(ubar, TensorParameter, allow_none=True)
- check_instance(formulation, PenalizationFormulation, allow_none=True)
- check_instance(obstacles, (tuple, dict), values=Field,
- keys=(ScalarParameter, float, type(lambda x: x)), check_kwds=False)
input_fields = {velocity: variables[velocity], }
output_fields = {velocity: variables[velocity]}
input_params = {dt.name: dt}
-
- if isinstance(coeff, ScalarParameter):
- self.coeff = lambda o: coeff()*o
- elif isinstance(coeff, type(lambda x: x)):
- self.coeff = coeff
- elif not coeff is None:
- c = ScalarParameter("penal_coeff", initial_value=coeff)
- self.coeff = lambda o: c()*o
-
- if isinstance(obstacles, dict):
- 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:
- obs[self.coeff] = o
- for o in obs.values():
- assert o.nb_components == 1
+ for o in (obstacles.values() if isinstance(obstacles, dict) else obstacles):
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.obstacles = obs
self.dt = dt
super(PythonPenalizeVelocity, self).__init__(
input_fields=input_fields, output_fields=output_fields,
- input_params=input_params, **kwds)
-
- @debug
- def handle_method(self, method):
- super(PythonPenalizeVelocity, self).handle_method(method)
- self.space_discretization = method.pop(SpaceDiscretization)
- csg = CenteredStencilGenerator()
- csg.configure(dtype=MPQ, dim=1)
- try:
- self.order = int(str(self.space_discretization)[3:])
- except:
- self.order = self.space_discretization
- self.stencil = tuple(csg.generate_exact_stencil(
- derivative=1,
- order=self.order) for _ in range(3))
-
- @debug
- def get_field_requirements(self):
- requirements = super(PythonPenalizeVelocity,
- 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():
- min_ghosts = (max(g, G) for g in req.min_ghosts.copy())
- req.min_ghosts = min_ghosts
- req.axes = (tuple(range(field.dim)), )
- return requirements
+ input_params=input_params,
+ obstacles=obstacles, dt=dt, coeff=coeff, ubar=ubar, formulation=formulation, **kwds)
@debug
def discretize(self):
if self.discretized:
return
super(PythonPenalizeVelocity, self).discretize()
- dim = self.velocity.dim
self.dvelocity = self.input_discrete_tensor_fields[self.velocity]
dv = self.dvelocity
- stencil = self.stencil[0]
- G = max(max(a, b) for a, b in zip(stencil.L, stencil.R))
- view = dv.local_slices(ghosts=(G, )*dim)
- V = tuple(Vi[view] for Vi in dv.buffers)
- self.V = V
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, )*dim)]
-
- for s, dx in zip(self.stencil, dv.space_step):
- s.replace_symbols({s.dx: dx})
+ self.dobstacles[c] = o_df.data[0]
self.v_ghost_exchanger = dv.build_ghost_exchanger()
@@ -386,14 +320,13 @@ class PythonPenalizeVelocity(HostOperator):
def get_work_properties(self):
requests = super(PythonPenalizeVelocity, self).get_work_properties()
buffers = MemoryRequest.empty_like(
- a=self.V[0], nb_components=4, backend=self.dvelocity.backend)
+ a=self.dvelocity.buffers[0], nb_components=1, backend=self.dvelocity.backend)
requests.push_mem_request('Wtmp', buffers)
return requests
def setup(self, work):
super(PythonPenalizeVelocity, self).setup(work=work)
Wtmp = work.get_buffer(self, 'Wtmp', handle=True)
- self.tmp_v = Wtmp[0:3]
self.tmp = Wtmp[-1]
@classmethod
@@ -404,19 +337,6 @@ class PythonPenalizeVelocity(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(PythonPenalizeVelocity, self).apply(**kwds)
@@ -424,6 +344,6 @@ class PythonPenalizeVelocity(HostOperator):
# Penalize velocity
self._compute_penalization()
- for ubar, v, tmp in zip(self._ubar(), self.V, self.tmp_v):
- tmp[...] = (v - ubar) * self.tmp
+ for ubar, v in zip(self._ubar()[::-1], self.dvelocity.buffers):
+ v[...] = (v - ubar) * self.tmp
self.v_ghost_exchanger()
diff --git a/hysop/problem.py b/hysop/problem.py
index 877ede4d7ec4f6e28dfa25a4efcd6acd99b7ce43..6e16c55258b27f5b2c48ad63e95364b1a92235fc 100644
--- a/hysop/problem.py
+++ b/hysop/problem.py
@@ -148,7 +148,8 @@ class Problem(ComputationalGraph):
if all((df in initialized) for df in dfield.discrete_fields()):
# all contained scalar fields were already initialized
continue
- elif mpi_params and mpi_params.task_id != dfield.topology.task_id:
+ elif mpi_params and not all([mpi_params.task_id == df.topology.task_id
+ for df in dfield.discrete_fields()]):
# Topology task does not matches given mpi_params task
continue
else: