From f599aed0746f9f8203d9b5e7e4a7b11901233389 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Thu, 18 Sep 2014 16:01:13 +0200
Subject: [PATCH] bye bye monitors (2)
---
HySoP/hysop/numerics/utils.py | 64 +++++
.../hysop/operator/discrete/drag_and_lift.py | 248 ++++++++++++++++++
HySoP/hysop/operator/drag_and_lift.py | 79 ++++++
3 files changed, 391 insertions(+)
create mode 100644 HySoP/hysop/numerics/utils.py
create mode 100644 HySoP/hysop/operator/discrete/drag_and_lift.py
create mode 100644 HySoP/hysop/operator/drag_and_lift.py
diff --git a/HySoP/hysop/numerics/utils.py b/HySoP/hysop/numerics/utils.py
new file mode 100644
index 000000000..b3d4bf44b
--- /dev/null
+++ b/HySoP/hysop/numerics/utils.py
@@ -0,0 +1,64 @@
+from parmepy.constants import XDIR, YDIR, ZDIR
+import parmepy.tools.numpywrappers as npw
+import numpy as np
+
+
+class Utils(object):
+
+ i1 = [YDIR, ZDIR, XDIR]
+ i2 = [ZDIR, XDIR, YDIR]
+ gen_indices = zip(i1, i2)
+
+ @staticmethod
+ def sum_cross_product(x, y, sl, work, subsets=None):
+ """
+ @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
+ res = npw.zeros(3)
+ 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 sum_cross_product_2(x, y, ind, work):
+ current_dir = 0
+ res = npw.zeros(3)
+ ilist = np.where(ind)
+ nb = len(ilist[0])
+ for (i, j) in Utils.gen_indices:
+ work.flat[:nb] = x[i].flat[ilist[i]] * y[j][ind]\
+ - x[j].flat[ilist[j]] * y[i][ind]
+ res[current_dir] = npw.real_sum(work.flat[:nb])
+ current_dir += 1
+ return res
+
+ @staticmethod
+ def sum_cross_product_3(x, y, ind):
+ """
+ Integrate over the control box using python loops.
+ ---> wrong way, seems to be really slower although
+ it costs less in memory.
+ Used only for tests (timing).
+ """
+ ilist = np.where(ind)
+ res = npw.zeros(3)
+ for(ix, iy, iz) in zip(ilist[0], ilist[YDIR], ilist[ZDIR]):
+ res[XDIR] += x[YDIR][0, iy, 0] * y[ZDIR][ix, iy, iz]\
+ - x[ZDIR][0, 0, iz] * y[YDIR][ix, iy, iz]
+ res[YDIR] += x[ZDIR][0, 0, iz] * y[XDIR][ix, iy, iz]\
+ - x[XDIR][ix, 0, 0] * y[ZDIR][ix, iy, iz]
+ res[ZDIR] += x[XDIR][ix, 0, 0] * y[YDIR][ix, iy, iz]\
+ - x[YDIR][0, iy, 0] * y[XDIR][ix, iy, iz]
+ return res
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
new file mode 100644
index 000000000..0d54f4284
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -0,0 +1,248 @@
+# -*- 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 sum_cross_product
+
+
+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, normalization,
+ volumeOfControl, 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)
+
+ ## 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)
+ self._voc.discretize(self._topology)
+ if obstacles is None:
+ self._obsinds = []
+ else:
+ for obs in obstacles:
+ obs.discretize(self._topology)
+ self._obsinds = [ind for ind in obs.ind[self._topology]
+ for obs in 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)
+
+ 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 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
+ 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 _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 = 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
new file mode 100644
index 000000000..c92be7f60
--- /dev/null
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -0,0 +1,79 @@
+# -*- 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, nu, normalization,
+ volumeOfControl, 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]
--
GitLab