From 5f219872f53e5882f39b5d16e928a6aa9bbb2590 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Mon, 23 May 2016 14:47:10 +0200
Subject: [PATCH] Update differential operations: internal work array
management, tests and docs
---
.../users_guide/differential_operations.rst | 8 +
hysop/numerics/differential_operations.py | 751 +++++++++---
hysop/numerics/tests/test_diffOp.py | 222 ----
.../tests/test_differential_operations.py | 1015 +++++++++++++++++
4 files changed, 1641 insertions(+), 355 deletions(-)
create mode 100644 docs/sphinx/users_guide/differential_operations.rst
delete mode 100755 hysop/numerics/tests/test_diffOp.py
create mode 100755 hysop/numerics/tests/test_differential_operations.py
diff --git a/docs/sphinx/users_guide/differential_operations.rst b/docs/sphinx/users_guide/differential_operations.rst
new file mode 100644
index 000000000..94552b4ce
--- /dev/null
+++ b/docs/sphinx/users_guide/differential_operations.rst
@@ -0,0 +1,8 @@
+.. _differential_operations:
+
+.. currentmodule:: hysop.numerics.differential_operations
+
+Vector calculous based on finite differences
+--------------------------------------------
+
+To be done ...
diff --git a/hysop/numerics/differential_operations.py b/hysop/numerics/differential_operations.py
index 67559774c..c6ba55445 100755
--- a/hysop/numerics/differential_operations.py
+++ b/hysop/numerics/differential_operations.py
@@ -5,11 +5,18 @@
* :class:`~hysop.numerics.differential_operations.Curl`,
* :class:`~hysop.numerics.differential_operations.DivRhoV`,
* :class:`~hysop.numerics.differential_operations.DivWV`,
+* :class:`~hysop.numerics.differential_operations.Laplacian`,
* :class:`~hysop.numerics.differential_operations.GradS`,
* :class:`~hysop.numerics.differential_operations.GradV`,
* :class:`~hysop.numerics.differential_operations.GradVxW`,
* :class:`~hysop.numerics.differential_operations.DivAdvection`,
-* :class:`~hysop.numerics.differential_operations.Laplacian`,
+* :class:`~hysop.numerics.differential_operations.Strain`,
+* :class:`~hysop.numerics.differential_operations.StrainCriteria`,
+* :class:`~hysop.numerics.differential_operations.MaxDiagGradV`,
+* :class:`~hysop.numerics.differential_operations.StretchLike`,
+* :class:`~hysop.numerics.differential_operations.DiagAndStretch`,
+* :class:`~hysop.numerics.differential_operations.StrainCriteria`,
+* :class:`~hysop.numerics.differential_operations.StrainAndStretch`,
* :class:`~hysop.numerics.differential_operations.DifferentialOperation`,
(abstract base class).
@@ -25,6 +32,8 @@ For example :
* a = DivRhoV(c,d), d a vector field, a and c scalar fields. d is a list of
3 arrays and a and c lists of 1 array.
+For practical examples of use, see test_differential_operations.py file.
+
"""
from hysop.constants import debug, XDIR, YDIR, ZDIR
from abc import ABCMeta
@@ -46,7 +55,7 @@ class DifferentialOperation(object):
# return object.__new__(cls, *args, **kw)
@debug
- def __init__(self, topo, indices=None, reduce_output_shape=None,
+ def __init__(self, topo, indices=None, reduce_output_shape=False,
method=None, work=None):
"""
Parameters
@@ -81,8 +90,6 @@ class DifferentialOperation(object):
* If work = None, some work arrays will be allocated internally.
Else, you must provide a list of lwk arrays of shape
topo.mesh.resolution,
- with lwk = ClassName.get_work_length(), ClassName being Curl, DivV ...
- or any other DifferentialOperation.
"""
self._dim = topo.domain.dimension
@@ -93,49 +100,68 @@ class DifferentialOperation(object):
self.method = method
if indices is None:
indices = topo.mesh.compute_index
- reduce_output_shape = False
- self._indices = indices
+ self.in_indices = indices
+ # True if work has to be done for this set
+ # of indices on the current proc.
+ self._on_proc = Utils.is_on_proc(indices)
self.fd_scheme = self._init_fd_method(topo, reduce_output_shape)
self.output_indices = self.fd_scheme.output_indices
- self._work = self._set_work_arrays(work, topo, reduce_output_shape)
+ self._work = self._set_work_arrays(topo, reduce_output_shape, work)
- def _set_work_arrays(self, work, topo, reduce_output_shape):
- """Check and allocate internal work buffers.
+ def _set_work_arrays(self, topo, reduce_output_shape, rwork=None):
+ """Check and/or allocate internal work buffers.
"""
- lwk = self.get_work_length()
+ wk_prop = self.get_work_properties(topo, self.in_indices)['rwork']
+ if wk_prop is None:
+ return []
if reduce_output_shape:
- shape = self.fd_scheme.result_shape
+ subshape = np.asarray([self.in_indices[i].stop -
+ self.in_indices[i].start
+ for i in xrange(self._dim)])
else:
- shape = tuple(topo.mesh.resolution)
- if work is None:
- work = []
- for _ in xrange(lwk):
- work.append(npw.asrealarray(shape))
- else:
- msg = 'Wrong input for work arrays.'
- assert isinstance(work, list), msg
- assert len(work) == lwk, msg
- for wk in work:
- assert wk.shape == shape
- return work
+ subshape = topo.mesh.resolution
+ subshape = tuple(subshape)
+ return WorkSpaceTools.check_work_array(len(wk_prop), subshape, rwork)
@staticmethod
- def get_work_length():
+ def get_work_properties(topo, indices=None):
+ """Default : no work vector
"""
- Compute the number of required work arrays for this method.
+ return DifferentialOperation._find_work(0, topo, indices)
+
+ @staticmethod
+ def _find_work(lwork, topo, indices=None):
+ """Internal function to find the shape of work vector
+ for a given topo and a set of indices
"""
- return 0
+ if lwork == 0:
+ return {'rwork': None, 'iwork': None}
+
+ if indices is not None:
+ if not Utils.is_on_proc(indices):
+ shape = (0, ) * topo.domain.dimension
+ else:
+ shape = np.prod(tuple([indices[i].stop - indices[i].start
+ for i in xrange(len(indices))]))
+ else:
+ shape = np.prod(topo.mesh.resolution)
+ return {'rwork': [(shape, ), ] * lwork, 'iwork': None}
def _init_fd_method(self, topo, reduce_output_shape):
"""Build the finite difference scheme
"""
msg = 'FD scheme Not yet implemented for this operation.'
assert self.method.__mro__[0] in self._authorized_methods, msg
+ if not self._on_proc:
+ empty_set = [slice(0, 0), ] * self._dim
+ return self.method(topo.mesh.space_step,
+ empty_set)
+
fd_scheme = self.method(topo.mesh.space_step,
- self._indices,
+ self.in_indices,
reduce_output_shape)
msg = 'Ghost layer is too small for the chosen FD scheme.'
- required_ghost_layer = fd_scheme.minimal_ghost_layer
+ required_ghost_layer = fd_scheme.ghosts_layer_size
assert (topo.ghosts() >= required_ghost_layer).all(), msg
return fd_scheme
@@ -151,19 +177,22 @@ class Curl(DifferentialOperation):
"""Curl of a vector field
"""
super(Curl, self).__init__(**kwds)
- assert len(self._indices) > 1
+ assert len(self.in_indices) > 1
# connect to fd function call
self.fcall = self._central
- if len(self._indices) == 2:
+ if len(self.in_indices) == 2:
# a 'fake' curl for 2D case
self.fcall = self._central_2d
@staticmethod
- def get_work_length():
- return 1
+ def get_work_properties(topo, indices=None):
+ return Curl._find_work(Curl._lwork, topo, indices)
def __call__(self, variable, result):
- return self.fcall(variable, result)
+ if self._on_proc:
+ return self.fcall(variable, result)
+ else:
+ return result
def _central(self, variable, result):
""" 3D Curl
@@ -239,13 +268,83 @@ class DivRhoV(DifferentialOperation):
Default method : FDC4
"""
+ self._in_gh = None
+ self._wk_indices = None
+ self._fd_work = None
super(DivRhoV, self).__init__(**kwds)
- #self.fcall = self._central4_caa
- self.fcall = self._central4
+
+ def _init_fd_method(self, topo, reduce_output_shape):
+ """Build the finite difference scheme
+ """
+ msg = 'FD scheme Not yet implemented for this operation.'
+ assert self.method.__mro__[0] in self._authorized_methods, msg
+ # for this operation, the fd scheme uses
+ # work array as input, which means that work
+ # array must include ghost points. This is done
+ # thanks to self.wk_indices
+ if not self._on_proc:
+ empty_set = [slice(0, 0), ] * self._dim
+ self._fd_work = self.method(topo.mesh.space_step,
+ empty_set)
+ self._in_gh = empty_set
+ self._wk_indices = empty_set
+ return self._fd_work
+
+ gh = self.method.ghosts_layer_size
+ ref = self.in_indices
+ self._in_gh = [slice(ref[i].start - gh, ref[i].stop + gh)
+ for i in xrange(self._dim)]
+ for sl in self._in_gh:
+ assert sl.start >= 0
+ assert sl.stop >= 0
+ wk_shape = tuple([self._in_gh[i].stop - self._in_gh[i].start
+ for i in xrange(self._dim)])
+ self._wk_indices = [slice(gh, wk_shape[i] - gh)
+ for i in xrange(self._dim)]
+ if reduce_output_shape:
+ iout = [slice(0, wk_shape[i] - 2 * gh)
+ for i in xrange(self._dim)]
+ else:
+ iout = None
+ fd_scheme = self.method(topo.mesh.space_step,
+ self._wk_indices,
+ indices_out=iout)
+ self._fd_work = self.method(topo.mesh.space_step,
+ self._wk_indices)
+ msg = 'Ghost layer is too small for the chosen FD scheme.'
+ required_ghost_layer = fd_scheme.ghosts_layer_size
+ assert (topo.ghosts() >= required_ghost_layer).all(), msg
+ return fd_scheme
+
+ def _set_work_arrays(self, topo, reduce_output_shape, rwork=None):
+ """Check and allocate internal work buffers.
+ """
+ wk_prop = self.get_work_properties(topo, self.in_indices)['rwork']
+ if wk_prop is None:
+ return []
+ subshape = np.asarray([self.in_indices[i].stop -
+ self.in_indices[i].start
+ for i in xrange(self._dim)])
+ if self._on_proc:
+ subshape += 2 * self.fd_scheme.ghosts_layer_size
+ subshape = tuple(subshape)
+
+ return WorkSpaceTools.check_work_array(len(wk_prop), subshape, rwork)
@staticmethod
- def get_work_length():
- return 2
+ def get_work_properties(topo, indices=None):
+ # fd scheme will be applied on work so it needs ghost points.
+ if indices is not None:
+ if not Utils.is_on_proc(indices):
+ shape = (0, ) * topo.domain.dimension
+ else:
+ shape = np.asarray([indices[i].stop - indices[i].start
+ for i in xrange(len(indices))])
+ shape += 2 * topo.ghosts()
+ else:
+ shape = np.asarray(topo.mesh.resolution).copy()
+ shape = np.prod(shape)
+ return {'rwork': [(shape,), ] * 2, 'iwork': None}
def __call__(self, var1, scal, result):
"""Apply operation
@@ -267,7 +366,10 @@ class DivRhoV(DifferentialOperation):
assert len(result) == len(scal)
for i in xrange(len(var1)):
assert var1[i] is not result
- return self.fcall(var1, scal, result)
+ if self._on_proc:
+ return self._central4(var1, scal, result)
+ else:
+ return result
def _central4(self, var1, scal, result):
"""
@@ -282,45 +384,23 @@ class DivRhoV(DifferentialOperation):
# Result does not need initialisation to zero.
# d/dx (scal * var1x), saved into result
- np.multiply(scal[0], var1[XDIR], self._work[0])
- self.fd_scheme.compute(self._work[0], XDIR, result[0])
- # other components (if any) saved into _work[0] and added into result
- # d/dy (scal * var1y), saved in work and added into result
- # d/dz(scal * var1z), saved in work and added into result (if 3D)
- for cdir in xrange(1, len(var1)):
- np.multiply(scal[0], var1[cdir], self._work[0])
- self.fd_scheme.compute(self._work[0], cdir, self._work[1])
- result[0][self.output_indices] += \
- self._work[1][self.output_indices]
-
- return result
-
- def _central4_caa(self, var1, scal, result):
- """
- Compute central finite difference scheme, order 4.
- Use fd_scheme.compute_and_add.
- """
- # _work[0] is used as temporary space
- # for computation
- # div computations are accumulate into result.
- # result does not need initialisation to zero.
-
- # d/dx (scal * var1x), saved into result
- np.multiply(scal[0][self._indices], var1[XDIR][self._indices],
- self._work[0][self.output_indices])
- self.fd_scheme.compute(self._work[0], XDIR, result[0])
+ np.multiply(scal[0][self._in_gh],
+ var1[XDIR][self._in_gh], self._work[0])
+ result[0] = self.fd_scheme.compute(self._work[0], XDIR, result[0])
# other components (if any) saved into _work[0] and added into result
# d/dy (scal * var1y), saved in work and added into result
# d/dz(scal * var1z), saved in work and added into result (if 3D)
for cdir in xrange(1, len(var1)):
- np.multiply(scal[0][self._indices], var1[cdir][self._indices],
- self._work[0][self.output_indices])
- self.fd_scheme.compute_and_add(self._work[0], cdir, result[0])
+ np.multiply(scal[0][self._in_gh],
+ var1[cdir][self._in_gh], self._work[0])
+ self._work[1] = self._fd_work.compute(self._work[0], cdir,
+ self._work[1])
+ result[0][self.output_indices] += self._work[1][self._wk_indices]
return result
-class DivWV(DifferentialOperation):
+class DivWV(DivRhoV):
"""
Computes nabla.(W.Vx, W.Vy, W.Vz), W and V some vector fields.
@@ -333,21 +413,13 @@ class DivWV(DifferentialOperation):
Parameters
----------
- fd_optim : bool, optional
- if 'CAA', compute result += div(rhoV) else
- result = div(rhoV), which is the default.
**kwds : parameters for base class
Default method : FDC4
"""
super(DivWV, self).__init__(**kwds)
- self.fcall = self._central4
- # set div(scal.var1) operator
- self.div = DivRhoV(**kwds)
-
- @staticmethod
- def get_work_length():
- return DivRhoV.get_work_length()
+ msg = 'Implemented only in 3D.'
+ assert self._dim == 3, msg
def __call__(self, var1, var2, result):
"""Apply operation
@@ -365,16 +437,11 @@ class DivWV(DifferentialOperation):
list of numpy arrays
"""
- assert len(result) == len(var1)
- return self.fcall(var1, var2, result)
-
- def _central4(self, var1, var2, result):
- # Note FP var1[dir] and var2[dir] must be different from result.
- # This is checked inside divV call.
-
- for cdir in xrange(len(var1)):
- result[cdir:cdir + 1] = self.div(var1, var2[cdir:cdir + 1],
- result=result[cdir:cdir + 1])
+ assert len(result) == len(var2)
+ if self._on_proc:
+ for cdir in xrange(len(var2)):
+ result[cdir:cdir + 1] = self._central4(
+ var1, var2[cdir:cdir + 1], result=result[cdir:cdir + 1])
return result
@@ -386,8 +453,8 @@ class Laplacian(DifferentialOperation):
_lwork = 1
@staticmethod
- def get_work_length():
- return 0
+ def get_work_properties(topo, indices=None):
+ return Laplacian._find_work(Laplacian._lwork, topo, indices)
def __call__(self, var, result):
""""Apply laplacian
@@ -403,10 +470,12 @@ class Laplacian(DifferentialOperation):
numpy array
"""
assert len(var) == len(result)
- for d in xrange(len(var)):
- self.fd_scheme.compute(var[d], 0, result[d])
- for cdir in xrange(1, self._dim):
- self.fd_scheme.compute_and_add(var[d], cdir, result[d])
+ if self._on_proc:
+ for d in xrange(len(var)):
+ self.fd_scheme.compute(var[d], 0, result[d])
+ for cdir in xrange(1, self._dim):
+ self.fd_scheme.compute_and_add(var[d], cdir, result[d],
+ self._work[0])
return result
@@ -426,9 +495,10 @@ class GradS(DifferentialOperation):
in/out result
"""
assert len(result) == self._dim
- for cdir in xrange(self._dim):
- # d/dcdir (scal), saved in data[cdir]
- self.fd_scheme.compute(scal[0], cdir, result[cdir])
+ if self._on_proc:
+ for cdir in xrange(self._dim):
+ # d/dcdir (scal), saved in data[cdir]
+ self.fd_scheme.compute(scal[0], cdir, result[cdir])
return result
@@ -438,10 +508,6 @@ class GradV(DifferentialOperation):
"""
_authorized_methods = [FDC4, FDC2]
- def __init__(self, **kwds):
- super(GradV, self).__init__(**kwds)
- self.grad = GradS(**kwds)
-
def __call__(self, var1, result):
"""Apply gradient
@@ -454,10 +520,13 @@ class GradV(DifferentialOperation):
"""
nbc = len(var1)
assert len(result) == nbc * self._dim
- for cdir in xrange(self._dim):
- i1 = cdir * nbc
- i2 = i1 + nbc
- result[i1:i2] = self.grad(var1[cdir:cdir + 1], result[i1:i2])
+ if self._on_proc:
+ for cdir in xrange(self._dim):
+ pos = nbc * cdir
+ for ddir in xrange(self._dim):
+ # d/dddir (v[cdir]), saved in result[pos]
+ self.fd_scheme.compute(var1[cdir], ddir, result[pos])
+ pos += 1
return result
@@ -470,8 +539,8 @@ class GradVxW(DifferentialOperation):
_lwork = 2
@staticmethod
- def get_work_length():
- return 2
+ def get_work_properties(topo, indices=None):
+ return DifferentialOperation._find_work(GradVxW._lwork, topo, indices)
def __call__(self, var1, var2, result, diagnostics):
"""Apply gradient, with central finite difference scheme.
@@ -487,29 +556,33 @@ class GradVxW(DifferentialOperation):
In/out param, overwritten
"""
assert len(result) == len(var1)
- nbc = len(var1)
- diagnostics[:] = 0.0
- for comp in xrange(nbc):
- result[comp][...] = 0.0
- self._work[1][...] = 0.0
- for cdir in xrange(self._dim):
- # self._work = d/dcdir (var1_comp)
- self.fd_scheme.compute(var1[comp], cdir, self._work[0])
- # some diagnostics ...
- if cdir == comp:
- tmp = np.max(abs(self._work[0]))
- diagnostics[0] = max(diagnostics[0], tmp)
- np.add(abs(self._work[0]), self._work[1], self._work[1])
-
- # compute self._work = self._work.var2[cdir]
- np.multiply(self._work[0][self.output_indices],
- var2[cdir][self._indices],
- self._work[0][self.output_indices])
- # sum to obtain nabla(var_comp) . var2, saved into result[comp]
- npw.add(result[comp][self.output_indices],
- self._work[0][self.output_indices],
- result[comp][self.output_indices])
- diagnostics[1] = max(diagnostics[1], np.max(self._work[1]))
+ if self._on_proc:
+ diagnostics[:] = 0.
+ nbc = len(var1)
+ for comp in xrange(nbc):
+ result[comp][...] = 0.0
+ self._work[1][...] = 0.0
+ for cdir in xrange(self._dim):
+ # self._work = d/dcdir (var1_comp)
+ self.fd_scheme.compute(var1[comp], cdir, self._work[0])
+ # some diagnostics ...
+ if cdir == comp:
+ tmp = np.max(abs(self._work[0][self.output_indices]))
+ diagnostics[0] = max(diagnostics[0], tmp)
+ np.add(abs(self._work[0]), self._work[1], self._work[1])
+
+ # compute self._work = self._work.var2[cdir]
+ np.multiply(self._work[0][self.output_indices],
+ var2[cdir][self.in_indices],
+ self._work[0][self.output_indices])
+ # sum to obtain nabla(var_comp) . var2,
+ # saved into result[comp]
+ npw.add(result[comp][self.output_indices],
+ self._work[0][self.output_indices],
+ result[comp][self.output_indices])
+ diagnostics[1] = max(
+ diagnostics[1], np.max(self._work[1][self.output_indices]))
+
return result, diagnostics
@@ -520,8 +593,9 @@ class DivAdvection(DifferentialOperation):
_lwork = 3
@staticmethod
- def get_work_length():
- return 3
+ def get_work_properties(topo, indices=None):
+ return DifferentialOperation._find_work(DivAdvection._lwork,
+ topo, indices)
def __call__(self, var1, result):
"""Apply divergence, with central finite difference scheme.
@@ -535,6 +609,8 @@ class DivAdvection(DifferentialOperation):
"""
assert len(result) == 1
nbc = len(var1)
+ if not self._on_proc:
+ return result
assert nbc == 3, 'Only 3D case is implemented.'
# Compute diff(var[dir], dir), saved in work[dir]
for d in xrange(nbc):
@@ -591,3 +667,412 @@ class DivAdvection(DifferentialOperation):
result[0][self.output_indices] *= 2.0
return result
+
+
+class Strain(DifferentialOperation):
+ """Compute 0.5(grad U + grad U^T), U a vector field
+ """
+ _authorized_methods = [FDC4, FDC2]
+
+ def __init__(self, **kwds):
+ super(Strain, self).__init__(**kwds)
+ if self._dim == 3:
+ self._i1 = [0, 0, 1]
+ self._i2 = [1, 2, 2]
+ elif self._dim == 2:
+ self._i1 = [0, ]
+ self._i2 = [1, ]
+
+ def __call__(self, var1, result):
+ """Compute strain
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : list of numpy arrays
+ in/out result
+
+ result = [diagonal components, extra-diag comp]
+
+ For example in 3D :
+
+ result = [Strain00, Strain11, Strain22, Strain01, Strain02, Strain12]
+ """
+ nbc = len(var1)
+ assert len(result) == self._dim * (self._dim + 1) / 2
+ if not self._on_proc:
+ return result
+
+ pos = nbc
+ # Extra diagonal terms. result[0] is used as temporary
+ # work array.
+ for i1, i2 in zip(self._i1, self._i2):
+ self.fd_scheme.compute(var1[i1], i2, result[pos])
+ self.fd_scheme.compute_and_add(var1[i2], i1, result[pos],
+ result[0])
+ np.multiply(result[pos], 0.5, result[pos])
+ pos += 1
+ # Then diagonal terms
+ for cdir in xrange(self._dim):
+ self.fd_scheme.compute(var1[cdir], cdir, result[cdir])
+
+ return result
+
+
+class StrainCriteria(DifferentialOperation):
+ """For Strain = 0.5(grad U + grad U^T), U a vector field
+ compute max(sum_j abs(Strain_i,j))
+
+ Warning : this is a local (mpi) max computation, since input array var1
+ is also a local array.
+ """
+ _authorized_methods = [FDC4, FDC2]
+
+ def __init__(self, **kwds):
+ super(StrainCriteria, self).__init__(**kwds)
+ if self._dim == 3:
+ self._i1 = [0, 0, 1]
+ self._i2 = [1, 2, 2]
+ self._pos = [0, 1, 2]
+ self._adds = [[0, 1], [0, 2], [1, 2]]
+ elif self._dim == 2:
+ self._i1 = [0, ]
+ self._i2 = [1, ]
+ self._pos = [0]
+ self._adds = [[0, ], [0, ]]
+
+ @staticmethod
+ def get_work_properties(topo, indices=None):
+ dime = topo.domain.dimension
+ if dime == 2:
+ lwork = 2
+ elif dime == 3:
+ lwork = 4
+
+ return DifferentialOperation._find_work(lwork, topo, indices)
+
+ def __call__(self, var1, result):
+ """Compute criteria
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : numpy array
+ criteria values for each direction.
+
+ """
+ work = self._work
+ assert len(result) == self._dim
+ if not self._on_proc:
+ return result
+
+ # First, compute extra-diagonal terms
+ for i1, i2, ipos in zip(self._i1, self._i2, self._pos):
+ self.fd_scheme.compute(var1[i1], i2, work[ipos])
+ self.fd_scheme.compute_and_add(var1[i2], i1, work[ipos],
+ work[-1])
+ np.multiply(work[ipos], 0.5, work[ipos])
+ np.abs(work[ipos], work[ipos])
+
+ # Then add those terms to diagonal term and compute max
+ for idiag in xrange(self._dim):
+ self.fd_scheme.compute(var1[idiag], idiag, work[-1])
+ np.abs(work[-1], work[-1])
+ for i1 in self._adds[idiag]:
+ np.add(work[i1], work[-1], work[-1])
+ if self._on_proc:
+ result[idiag] = np.max(work[-1][self.output_indices])
+ else:
+ result[idiag] = 0.0
+ return result
+
+
+class MaxDiagGradV(DifferentialOperation):
+ """Maximum of each diagonal term (abs) of the Gradient of a vector field
+ """
+ _authorized_methods = [FDC4, FDC2]
+ _lwork = 1
+
+ @staticmethod
+ def get_work_properties(topo, indices=None):
+ return DifferentialOperation._find_work(MaxDiagGradV._lwork,
+ topo, indices)
+
+ def __call__(self, var1, result):
+ """Compute max
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : numpy array
+ in/out result
+ """
+ assert len(result) == self._dim
+ if not self._on_proc:
+ return result
+
+ for cdir in xrange(self._dim):
+ # d/dcdir (var1[cdir]) saved in work
+ self._work[0] = self.fd_scheme.compute(var1[cdir],
+ cdir, self._work[0])
+ np.abs(self._work[0], self._work[0])
+ if self._on_proc:
+ result[cdir] = np.max(self._work[0][self.output_indices])
+
+ return result
+
+
+class StretchLike(DifferentialOperation):
+ """Maximum of the sum of each 'line' (abs) of the Gradient of a vector field
+
+ Warning: local (mpi) maximum.
+ """
+ _authorized_methods = [FDC4, FDC2]
+ _lwork = 2
+
+ @staticmethod
+ def get_work_properties(topo, indices=None):
+ return DifferentialOperation._find_work(StretchLike._lwork,
+ topo, indices)
+
+ def __call__(self, var1, result):
+ """Compute max
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : numpy array
+ in/out result
+
+ result[i] = max(sum_j abs(diff(u_i, j)) for i = 0..dim
+ """
+ assert len(result) == self._dim
+ result[...] = 0.
+ if not self._on_proc:
+ return result
+ for cdir in xrange(self._dim):
+ self._work[1][...] = 0.0
+ for ddir in xrange(self._dim):
+ # d/dddir (var1[cdir]) saved in work
+ self._work[0] = self.fd_scheme.compute(var1[cdir],
+ ddir, self._work[0])
+ np.abs(self._work[0], self._work[0])
+ np.add(self._work[0], self._work[1], self._work[1])
+ if self._on_proc:
+ result[cdir] = np.max(self._work[1][self.output_indices])
+
+ return result
+
+
+class DiagAndStretch(DifferentialOperation):
+ """Maximum of the sum of each 'line' of abs of the Gradient of a vector field
+ and max of each abs of 'diagonal' block of the gradient.
+
+ Warning : this is a local (mpi) max computation, since input array var1
+ is also a local array.
+ """
+ _authorized_methods = [FDC4, FDC2]
+ _lwork = 2
+
+ @staticmethod
+ def get_work_properties(topo, indices=None):
+ return DifferentialOperation._find_work(DiagAndStretch._lwork,
+ topo, indices)
+
+ def __call__(self, var1, result):
+ """Compute max
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : numpy array
+ in/out result
+
+ result[i] = max(abs(diff(u_i,i))) for i = 0..dim
+ result[dim:...] = max(sum_j abs(diff(u_i, j)) for i = 0..dim
+ """
+ assert len(result) == 2 * self._dim
+ result[...] = 0.
+
+ if not self._on_proc:
+ return result
+
+ pos = 0
+ for cdir in xrange(self._dim):
+ self._work[1][...] = 0.0
+ for ddir in xrange(self._dim):
+ # d/dcdir (var1[cdir]) saved in work
+ self._work[0] = self.fd_scheme.compute(var1[cdir],
+ ddir, self._work[0])
+ np.abs(self._work[0], self._work[0])
+ if cdir == ddir and self._on_proc:
+ result[pos] = np.max(self._work[0][self.output_indices])
+ pos += 1
+ np.add(self._work[0], self._work[1], self._work[1])
+ if self._on_proc:
+ result[self._dim + cdir] = np.max(
+ self._work[1][self.output_indices])
+
+ return result
+
+
+class StrainAndStretch(DifferentialOperation):
+ """For Strain = 0.5(grad U + grad U^T), U a vector field
+ compute max(sum_j abs(Strain_i,j)) and maximum of the sum
+ of each 'line' (abs) of the Gradient of U.
+
+ Warning : this is a local (mpi) max computation, since input array var1
+ is also a local array.
+
+ result = [ max(sum_j(|strain_0,j|), ...],
+ max(sum_j(|u0,j|), max(sum_j(|u1,j|)], ...]
+ """
+ _authorized_methods = [FDC4, FDC2]
+
+ def __init__(self, **kwds):
+ super(StrainAndStretch, self).__init__(**kwds)
+ if self._dim == 3:
+ self._fcall = self._call_3d
+ elif self._dim == 2:
+ self._fcall = self._call_2d
+
+ @staticmethod
+ def get_work_properties(topo, indices=None):
+ dime = topo.domain.dimension
+ if dime == 2:
+ lwork = 4
+ elif dime == 3:
+ lwork = 6
+
+ return DifferentialOperation._find_work(lwork, topo, indices)
+
+ def __call__(self, var1, result):
+ """Compute criteria
+
+ Parameters
+ ----------
+ var1 : list of numpy arrays
+ the input vector field
+ result : numpy array
+ criteria values for each direction.
+
+ """
+ if self._on_proc:
+ return self._fcall(var1, result)
+ else:
+ return result
+
+ def _call_3d(self, var1, result):
+ work = self._work
+ assert len(result) == 2 * self._dim
+ # work[0] = abs(eps_x,x)
+ self.fd_scheme.compute(var1[XDIR], XDIR, work[0])
+ np.abs(work[0], work[0])
+
+ # work[1] = diff(u(x), i2)
+ # work[2] = abs(work[1])
+ # work[3] = diff(u(x), i3)
+ # work[4] = abs(work[3])
+ self.fd_scheme.compute(var1[XDIR], YDIR, work[1])
+ np.abs(work[1], work[2])
+ self.fd_scheme.compute(var1[XDIR], ZDIR, work[3])
+ np.abs(work[3], work[4])
+ # max(sum_j(abs(diff(u(x),j))))
+ np.add(work[2], work[4], work[2])
+ np.add(work[2], work[0], work[2])
+ result[3] = np.max(work[2][self.output_indices])
+ # free: p[2, 4, 5]
+ # work[5] = diff(u(y), x)
+ self.fd_scheme.compute(var1[YDIR], XDIR, work[5])
+ # ... used to compute work[1] = abs(eps_x,y)
+ np.add(work[1], work[5], work[1])
+ np.multiply(work[1], 0.5, work[1])
+ np.abs(work[1], work[1])
+ # work[4] = diff(u(z), x)
+ self.fd_scheme.compute(var1[ZDIR], XDIR, work[4])
+ # ... used to compute work[2] = abs(eps_x,z)
+ np.add(work[4], work[3], work[2])
+ np.multiply(work[2], 0.5, work[2])
+ np.abs(work[2], work[2])
+ # At this point,
+ # work[p[0, 1, 2] = abs([epsxx, epsxy, epsxz])
+ # compute max(epsxx + epsxy + epsxz)
+ np.add(work[1], work[0], work[0])
+ np.add(work[2], work[0], work[0])
+ result[0] = np.max(work[0][self.output_indices])
+ # work[4] = abs(diff(u(z), x))
+ # work[5] = abs(diff(u(y), x))
+ np.abs(work[4], work[4])
+ np.abs(work[5], work[5])
+ # free : p[0, 3]
+ # work[0] = abs(eps_y,y)
+ self.fd_scheme.compute(var1[YDIR], YDIR, work[0])
+ np.abs(work[0], work[0])
+ np.add(work[0], work[1], work[1])
+ np.add(work[0], work[5], work[5])
+ # work[0] = abs(diff(u(y), z))
+ self.fd_scheme.compute(var1[YDIR], ZDIR, work[3])
+ np.abs(work[3], work[0])
+ np.add(work[0], work[5], work[5])
+ result[4] = np.max(work[5][self.output_indices])
+ # free : p[0, 5]
+ # work[0] = diff(u(i1), i2)
+ # work[5] = abs(diff(u(i1), i2))
+ self.fd_scheme.compute(var1[ZDIR], YDIR, work[0])
+ np.abs(work[0], work[5])
+ # ... used to compute work[0] = abs(eps_yz)
+ np.add(work[0], work[3], work[0])
+ np.multiply(work[0], 0.5, work[0])
+ np.abs(work[0], work[0])
+ # compute max(epsyy + epsxy + epsyz)
+ np.add(work[1], work[0], work[1])
+ result[1] = np.max(work[1][self.output_indices])
+ np.add(work[2], work[0], work[2])
+ # work[0] = abs(diff(uzz)
+ self.fd_scheme.compute(var1[2], 2, work[0])
+ np.abs(work[0], work[0])
+ np.add(work[2], work[0], work[2])
+ result[2] = np.max(work[2][self.output_indices])
+ np.add(work[0], work[4], work[0])
+ np.add(work[0], work[5], work[0])
+ result[5] = np.max(work[0][self.output_indices])
+ # ouf!
+
+ return result
+
+ def _call_2d(self, var1, result):
+ work = self._work
+ assert len(result) == 2 * self._dim
+ # work[0] = abs(eps_x,x)
+ self.fd_scheme.compute(var1[XDIR], XDIR, work[0])
+ np.abs(work[0], work[0])
+
+ # work[1] = diff(u(x), y)
+ # work[2] = abs(diff(u(x), y))
+ # work[3] = diff(u(y), x)
+ self.fd_scheme.compute(var1[XDIR], YDIR, work[1])
+ np.abs(work[1], work[2])
+ self.fd_scheme.compute(var1[YDIR], XDIR, work[3])
+ np.add(work[0], work[2], work[2])
+ result[2] = np.max(work[2][self.output_indices])
+ # work[1] = epsxy
+ np.add(work[1], work[3], work[1])
+ np.multiply(work[1], 0.5, work[1])
+ np.abs(work[1], work[1])
+ np.add(work[0], work[1], work[0])
+ result[0] = np.max(work[0][self.output_indices])
+ # work[0] = abs(eps_y,y)
+ self.fd_scheme.compute(var1[YDIR], YDIR, work[0])
+ np.abs(work[0], work[0])
+ # work[2] = abs(uy,x)
+ np.abs(work[3], work[2])
+ np.add(work[2], work[0], work[2])
+ result[3] = np.max(work[2][self.output_indices])
+ np.add(work[0], work[1], work[0])
+ result[0] = np.max(work[0][self.output_indices])
+ return result
diff --git a/hysop/numerics/tests/test_diffOp.py b/hysop/numerics/tests/test_diffOp.py
deleted file mode 100755
index 555ea48a3..000000000
--- a/hysop/numerics/tests/test_diffOp.py
+++ /dev/null
@@ -1,222 +0,0 @@
-# -*- coding: utf-8 -*-
-from hysop import Field, Box
-import numpy as np
-import hysop.numerics.differential_operations as diffop
-import hysop.tools.numpywrappers as npw
-from hysop.tools.parameters import Discretization
-import math as m
-pi = m.pi
-cos = np.cos
-sin = np.sin
-
-
-def computeVel(res, x, y, z, t):
- res[0][...] = sin(x) * cos(y) * cos(z)
- res[1][...] = - cos(x) * sin(y) * cos(z)
- res[2][...] = 0.
- return res
-
-
-def computeVort(res, x, y, z, t):
- res[0][...] = - cos(x) * sin(y) * sin(z)
- res[1][...] = - sin(x) * cos(y) * sin(z)
- res[2][...] = 2. * sin(x) * sin(y) * cos(z)
- return res
-
-
-def laplacian_func(res, x, y, z, t):
- res[0][...] = 3 * cos(x) * sin(y) * sin(z)
- res[1][...] = 3 * sin(x) * cos(y) * sin(z)
- res[2][...] = -6. * sin(x) * sin(y) * cos(z)
- return res
-
-
-def computeVel2(res, x, y, t):
- res[0][...] = sin(x) * cos(y)
- res[1][...] = - cos(x) * sin(y)
- return res
-
-
-def computeVort2(res, x, y, t):
- res[0][...] = 2. * sin(x) * sin(y)
- return res
-
-
-def analyticalDivWV(res, x, y, z, t):
- res[0][...] = - sin(y) * cos(y) * sin(z) * cos(z) * \
- (- sin(x) * sin(x) + cos(x) * cos(x)) + \
- 2. * cos(x) * cos(x) * sin(z) * cos(z) * sin(y) * cos(y)
- res[1][...] = - 2. * cos(y) * cos(y) * sin(z) * cos(z) * sin(x) * cos(x) +\
- sin(x) * cos(x) * sin(z) * cos(z) * (cos(y) * cos(y) - sin(y) * sin(y))
- res[2][...] = 0.
- return res
-
-
-def analyticalDivWV2D(res, x, y, t):
- res[0][...] = 0.
- return res
-
-def analyticalGradVxW(res, x, y, z, t):
- res[0][...] = - sin(y) * cos(y) * sin(z) * cos(z)
- res[1][...] = sin(x) * cos(x) * sin(z) * cos(z)
- res[2][...] = 0.
- return res
-
-
-def analyticalDivStressTensor(res, x, y, z, t):
- res[0][...] = - 3. * sin(x) * cos(y) * cos(z)
- res[1][...] = 3. * cos(x) * sin(y) * cos(z)
- res[2][...] = 0.
- return res
-
-def analyticalDivAdvection(res, x, y, z, t):
- res[0][...] = - cos(z) * cos(z) * (cos(x) * cos(x) - sin(x) * sin(x)) - \
- cos(z) * cos(z) * (cos(y) * cos(y) - sin(y) * sin(y))
- return res
-
-
-Nx = 128
-Ny = 110
-Nz = 162
-g = 2
-discr3D = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
-discr2D = Discretization([Nx + 1, Ny + 1], [g, g])
-ldom = npw.asrealarray([pi * 2., ] * 3)
-xdom = [0., 0., 0.]
-
-
-def init(discr, vform, wform):
- dim = len(discr.resolution)
- dom = Box(dimension=dim, length=ldom[:dim],
- origin=xdom[:dim])
- topo = dom.create_topology(discr)
- work = []
- shape = tuple(topo.mesh.resolution)
- for _ in xrange(3):
- work.append(npw.zeros(shape))
- velo = Field(domain=dom, formula=vform,
- name='Velocity', is_vector=True)
- vorti = Field(domain=dom, formula=wform,
- name='Vorticity', is_vector=dim == 3)
- vd = velo.discretize(topo)
- wd = vorti.discretize(topo)
- velo.initialize(topo=topo)
- vorti.initialize(topo=topo)
- return vd, wd, topo, work
-
-
-# Init 2D and 3D
-v3d, w3d, topo3, work3 = init(discr3D, computeVel, computeVort)
-v2d, w2d, topo2, work2 = init(discr2D, computeVel2, computeVort2)
-ic3 = topo3.mesh.compute_index
-ic2 = topo2.mesh.compute_index
-
-
-def build_op(op_class, len_result, topo, work):
- lwk = op_class.get_work_length()
- op = op_class(topo=topo, work=work[:lwk])
- memshape = tuple(topo.mesh.resolution)
- result = [npw.zeros(memshape) for _ in xrange(len_result)]
- return op, result
-
-
-def test_curl():
- op, result = build_op(diffop.Curl, 3, topo3, work3)
- result = op(v3d.data, result)
- rtol = np.max(topo3.mesh.space_step ** 2)
- for i in xrange(3):
- assert np.allclose(w3d.data[i][ic3], result[i][ic3], rtol=rtol)
-
-
-def test_curl_2d():
- op, result = build_op(diffop.Curl, 1, topo2, work2)
- result = op(v2d.data, result)
- rtol = np.max(topo3.mesh.space_step ** 2)
- assert np.allclose(w2d.data[0][ic2], result[0][ic2], rtol=rtol)
-
-
-def test_laplacian():
- ref = Field(domain=topo3.domain, formula=laplacian_func,
- name='Analytical', is_vector=True)
- rd = ref.discretize(topo3)
- ref.initialize(topo=topo3)
-
- op, result = build_op(diffop.Laplacian, 3, topo3, work3)
- result = op(w3d.data, result)
- tol = np.max(topo3.mesh.space_step ** 2)
- for i in xrange(3):
- assert np.allclose(rd.data[i][ic3], result[i][ic3], rtol=tol)
-
-
-def test_div_rho_v_2d():
- # Reference field
- ref = Field(domain=topo2.domain, formula=analyticalDivWV2D,
- name='Analytical', is_vector=False)
- rd = ref.discretize(topo2)
- ref.initialize(topo=topo2)
- op, result = build_op(diffop.DivRhoV, 1, topo2, work2)
- tol = np.max(topo2.mesh.space_step ** 4)
- result = op(v2d.data, w2d.data[0:1], result)
- # Numerical VS analytical
- assert np.allclose(rd[0][ic2], result[0][ic2], atol=tol)
-
-
-def test_div_rho_v():
- # Reference field
- ref = Field(domain=topo3.domain, formula=analyticalDivWV,
- name='Analytical', is_vector=True)
- rd = ref.discretize(topo3)
- ref.initialize(topo=topo3)
- op, result = build_op(diffop.DivRhoV, 3, topo3, work3)
- tol = np.max(topo3.mesh.space_step ** 4)
- for i in xrange(3):
- result[i:i + 1] = op(v3d.data, w3d.data[i:i + 1], result[i:i + 1])
- # Numerical VS analytical
- assert np.allclose(rd[i][ic3], result[i][ic3], atol=tol)
-
-
-def test_divwv():
- # Reference field
- ref = Field(domain=topo3.domain, formula=analyticalDivWV,
- name='Analytical', is_vector=True)
- rd = ref.discretize(topo3)
- ref.initialize(topo=topo3)
- op, result = build_op(diffop.DivWV, 3, topo3, work3)
- result = op(v3d.data, w3d.data, result)
-
- # Numerical VS analytical
- tol = np.max(topo3.mesh.space_step ** 2)
- for i in xrange(3):
- assert np.allclose(rd[i][ic3], result[i][ic3], atol=tol)
-
-
-def test_grad_vxw():
- # Reference field
- ref = Field(domain=topo3.domain, formula=analyticalGradVxW,
- name='Analytical', is_vector=True)
- rd = ref.discretize(topo3)
- ref.initialize(topo=topo3)
- op, result = build_op(diffop.GradVxW, 3, topo3, work3)
- diag = npw.zeros(2)
- result, diag = op(v3d.data, w3d.data, result, diag)
-
- # Numerical VS analytical
- tol = np.max(topo3.mesh.space_step ** 2)
- for i in xrange(3):
- assert np.allclose(rd[i][ic3], result[i][ic3], atol=tol)
-
-
-# test for RHS of Pressure-Poisson equation
-def test_div_advection():
- # Reference scalar field
- ref = Field(domain=topo3.domain, formula=analyticalDivAdvection,
- name='Analytical')
- rd = ref.discretize(topo3)
- ref.initialize(topo=topo3)
- op, result = build_op(diffop.DivAdvection, 1, topo3, work3)
- result = op(v3d.data, result)
-
- # Numerical VS analytical
- errx = (topo3.domain.length[0] / (Nx - 1)) ** 4
- assert np.allclose(rd[0][ic3], result[0][ic3], rtol=errx)
-
diff --git a/hysop/numerics/tests/test_differential_operations.py b/hysop/numerics/tests/test_differential_operations.py
new file mode 100755
index 000000000..35139d349
--- /dev/null
+++ b/hysop/numerics/tests/test_differential_operations.py
@@ -0,0 +1,1015 @@
+# -*- coding: utf-8 -*-
+from hysop import Field, Box
+import numpy as np
+import hysop.numerics.differential_operations as diffop
+import hysop.tools.numpywrappers as npw
+from hysop.tools.parameters import Discretization
+from hysop.domain.subsets import SubBox
+import math as m
+pi = m.pi
+cos = np.cos
+sin = np.sin
+
+
+def compute_vel(res, x, y, z, t):
+ res[0][...] = sin(x) * cos(y) * cos(z)
+ res[1][...] = - cos(x) * sin(y) * cos(z)
+ res[2][...] = 0.
+ return res
+
+
+def compute_vort(res, x, y, z, t):
+ res[0][...] = - cos(x) * sin(y) * sin(z)
+ res[1][...] = - sin(x) * cos(y) * sin(z)
+ res[2][...] = 2. * sin(x) * sin(y) * cos(z)
+ return res
+
+
+def compute_vel2(res, x, y, t):
+ res[0][...] = sin(x) * cos(y)
+ res[1][...] = - cos(x) * sin(y)
+ return res
+
+
+def compute_vort2(res, x, y, t):
+ res[0][...] = 2. * sin(x) * sin(y)
+ return res
+
+
+def analyticalDivRhoV(res, x, y, z, t):
+ res[0][...] = cos(y) * cos(z) * sin(y) * sin(z)
+ return res
+
+
+def analyticalDivRhoV2D(res, x, y, t):
+ res[0][...] = 0.
+ return res
+
+
+def analyticalDivWV(res, x, y, z, t):
+ res[0][...] = - cos(y) * cos(z) * sin(y) * sin(z)
+ res[1][...] = cos(x) * cos(z) * sin(z) * sin(x)
+ res[2][...] = 0.
+ return res
+
+
+def analyticalGradVxW(res, x, y, z, t):
+ res[0][...] = - sin(y) * cos(y) * sin(z) * cos(z)
+ res[1][...] = sin(x) * cos(x) * sin(z) * cos(z)
+ res[2][...] = 0.
+ return res
+
+
+def analyticalDivStressTensor(res, x, y, z, t):
+ res[0][...] = - 3. * sin(x) * cos(y) * cos(z)
+ res[1][...] = 3. * cos(x) * sin(y) * cos(z)
+ res[2][...] = 0.
+ return res
+
+def analyticalDivAdvection(res, x, y, z, t):
+ res[0][...] = - cos(z) * cos(z) * (cos(x) * cos(x) - sin(x) * sin(x)) - \
+ cos(z) * cos(z) * (cos(y) * cos(y) - sin(y) * sin(y))
+ return res
+
+
+Nx = 88
+Ny = 50
+Nz = 64
+g = 2
+discr3D = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+discr2D = Discretization([Nx + 1, Ny + 1], [g, g])
+ldom = npw.asrealarray([2., pi, 3.5])
+xdom = [0., -0.3, 0.4]
+
+
+def init(discr, vform, wform):
+ dim = len(discr.resolution)
+ dom = Box(dimension=dim, length=ldom[:dim],
+ origin=xdom[:dim])
+ topo = dom.create_topology(discr)
+ work = []
+ shape = tuple(topo.mesh.resolution)
+ subsize = np.prod(shape)
+ for _ in xrange(6):
+ work.append(npw.zeros(subsize,))
+ velo = Field(domain=dom, formula=vform,
+ name='Velocity', is_vector=True)
+ vorti = Field(domain=dom, formula=wform,
+ name='Vorticity', is_vector=dim == 3)
+ vd = velo.discretize(topo)
+ wd = vorti.discretize(topo)
+ velo.initialize(topo=topo)
+ vorti.initialize(topo=topo)
+ return vd, wd, topo, work
+
+
+# Init 2D and 3D
+v3d, w3d, topo3, work3 = init(discr3D, compute_vel, compute_vort)
+v2d, w2d, topo2, work2 = init(discr2D, compute_vel2, compute_vort2)
+
+
+def build_op(op_class, len_result, topo, work=None, reduced=False,
+ reduce_output=False):
+ if reduced:
+ # computations will be done on a subset of the domain
+ sl = topo.domain.length * 0.5
+ orig = topo.domain.origin + 0.1 * topo.domain.length
+ subbox = SubBox(parent=topo.domain, origin=orig, length=sl)
+ indices = subbox.discretize(topo)[0]
+ if reduce_output:
+ outputshape = subbox.mesh[topo].resolution
+ else:
+ outputshape = np.asarray(topo.mesh.resolution).copy()
+ else:
+ indices = None
+ if reduce_output:
+ outputshape = np.asarray(topo.mesh.resolution).copy()
+ gh = topo.ghosts()
+ outputshape -= 2 * gh
+ outputshape = tuple(outputshape)
+ else:
+ outputshape = tuple(topo.mesh.resolution)
+
+ if work is not None:
+ wk_prop = op_class.get_work_properties(topo, indices)
+ if wk_prop['rwork'] is not None:
+ lwork = len(wk_prop['rwork'])
+ else:
+ lwork = 0
+ op = op_class(topo=topo, indices=indices, work=work[:lwork],
+ reduce_output_shape=reduce_output)
+ else:
+ op = op_class(topo=topo, indices=indices,
+ reduce_output_shape=reduce_output)
+ result = [npw.zeros(outputshape) for _ in xrange(len_result)]
+ return op, result
+
+
+def assert_curl(wdim, topo, velo, vorti, work=None, reduced=False,
+ reduce_output=False):
+ op, result = build_op(diffop.Curl, wdim, topo, work, reduced,
+ reduce_output)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ result = op(velo.data, result)
+ atol = np.max(topo.mesh.space_step ** 2)
+ for i in xrange(wdim):
+ assert np.allclose(vorti.data[i][ic_in], result[i][ic_out], atol=atol)
+
+
+def test_curl():
+ # 2d
+ assert_curl(1, topo2, v2d, w2d)
+ # 2d work
+ assert_curl(1, topo2, v2d, w2d, work2)
+ # 2d reduced
+ assert_curl(1, topo2, v2d, w2d, reduced=True)
+ # 2d reduced and work
+ assert_curl(1, topo2, v2d, w2d, work2, reduced=True)
+ # 2d reduce output
+ assert_curl(1, topo2, v2d, w2d, reduce_output=True)
+ # 2d reduce output and work
+ assert_curl(1, topo2, v2d, w2d, work2, reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_curl(1, topo2, v2d, w2d, work2, reduced=True, reduce_output=True)
+ # 3d
+ assert_curl(3, topo3, v3d, w3d)
+ # 3d work
+ assert_curl(3, topo3, v3d, w3d, work3)
+ # 3d reduced
+ assert_curl(3, topo3, v3d, w3d, reduced=True)
+ # 3d reduced and work
+ assert_curl(3, topo3, v3d, w3d, work3, reduced=True)
+ # 3d reduce output
+ assert_curl(3, topo3, v3d, w3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_curl(3, topo3, v3d, w3d, work3, reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_curl(3, topo3, v3d, w3d, work3, reduced=True, reduce_output=True)
+
+
+def assert_div_rho_v(topo, velo, scal, formula, work=None, reduced=False,
+ reduce_output=False):
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', is_vector=False)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.DivRhoV, 1, topo, work, reduced,
+ reduce_output)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ tol = np.max(topo.mesh.space_step ** 4)
+ result = op(velo.data, [scal], result)
+ assert np.allclose(rd[0][ic_in], result[0][ic_out], atol=tol)
+
+
+def test_div_rho_v():
+ # 2d
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D)
+ # 2d work
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D, work2)
+ # 2d reduced
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D, reduced=True)
+ # 2d reduced and work
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D,
+ reduce_output=False)
+ # 2d reduce output and work
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D, work2,
+ reduce_output=False)
+ # 2d reduce output, on subset, with work
+ assert_div_rho_v(topo2, v2d, w2d[0], analyticalDivRhoV2D, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV)
+ # 3d work
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV, work3)
+ # 3d reduced
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV, reduced=True)
+ # 3d reduced and work
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV,
+ reduce_output=False)
+ # 3d reduce output and work
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV, work3,
+ reduce_output=False)
+ # 3d, work, reduce output, on subset
+ assert_div_rho_v(topo3, v3d, w3d[0], analyticalDivRhoV, work3,
+ reduce_output=False, reduced=True)
+
+
+def assert_div_w_v(wdim, topo, velo, vorti, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference field
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.DivWV, wdim, topo, work, reduced,
+ reduce_output)
+ result = op(vorti.data, velo.data, result)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+
+ # Numerical VS analytical
+ tol = np.max(topo.mesh.space_step ** 2)
+ for i in xrange(wdim):
+ assert np.allclose(rd[i][ic_in], result[i][ic_out], atol=tol)
+
+
+def test_div_w_v():
+ # 3d
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV)
+ # 3d work
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, work3)
+ # 3d reduced
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, reduced=True)
+ # 3d reduced and work
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, reduce_output=True)
+ # 3d reduce output and work
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_div_w_v(3, topo3, v3d, w3d, analyticalDivWV, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_laplacian(wdim, topo, vorti, formula, work=None, reduced=False,
+ reduce_output=False):
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.Laplacian, wdim, topo, work, reduced,
+ reduce_output)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ result = op(vorti.data, result)
+ tol = np.max(topo.mesh.space_step ** 2)
+ for i in xrange(wdim):
+ assert np.allclose(rd.data[i][ic_in], result[i][ic_out], rtol=tol)
+
+
+def laplacian_func(res, x, y, z, t):
+ res[0][...] = 3 * cos(x) * sin(y) * sin(z)
+ res[1][...] = 3 * sin(x) * cos(y) * sin(z)
+ res[2][...] = -6. * sin(x) * sin(y) * cos(z)
+ return res
+
+
+def laplacian_func_2d(res, x, y, t):
+ res[0][...] = -4. * sin(x) * sin(y)
+ return res
+
+
+def test_laplacian():
+ # 2d
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d)
+ # 2d work
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, work2)
+ # 2d reduced
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, reduce_output=True)
+ # 2d reduce output and work
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, work2,
+ reduce_output=True)
+ # 2d, work, reduce output, on subset
+ assert_laplacian(1, topo2, w2d, laplacian_func_2d, work2,
+ reduce_output=True, reduced=True)
+
+ # 3d
+ assert_laplacian(3, topo3, w3d, laplacian_func)
+ # 3d work
+ assert_laplacian(3, topo3, w3d, laplacian_func, work3)
+
+ # 3d reduced
+ assert_laplacian(3, topo3, w3d, laplacian_func, reduced=True)
+ # 3d reduced and work
+ assert_laplacian(3, topo3, w3d, laplacian_func, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_laplacian(3, topo3, w3d, laplacian_func, reduce_output=True)
+ # 3d reduce output and work
+ assert_laplacian(3, topo3, w3d, laplacian_func, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_laplacian(3, topo3, w3d, laplacian_func, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_grad_s(topo, scal, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference field
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.GradS, topo.domain.dimension,
+ topo, work, reduced, reduce_output)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ result = op([scal], result)
+ # Numerical VS analytical
+ tol = np.max(topo.mesh.space_step ** 4)
+ for i in xrange(topo.domain.dimension):
+ assert np.allclose(rd[i][ic_in], result[i][ic_out], atol=tol)
+
+
+def grad_s_func_2d(res, x, y, t):
+ res[0][...] = cos(x) * cos(y)
+ res[1][...] = -sin(x) * sin(y)
+ return res
+
+
+def grad_s_func_3d(res, x, y, z, t):
+ res[0][...] = cos(x) * cos(y) * cos(z)
+ res[1][...] = -sin(x) * sin(y) * cos(z)
+ res[2][...] = -sin(x) * cos(y) * sin(z)
+ return res
+
+
+def test_grad_s():
+ # 2d
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d)
+ # 2d work
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d, work2)
+ # 2d reduced
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_grad_s(topo2, v2d[0], grad_s_func_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d)
+ # 3d work
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, work3)
+ # 3d reduced
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, reduced=True)
+ # 3d reduced and work
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_grad_s(topo3, v3d[0], grad_s_func_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_grad_v(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference field
+ resdim = topo.domain.dimension ** 2
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.GradV, resdim, topo, work, reduced,
+ reduce_output)
+ result = op(velo.data, result)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ # Numerical VS analytical
+ tol = np.max(topo.mesh.space_step ** 4)
+ for i in xrange(resdim):
+ assert np.allclose(rd[i][ic_in], result[i][ic_out], atol=tol)
+
+
+def grad_v_func_2d(res, x, y, t):
+ res[0][...] = cos(x) * cos(y)
+ res[1][...] = -sin(x) * sin(y)
+ res[2][...] = sin(x) * sin(y)
+ res[3][...] = -cos(x) * cos(y)
+ return res
+
+
+def grad_v_func_3d(res, x, y, z, t):
+ res[0][...] = cos(x) * cos(y) * cos(z)
+ res[1][...] = -sin(x) * sin(y) * cos(z)
+ res[2][...] = -sin(x) * cos(y) * sin(z)
+ res[3][...] = sin(x) * sin(y) * cos(z)
+ res[4][...] = - cos(x) * cos(y) * cos(z)
+ res[5][...] = cos(x) * sin(y) * sin(z)
+ res[6][...] = 0.
+ res[7][...] = 0.
+ res[8][...] = 0.
+ return res
+
+
+def test_grad_v():
+ # 2d
+ assert_grad_v(topo2, v2d, grad_v_func_2d)
+ # 2d work
+ assert_grad_v(topo2, v2d, grad_v_func_2d, work2)
+ # 2d reduced
+ assert_grad_v(topo2, v2d, grad_v_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_grad_v(topo2, v2d, grad_v_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_grad_v(topo2, v2d, grad_v_func_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_grad_v(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_grad_v(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_grad_v(topo3, v3d, grad_v_func_3d)
+ # 3d work
+ assert_grad_v(topo3, v3d, grad_v_func_3d, work3)
+ # 3d reduced
+ assert_grad_v(topo3, v3d, grad_v_func_3d, reduced=True)
+ # 3d reduced and work
+ assert_grad_v(topo3, v3d, grad_v_func_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_grad_v(topo3, v3d, grad_v_func_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_grad_v(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_grad_v(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_grad_vxw(topo, velo, vorti, work=None, reduced=False,
+ reduce_output=False):
+ # Reference field
+ ref = Field(domain=topo.domain, formula=analyticalGradVxW,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.GradVxW, velo.nb_components, topo,
+ work, reduced, reduce_output)
+ diag = npw.zeros(2)
+ result, diag = op(velo.data, vorti.data, result, diag)
+ # Numerical VS analytical
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ tol = np.max(topo3.mesh.space_step ** 2)
+ for i in xrange(topo.domain.dimension):
+ assert np.allclose(rd[i][ic_in], result[i][ic_out], atol=tol)
+
+
+def test_grad_vxw():
+ # 3d
+ assert_grad_vxw(topo3, v3d, w3d)
+ # 3d work
+ assert_grad_vxw(topo3, v3d, w3d, work3)
+ # 3d reduced
+ assert_grad_vxw(topo3, v3d, w3d, reduced=True)
+ # 3d reduced and work
+ assert_grad_vxw(topo3, v3d, w3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_grad_vxw(topo3, v3d, w3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_grad_vxw(topo3, v3d, w3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_grad_vxw(topo3, v3d, w3d, work3,
+ reduce_output=True, reduced=True)
+
+
+# test for RHS of Pressure-Poisson equation
+def assert_div_advection(topo, velo, work=None, reduced=False,
+ reduce_output=False):
+ # Reference scalar field
+ ref = Field(domain=topo.domain, formula=analyticalDivAdvection,
+ name='Analytical')
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.DivAdvection, 1, topo, work, reduced,
+ reduce_output)
+ result = op(velo.data, result)
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+
+ # Numerical VS analytical
+ tol = np.max(topo.mesh.space_step ** 4)
+ assert np.allclose(rd[0][ic_in], result[0][ic_out], rtol=tol)
+
+
+def test_div_advection():
+ # 3d
+ assert_div_advection(topo3, v3d)
+ # 3d work
+ assert_div_advection(topo3, v3d, work3)
+ # 3d reduced
+ assert_div_advection(topo3, v3d, reduced=True)
+ # 3d reduced and work
+ assert_div_advection(topo3, v3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_div_advection(topo3, v3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_div_advection(topo3, v3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_div_advection(topo3, v3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_strain(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference scalar field
+ resdim = topo.domain.dimension * (topo.domain.dimension + 1) / 2
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.Strain, resdim, topo, work, reduced,
+ reduce_output)
+ result = op(velo.data, result)
+
+ ic_out = op.output_indices
+ ic_in = op.in_indices
+ # Numerical VS analytical
+ tol = np.max(topo.mesh.space_step ** 4)
+ for d in xrange(resdim):
+ assert np.allclose(rd[d][ic_in], result[d][ic_out], atol=tol)
+
+
+def strain_f_2d(res, x, y, t):
+ res[0][...] = cos(x) * cos(y)
+ res[1][...] = - cos(x) * cos(y)
+ res[2][...] = 0.
+ return res
+
+
+def strain_f_3d(res, x, y, z, t):
+ res[0][...] = cos(x) * cos(y) * cos(z)
+ res[1][...] = - cos(x) * cos(y) * cos(z)
+ res[2][...] = 0.
+ res[3][...] = 0.
+ res[4][...] = -0.5 * sin(x) * cos(y) * sin(z)
+ res[5][...] = 0.5 * cos(x) * sin(y) * sin(z)
+ return res
+
+
+def test_strain():
+ # 2d
+ assert_strain(topo2, v2d, strain_f_2d)
+ # 2d work
+ assert_strain(topo2, v2d, strain_f_2d, work2)
+ # 2d reduced
+ assert_strain(topo2, v2d, strain_f_2d, reduced=True)
+ # 2d reduced and work
+ assert_strain(topo2, v2d, strain_f_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_strain(topo2, v2d, strain_f_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_strain(topo2, v2d, strain_f_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_strain(topo2, v2d, strain_f_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_strain(topo3, v3d, strain_f_3d)
+ # 3d work
+ assert_strain(topo3, v3d, strain_f_3d, work3)
+ # 3d reduced
+ assert_strain(topo3, v3d, strain_f_3d, reduced=True)
+ # 3d reduced and work
+ assert_strain(topo3, v3d, strain_f_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_strain(topo3, v3d, strain_f_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_strain(topo3, v3d, strain_f_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_strain(topo3, v3d, strain_f_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def compute_max(field, ic):
+ """Compute max of sum_i (abs(field_i[ic])) for i in len(field)
+ """
+ temp = npw.zeros_like(field[0])
+ for d in xrange(len(field)):
+ temp[ic] += np.abs(field[d][ic])
+ return np.max(temp)
+
+
+def assert_strain_criteria(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference scalar field
+ resdim = topo.domain.dimension * (topo.domain.dimension + 1) / 2
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.StrainCriteria, 0, topo, work, reduced,
+ reduce_output)
+ result = npw.zeros(topo.domain.dimension)
+ result = op(velo.data, result)
+ ic_in = op.in_indices
+
+ dd = topo.domain.dimension
+ if dd == 3:
+ ll = [None, ] * dd
+ ll[0] = [rd[0], rd[3], rd[4]]
+ ll[1] = [rd[1], rd[3], rd[5]]
+ ll[2] = [rd[2], rd[4], rd[5]]
+ elif dd == 2:
+ ll = [None, ] * dd
+ ll[0] = [rd[0], rd[2]]
+ ll[1] = [rd[1], rd[2]]
+
+ for d in xrange(dd):
+ refmax = compute_max(ll[d], ic_in)
+ assert np.allclose(refmax, result[d])
+
+
+def test_strain_criteria():
+ # 2d
+ assert_strain_criteria(topo2, v2d, strain_f_2d)
+ # 2d work
+ assert_strain_criteria(topo2, v2d, strain_f_2d, work2)
+ # 2d reduced
+ assert_strain_criteria(topo2, v2d, strain_f_2d, reduced=True)
+ # 2d reduced and work
+ assert_strain_criteria(topo2, v2d, strain_f_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_strain_criteria(topo2, v2d, strain_f_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_strain_criteria(topo2, v2d, strain_f_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_strain_criteria(topo2, v2d, strain_f_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_strain_criteria(topo3, v3d, strain_f_3d)
+ # 3d work
+ assert_strain_criteria(topo3, v3d, strain_f_3d, work3)
+ # 3d reduced
+ assert_strain_criteria(topo3, v3d, strain_f_3d, reduced=True)
+ # 3d reduced and work
+ assert_strain_criteria(topo3, v3d, strain_f_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_strain_criteria(topo3, v3d, strain_f_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_strain_criteria(topo3, v3d, strain_f_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_strain_criteria(topo3, v3d, strain_f_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_max_diag_grad(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference field
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.MaxDiagGradV, 0, topo, work, reduced,
+ reduce_output)
+ result = npw.zeros(topo.domain.dimension)
+ result = op(velo.data, result)
+ ic_in = op.in_indices
+ for d in xrange(topo.domain.dimension):
+ refmax = compute_max([rd[d]], ic_in)
+ assert np.allclose(refmax, result[d])
+
+
+def diag_grad_func_2d(res, x, y, t):
+ res[0][...] = cos(x) * cos(y)
+ res[1][...] = -cos(x) * cos(y)
+ return res
+
+
+def diag_grad_func_3d(res, x, y, z, t):
+ res[0][...] = cos(x) * cos(y) * cos(z)
+ res[1][...] = - cos(x) * cos(y) * cos(z)
+ res[2][...] = 0.
+ return res
+
+
+def test_max_diag_grad():
+ # 2d
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d)
+ # 2d work
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d, work2)
+ # 2d reduced
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_max_diag_grad(topo2, v2d, diag_grad_func_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d)
+ # 3d work
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, work3)
+ # 3d reduced
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, reduced=True)
+ # 3d reduced and work
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_max_diag_grad(topo3, v3d, diag_grad_func_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_stretch_like(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference scalar field
+ resdim = topo.domain.dimension ** 2
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.StretchLike, 0, topo, work, reduced,
+ reduce_output)
+ result = npw.zeros(topo.domain.dimension)
+ result = op(velo.data, result)
+ ic_in = op.in_indices
+
+ dd = topo.domain.dimension
+ for d in xrange(dd):
+ pos = d * dd
+ refmax = compute_max(rd[pos: pos + dd], ic_in)
+ assert np.allclose(refmax, result[d])
+
+
+def test_stretch_like():
+ # 2d
+ assert_stretch_like(topo2, v2d, grad_v_func_2d)
+ # 2d work
+ assert_stretch_like(topo2, v2d, grad_v_func_2d, work2)
+ # 2d reduced
+ assert_stretch_like(topo2, v2d, grad_v_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_stretch_like(topo2, v2d, grad_v_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_stretch_like(topo2, v2d, grad_v_func_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_stretch_like(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_stretch_like(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_stretch_like(topo3, v3d, grad_v_func_3d)
+ # 3d work
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, work3)
+ # 3d reduced
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, reduced=True)
+ # 3d reduced and work
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_stretch_like(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_diag_and_stretch(topo, velo, formula, work=None, reduced=False,
+ reduce_output=False):
+ # Reference scalar field
+ resdim = topo.domain.dimension ** 2
+ ref = Field(domain=topo.domain, formula=formula,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(topo=topo)
+ op, result = build_op(diffop.DiagAndStretch, 0, topo, work, reduced,
+ reduce_output)
+ result = npw.zeros(2 * topo.domain.dimension)
+ result = op(velo.data, result)
+ ic_in = op.in_indices
+
+ dd = topo.domain.dimension
+ ll = [rd[i * (dd + 1)] for i in xrange(dd)]
+
+ for d in xrange(dd):
+ refmax = compute_max([ll[d]], ic_in)
+ assert np.allclose(refmax, result[d])
+ pos = d * dd
+ refmax = compute_max(rd[pos: pos + dd], ic_in)
+ assert np.allclose(refmax, result[d + dd])
+
+
+def test_diag_and_stretch():
+ # 2d
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d)
+ # 2d work
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d, work2)
+ # 2d reduced
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d, reduced=True)
+ # 2d reduced and work
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d, work2,
+ reduced=True)
+ # 2d reduce output
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_diag_and_stretch(topo2, v2d, grad_v_func_2d, work2,
+ reduce_output=True, reduced=True)
+ # 3d
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d)
+ # # 3d work
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, work3)
+ # 3d reduced
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, reduced=True)
+ # 3d reduced and work
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, work3,
+ reduced=True)
+ # 3d reduce output
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, reduce_output=True)
+ # 3d reduce output and work
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_diag_and_stretch(topo3, v3d, grad_v_func_3d, work3,
+ reduce_output=True, reduced=True)
+
+
+def assert_strain_stretch_criteria(topo, velo, formula1, formula2, work=None,
+ reduced=False, reduce_output=False):
+ # Reference scalar field
+ resdim = topo.domain.dimension * (topo.domain.dimension + 1) / 2
+ graddim = topo.domain.dimension ** 2
+ refstrain = Field(domain=topo.domain, formula=formula2,
+ name='Strain', nb_components=resdim)
+ refgrad = Field(domain=topo.domain, formula=formula1,
+ name='Gradl', nb_components=graddim)
+
+ rd = refstrain.discretize(topo)
+ gd = refgrad.discretize(topo)
+ refstrain.initialize(topo)
+ refgrad.initialize(topo=topo)
+ op, result = build_op(diffop.StrainAndStretch, 0, topo, work, reduced,
+ reduce_output)
+ result = npw.zeros(2 * topo.domain.dimension)
+ result = op(velo.data, result)
+ ic_in = op.in_indices
+
+ dd = topo.domain.dimension
+ if dd == 3:
+ ll = [None, ] * dd
+ ll[0] = [rd[0], rd[3], rd[4]]
+ ll[1] = [rd[1], rd[3], rd[5]]
+ ll[2] = [rd[2], rd[4], rd[5]]
+ elif dd == 2:
+ ll = [None, ] * dd
+ ll[0] = [rd[0], rd[2]]
+ ll[1] = [rd[1], rd[2]]
+
+ if dd == 3:
+ for d in xrange(dd):
+ refmaxstr = compute_max(ll[d], ic_in)
+ assert np.allclose(refmaxstr, result[d])
+ pos = d * dd
+ refmaxgrad = compute_max(gd[pos: pos + dd], ic_in)
+ assert np.allclose(refmaxgrad, result[d + dd])
+
+
+def test_strain_stretch():
+ # 2d
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d)
+ # 2d work
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ work2)
+ # 2d reduced
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ reduced=True)
+ # 2d reduced and work
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ work2, reduced=True)
+ # 2d reduce output
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ reduce_output=True)
+ # 2d reduce output and work
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ work2, reduce_output=True)
+ # 2d reduce output, on subset, with work
+ assert_strain_stretch_criteria(topo2, v2d, grad_v_func_2d, strain_f_2d,
+ work2, reduce_output=True, reduced=True)
+ # 3d
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d)
+ # 3d work
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ work3)
+ # 3d reduced
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ reduced=True)
+ # 3d reduced and work
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ work3, reduced=True)
+ # 3d reduce output
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ reduce_output=True)
+ # 3d reduce output and work
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ work3, reduce_output=True)
+ # 3d, work, reduce output, on subset
+ assert_strain_stretch_criteria(topo3, v3d, grad_v_func_3d, strain_f_3d,
+ work3, reduce_output=True, reduced=True)
+
+if __name__ == "__main__":
+ test_curl()
+ test_div_rho_v()
+ test_div_w_v()
+ test_laplacian()
+ test_grad_s()
+ test_grad_v()
+ test_grad_vxw()
+ test_div_advection()
+ test_strain()
+ test_strain_criteria()
+ test_max_diag_grad()
+ test_stretch_like()
+ test_diag_and_stretch()
+ test_strain_criteria()
+ test_strain_stretch()
--
GitLab