temp drag and lift. Untest, unstable

HySoP/hysop/domain/subsets/boxes.py

@@ -28,7 +28,10 @@ class SubBox(Subset):
         self.length = npw.asrealarray(length).copy()
         ## coordinates of the 'highest' point of this subset
         self.end = self.origin + self.length
+        ## directions of coordinate axes belonging to the subset
         self.s_dir = np.where(self.length != 0)[0]
+        ## directions of coordinate axes normal to the subset
+        self.n_dir = np.where(self.length == 0)[0]
         msg = 'Subset error, the origin is outside of the domain.'
         assert ((self.origin - self._parent.origin) >= 0).all(), msg
         msg = 'Subset error, the subset is too large for the domain.'
@@ -62,6 +65,7 @@ class SubBox(Subset):
         self.real_orig[topo] = self.mesh[topo].global_origin
         self.on_proc[topo] = self.mesh[topo].on_proc
         self.ind[topo] = self.mesh[topo].iCompute
+        self._reduce_topology(topo)
         return self.ind[topo]
 
     def integrate_field_on_proc(self, field, topo, component=0):

HySoP/hysop/domain/subsets/control_box.py

@@ -48,6 +48,7 @@ class ControlBox(SubBox):
 
         for s in self.surf:
             s.discretize(topo)
+        return self.ind[topo]
 
     def _check_boundaries(self, surf, topo):
         """

HySoP/hysop/domain/subsets/cylinder.py

@@ -59,6 +59,7 @@ class Cylinder(Subset):
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
+            self._reduce_topology(topo)
         return self.ind[topo]
 
     def chi(self, *args):
@@ -110,4 +111,6 @@ class HemiCylinder(Cylinder):
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
+            self._reduce_topology(topo)
+
         return self.ind[topo]

HySoP/hysop/domain/subsets/sphere.py

@@ -49,6 +49,8 @@ class Sphere(Subset):
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
+            self._reduce_topology(topo)
+
         return self.ind[topo]
 
     def __str__(self):
@@ -100,6 +102,8 @@ class HemiSphere(Sphere):
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
+            self._reduce_topology(topo)
+
         return self.ind[topo]
 
 
@@ -146,4 +150,5 @@ class Porous(Subset):
         self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
                                           for i in xrange(dim)])
                               != 0).all()
+        self._reduce_topology(topo)
         return self.ind[topo]

HySoP/hysop/domain/subsets/subset.py

@@ -62,6 +62,11 @@ class Subset(object):
         self.is_porous = False
         # list of space direction, used for integration.
         self.s_dir = [d for d in xrange(parent.dimension)]
+        ## dict of mpi communicators, keys=topo
+        ## such that self.subcomm[topo] represents the mpi processes
+        ## where on_proc[topo] is true. Useful to reduce collective comm
+        ## on the subset
+        self.subcomm = {}
 
     def chi(self, *args):
         return self._is_inside(*args) <= 0
@@ -79,6 +84,16 @@ class Subset(object):
             self.on_proc[topo] = (np.asarray([len(self.ind[topo][i])
                                               for i in xrange(dim)])
                                   != 0).all()
+            self._reduce_topology(topo)
+        return self.ind[topo]
+
+    def _reduce_topology(self, topo):
+        plist = np.asarray(topo.comm.allgather(self.on_proc[topo]),
+                           dtype=np.bool)
+        gtopo = topo.comm.Get_group()
+        rks = np.where(plist)[0]
+        subgroup = gtopo.Incl(rks)
+        self.subcomm[topo] = topo.comm.Create(subgroup)
 
     def apply(self, field):
         """

HySoP/hysop/numerics/utils.py

@@ -10,23 +10,40 @@ class Utils(object):
     gen_indices = zip(i1, i2)
 
     @staticmethod
-    def sum_cross_product(x, y, sl, work, subsets=None):
+    def sum_cross_product(x, y, sl, work):
         """
         @param x : a tuple representing grid coordinates (like coords in mesh)
         @param y : a discrete field (i.e. list of numpy arrays)
         @param work : numpy array
         Sum on a volume (defined by sl) of cross products of
-        x with y at each grid point
-        . Result in work.
+        x with y at each grid point. Result in work.
         """
         current_dir = 0
-        res = npw.zeros(3)
+        dim = len(y)
+        res = npw.zeros(dim)
         for (i, j) in Utils.gen_indices:
             work[sl] = x[i] * y[j][sl]
             work[sl] -= x[j] * y[i][sl]
-            for sub in subsets:
-                work[sub] = 0.0
+            res[current_dir] = npw.real_sum(work[sl])
+            current_dir += 1
+        return res
 
+    @staticmethod
+    def noca_1(x, v, w, sl, work):
+        """
+        @param x : a tuple representing grid coordinates (like coords in mesh)
+        @param y : a discrete field (i.e. list of numpy arrays)
+        @param work : numpy array
+        Sum on a volume (defined by sl) of cross products of
+        x with y at each grid point. Result in work.
+        """
+        current_dir = 0
+        dim = len(w)
+        res = npw.zeros(dim)
+        for (i, j) in Utils.gen_indices:
+            work[sl] = x[i] * w[j][sl]
+            work[sl] -= x[j] * w[i][sl]
+            work[sl] *= v[sl]
             res[current_dir] = npw.real_sum(work[sl])
             current_dir += 1
         return res

HySoP/hysop/operator/discrete/drag_and_lift.py

@@ -3,24 +3,23 @@
 @file compute_forces.py
 Compute forces
 """
-from parmepy.constants import np, XDIR, YDIR, ZDIR, debug
 from parmepy.numerics.update_ghosts import UpdateGhosts
-from parmepy.numerics.differential_operations import Laplacian
-from parmepy.numerics.finite_differences import FD_C_2
 from parmepy.operator.discrete.discrete import DiscreteOperator
 import parmepy.tools.numpywrappers as npw
-from parmepy.tools.profiler import profile
+from abc import ABCMeta, abstractmethod
 from parmepy.numerics.utils import Utils
 
 
-class DragAndLift(DiscreteOperator):
+class Forces(DiscreteOperator):
     """
     Compute drag and lift using Noca's formula.
     See Noca99 or Plouhmans, 2002, Journal of Computational Physics
     The present class implements formula (52) of Plouhmans2002.
     Integral inside the obstacle is not taken into account.
     """
-    def __init__(self, velocity, subsets, **kwds):
+    __metaclass__ = ABCMeta
+
+    def __init__(self, obstacles, **kwds):
         """
         @param velocity field
         @param vorticity field
@@ -32,43 +31,65 @@ class DragAndLift(DiscreteOperator):
         @param obstacles a list of parmepy.domain.subsets inside the volume of
         control
         """
-        assert 'variables' not in kwds, 'variables parameter is useless.'
-        ## discrete velocity field
-        self.velocity = velocity
-        super(DragAndLift, self).__init__(variables=[velocity], **kwds)
-        self._voc = subsets
+        super(Forces, self).__init__(**kwds)
+        # deal with obstacles, volume of control ...
+        self._indices = self._init_indices(obstacles)
+
         self._dim = self.domain.dimension
         msg = 'Force computation undefined for domain of dimension 1.'
         assert self._dim > 1, msg
+
         # Local buffers, used for time-derivative computation
         self._previous = npw.zeros(self._dim)
         self._buffer = npw.zeros(self._dim)
-        ## The computed forces
+        ## The force we want to compute (lift(s) and drag)
         self.force = npw.zeros(self._dim)
-        # Ghost points synchronizer
-        self._synchronize = None
-        self._topology = self.velocity.topology
+
+        # topology common to all variables
+        self._topology = self.input[0].topology
         assert (self._topology.ghosts() >= 1).all()
-        # Get the list of indices of points used to compute the
-        # integrals over the volume.
-        self._indices = self._voc.ind[self._topology]
-        nbc = len(self.velocity.data)
+
+        # list of np arrays to be synchronized
+        self._datalist = []
+        for v in self.input:
+            self._datalist += v.data
+        nbc = len(self._datalist)
+        # Ghost points synchronizer
         self._synchronize = UpdateGhosts(self._topology, nbc)
-        # Which formula must be used to compute the forces ...
-        self._formula = self._switch_method()
-        self._on_proc = self._voc.on_proc[self._topology]
 
-    def _switch_method(self):
+        # Which formula must be used to compute the forces.
+        # Must be set in derived class.
+        self._formula = lambda dt: 0
+
+        # true if the operator needs to work on the current process.
+        # Updated in derived class.
+        self._on_proc = True
+        self.input = self.variables
+        # Set how reduction will be performed
+        # Default = reduction over all process.
+        # \todo : add param to choose this option
+        self.mpi_sum = self._mpi_allsum
+        # A 'reduced' communicator used to mpi-reduce the force.
+        # Set in derived class
+        self._subcomm = None
+
+    @abstractmethod
+    def _init_indices(self, obstacles):
         """
-        select the proper method to compute drag and lift
+        @param obstacles: a list of parmepy.domain.subsets.subset.Subset
+        @return : a list of np arrays representing points indices
+        (like result from np.where)
+        Discretize obstacles, volume of control ... and
+        compute a list of points representing these sets.
+        What is inside indices depends on the chosen method.
+        See derived class 'init_indices' function for details.
         """
-        # Heloise formula
-        return self._momentum
+        pass
 
     def _set_work_arrays(self, rwork=None, iwork=None):
 
-        v_ind = self.velocity.topology.mesh.iCompute
-        shape_v = self.velocity.data[0][v_ind].shape
+        v_ind = self.input[0].topology.mesh.iCompute
+        shape_v = self.input[0].data[0][v_ind].shape
         # setup for rwork, iwork is useless.
         if rwork is None:
             # ---  Local allocation ---
@@ -85,7 +106,7 @@ class DragAndLift(DiscreteOperator):
         Performs MPI reduction (sum result value over all process)
         All process get the result of the sum.
         """
-        self.force = self._topology.comm.allreduce(self.force)
+        self.force = self._subcomm.allreduce(self.force)
 
     def _mpi_sum(self, root=0):
         """
@@ -93,31 +114,270 @@ class DragAndLift(DiscreteOperator):
         Result send only to 'root' process.
         @param root : number of the process which get the result.
         """
-        self.force = self._topology.comm.reduce(self.force, root=root)
+        self.force = self._subcomm.reduce(self.force, root=root)
 
-    def _momentum(self, dt):
+    def apply(self, simulation=None):
+        """
+        Perform integrals on volume and surfaces of the control box
+        @param parmepy.problem.simulation : object describing
+        simulation parameters
+        """
+        assert simulation is not None,\
+            "Simulation parameter is required for Forces apply."
         # Synchro of ghost points is required for fd schemes
-        self._synchronize(self.velocity.data)
-
-        # -- Integration over the volume of control --
+        self._synchronize(self._datalist)
+        # Compute forces locally
+        dt = simulation.timeStep
         if not self._on_proc:
             self._buffer[...] = 0.0
             self.force[...] = 0.0
         else:
-            self._buffer = self._voc.integrate_dfield_allc(self.velocity)
-            print self._buffer
-            self.force[...] = 1. / dt * (self._buffer - self._previous)
+            self._formula(dt)
+        # Reduce results over MPI processes
+        self.mpi_sum()
+        # Print results, if required
+        ite = simulation.currentIteration
+        if self._writer is not None and self._writer.do_write(ite):
+            self._writer.buffer[0, 0] = simulation.time
+            self._writer.buffer[0, 1:] = self.force
+            self._writer.write()
+
+
+class MomentumForces(Forces):
+    """
+    Compute drag and lift using Noca's formula.
+    See Noca99 or Plouhmans, 2002, Journal of Computational Physics
+    The present class implements formula (52) of Plouhmans2002.
+    Integral inside the obstacle is not taken into account.
+    """
+    def __init__(self, velocity, **kwds):
+        """
+        @param velocity field
+        @param vorticity field
+        @param nu viscosity
+        @param a volume of control
+        @param normalization : a normalization coefficient applied to the force
+        (default = 1.)
+        (parmepy.domain.subset.boxes.SubBox)
+        @param obstacles a list of parmepy.domain.subsets inside the volume of
+        control
+        """
+        assert 'variables' not in kwds, 'variables parameter is useless.'
+        ## discrete velocity field
+        self.velocity = velocity
+        #
+        self._obstacle = None
+        super(MomentumForces, self).__init__(variables=[velocity], **kwds)
+        self._formula = self._momentum
+        self._subcomm = self._obstacle.subcomm[self._topology]
+
+    def _init_indices(self, obstacles):
+        msg = 'obstacles arg must be a list.'
+        assert isinstance(obstacles, list), msg
+        # only one obstacle allowed for the moment
+        assert len(obstacles) == 1
+        self._obstacle = obstacles[0]
+        self._obstacle.discretize(self._topology)
+        self._on_proc = self._obstacle.on_proc[self._topology]
+        return self._obstacle.ind[self._topology]
+
+    def _momentum(self, dt):
+        # -- Integration over the volume of control --
+        nbc = self.velocity.nbComponents
+        self._buffer[...] = \
+            [self._obstacle.integrate_dfield_on_proc(self.velocity,
+                                                     component=d)
+             for d in xrange(nbc)]
+        self.force[...] = 1. / dt * (self._buffer - self._previous)
 
         # Update previous for next time step ...
         self._previous[...] = self._buffer[...]
 
-    def apply(self, simulation=None):
+
+class NocaForces(Forces):
+    """
+    Compute drag and lift using Noca's formula.
+    See Noca99 or Plouhmans, 2002, Journal of Computational Physics
+    The present class implements formula (52) of Plouhmans2002.
+    Integral inside the obstacle is not taken into account.
+    """
+
+    def __init__(self, velocity, vorticity, nu, volume_of_control,
+                 normalization=1., **kwds):
+        """
+        @param vorticity field
+        @param nu viscosity
+        @param a volume of control
+        @param normalization : a normalization coefficient applied to the force
+        (default = 1.)
+        (parmepy.domain.subset.boxes.SubBox)
+        @param obstacles a list of parmepy.domain.subsets inside the volume of
+        control
+        """
+        # A volume of control, in which forces are computed
+        self._voc = volume_of_control
+        ## discrete velocity field
+        self.velocity = velocity
+        ## discrete vorticity field
+        self.vorticity = vorticity
+
+        assert 'variables' not in kwds, 'variables parameter is useless.'
+        super(NocaForces, self).__init__(variables=[velocity, vorticity],
+                                         **kwds)
+        # Coef in the Noca formula
+        self._coeff = 1. / (self._dim - 1)
+        ## viscosity
+        self.nu = nu
+        # Normalizing coefficient for forces
+        # (based on the physics of the flow)
+        self._normalization = normalization
+        #
+        self._dvol = npw.prod(self._topology.mesh.space_step)
+        # function to compute the laplacian of a
+        # scalar field. Default fd scheme. (See Laplacian)
+        from parmepy.numerics.differential_operations import Laplacian
+        self._laplacian = Laplacian(self._topology)
+        # function used to compute first derivative of
+        # a scalar field in a given direction.
+        # Default = FD_C_2. Todo : set this as an input method value.
+        from parmepy.numerics.finite_differences import FD_C_2
+        self._fd_scheme = FD_C_2(self._topology.mesh.space_step)
+
+        self._formula = self._noca
+        self._subcomm = self._voc.subcomm[self._topology]
+
+    def _init_indices(self, obstacles):
+        """
+        Compute a list of indices corresponding to points inside
+        the volume of control minus those inside the obstacles
+        """
+        from parmepy.domain.subsets.control_box import ControlBox
+        assert isinstance(self._voc, ControlBox)
+        self._on_proc = self._voc.on_proc[self._topology]
+        msg = 'obstacles arg must be a list.'
+        assert isinstance(obstacles, list), msg
+        from parmepy.domain.subsets.subset import Subset
+        # no obstacle in the box, just for test purpose.
+        if obstacles is None or len(obstacles) == 0:
+            return self._voc.ind[self._topology]
+        else:
+            for obs in obstacles:
+                obs.discretize(self._topology)
+            # return indices of points inside the box, excluding points
+            # inside the obstacles.
+            return Subset.substract_list_of_sets([self._voc], obstacles,
+                                                 self._topology)
+
+    def _set_work_arrays(self, rwork=None, iwork=None):
+
+        v_ind = self.velocity.topology.mesh.iCompute
+        shape_v = self.velocity.data[0][v_ind].shape
+        # setup for rwork, iwork is useless.
+        if rwork is None:
+            # ---  Local allocation ---
+            self._rwork = [npw.zeros(shape_v)]
+        else:
+            assert isinstance(rwork, list), 'rwork must be a list.'
+            # --- External rwork ---
+            self._rwork = rwork
+            assert len(self._rwork) == 1
+            assert self._rwork[0].shape == shape_v
+
+    def _noca(self, dt):
         """
         Perform integrals on volume and surfaces of the control box
-        @param parmepy.problem.simulation : object describing
-        simulation parameters
+        @param dt : current time step
+        @return the local (i.e. current mpi process) force
         """
-        assert simulation is not None,\
-            "Simulation parameter is required for DragAndLift apply."
-        dt = simulation.timeStep
-        self._formula(dt)
+        # -- Integration over the volume of control --
+        # -1/(N-1) . d/dt int(x ^ w)
+        self._buffer = self._integrateOnBox(self._buffer)
+        self.force[...] = -1. / dt * self._coeff * (self._buffer
+                                                    - self._previous)
+
+        # Update previous for next time step ...
+        self._previous[...] = self._buffer[...]
+        # -- Integrals on surfaces --
+        # Only on surf. normal to XDIR
+        s_x0 = self._voc.surf[0]
+        s_xn = self._voc.surf[1]
+        self._buffer = self._integrateOnSurface(s_x0, self._buffer, 1)
+        self.force += self._buffer
+        self._buffer = self._integrateOnSurface(s_xn, self._buffer, -1)
+        self.force += self._buffer
+
+        self.force *= self._normalization
+
+    def _integrateOnSurface(self, surf, res, normal):
+        res[...] = 0.0
+        if not surf.on_proc[self._topology]:
+            return res
+
+        # Get indices of the surface
+        sl = surf.slices[self._topology]
+        coords = surf.coords[self._topology]
+        vd = self.velocity.data
+        # i_n : normal dir
+        # i2 : other dirs
+        i_n = surf.n_dir
+        i_t = surf.s_dir
+        dsurf = npw.prod(self._topology.mesh.space_step[i_t])
+        # Indices used for cross-product
+        j1 = [YDIR, ZDIR, XDIR]
+        j2 = [ZDIR, XDIR, YDIR]
+
+        # i_n component
+        res[i_n] = 0.5 * normal * npw.real_sum(
+            sum([vd[j][sl] ** 2 for j in i_t]) - vd[i_n][sl] ** 2)
+        # other components
+        for j in i_t:
+            res[j] = - normal * npw.real_sum(vd[i_n][sl] * vd[j][sl])
+
+        # Second part of integral on surface ...
+        wd = self.vorticity.data
+        buff = self._rwork[i_n]
+        for j in i_t:
+            buff[...] = vd[j1[j]][sl] * wd[j2[j]][sl]\
+                - vd[j2[j]][sl] * wd[j1[j]][sl]
+            res[i_n] -= self._coeff * normal[i_n] * npw.real_sum(coords[j]
+                                                               * buff)
+            res[j] -= self._coeff * normal[i_n] * coords[i_n] *\
+                npw.real_sum(buff)
+
+        # Last part
+        # Update fd schemes in order to compute laplacian and other derivatives
+        # only on the surface (i.e. for list of indices in sl)
+        self._laplacian.fd_scheme.computeIndices(sl)
+        for j in i_t:
+            self._rwork[...] = self._laplacian(vd[j], self._rwork)
+            res[i_n] += self._coeff * self.nu * normal[i_n]\
+                * npw.real_sum(coords[j] * self._rwork[sl])
+            res[j] -= self._coeff * self.nu * normal[i_n] * coords[i_n] * \
+                npw.real_sum(self._rwork[sl])
+        self._fd_scheme.computeIndices(sl)
+        self._fd_scheme.compute(vd[i_n], i_n, self._rwork)
+        res[i_n] += 2.0 * normal[i_n] * self.nu * npw.real_sum(self._rwork[sl])
+        for j in i_t:
+            self._fd_scheme.compute(vd[i_n], j, self._rwork)
+            res[j] += normal[i_n] * self.nu * npw.real_sum(self._rwork[sl])
+            self._fd_scheme.compute(vd[j], i_n, self._rwork)
+            res[j] += normal[i_n] * self.nu * npw.real_sum(self._rwork[sl])
+
+        res *= dsurf
+        return res
+
+    def _integrateOnBox(self, res):
+        """
+        Compute the integral over the volume of control
+        of \f$ x ^ \omega \f$
+        """
+        assert self._dim == 3, 'Not defined for dim < 3'
+        # coords of the points in the volume of control
+        coords = self._voc.coords[self._topology]
+        # slices representing the points inside the subset
+        # (global indices relative to the global mesh of the current topology)
+        sl = self._voc.ind[self._topology]
+        res = Utils.sum_cross_product(coords, self.vorticity.data,
+                                      sl, self._rwork)
+        res *= self._dvol
+        return res

HySoP/hysop/operator/discrete/drag_and_lift_noca.py

-# -*- coding: utf-8 -*-
-"""
-@file compute_forces.py
-Compute forces
-"""
-from parmepy.constants import np, XDIR, YDIR, ZDIR, debug
-from parmepy.numerics.update_ghosts import UpdateGhosts
-from parmepy.numerics.differential_operations import Laplacian
-from parmepy.numerics.finite_differences import FD_C_2
-from parmepy.operator.discrete.discrete import DiscreteOperator
-import parmepy.tools.numpywrappers as npw
-from parmepy.tools.profiler import profile
-from parmepy.numerics.utils.Utils import Utils
-
-
-class DragAndLift(DiscreteOperator):
-    """
-    Compute drag and lift using Noca's formula.
-    See Noca99 or Plouhmans, 2002, Journal of Computational Physics
-    The present class implements formula (52) of Plouhmans2002.
-    Integral inside the obstacle is not taken into account.
-    """
-    def __init__(self, velocity, vorticity, nu, volumeOfControl,
-                 normalization=1., obstacles=None, **kwds):
-        """
-        @param velocity field
-        @param vorticity field
-        @param nu viscosity
-        @param a volume of control
-        @param normalization : a normalization coefficient applied to the force
-        (default = 1.)
-        (parmepy.domain.subset.boxes.SubBox)
-        @param obstacles a list of parmepy.domain.subsets inside the volume of
-        control
-        """
-        assert 'variables' not in kwds, 'variables parameter is useless.'
-        super(DragAndLift, self).__init__(variables=[velocity, vorticity],
-                                          **kwds)
-
-        ## discrete velocity field
-        self.velocity = velocity
-        ## discrete vorticity field
-        self.vorticity = vorticity
-        # volume on which forces are computed
-        self._voc = volumeOfControl
-        self._dim = self.domain.dimension
-        msg = 'Force computation undefined for domain of dimension 1.'
-        assert self._dim > 1, msg
-        # Local buffers, used for time-derivative computation
-        self._previous = npw.zeros(self._dim)
-        self._buffer = npw.zeros(self._dim)
-        ## The computed forces
-        self.force = npw.zeros(self._dim)
-        # Coef in the Noca formula
-        self._coeff = 1. / (self._dim - 1)
-        # Set how reduction will be performed
-        # Default = reduction over all process.
-        # \todo : add param to choose this option
-        self.mpi_sum = self._mpi_allsum
-        # Ghost points synchronizer
-        self._synchronize = None
-        ## viscosity
-        self.nu = nu
-        # Normalizing coefficient for forces
-        # (based on the physics of the flow)
-        self._normalization = normalization
-        self._topology = self.velocity.topology
-        self._dvol = npw.prod(self._topology.mesh.space_step)
-        assert (self._topology.ghosts() >= 1).all()
-        # function to compute the laplacian of a
-        # scalar field. Default fd scheme. (See Laplacian)
-        self._laplacian = Laplacian(self._topology)
-        # function used to compute first derivative of
-        # a scalar field in a given direction.
-        # Default = FD_C_2. Todo : set this as an input method value.
-        self._fd_scheme = FD_C_2(self._topology.mesh.space_step)
-
-        # Get the list of indices of points used to compute the
-        # integrals over the volume.
-        self._indices = self._init_indices(obstacles)
-
-        # prepare ghost points synchro for velocity and vorticity used
-        # in fd schemes
-        nbc = self.velocity.nbComponents + self.vorticity.nbComponents
-        self._synchronize = UpdateGhosts(self._topology, nbc)
-
-        # Which formula must be used to compute the forces ...
-        self._formula = self._switch_method()
-
-        self._on_proc = self._voc.on_proc[self._topology]
-
-    def _init_indices(self, obstacles):
-        """
-        Compute a list of indices corresponding to points inside
-        the volume of control minus those inside the obstacles
-        """
-        from parmepy.domain.subsets.subset import Subset
-        self._voc.discretize(self._topology)
-        if obstacles is None:
-            return self._voc.ind[self._topology]
-        else:
-            for obs in obstacles:
-                obs.discretize(self._topology)
-            return Subset.substract_list_of_sets([self._voc], obstacles,
-                                                 self._topology)
-
-    def _switch_method(self):
-        """
-        select the proper method to compute drag and lift
-        """
-        # Heloise formula
-        # assert self._dim == 3, 'Not defined for dim < 3'
-        return self._momentum
-
-    def _set_work_arrays(self, rwork=None, iwork=None):
-
-        v_ind = self.velocity.topology.mesh.iCompute
-        shape_v = self.velocity.data[0][v_ind].shape
-        # setup for rwork, iwork is useless.
-        if rwork is None:
-            # ---  Local allocation ---
-            self._rwork = [npw.zeros(shape_v)]
-        else:
-            assert isinstance(rwork, list), 'rwork must be a list.'
-            # --- External rwork ---
-            self._rwork = rwork
-            assert len(self._rwork) == 1
-            assert self._rwork[0].shape == shape_v
-
-    def _mpi_allsum(self):
-        """
-        Performs MPI reduction (sum result value over all process)
-        All process get the result of the sum.
-        """
-        self.force = self._topology.comm.allreduce(self.force)
-
-    def _mpi_sum(self, root=0):
-        """
-        Performs MPI reduction (sum result value over all process)
-        Result send only to 'root' process.
-        @param root : number of the process which get the result.
-        """
-        self.force = self._topology.comm.reduce(self.force, root=root)
-
-    @debug
-    @profile
-    def _noca(self, simulation=None):
-        """
-        Perform integrals on volume and surfaces of the control box
-        @param parmepy.problem.simulation : object describing
-        simulation parameters
-        """
-        assert simulation is not None,\
-            "Simulation parameter is required for DragAndLift apply."
-        dt = simulation.timeStep
-        ite = simulation.currentIteration
-        # Synchro of ghost points is required for fd schemes
-        self._synchronize(self.velocity.data + self.vorticity.data)
-
-        # -- Integration over the volume of control --
-        # -1/(N-1) . d/dt int(x ^ w)
-        if self._voc.isEmpty[self._topology]:
-            self._buffer[...] = 0.0
-            self.force[...] = 0.0
-        else:
-            self._buffer = self._integrateOnBox(self._buffer)
-            self.force[...] = -1. / dt * self._coeff * (self._buffer
-                                                        - self._previous)
-
-        # Update previous for next time step ...
-        self._previous[...] = self._buffer[...]
-        # -- Integrals on surfaces --
-        # Only on surf. normal to XDIR
-        self._buffer = self._integrateOnSurface(self._voc.upperS[XDIR],
-                                                self._buffer)
-        self.force += self._buffer
-        self._buffer = self._integrateOnSurface(self._voc.lowerS[XDIR],
-                                                self._buffer)
-        self.force += self._buffer
-
-        # Reduce results over all MPI process in topo
-        self.mpi_sum()
-        self.force *= self._normalization
-
-        # Print results, if required
-        if self._writer is not None and self._writer.do_write(ite):
-            self._writer.buffer[0, 0] = simulation.time
-            self._writer.buffer[0, 1:] = self.force
-            self._writer.write()
-
-        return self.force
-
-    def _momentum(self, dt):
-        # Synchro of ghost points is required for fd schemes
-        self._synchronize(self.velocity.data + self.vorticity.data)
-
-        # -- Integration over the volume of control --
-        if not self._on_proc:
-            self._buffer[...] = 0.0
-            self.force[...] = 0.0
-        else:
-            self._buffer = self._voc.integrate_dfield_allc(self.velocity)
-            self.force[...] = 1. / dt * (self._buffer - self._previous)
-
-        # Update previous for next time step ...
-        self._previous[...] = self._buffer[...]
-
-    def apply(self, simulation=None):
-        """
-        Perform integrals on volume and surfaces of the control box
-        @param parmepy.problem.simulation : object describing
-        simulation parameters
-        """
-        assert simulation is not None,\
-            "Simulation parameter is required for DragAndLift apply."
-        dt = simulation.timeStep
-        self._formula(dt)
-
-    def _integrateOnSurface(self, surf, res):
-
-        res[...] = 0.0
-        if surf.isEmpty[self._topology]:
-            return res
-
-        # Get normal of the surface
-        normal = surf.normal
-        # Get indices of the surface
-        sl = surf.slices[self._topology]
-        coords = surf.coords[self._topology]
-        vdata = self.velocity.data
-        wdata = self.vorticity.data
-        # i1 : normal dir
-        # i2 : other dirs
-        i1 = np.where(normal)[0][0]
-        i2 = np.where(normal == 0)[0]
-        dsurf = npw.prod(self._topology.mesh.space_step[i2])
-        # Indices used for cross-product
-        j1 = [YDIR, ZDIR, XDIR]
-        j2 = [ZDIR, XDIR, YDIR]
-
-        # i1 component
-        res[i1] = normal[i1] * 0.5 * npw.real_sum((-vdata[i1][sl] ** 2
-                                                   + sum([vdata[j][sl] ** 2
-                                                          for j in i2])))
-        # other components
-        for j in i2:
-            res[j] = -normal[i1] * npw.real_sum(vdata[i1][sl] * vdata[j][sl])
-
-        # Second part of integral on surface ...
-        buff = npw.zeros(vdata[0][sl].shape)
-        for j in i2:
-            buff[...] = vdata[j1[j]][sl] * wdata[j2[j]][sl]\
-                - vdata[j2[j]][sl] * wdata[j1[j]][sl]
-            res[i1] -= self._coeff * normal[i1] * npw.real_sum(coords[j]
-                                                               * buff)
-            res[j] -= self._coeff * normal[i1] * coords[i1] *\
-                npw.real_sum(buff)
-
-        # Last part
-        # Update fd schemes in order to compute laplacian and other derivatives
-        # only on the surface (i.e. for list of indices in sl)
-        self._laplacian.fd_scheme.computeIndices(sl)
-        for j in i2:
-            self._rwork[...] = self._laplacian(vdata[j], self._rwork)
-            res[i1] += self._coeff * self.nu * normal[i1]\
-                * npw.real_sum(coords[j] * self._rwork[sl])
-            res[j] -= self._coeff * self.nu * normal[i1] * coords[i1] * \
-                npw.real_sum(self._rwork[sl])
-        self._fd_scheme.computeIndices(sl)
-        self._fd_scheme.compute(vdata[i1], i1, self._rwork)
-        res[i1] += 2.0 * normal[i1] * self.nu * npw.real_sum(self._rwork[sl])
-        for j in i2:
-            self._fd_scheme.compute(vdata[i1], j, self._rwork)
-            res[j] += normal[i1] * self.nu * npw.real_sum(self._rwork[sl])
-            self._fd_scheme.compute(vdata[j], i1, self._rwork)
-            res[j] += normal[i1] * self.nu * npw.real_sum(self._rwork[sl])
-
-        res *= dsurf
-        return res
-
-    def _integrateOnBox(self, res):
-        """
-        Compute the integral over the volume of control
-        of \f$ x ^ \omega \f$
-        """
-        assert self._dim == 3, 'Not defined for dim < 3'
-        # coords of the points in the volume of control
-        coords = self._voc.coords[self._topology]
-        # slices representing the points inside the subset
-        # (global indices relative to the global mesh of the current topology)
-        sl = self._voc.slices[self._topology]
-        res = Utils.sum_cross_product(coords, self.vorticity.data,
-                                      sl, self._rwork, self._obsinds)
-        res *= self._dvol
-        return res

HySoP/hysop/operator/drag_and_lift.py

@@ -5,50 +5,30 @@ Compute forces
 """
 from parmepy.operator.computational import Computational
 from parmepy.operator.continuous import opsetup
+from abc import ABCMeta, abstractmethod
 
 
-class DragAndLift(Computational):
-
-    Momentum = 1
-    Noca = 2
-
+class Forces(Computational):
     """
-    Compute drag and lift using Noca's formula.
-    See Noca99 or Plouhmans, 2002, Journal of Computational Physics
-    The present class implements formula (52) of Plouhmans2002.
-    Integral inside the obstacle is not taken into account.
+    Abstract interface to classes dedicated to drag/lift computation
+    for a flow around a predefined obstacle.
     """
-    def __init__(self, velocity, subsets, formula=None,
-                 vorticity=None, **kwds):
+
+    __metaclass__ = ABCMeta
+
+    def __init__(self, obstacles, **kwds):
         """
-        @param velocity field
-        @param vorticity field
-        @param nu viscosity
-        @param the topology on which forces will be computed
-        @param a volume of control
-        (parmepy.domain.obstacle.controlBox.ControlBox object)
         @param obstacles a list of parmepy.domain.obstacles inside
         the box
-        @param io_params : parameters (dict) to set file output.
-        If  None, no output. Set io_params = {} if you want output,
-        with default parameters values. Default file name = 'drag_and_lift'
-        See parmepy.tools.io_utils.Writer for details
         """
-        assert 'variables' not in kwds, 'variables parameter is useless.'
-        variables = [velocity]
-        if vorticity is not None:
-            variables.append(vorticity)
-        super(DragAndLift, self).__init__(variables=variables,
-                                          **kwds)
-        self.vorticity = vorticity
-        self.velocity = velocity
+        super(Forces, self).__init__(**kwds)
         self.input = self.variables
-        if isinstance(subsets, list):
-            self._voc = subsets[0]
-        else:
-            self._voc = subsets
-        self.formula = formula
+        ## List of parmepy.domain.subsets, obstacles to the flow
+        self.obstacles = obstacles
+        # Minimal number of ghost points required. Defined in derived class
+        self._min_ghosts = 0
 
+    @abstractmethod
     def get_work_properties(self):
         if not self._is_discretized:
             msg = 'The operator must be discretized '
@@ -60,27 +40,14 @@ class DragAndLift(Computational):
         return {'rwork': [shape_v], 'iwork': None}
 
     def discretize(self):
-        super(DragAndLift, self)._standard_discretize()
+        super(Forces, self)._standard_discretize(self._min_ghosts)
         # all variables must have the same resolution
         assert self._single_topo, 'multi-resolution case not allowed.'
-        topo = self.variables.values()[0]
-        self._voc.discretize(topo)
 
+    @abstractmethod
     @opsetup
     def setup(self, rwork=None, iwork=None):
-        msg = 'Multi-resolution is not allowed.'
-        assert self._single_topo, msg
-        if not self._is_uptodate:
-            from parmepy.operator.discrete.drag_and_lift \
-                import DragAndLift as DAL
-            self.discreteOperator = DAL(
-                velocity=self.discreteFields[self.velocity],
-                subsets=self._voc)
-
-            # output setup
-            self._set_io('drag_and_lift', (1, 4))
-            self.discreteOperator.setWriter(self._writer)
-            self._is_uptodate = True
+        pass
 
     def drag(self):
         """
@@ -93,3 +60,126 @@ class DragAndLift(Computational):
         @return the last computed value for lift
         """
         return self.discreteOperator.force[:self.domain.dimension - 1]
+
+
+class MomentumForces(Forces):
+    """
+    Computation of forces (drag and lift) around an obstacle using
+    Momentum formula
+    \f{eqnarray}
+    force = \int_{obst} \frac{dv}{dt} \\
+    v = velocity, \ obst = volume \ of \ the \ obstacle
+    \f}
+    """
+
+    __metaclass__ = ABCMeta
+
+    def __init__(self, velocity, **kwds):
+        """
+        @param velocity : velocity continuous vector field
+        """
+        assert 'variables' not in kwds, 'variables parameter is useless.'
+        self.velocity = velocity
+        super(MomentumForces, self).__init__(variables=[velocity], **kwds)
+
+    def get_work_properties(self):
+        return {'rwork': None, 'iwork': None}
+
+    @opsetup
+    def setup(self, rwork=None, iwork=None):
+        if not self._is_uptodate:
+            from parmepy.operator.discrete.drag_and_lift \
+                import MomentumForces as DMF
+            self.discreteOperator = DMF(
+                velocity=self.discreteFields[self.velocity],
+                obstacles=self.obstacles)
+
+            # output setup
+            self._set_io('drag_and_lift', (1, 4))
+            self.discreteOperator.setWriter(self._writer)
+            self._is_uptodate = True
+
+
+class NocaForces(Forces):
+    """
+    Computation of forces (drag and lift) around an obstacle using
+    Noca's formula
+    (See Noca99 or Plouhmans, 2002, Journal of Computational Physics)
+    """
+
+    __metaclass__ = ABCMeta
+
+    def __init__(self, velocity, vorticity, nu, volume_of_control=None,
+                 normalization=1., **kwds):
+        """
+        @param velocity field
+        @param vorticity field
+        @param nu viscosity
+        @param a volume of control
+        (parmepy.domain.obstacle.controlBox.ControlBox object)
+        @param normalization : a normalization coefficient applied to the force
+        (default = 1.)
+        """
+        assert 'variables' not in kwds, 'variables parameter is useless.'
+        ## velocity field
+        self.velocity = velocity
+        ## vorticity field
+        self.vorticity = vorticity
+        super(NocaForces, self).__init__(variables=[velocity, vorticity],
+                                         **kwds)
+        ## viscosity
+        self.nu = nu
+        ## Normalizing coefficient for forces
+        ## (based on the physics of the flow)
+        self.normalization = normalization
+
+        # setup for finite differences
+        if self.method is None:
+            import parmepy.default_methods as default
+            self.method = default.FORCES
+        from parmepy.methods_keys import SpaceDiscretisation
+        from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
+        assert SpaceDiscretisation in self.method.keys()
+        if SpaceDiscretisation is FD_C_2:
+            self._min_ghosts = 1
+        elif SpaceDiscretisation is FD_C_4:
+            self._min_ghosts = 2
+
+        # The volume of control, in which forces are computed
+        if volume_of_control is None:
+            xr = self.domain.origin
+            lr = self.domain.length
+        from parmepy.domain.subsets.control_box import ControlBox
+        volume_of_control = ControlBox(parent=self.domain,
+                                       origin=xr, length=lr)
+        self._voc = volume_of_control
+
+    def get_work_properties(self):
+        if not self._is_discretized:
+            msg = 'The operator must be discretized '
+            msg += 'before any call to this function.'
+            raise RuntimeError(msg)
+        vd = self.discreteFields[self.input[0]]
+        v_ind = vd.topology.mesh.iCompute
+        shape_v = vd[0][v_ind].shape
+        return {'rwork': [shape_v], 'iwork': None}
+
+    @opsetup
+    def setup(self, rwork=None, iwork=None):
+        if not self._is_uptodate:
+            from parmepy.operator.discrete.drag_and_lift \
+                import NocaForces as DNF
+            topo = self.discreteFields[self.velocity].topology
+            self._voc.discretize(topo)
+            self.discreteOperator = DNF(
+                velocity=self.discreteFields[self.velocity],
+                vorticity=self.discreteFields[self.vorticity],
+                nu=self.nu,
+                volume_of_control=self._voc,
+                normalization=self.normalization,
+                obstacles=self.obstacles)
+
+            # output setup
+            self._set_io('drag_and_lift', (1, 4))
+            self.discreteOperator.setWriter(self._writer)
+            self._is_uptodate = True

HySoP/hysop/operator/drag_and_lift_noca.py

-# -*- coding: utf-8 -*-
-"""
-@file compute_forces.py
-Compute forces
-"""
-from parmepy.operator.computational import Computational
-from parmepy.operator.continuous import opsetup
-
-
-class DragAndLift(Computational):
-    """
-    Compute drag and lift using Noca's formula.
-    See Noca99 or Plouhmans, 2002, Journal of Computational Physics
-    The present class implements formula (52) of Plouhmans2002.
-    Integral inside the obstacle is not taken into account.
-    """
-    def __init__(self, velocity, vorticity, volumeOfControl,
-                 nu=0., normalization=1., obstacles=None, **kwds):
-        """
-        @param velocity field
-        @param vorticity field
-        @param nu viscosity
-        @param the topology on which forces will be computed
-        @param a volume of control
-        (parmepy.domain.obstacle.controlBox.ControlBox object)
-        @param obstacles a list of parmepy.domain.obstacles inside
-        the box
-        @param io_params : parameters (dict) to set file output.
-        If  None, no output. Set io_params = {} if you want output,
-        with default parameters values. Default file name = 'drag_and_lift'
-        See parmepy.tools.io_utils.Writer for details
-        """
-        assert 'variables' not in kwds, 'variables parameter is useless.'
-        super(DragAndLift, self).__init__(variables=[velocity, vorticity],
-                                          **kwds)
-        if obstacles is None:
-            obstacles = []
-        self.input = [velocity, vorticity]
-        self._discr_args = {velocity: velocity,
-                            vorticity: vorticity, nu: nu,
-                            normalization: normalization,
-                            volumeOfControl: volumeOfControl,
-                            obstacles: obstacles}
-
-    def get_work_properties(self):
-        if not self._is_discretized:
-            msg = 'The operator must be discretized '
-            msg += 'before any call to this function.'
-            raise RuntimeError(msg)
-        vd = self.discreteFields[self.input[0]]
-        v_ind = vd.topology.mesh.iCompute
-        shape_v = vd[0][v_ind].shape
-        return {'rwork': [shape_v], 'iwork': None}
-
-    @opsetup
-    def setup(self, rwork=None, iwork=None):
-        msg = 'Multi-resolution is not allowed.'
-        assert self._single_topo, msg
-        if not self._is_uptodate:
-            from parmepy.operator.discrete.drag_and_lift \
-                import DragAndLift as DAL
-            self.discreteOperator = DAL(**self._discr_args)
-
-            # output setup
-            self._set_io('drag_and_lift', (1, 4))
-            self.discreteOperator.setWriter(self._writer)
-            self._is_uptodate = True
-
-    def drag(self):
-        """
-        @return the last computed value for drag force
-        """
-        return self.discreteOperator.force[2]
-
-    def lift(self):
-        """
-        @return the last computed value for lift
-        """
-        return self.discreteOperator.force[:1]

HySoP/hysop/operator/tests/test_drag_and_lift.py

@@ -9,7 +9,7 @@ import numpy as np
 import os
 from parmepy import Field, Box
 from parmepy.operator.hdf_io import HDF_Reader
-from parmepy.operator.drag_and_lift import DragAndLift
+from parmepy.operator.drag_and_lift import MomentumForces, NocaForces
 
 
 def v2d(res, x, y, t):
@@ -92,7 +92,7 @@ def check_drag(penal, wref, vref, vorti, velo):
         assert np.allclose(wd.data[d][ind], wref.data[d][ind])
 
 
-def test_penal_vort_3D():
+def test_momentum_forces_3D():
     """
     Compute drag/lift in 3D, flow around a sphere
     """
@@ -109,7 +109,7 @@ def test_penal_vort_3D():
 
     penal.discretize()
     penal.setup()
-    dg = DragAndLift(velocity=velo, discretization=topo, subsets=[hsphere])
+    dg = MomentumForces(velocity=velo, discretization=topo, subsets=[hsphere])
     dg.discretize()
     dg.setup()
     simu = Simulation(nbIter=3)