diff --git a/HySoP/hysop/domain/subsets/boxes.py b/HySoP/hysop/domain/subsets/boxes.py
index 1505091dfc9cb398f185ad5586443a5ca7c731be..0beca1b80c6e692bd84f41725e06705317da82a7 100644
--- a/HySoP/hysop/domain/subsets/boxes.py
+++ b/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):
diff --git a/HySoP/hysop/domain/subsets/control_box.py b/HySoP/hysop/domain/subsets/control_box.py
index 1821ac8a886f6ac08ac164c66db03e32687f094f..5d162b561821b08437c65a334f277e98e16130d1 100644
--- a/HySoP/hysop/domain/subsets/control_box.py
+++ b/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):
"""
diff --git a/HySoP/hysop/domain/subsets/cylinder.py b/HySoP/hysop/domain/subsets/cylinder.py
index a4dd309c9a152f8658a0fa450f7c17fa76985df5..1459d7db1d12947f6635a85134b9eab35b45b762 100644
--- a/HySoP/hysop/domain/subsets/cylinder.py
+++ b/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]
diff --git a/HySoP/hysop/domain/subsets/sphere.py b/HySoP/hysop/domain/subsets/sphere.py
index 469b098546ced68fc202ea17c18b1467e473af81..cc48610b43582e12616e195bdcacafbd9e7ca0a7 100644
--- a/HySoP/hysop/domain/subsets/sphere.py
+++ b/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]
diff --git a/HySoP/hysop/domain/subsets/subset.py b/HySoP/hysop/domain/subsets/subset.py
index 2e8aea2dd353933924c07799868e305aa0545dac..83d74a1ace89e139476ae726ad663c6e81df95f6 100644
--- a/HySoP/hysop/domain/subsets/subset.py
+++ b/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):
"""
diff --git a/HySoP/hysop/numerics/utils.py b/HySoP/hysop/numerics/utils.py
index b3d4bf44b0782cd5fe8a7973782b7c1f09f5976d..c07d0c8def6a152f14183333f80abc62953de6c6 100644
--- a/HySoP/hysop/numerics/utils.py
+++ b/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
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index 29ef4108976c07a10d32aa5e6ba028572fe28959..efe88ecb1fcc4b9e48d559991163cb2c5e0961c0 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/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
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift_noca.py b/HySoP/hysop/operator/discrete/drag_and_lift_noca.py
deleted file mode 100644
index ad272ceec21dce3a1270ee622099e3f4c7aef95c..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/discrete/drag_and_lift_noca.py
+++ /dev/null
@@ -1,295 +0,0 @@
-# -*- 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
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
index 66aa09314b92fbd13e7d868317fe6bbfd059b879..99ab378dca46f3859b8802f3f9655b1cbe6ff62a 100644
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/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
diff --git a/HySoP/hysop/operator/drag_and_lift_noca.py b/HySoP/hysop/operator/drag_and_lift_noca.py
deleted file mode 100644
index 4d72c35c34811d015381b732900c2c7963555704..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/drag_and_lift_noca.py
+++ /dev/null
@@ -1,79 +0,0 @@
-# -*- 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]
diff --git a/HySoP/hysop/operator/tests/test_drag_and_lift.py b/HySoP/hysop/operator/tests/test_drag_and_lift.py
index c598e83d676d25638068a01cff1d0007f75a07ea..150c1019934fd5c90d9292d0517e1be2a4da35ea 100644
--- a/HySoP/hysop/operator/tests/test_drag_and_lift.py
+++ b/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)