From 8e1972311db90d2b88ab2a500ba3516ccfdb93f5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Mon, 17 Nov 2014 17:19:40 +0100
Subject: [PATCH] pre-merge v0
---
.../hysop/operator/monitors/compute_forces.py | 296 ------------------
1 file changed, 296 deletions(-)
delete mode 100644 HySoP/hysop/operator/monitors/compute_forces.py
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
deleted file mode 100644
index 701f711a8..000000000
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ /dev/null
@@ -1,296 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@file compute_forces.py
-Compute forces
-"""
-from parmepy.constants import np, XDIR, YDIR, ZDIR
-from parmepy.numerics.updateGhosts import UpdateGhosts
-from parmepy.numerics.differential_operations import Laplacian
-from parmepy.numerics.finite_differences import FD_C_2
-from parmepy.operator.monitors.monitoring import Monitoring
-import parmepy.tools.numpywrappers as npw
-
-
-class DragAndLift(Monitoring):
- """
- 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, coefForce,
- 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
- """
- if 'io_params' in kwds:
- params = kwds['io_params']
- if not "filename" in params:
- params["filename"] = "drag_and_lift"
- # Set output buffer shape
- params["writebuffshape"] = (1, 4)
- super(DragAndLift, self).__init__(variables=[velocity, vorticity],
- **kwds)
-
- self.velocity = velocity
- self.vorticity = vorticity
- self._voc = volumeOfControl
- self._dim = self.domain.dimension
- msg = 'Force computation undefined for domain of dimension 1.'
- assert self._dim > 1, msg
- if obstacles is None:
- obstacles = []
- ## A list of obstacles (rigid bodies) in the control box
- self.obstacles = obstacles
- # 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
- # discrete velocity field
- self.vd = None
- # discrete vorticity field
- self.wd = None
- ## viscosity
- self.nu = nu
- # Normalizing coefficient for forces
- # (based on the physics of the flow)
- self.coefForce = coefForce
-
- 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)
-
- def setUp(self):
- """
- """
- if not self._isUpToDate:
- self._step = np.asarray(self.topology.mesh.space_step)
- self._dvol = np.prod(self._step)
- self._work = npw.zeros(self.topology.mesh.resolution)
- 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)
-
- for v in self.variables:
- # the discrete fields
- self.discreteFields[v] = v.discretize(self.topology)
- self.vd = self.discreteFields[self.velocity]
- self.wd = self.discreteFields[self.vorticity]
- self._voc.discretize(self.topology)
- for obs in self.obstacles:
- obs.discretize(self.topology)
- # prepare ghost points synchro for velocity and vorticity used
- # in fd schemes
- self._synchronize = UpdateGhosts(self.topology,
- self.vd.nbComponents
- + self.wd.nbComponents)
- self._isUpToDate = True
-
- 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.vd.data + self.wd.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 --
- ## for s in self._voc.upperS:
- ## self._buffer = self._integrateOnSurface(s, self._buffer)
- ## self.force += self._buffer
- ## for s in self._voc.lowerS:
- ## self._buffer = self._integrateOnSurface(s, self._buffer)
- ## self.force += self._buffer
- self._buffer = self._integrateOnSurface(self._voc.upperS[0], self._buffer)
- self.force += self._buffer
- self._buffer = self._integrateOnSurface(self._voc.lowerS[0], self._buffer)
- self.force += self._buffer
-
- # Reduce results over all MPI process in topo
- self.mpi_sum()
-
- self.force *= self.coefForce
-
- # Print results, if required
- if self._writer is not None and self._writer.doWrite(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.vd.data
- wdata = self.wd.data
- # i1 : normal dir
- # i2 : other dirs
- i1 = np.where(normal)[0][0]
- i2 = np.where(normal == 0)[0]
- dsurf = np.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 * np.sum((-vdata[i1][sl] ** 2
- + sum([vdata[j][sl] ** 2
- for j in i2])))
- # other components
- for j in i2:
- res[j] = -normal[i1] * np.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] * np.sum(coords[j] * buff)
- res[j] -= self._coeff * normal[i1] * coords[i1] * np.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._work[...] = self._laplacian(vdata[j], self._work)
- res[i1] += self._coeff * self.nu * normal[i1]\
- * np.sum(coords[j] * self._work[sl])
- res[j] -= self._coeff * self.nu * normal[i1] * coords[i1] * \
- np.sum(self._work[sl])
- self._fd_scheme.computeIndices(sl)
- self._fd_scheme.compute(vdata[i1], i1, self._work)
- res[i1] += 2.0 * normal[i1] * self.nu * np.sum(self._work[sl])
- for j in i2:
- self._fd_scheme.compute(vdata[i1], j, self._work)
- res[j] += normal[i1] * self.nu * np.sum(self._work[sl])
- self._fd_scheme.compute(vdata[j], i1, self._work)
- res[j] += normal[i1] * self.nu * np.sum(self._work[sl])
-
- res *= dsurf
- return res
-
- def _integrateOnBox(self, res):
- assert self._dim == 3, 'Not defined for dim < 3'
- coords = self._voc.coords[self.topology]
- wdata = self.wd.data
- i1 = [YDIR, ZDIR, XDIR]
- i2 = [ZDIR, XDIR, YDIR]
- direc = 0
- sl = self._voc.slices[self.topology]
- for (i, j) in zip(i1, i2):
- self._work[sl] = coords[i] * wdata[j][sl]
- self._work[sl] -= coords[j] * wdata[i][sl]
- for obs in self.obstacles:
- for inds in obs.ind[self.topology]:
- self._work[inds] = 0.0
- res[direc] = np.sum(self._work[sl])
- direc += 1
- res *= self._dvol
- return res
-
- def _integrateOnBox2(self, res):
- assert self._dim == 3, 'Not defined for dim < 3'
- coords = self.topology.mesh.coords
- wdata = self.wd.data
- i1 = [YDIR, ZDIR, XDIR]
- i2 = [ZDIR, XDIR, YDIR]
- direc = 0
- ind = self._voc.ind[self.topology][0]
- ilist = np.where(ind)
- nb = len(ilist[0])
- ind = self._voc.ind[self.topology][0]
- for (i, j) in zip(i1, i2):
- self._work.flat[:nb] = coords[i].flat[ilist[i]] * wdata[j][ind]\
- - coords[j].flat[ilist[j]] * wdata[i][ind]
- res[direc] = np.sum(self._work.flat[:nb])
- direc += 1
- res *= self._dvol
- return res
-
- def _integrateOnBoxLoop(self, res):
- """
- 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).
- """
- assert self._dim == 3, 'Not defined for dim < 3'
- coords = self.topology.mesh.coords
- ind = self._voc.ind[self.topology][0]
- ilist = np.where(ind)
- wdata = self.wd.data
- for(ix, iy, iz) in zip(ilist[0], ilist[YDIR], ilist[ZDIR]):
- res[XDIR] += coords[YDIR][0, iy, 0] * wdata[ZDIR][ix, iy, iz]\
- - coords[ZDIR][0, 0, iz] * wdata[YDIR][ix, iy, iz]
- res[YDIR] += coords[ZDIR][0, 0, iz] * wdata[XDIR][ix, iy, iz]\
- - coords[XDIR][ix, 0, 0] * wdata[ZDIR][ix, iy, iz]
- res[ZDIR] += coords[XDIR][ix, 0, 0] * wdata[YDIR][ix, iy, iz]\
- - coords[YDIR][0, iy, 0] * wdata[XDIR][ix, iy, iz]
-
- res *= self._dvol
- return res
--
GitLab