Skip to content
Snippets Groups Projects
Commit 8e197231 authored by Franck Pérignon's avatar Franck Pérignon
Browse files

pre-merge v0

parent 426249b5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
HySoP/hysop/operator/monitors/compute_forces.py deleted 100644 → 0
+ 0
− 296
View file @ 426249b5
# -*- 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment