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

Fix differential op tests

parent 5f219872
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
hysop/numerics/update_ghosts.py
+ 78
− 95
View file @ 5ffe076f
"""
@file numerics/update_ghosts.py
Update ghost points for a list of numpy arrays
"""Update ghost points of a list of numpy arrays
for a given topology.
Setup for send/recv process of ghosts points for a list
of numpy arrays, for a given topology.
.. currentmodule:: hysop.numerics.update_ghosts
* :class:`~UpdateGhosts` : udpate ghosts layers in a 'classical' way
* :class:`~UpdateGhostsFull` : update also points in the corners of the grid,
i.e. at the intersection of ghosts points lines.
Usage :
.. code::
# init
up = UpdateGhosts(topo, 2)
...
# update
up([field1, field2])
"""
from hysop.constants import debug, PERIODIC, HYSOP_MPI_REAL
import hysop.tools.numpywrappers as npw
......@@ -10,25 +27,28 @@ import numpy as np
class UpdateGhosts(object):
"""
Ghost points synchronization for a list of numpy arrays
"""Ghost points synchronization for a list of numpy arrays
"""
@debug
def __init__(self, topo, nbElements):
def __init__(self, topo, nb_elements):
"""
Setup for send/recv process of ghosts points for a list
of numpy arrays, for a given topology.
@param topology : the topology common to all fields.
@param nbElements : max number of arrays that will be update
at each call.
nbElements and memshape will be used to allocate memory for local
buffers used for send-recv process.
Parameters
----------
topo : :class:`hysop.mpi.topology.Cartesian`
data and mpi process distribution
nb_elements : int
number of arrays that will be update
at each call.
Notes
------
* nb_elements and topo.mesh.resolution will be used to
allocate memory for local buffers used for send-recv process.
"""
## The mpi topology and mesh distribution
# The mpi topology and mesh distribution
self.topology = topo
## Ghost layer
# Ghost layer
self.ghosts = self.topology.mesh.discretization.ghosts
# Indices of points to be filled from previous neighbour
# Each component is a slice and corresponds to a direction in
......@@ -49,12 +69,12 @@ class UpdateGhosts(object):
self._setup_slices()
## shape of numpy arrays to be updated.
# shape of numpy arrays to be updated.
self.memshape = tuple(self.topology.mesh.resolution)
# length of memory required to save one numpy array
self._memoryblocksize = np.prod(self.memshape)
## Number of numpy arrays that will be updated
self.nbElements = nbElements
# Number of numpy arrays that will be updated
self.nb_elements = nb_elements
if self.topology.size > 1: # else no need to set buffers ...
# Size computation below assumes that what we send in one
# dir has the same size as what we receive from process in the
......@@ -66,13 +86,14 @@ class UpdateGhosts(object):
for d in exchange_dir:
buffsize = 0
buffsize += temp[self._g_tonext[d]].size * self.nbElements
buffsize += temp[self._g_tonext[d]].size * self.nb_elements
self._sendbuffer.append(npw.zeros((buffsize)))
self._recvbuffer.append(npw.zeros((buffsize)))
def _setup_slices(self):
"""
Compute slices to send and recieve ghosts values.
"""Compute slices used to describe what to send
and receive in ghosts points.
"""
defslice = [slice(None, None, None)] * self._dim
nogh_slice = [slice(0)] * self._dim
......@@ -89,8 +110,8 @@ class UpdateGhosts(object):
self._g_tonext[d][d] = slice(-2 * self.ghosts[d],
-self.ghosts[d])
otherDim = [x for x in xrange(self._dim) if x != d]
for d2 in otherDim:
other_dim = self._other_dim(d)
for d2 in other_dim:
self._g_fromprevious[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
self._g_fromnext[d][d2] = slice(self.ghosts[d2],
......@@ -105,26 +126,31 @@ class UpdateGhosts(object):
self._g_toprevious.append(list(nogh_slice))
self._g_tonext.append(list(nogh_slice))
def _other_dim(self, d):
return [x for x in xrange(self._dim) if x != d]
def __call__(self, variables):
return self.apply(variables)
def applyBC(self, variables):
"""
Apply boundary conditions for non-distributed directions (only).
@param variables : a list of ndarrays
def apply_bc(self, variables):
"""Apply boundary conditions for non-distributed directions (only).
Parameters
----------
variables : a list of ndarrays
Note that in distributed directions, BC are automatically set
during ghosts update (apply).
"""
"""
assert (self.topology.domain.boundaries == PERIODIC).all(),\
'Only implemented for periodic boundary conditions.'
assert isinstance(variables, list)
dirs = [d for d in xrange(self._dim)
if self.topology.shape[d] == 1]
for d in dirs:
self._applyBC_in_dir(variables, d)
self._apply_bc_in_dir(variables, d)
def _applyBC_in_dir(self, variables, d):
def _apply_bc_in_dir(self, variables, d):
"""Apply periodic boundary condition in direction d."""
for v in variables:
assert v.shape == self.memshape
......@@ -148,7 +174,7 @@ class UpdateGhosts(object):
i += 1
# End of loop through send/recv directions.
# Apply boundary conditions for non-distributed directions
self.applyBC(variables)
self.apply_bc(variables)
def _apply_in_dir(self, variables, d, i):
"""Communicate ghosts values in direction d for neighbour
......@@ -174,7 +200,7 @@ class UpdateGhosts(object):
recvbuf=self._recvbuffer[i],
source=from_rk, recvtag=from_rk)
# 3 - Print recvbuffer back to variables and update sendbuffer
# 3 - fill variables with recvbuffer and update sendbuffer
# for next send
pos = 0
nextpos = 0
......@@ -193,7 +219,7 @@ class UpdateGhosts(object):
dest=dest_rk, sendtag=rank,
recvbuf=self._recvbuffer[i],
source=from_rk, recvtag=from_rk)
# 5 - Print recvbuffer back to variables.
# 5 - fill variables with recvbuffer
pos = 0
nextpos = 0
for v in variables:
......@@ -204,68 +230,25 @@ class UpdateGhosts(object):
class UpdateGhostsFull(UpdateGhosts):
"""
Ghost points synchronization for a list of numpy arrays
"""
"""Ghost points synchronization for a list of numpy arrays
@debug
def __init__(self, topo, nbElements):
"""
Setup for send/recv process of ghosts points for a list
of numpy arrays, for a given topology.
@param topology : the topology common to all fields.
@param nbElements : max number of arrays that will be update
at each call.
nbElements and memshape will be used to allocate memory for local
buffers used for send-recv process.
This version differs from UpdateGhosts by computing also ghosts values
in edges and corners of the domain. The directions are computed in reversed order.
"""
super(UpdateGhostsFull, self).__init__(topo, nbElements)
Notes
-----
* This version differs from UpdateGhosts
by computing also ghosts values
in edges and corners of the domain.
The directions are computed in reversed order.
def _setup_slices(self):
"""
Computes slices to send and recieve ghosts values.
It assumes that directions are computed from an xrange() loop so that
ghosts in previous directions are completed.
"""
defslice = [slice(None, None, None)] * self._dim
nogh_slice = [slice(0)] * self._dim
for d in xrange(self._dim):
if self.ghosts[d] > 0:
self._g_fromprevious.append(list(defslice))
self._g_fromprevious[d][d] = slice(self.ghosts[d])
self._g_fromnext.append(list(defslice))
self._g_fromnext[d][d] = slice(-self.ghosts[d], None, None)
self._g_toprevious.append(list(defslice))
self._g_toprevious[d][d] = slice(self.ghosts[d],
2 * self.ghosts[d], None)
self._g_tonext.append(list(defslice))
self._g_tonext[d][d] = slice(-2 * self.ghosts[d],
-self.ghosts[d])
"""
## Slices for other directions corresponding to directions not
## yet exchanged : x > d. For directions x < d, the slices is a
## full slice that includes ghosts. This assumes that directions
## x < d have already been computed (by communications or local
## exchanges)
# if d == 0:
otherDim = [x for x in xrange(self._dim) if x > d]
for d2 in otherDim:
self._g_fromprevious[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
self._g_fromnext[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
self._g_toprevious[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
self._g_tonext[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
else:
self._g_fromprevious.append(list(nogh_slice))
self._g_fromnext.append(list(nogh_slice))
self._g_toprevious.append(list(nogh_slice))
self._g_tonext.append(list(nogh_slice))
def _other_dim(self, d):
# Slices for other directions corresponding to directions not
# yet exchanged : x > d. For directions x < d, the slices is a
# full slice that includes ghosts. This assumes that directions
# x < d have already been computed (by communications or local
# exchanges)
# if d == 0:
return [x for x in xrange(self._dim) if x > d]
@debug
def apply(self, variables):
......@@ -285,7 +268,7 @@ class UpdateGhostsFull(UpdateGhosts):
i = 0
for d in xrange(self._dim):
if d in local_bc_dir:
self._applyBC_in_dir(variables, d)
self._apply_bc_in_dir(variables, d)
elif d in exchange_dir:
self._apply_in_dir(variables, d, i)
# update index in neighbours list
......
hysop/operator/differential.py
+ 21
− 22
View file @ 5ffe076f
......@@ -19,6 +19,8 @@ from hysop.numerics.finite_differences import FiniteDifference
import hysop.default_methods as default
from abc import ABCMeta, abstractmethod
from hysop.numerics.differential_operations import Curl as NumCurl
from hysop.numerics.differential_operations import DivAdvection as NumDA
from hysop import __FFTW_ENABLED__
class Differential(Computational):
......@@ -105,38 +107,23 @@ class Curl(Differential):
"""
def _init_space_discr_method(self):
if self.method[SpaceDiscretisation] is 'fftw':
if self.method[SpaceDiscretisation] is 'fftw' and __FFTW_ENABLED__:
op_class = CurlFFT
elif self.method[SpaceDiscretisation].mro()[1] is FiniteDifference:
op_class = CurlFD
elif not isinstance(self.method[SpaceDiscretisation], str):
if self.method[SpaceDiscretisation].mro()[1] is FiniteDifference:
op_class = CurlFD
else:
raise ValueError("The required Space Discretisation is\
not available for Curl.")
return op_class
def get_work_properties(self):
"""Get properties of internal work arrays. Must be call after discretize
but before setup.
Returns
-------
dictionnary
keys = 'rwork' and 'iwork', values = list of shapes of internal
arrays required by this operator (real arrays for rwork, integer
arrays for iwork).
"""
if not self._is_discretized:
msg = 'The operator must be discretized '
msg += 'before any call to this function.'
raise RuntimeError(msg)
super(Curl, self).get_work_properties()
res = {'rwork': None, 'iwork': None}
# Only FD methods need internal work space
if self.method[SpaceDiscretisation].mro()[1] is FiniteDifference:
work_length = NumCurl.get_work_length()
shape = self.discreteFields[self.invar].data[0].shape
res['rwork'] = []
for _ in xrange(work_length):
res['rwork'].append(shape)
toporef = self.variables[self.invar]
res = NumCurl.get_work_properties(toporef)
return res
......@@ -152,6 +139,9 @@ class Grad(Differential):
not available for Grad.")
return op_class
def get_work_properties(self):
return {'rwork': None, 'iwork': None}
class DivAdvection(Differential):
"""Computes outVar = -nabla .(invar . nabla(nvar))
......@@ -163,3 +153,12 @@ class DivAdvection(Differential):
raise ValueError("The required Space Discretisation is\
not available for DivAdvection.")
return op_class
def get_work_properties(self):
super(DivAdvection, self).get_work_properties()
res = {'rwork': None, 'iwork': None}
# Only FD methods need internal work space
if self.method[SpaceDiscretisation].mro()[1] is FiniteDifference:
toporef = self.variables[self.invar]
res = NumDA.get_work_properties(toporef)
return res
hysop/operator/discrete/differential.py
+ 14
− 26
View file @ 5ffe076f
......@@ -13,7 +13,6 @@ from hysop.constants import debug
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.numerics.differential_operations import Curl, GradV,\
DivAdvection
import hysop.tools.numpywrappers as npw
from abc import ABCMeta, abstractmethod
from hysop.numerics.update_ghosts import UpdateGhosts
from hysop.methods_keys import SpaceDiscretisation
......@@ -23,13 +22,14 @@ except ImportError:
from hysop.fakef2py import fftw2py
import hysop.default_methods as default
from hysop.tools.profiler import profile
from hysop.tools.misc import WorkSpaceTools
class Differential(DiscreteOperator):
"""Abstract base class for discrete differential operators
"""
__metaclass__ = ABCMeta
# @debug
# def __new__(cls, *args, **kw):
# return object.__new__(cls, *args, **kw)
......@@ -118,18 +118,12 @@ class CurlFD(Differential):
self._function = Curl(topo=self.invar.topology, work=self._rwork,
method=self.method[SpaceDiscretisation])
def _set_work_arrays(self, rwork, iwork):
worklength = Curl.get_work_length()
memshape = self.invar.data[0].shape
if rwork is None:
self._rwork = [npw.zeros(memshape) for _ in xrange(worklength)]
else:
assert isinstance(rwork, list), 'rwork must be a list.'
# --- External rwork ---
self._rwork = rwork
assert len(self._rwork) == worklength
for wk in self._rwork:
assert wk.shape == memshape
def _set_work_arrays(self, rwork=None, iwork=None):
work_prop = Curl.get_work_properties(self.invar.topology)
lwork = len(work_prop['rwork'])
memshape = work_prop['rwork'][0]
self._rwork = WorkSpaceTools.check_work_array(lwork, memshape,
rwork)
@debug
@profile
......@@ -175,18 +169,12 @@ class DivAdvectionFD(Differential):
method=self.method[SpaceDiscretisation],
work=self._rwork)
def _set_work_arrays(self, rwork, iwork):
worklength = DivAdvection.get_work_length()
memshape = self.invar.data[0].shape
if rwork is None:
self._rwork = [npw.zeros(memshape) for _ in xrange(worklength)]
else:
assert isinstance(rwork, list), 'rwork must be a list.'
# --- External rwork ---
self._rwork = rwork
assert len(self._rwork) == worklength
for wk in self._rwork:
assert wk.shape == memshape
def _set_work_arrays(self, rwork=None, iwork=None):
work_prop = DivAdvection.get_work_properties(self.invar.topology)
lwork = len(work_prop['rwork'])
memshape = work_prop['rwork'][0]
self._rwork = WorkSpaceTools.check_work_array(lwork, memshape,
rwork)
@debug
@profile
......
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