diff --git a/docs/sphinx/users_guide/finite_differences.rst b/docs/sphinx/users_guide/finite_differences.rst
new file mode 100644
index 0000000000000000000000000000000000000000..81d277931d12716dff1f41adee43ac36cf2b6e3f
--- /dev/null
+++ b/docs/sphinx/users_guide/finite_differences.rst
@@ -0,0 +1,199 @@
+.. _finite_differences:
+
+.. currentmodule:: hysop.numerics.finite_differences
+
+Finite differences schemes
+--------------------------
+
+Differentiate some fields in one direction using finite differences.
+
+So, to compute
+
+.. math::
+ \begin{eqnarray}
+ result &=& \frac{\partial \rho}{\partial y}
+ \end{eqnarray}
+
+if tab is a numpy array representing the discrete values of the scalar field :math:`\rho`
+on a grid, then the basic usage of such schemes is :
+
+.. code::
+
+ # Build/declare the scheme for a given space discretization
+ scheme = FDC2(step, indices)
+ # Apply scheme on the array
+ dir = 1 # y direction
+ result = scheme(tab, dir, result)
+
+This will compute :
+
+result = diff(tab[indices], dir), diff depending on
+the chosen scheme.
+
+
+Available schemes
+^^^^^^^^^^^^^^^^^
+
+* :class:`~FDC2` : first derivative, 2nd order centered scheme
+* :class:`~FD2C2`: second derivative, 2nd order centered scheme
+* :class:`~FDC4`: first derivative, 4th order centered scheme
+
+
+A few important remarks
+^^^^^^^^^^^^^^^^^^^^^^^
+
+* step is the space discretization step size in each direction, i.e. a list or numpy array
+ with d values, d being the dimension of the domain. A common case is :code::
+
+ step = topo.mesh.space_step
+
+* indices represent the set of points on which the scheme must be applied. This is usually
+ a list of slices, for example, :code::
+
+ indices = topo.mesh.iCompute
+
+* result must be a predefined numpy array of the 'right' size, here the same size/shape as tab.
+
+* To optimize memory usage and computational time, it's possible to reduce the size of
+ the output and/or to apply the scheme on a subset of the domain. All available possibilities
+ are summarized through the examples below.
+
+* the size of the ghost layer depends on the scheme but is not checked! You must ensure
+ that topo.ghosts() >= scheme.minimal_ghost_layer.
+
+* for some schemes a work array may be provided during call. It must be a numpy array of
+ the same size of the result. It's shape is not really important, since during call
+ work will be reshaped to be of the same shape as result. This allows us to provide 1D array
+ that may be shared between different operators/methods whatever the internal required shapes are.
+
+* Notice that most of the time, finite-differences are defined as internal methods of operators and work arrays management, indices list or ghosts layers are set/checked internally.
+
+Default case : apply scheme on all the domain with a full-size result
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ from hysop.numerics.finite_differences import FDC2
+ box = Box(length=[1., 1.])
+ d2d = Discretization([33, 33], [1, 1])
+ topo = box.create_topology(d2d)
+ rho = Field(box, name='rho')
+ # discretize and initialize rho
+ rd = rho.randomize(topo)
+ # Get 'computational' points, i.e. grid points excluding ghosts
+ ic = topo.mesh.iCompute
+ # space step
+ step = topo.mesh.space_step
+ # field resolution on the grid defined by topo
+ shape = topo.mesh.resolution
+ result = npw.zeros(shape)
+ scheme = FDC2(step, ic)
+ assert (topo.ghosts() >= scheme.minimal_ghost_layer).all()
+ result = scheme(rd.data[0], 1, result)
+
+In that case:
+
+* result.shape = (34, 34)
+* scheme.indices = [slice(1,33), slice(1,33)]
+* scheme.output_indices = [slice(0,32), slice(0,32)]
+
+
+Apply scheme on all the domain with a reduced result
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If you do not want to allocate ghosts points for the result, then use:
+
+.. code::
+
+ shape = np.asarray(topo.mesh.resolution).copy()
+ shape -= 2 * topo.ghosts()
+ shape = tuple(shape)
+ result = npw.zeros(shape)
+ scheme = FDC2(step, ic, indices_out=True)
+ assert (topo.ghosts() >= scheme.minimal_ghost_layer).all()
+ result = scheme(rd.data[0], 1, result)
+
+In that case:
+
+* result.shape = (32,32)
+* scheme.indices = [slice(1,33), slice(1,33)]
+* scheme.output_indices = [slice(0,32), slice(0,32)]
+
+
+Apply scheme on a subset of the domain with a full-size result
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ # We define a subset of the domain,
+ # a box of size half of the domain size.
+ from hysop.domain.subsets import SubBox
+ 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]
+ scheme = FDC2(step, indices)
+ result = npw.zeros_like(rd.data[0])
+ result = scheme(rd.data[0], 1, result)
+
+In that case:
+
+* result.shape = (34,34)
+* scheme.indices = [slice(4,21), slice(4,21)]
+* scheme.output_indices = scheme.indices
+
+
+Apply scheme on a subset of the domain with a reduced-size result
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ # We define a subset of the domain,
+ # a box of size half of the domain size.
+ from hysop.domain.subsets import SubBox
+ 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]
+ scheme = FDC2(step, indices, indices_out=True)
+ shape = subbox.mesh[topo].resolution
+ result = npw.zeros(shape)
+ result = scheme(rd.data[0], 1, result)
+
+In that case:
+
+* result.shape = (17,17)
+* scheme.indices = [slice(4,21), slice(4,21)]
+* scheme.output_indices = [slice(0,17), slice(0,17)]scheme.indices
+
+
+Apply scheme on a subset of the domain with a predifined-size result
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Here you can choose explicitely where in result you want to put
+the values of the computed derivative.
+
+
+.. code::
+
+ # We define a subset of the domain,
+ # a box of size half of the domain size.
+ from hysop.domain.subsets import SubBox
+ 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]
+ shape = subbox.mesh[topo].resolution
+ result = npw.zeros_like(rd.data[0)
+ iout = [slice(2, 19), slice(3, 20)]
+ scheme = FDC2(step, indices, indices_out=iout)
+ result = scheme(rd.data[0], 1, result)
+
+In that case:
+
+* result.shape = (17,17)
+* scheme.indices = [slice(4,21), slice(4,21)]
+* scheme.output_indices = iout
+
+This option is usefull when input tab is a work array, with a size different
+from the topology resolution, and the result a field defined on the whole grid.
+
diff --git a/docs/sphinx/users_guide/numerics.rst b/docs/sphinx/users_guide/numerics.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3f420a06f998b8dcc8e709db9f5585ed74a5072a
--- /dev/null
+++ b/docs/sphinx/users_guide/numerics.rst
@@ -0,0 +1,15 @@
+ .. _numerics:
+
+Numerical operations
+====================
+
+This module handles all the 'low-level' operations mostly applied on numpy arrays.
+
+.. toctree::
+ :maxdepth: 2
+
+ finite_differences
+ differential_operations
+ remeshing_and_interpolation
+ odesolvers
+
diff --git a/hysop/default_methods.py b/hysop/default_methods.py
index 086c396101650428a7ff3b7d7858f97fe7eba9d4..425cafd5f9e8e7f7886dead7cf75c825aa6c5fd8 100644
--- a/hysop/default_methods.py
+++ b/hysop/default_methods.py
@@ -1,6 +1,4 @@
-"""
-@file default_methods.py
-Default parameter values for methods in operators.
+"""Default parameters values for numerical methods in operators.
"""
from hysop.methods_keys import TimeIntegrator, Interpolation, GhostUpdate,\
Remesh, Support, Splitting, MultiScale, Formulation, SpaceDiscretisation, \
@@ -8,6 +6,7 @@ from hysop.methods_keys import TimeIntegrator, Interpolation, GhostUpdate,\
from hysop.constants import HYSOP_REAL
from hysop.numerics.integrators.runge_kutta2 import RK2
from hysop.numerics.integrators.runge_kutta3 import RK3
+from hysop.numerics.finite_differences import FDC4, FDC2
from hysop.numerics.interpolation import Linear
from hysop.numerics.remeshing import L2_1, Linear as Rmsh_Linear
#from hysop.operator.discrete.stretching import Conservative
@@ -17,18 +16,18 @@ ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear,
Remesh: L2_1, Support: '', Splitting: 'o2', MultiScale: L2_1,
Precision: HYSOP_REAL}
-from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+from hysop.numerics.finite_differences import FDC4, FDC2
-DIFFERENTIAL = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
+DIFFERENTIAL = {SpaceDiscretisation: FDC4, GhostUpdate: True}
-ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4,
+ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FDC4,
dtCrit: 'vort'}
-BAROCLINIC = {SpaceDiscretisation: FD_C_4}
+BAROCLINIC = {SpaceDiscretisation: FDC4}
-MULTIPHASEBAROCLINIC = {SpaceDiscretisation: FD_C_4}
+MULTIPHASEBAROCLINIC = {SpaceDiscretisation: FDC4}
-MULTIPHASEGRADP = {SpaceDiscretisation: FD_C_4}
+MULTIPHASEGRADP = {SpaceDiscretisation: FDC4}
DIFFUSION = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
@@ -36,8 +35,8 @@ POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True,
Formulation: 'velocity'}
STRETCHING = {TimeIntegrator: RK3, Formulation: "Conservative",
- SpaceDiscretisation: FD_C_4}
+ SpaceDiscretisation: FDC4}
-FORCES = {SpaceDiscretisation: FD_C_2}
+FORCES = {SpaceDiscretisation: FDC2}
MULTIRESOLUTION_FILER = {Remesh: Rmsh_Linear, }
diff --git a/hysop/gpu/gpu_multiphase_baroclinic_rhs.py b/hysop/gpu/gpu_multiphase_baroclinic_rhs.py
index 25500f61db308ee81aa5401590a6cebb7905e6db..32ce524cb79d6d3d97db0945f2bc04d5495a1ed4 100644
--- a/hysop/gpu/gpu_multiphase_baroclinic_rhs.py
+++ b/hysop/gpu/gpu_multiphase_baroclinic_rhs.py
@@ -15,7 +15,7 @@ from hysop.gpu.gpu_discrete import GPUDiscreteField
from hysop.tools.profiler import FProfiler
from hysop.mpi import Wtime
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
class BaroclinicRHS(DiscreteOperator, GPUOperator):
@@ -103,7 +103,7 @@ class BaroclinicRHS(DiscreteOperator, GPUOperator):
if 0 in self._cutdir_list:
raise ValueError("Not yet implemented with comm in X dir")
for d in self._cutdir_list:
- if self.method[SpaceDiscretisation] == FD_C_4:
+ if self.method[SpaceDiscretisation] == FDC4:
gh = 2
else:
gh = 1
diff --git a/hysop/gpu/tests/test_multiphase_baroclinic.py b/hysop/gpu/tests/test_multiphase_baroclinic.py
index e4d45594ae55d5e2b2eeda042fb0f972f2986d72..5c3aa47f847fb130ee077a32f1ee62a7efb58941 100755
--- a/hysop/gpu/tests/test_multiphase_baroclinic.py
+++ b/hysop/gpu/tests/test_multiphase_baroclinic.py
@@ -5,7 +5,7 @@ from hysop.tools.parameters import Discretization, MPIParams
from hysop.domain.box import Box
from hysop.fields.continuous import Field
from hysop.operator.multiphase_baroclinic_rhs import MultiphaseBaroclinicRHS
-from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+from hysop.numerics.finite_differences import FDC4, FDC2
import hysop.tools.numpywrappers as npw
import numpy as np
from hysop.methods_keys import Support, SpaceDiscretisation, ExtraArgs
@@ -91,7 +91,7 @@ def call_operator(func, grad_func, vfunc):
gradp: d_coarse,
rho: d_fine},
method={Support: 'gpu',
- SpaceDiscretisation: FD_C_4,
+ SpaceDiscretisation: FDC4,
ExtraArgs: {'density_func': 'x', }})
op.discretize()
op.setup()
diff --git a/hysop/methods.py b/hysop/methods.py
index 88858e611ecf5aa16661735798696685021c7006..3b971ad3f8ebd1487daa6a35c7baf04c49a8b812 100644
--- a/hysop/methods.py
+++ b/hysop/methods.py
@@ -1,19 +1,24 @@
-"""
-@file methods.py
-A list of numerical methods available in HySoP that may be used to set methods
-in operators.
+"""Some shortcuts to numerical methods available in HySoP that may be used
+to set methods attribute in operators.
+
Usage:
+
method = {key: value, ...}
+
Keys must be one of the constants given in methods_keys.py.
+
Value is usually a class name
and sometimes a string. See details in each operator.
-Example, the stretching case :
-method = {TimeIntegrator: RK3, Formulation: Conservative,
- SpaceDiscretisation: FD_C_4}
+For example, consider the stretching case where time-integration
+is done with Runge-Kutta 3, order 4 finite-difference scheme is used
+for space discretisation and the default formulation is conservative::
+
+ method = {TimeIntegrator: RK3, Formulation: Conservative,
+ SpaceDiscretisation: FDC4}
Note FP: to avoid cycling, this file must never be imported
-inside a HySoP module. It's only a review of all the methods
+inside an HySoP module. It's only a review of all the methods
that can be imported by final user.
"""
@@ -50,9 +55,9 @@ Linear = interpolation.Linear
# Finite differences
import hysop.numerics.finite_differences as fd
-FD_C_4 = fd.FD_C_4
-FD_C_2 = fd.FD_C_2
-FD2_C_2 = fd.FD2_C_2
+FDC4 = fd.FDC4
+FDC2 = fd.FDC2
+FD2C2 = fd.FD2C2
# Stretching formulations
import hysop.operator.discrete.stretching as strd
diff --git a/hysop/methods_keys.py b/hysop/methods_keys.py
index a4505dee4b5301303b6c6941c356bd9583a6cc11..ec4601ed3d02e5d1432c25dd7f519beeee090b69 100644
--- a/hysop/methods_keys.py
+++ b/hysop/methods_keys.py
@@ -9,7 +9,7 @@ and sometimes a string. See details in each operator.
Example, the stretching case :
method = {TimeIntegrator: RK3, Formulation: Conservative,
- SpaceDiscretisation: FD_C_4}
+ SpaceDiscretisation: FDC4}
"""
diff --git a/hysop/numerics/__init__.py b/hysop/numerics/__init__.py
index c991ed075ff1f03868bc14e421f731f1f9f6f649..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/hysop/numerics/__init__.py
+++ b/hysop/numerics/__init__.py
@@ -1,17 +0,0 @@
-## @package hysop.numerics
-#
-# \todo write a proper doc for numerical methods
-#
-# All functions in this package are supposed to work with numpy arrays,
-# not with hysop fields.
-#
-#
-#
-#
-#
-# \section perfandmem About performances and optimisations
-#
-# \subsection Finite Differences
-# Two ways :
-# - fd.compute
-# - fd.compute_and_add
diff --git a/hysop/numerics/differential_operations.py b/hysop/numerics/differential_operations.py
index 3f78337a60ee3fdc12309aa4c4406d0f26fc759c..67559774c23e9d897297f697a90c0c39c0466159 100755
--- a/hysop/numerics/differential_operations.py
+++ b/hysop/numerics/differential_operations.py
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
"""Library of functions used to perform classical vector calculus
-(diff operations like grad, curl ...)
+(diff operations like grad, curl ...) based on finite differences schemes.
* :class:`~hysop.numerics.differential_operations.Curl`,
* :class:`~hysop.numerics.differential_operations.DivRhoV`,
@@ -28,12 +28,16 @@ For example :
"""
from hysop.constants import debug, XDIR, YDIR, ZDIR
from abc import ABCMeta
-from hysop.numerics.finite_differences import FD_C_4, FD_C_2, FD2_C_2
+from hysop.numerics.finite_differences import FDC4, FDC2, FD2C2
import numpy as np
import hysop.tools.numpywrappers as npw
+from hysop.tools.misc import WorkSpaceTools, Utils
class DifferentialOperation(object):
+ """Abstract base class for all operations
+ based on finite differences.
+ """
__metaclass__ = ABCMeta
_authorized_methods = []
@@ -140,7 +144,8 @@ class Curl(DifferentialOperation):
"""
Computes nabla X V, V being a vector field.
"""
- _authorized_methods = [FD_C_2, FD_C_4]
+ _authorized_methods = [FDC2, FDC4]
+ _lwork = 1
def __init__(self, **kwds):
"""Curl of a vector field
@@ -217,18 +222,10 @@ class DivRhoV(DifferentialOperation):
V a vector field.
Works for any dimension.
- Methods : FD_C_4
- 1 - default : based on fd.compute.
- 2 - opt : based on fd.compute_and_add.
- To use this one, call self._central4_caa.
- shape is a tuple (like np.ndarray.shape), e.g. (nx,ny,nz).
-
- Note FP:
- - 1 seems to be faster than 2 for large systems (>200*3) and
- less memory consumming. For smaller systems, well, it depends ...
+ Methods : FDC4
"""
- _authorized_methods = [FD_C_4]
+ _authorized_methods = [FDC4]
def __init__(self, **kwds):
"""Divergence of rho V, rho a scalar field, V a vector field
@@ -240,7 +237,7 @@ class DivRhoV(DifferentialOperation):
result = div(rhoV), which is the default.
**kwds : parameters for base class
- Default method : FD_C_4
+ Default method : FDC4
"""
super(DivRhoV, self).__init__(**kwds)
#self.fcall = self._central4_caa
@@ -329,7 +326,7 @@ class DivWV(DifferentialOperation):
"""
- _authorized_methods = [FD_C_4]
+ _authorized_methods = [FDC4]
def __init__(self, **kwds):
"""Divergence of (W.Vx, W.Vy, W.Vz), W, V two vector fields
@@ -341,7 +338,7 @@ class DivWV(DifferentialOperation):
result = div(rhoV), which is the default.
**kwds : parameters for base class
- Default method : FD_C_4
+ Default method : FDC4
"""
super(DivWV, self).__init__(**kwds)
self.fcall = self._central4
@@ -385,7 +382,8 @@ class DivWV(DifferentialOperation):
class Laplacian(DifferentialOperation):
"""Computes the laplacian of a field.
"""
- _authorized_methods = [FD2_C_2]
+ _authorized_methods = [FD2C2]
+ _lwork = 1
@staticmethod
def get_work_length():
@@ -415,7 +413,7 @@ class Laplacian(DifferentialOperation):
class GradS(DifferentialOperation):
"""Gradient of a scalar field
"""
- _authorized_methods = [FD_C_4, FD_C_2]
+ _authorized_methods = [FDC4, FDC2]
def __call__(self, scal, result):
"""Apply gradient, with central finite difference scheme.
@@ -438,7 +436,7 @@ class GradS(DifferentialOperation):
class GradV(DifferentialOperation):
"""Gradient of a vector field
"""
- _authorized_methods = [FD_C_4, FD_C_2]
+ _authorized_methods = [FDC4, FDC2]
def __init__(self, **kwds):
super(GradV, self).__init__(**kwds)
@@ -468,7 +466,8 @@ class GradVxW(DifferentialOperation):
"""Computes [nabla(V)][W] with
V and W some vector fields.
"""
- _authorized_methods = [FD_C_4, FD_C_2]
+ _authorized_methods = [FDC4, FDC2]
+ _lwork = 2
@staticmethod
def get_work_length():
@@ -517,7 +516,8 @@ class GradVxW(DifferentialOperation):
class DivAdvection(DifferentialOperation):
""" Computes -nabla .(V . nabla V) with V a vector field.
"""
- _authorized_methods = [FD_C_4, FD_C_2]
+ _authorized_methods = [FDC4, FDC2]
+ _lwork = 3
@staticmethod
def get_work_length():
diff --git a/hysop/numerics/finite_differences.py b/hysop/numerics/finite_differences.py
index cecf245794feca8d17fd730f4395b4a4bf5cb341..dbf51487dd75011a5000718a51668e13639fe095 100644
--- a/hysop/numerics/finite_differences.py
+++ b/hysop/numerics/finite_differences.py
@@ -1,4 +1,49 @@
"""Finite difference schemes
+
+.. currentmodule hysop.numerics.finite_differences
+
+* :class:`~FDC2` : first derivative, 2nd order centered scheme
+* :class:`~FD2C2`: second derivative, 2nd order centered scheme
+* :class:`~FDC4`: first derivative, 4th order centered scheme
+* :class:`~FiniteDifference`, abstract base class.
+
+
+Usage:
+
+.. code::
+
+ # indices and space step from a predefined mesh
+ indices = topo.mesh.compute_index
+ step = topo.mesh.space_step
+ # init scheme
+ scheme = FDC4(step, indices)
+ # result = some array of the same shape as field.
+ # compute dfield/d direction in result
+ result = scheme(field, direction, result)
+ # or (with more assert and checks on inputs)
+ scheme.compute(field, direction, result)
+ # or result += dfield/d direction
+ scheme.compute_and_add(field, direction, result)
+
+
+Note
+----
+* fd are computed only on points described by indices :
+ result[ind] = diff(field[ind]),
+ which means that result and field must have the same shape.
+* You can also used a 'reduced' output result to minimize memory.
+ Use 'reduce_output_shape=True' in scheme init. In that case,
+ result must be of a shape corresponding to indices 'shape'.
+
+ Example:
+
+.. code::
+
+ scheme = FDC4(step, indices, reduce_output_shape=True)
+ # if indices is slice(4,10) in each direction
+ result = npw.zeros((6,6,6))
+ ...
+
"""
from abc import ABCMeta, abstractmethod
from hysop.constants import debug
@@ -13,7 +58,7 @@ class FiniteDifference(object):
Usage :
>> step = topo.mesh_space_step
- >> scheme = FD_C_4(step, topo.mesh.compute_index)
+ >> scheme = FDC4(step, topo.mesh.compute_index)
For a given numpy array (obviously discretized on the topo
used to compute indices), to compute result as the derivative
@@ -39,7 +84,7 @@ class FiniteDifference(object):
def __new__(cls, *args, **kw):
return object.__new__(cls, *args, **kw)
- def __init__(self, step, indices, reduce_output_shape=False):
+ def __init__(self, step, indices, indices_out=None):
"""
Parameters
@@ -48,26 +93,20 @@ class FiniteDifference(object):
resolution of the mesh
indices : list of slices
Represents the local mesh on which finite-differences
- will be applied, like compute_index in :class:`~hysop.domain.mesh.Mesh`.
- reduce_output_shape : boolean, optional
- True to return the result in a reduced array. See notes below.
+ will be applied, like compute_index
+ in :class:`~hysop.domain.mesh.Mesh`.
+ indices_out : list of slices, optional
+ set which indices of result will be modified. By default,
+ indices_out = indices
Attributes
----------
- minimal_ghost_layer : int
+ ghosts_layer_size : int, static
minimal number of points required inside the ghost layer
for this scheme.
- Notes
- -----
- Two ways to compute result = scheme(input)
- * 1 - either input and result are arrays of the same shape, and
- result[indices] = scheme(input[indices]) will be computed
- * 2 - or result is a smaller array than input, and
- result[...] = scheme(input(indices))
- To use case 2, set reduce_output_shape = True.
"""
-
+
step = np.asarray(step)
# dim of the field on which scheme will be applied
# (i.e dim of the domain)
@@ -78,18 +117,36 @@ class FiniteDifference(object):
self._a2 = None
self._coeff = None
# List of slices representing the mesh on which fd scheme is applied
- self._indices = indices
- if not reduce_output_shape:
- self.result_shape = None
- self.output_indices = self._indices
- else:
- self.result_shape = tuple([indices[i].stop - indices[i].start
- for i in xrange(len(indices))])
- self.output_indices = [slice(0, self.result_shape[i])
+ self.indices = indices
+ if indices_out is None or indices_out is False:
+ # This is the default case. It means that
+ # result shape is unknown. It will
+ # be the same as input tab shape during call
+ # and input/output indices are the same:
+ # we compute result[indices] = scheme(tab[indices])
+ self.output_indices = self.indices
+ elif isinstance(indices_out, list):
+ # In that case, result shape remains unknown but
+ # the position where result is written is not the
+ # same as input indices.
+ # We compute:
+ # result[indices_out] = scheme(tab[indices])
+ # Warning : result shape/indices_out compliance
+ # is not checked!
+ self.output_indices = indices_out
+
+ elif indices_out is True:
+ # Here, output shape is fixed to the minimal required
+ # shape, computed with indices.
+ # So, we compute:
+ # result[...] = scheme(tab[indices])
+ result_shape = tuple([indices[i].stop - indices[i].start
+ for i in xrange(len(indices))])
+ self.output_indices = [slice(0, result_shape[i])
for i in xrange(len(indices))]
+ else:
+ assert False
- # Minimal size of the ghost layer.
- self.minimal_ghost_layer = 1
self._step = step
self._compute_indices(step)
@@ -118,13 +175,30 @@ class FiniteDifference(object):
"""
assert result is not tab
assert result.__class__ is np.ndarray
- assert result.shape == self.result_shape or result.shape == tab.shape
assert tab.__class__ is np.ndarray
- self.__call__(tab, cdir, result)
+ result = self.__call__(tab, cdir, result)
+ return result
@abstractmethod
- def compute_and_add(self, tab, cdir, result):
- """Apply FD scheme and add the result inplace.
+ def compute_and_add(self, tab, cdir, result, work):
+ """Apply FD scheme and add the abs of the result inplace.
+
+ Parameters
+ ----------
+ tab : numpy array
+ input field
+ cdir : int
+ direction of differentiation
+ result : numpy array
+ in/out, abs(derivative of tab) + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+ """
+
+ @abstractmethod
+ def compute_and_add_abs(self, tab, cdir, result, work):
+ """Apply FD scheme and add np.abs of the derivative inplace.
Parameters
----------
@@ -134,10 +208,13 @@ class FiniteDifference(object):
direction of differentiation
result : numpy array
in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
"""
-class FD_C_2(FiniteDifference):
+class FDC2(FiniteDifference):
"""
1st derivative, centered scheme, 2nd order.
@@ -147,19 +224,18 @@ class FD_C_2(FiniteDifference):
def _compute_indices(self, step):
- self.minimal_ghost_layer = 1
self._coeff = npw.asarray(1. / (2. * step))
self._m1 = []
self._a1 = []
for dim in xrange(self._dim):
- self._m1.append(list(self._indices))
- self._m1[dim][dim] = slice(self._indices[dim].start - 1,
- self._indices[dim].stop - 1,
- self._indices[dim].step)
- self._a1.append(list(self._indices))
- self._a1[dim][dim] = slice(self._indices[dim].start + 1,
- self._indices[dim].stop + 1,
- self._indices[dim].step)
+ self._m1.append(list(self.indices))
+ self._m1[dim][dim] = slice(self.indices[dim].start - 1,
+ self.indices[dim].stop - 1,
+ self.indices[dim].step)
+ self._a1.append(list(self.indices))
+ self._a1[dim][dim] = slice(self.indices[dim].start + 1,
+ self.indices[dim].stop + 1,
+ self.indices[dim].step)
def __call__(self, tab, cdir, result):
result[self.output_indices] = tab[self._a1[cdir]]
@@ -167,8 +243,8 @@ class FD_C_2(FiniteDifference):
result[self.output_indices] *= self._coeff[cdir]
return result
- def compute_and_add(self, tab, cdir, result):
- """Apply FD scheme and add the result inplace.
+ def compute_and_add(self, tab, cdir, result, work):
+ """Apply FD scheme and add np.abs of the derivative inplace.
Parameters
----------
@@ -178,14 +254,55 @@ class FD_C_2(FiniteDifference):
direction of differentiation
result : numpy array
in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
"""
assert result.__class__ is np.ndarray
assert tab.__class__ is np.ndarray
- result[self.output_indices] += \
- self._coeff[cdir] * (tab[self._a1[cdir]] - tab[self._m1[cdir]])
+ assert work.size == result.size
+ wk = work.reshape(result.shape)
+ np.subtract(tab[self._a1[cdir]], tab[self._m1[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+ return result
+ def compute_and_add_abs(self, tab, cdir, result, work):
+ """Apply FD scheme and add np.abs of the derivative inplace.
-class FD2_C_2(FiniteDifference):
+ Parameters
+ ----------
+ tab : numpy array
+ input field
+ cdir : int
+ direction of differentiation
+ result : numpy array
+ in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
+ """
+ assert result.__class__ is np.ndarray
+ assert tab.__class__ is np.ndarray
+ assert work.__class__ is np.ndarray
+ assert work is not result
+ wk = work.reshape(result.shape)
+ np.subtract(tab[self._a1[cdir]], tab[self._m1[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.abs(wk[self.output_indices], wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+ return result
+
+
+class FD2C2(FiniteDifference):
"""
Second derivative, centered scheme, 2nd order.
"""
@@ -194,29 +311,28 @@ class FD2_C_2(FiniteDifference):
def _compute_indices(self, step):
- self.minimal_ghost_layer = 1
self._m1 = []
self._a1 = []
self._coeff = npw.asarray(1. / (step * step))
for dim in xrange(self._dim):
- self._m1.append(list(self._indices))
- self._m1[dim][dim] = slice(self._indices[dim].start - 1,
- self._indices[dim].stop - 1,
- self._indices[dim].step)
- self._a1.append(list(self._indices))
- self._a1[dim][dim] = slice(self._indices[dim].start + 1,
- self._indices[dim].stop + 1,
- self._indices[dim].step)
+ self._m1.append(list(self.indices))
+ self._m1[dim][dim] = slice(self.indices[dim].start - 1,
+ self.indices[dim].stop - 1,
+ self.indices[dim].step)
+ self._a1.append(list(self.indices))
+ self._a1[dim][dim] = slice(self.indices[dim].start + 1,
+ self.indices[dim].stop + 1,
+ self.indices[dim].step)
def __call__(self, tab, cdir, result):
- result[self.output_indices] = tab[self._indices]
+ result[self.output_indices] = tab[self.indices]
result[self.output_indices] *= -2
result[self.output_indices] += tab[self._a1[cdir]]
result[self.output_indices] += tab[self._m1[cdir]]
result[self.output_indices] *= self._coeff[cdir]
return result
- def compute_and_add(self, tab, cdir, result):
+ def compute_and_add(self, tab, cdir, result, work):
"""Apply FD scheme and add the result inplace.
Parameters
@@ -227,15 +343,59 @@ class FD2_C_2(FiniteDifference):
direction of differentiation
result : numpy array
in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
+ """
+ assert result.__class__ is np.ndarray
+ assert tab.__class__ is np.ndarray
+ wk = work.reshape(result.shape)
+ np.multiply(tab[self.indices], -2., wk[self.output_indices])
+ np.add(tab[self._a1[cdir]], wk[self.output_indices],
+ wk[self.output_indices])
+ np.add(tab[self._m1[cdir]], wk[self.output_indices],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+ return result
+
+ def compute_and_add_abs(self, tab, cdir, result, work):
+ """Apply FD scheme and add np.abs of the derivative inplace.
+
+ Parameters
+ ----------
+ tab : numpy array
+ input field
+ cdir : int
+ direction of differentiation
+ result : numpy array
+ in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
"""
assert result.__class__ is np.ndarray
assert tab.__class__ is np.ndarray
- result[self.output_indices] += self._coeff[cdir] * (
- -2 * tab[self._indices] + tab[self._m1[cdir]] + tab[self._a1[cdir]]
- )
+ wk = work.reshape(result.shape)
+ np.multiply(tab[self.indices], -2., wk[self.output_indices])
+ np.add(tab[self._a1[cdir]], wk[self.output_indices],
+ wk[self.output_indices])
+ np.add(tab[self._m1[cdir]], wk[self.output_indices],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.abs(wk[self.output_indices], wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+
+ return result
-class FD_C_4(FiniteDifference):
+class FDC4(FiniteDifference):
"""
1st derivative, centered scheme, 4th order.
"""
@@ -244,7 +404,6 @@ class FD_C_4(FiniteDifference):
def _compute_indices(self, step):
- self.minimal_ghost_layer = 2
self._m1 = []
self._m2 = []
self._a1 = []
@@ -252,22 +411,22 @@ class FD_C_4(FiniteDifference):
# FD scheme coefficients
self._coeff = npw.asarray(1. / (12. * step))
for dim in xrange(self._dim):
- self._m1.append(list(self._indices))
- self._m1[dim][dim] = slice(self._indices[dim].start - 1,
- self._indices[dim].stop - 1,
- self._indices[dim].step)
- self._m2.append(list(self._indices))
- self._m2[dim][dim] = slice(self._indices[dim].start - 2,
- self._indices[dim].stop - 2,
- self._indices[dim].step)
- self._a1.append(list(self._indices))
- self._a1[dim][dim] = slice(self._indices[dim].start + 1,
- self._indices[dim].stop + 1,
- self._indices[dim].step)
- self._a2.append(list(self._indices))
- self._a2[dim][dim] = slice(self._indices[dim].start + 2,
- self._indices[dim].stop + 2,
- self._indices[dim].step)
+ self._m1.append(list(self.indices))
+ self._m1[dim][dim] = slice(self.indices[dim].start - 1,
+ self.indices[dim].stop - 1,
+ self.indices[dim].step)
+ self._m2.append(list(self.indices))
+ self._m2[dim][dim] = slice(self.indices[dim].start - 2,
+ self.indices[dim].stop - 2,
+ self.indices[dim].step)
+ self._a1.append(list(self.indices))
+ self._a1[dim][dim] = slice(self.indices[dim].start + 1,
+ self.indices[dim].stop + 1,
+ self.indices[dim].step)
+ self._a2.append(list(self.indices))
+ self._a2[dim][dim] = slice(self.indices[dim].start + 2,
+ self.indices[dim].stop + 2,
+ self.indices[dim].step)
def __call__(self, tab, cdir, result):
result[self.output_indices] = tab[self._a1[cdir]]
@@ -278,7 +437,7 @@ class FD_C_4(FiniteDifference):
result[self.output_indices] *= self._coeff[cdir]
return result
- def compute_and_add(self, tab, cdir, result):
+ def compute_and_add(self, tab, cdir, result, work):
"""Apply FD scheme and add the result inplace.
Parameters
@@ -289,9 +448,62 @@ class FD_C_4(FiniteDifference):
direction of differentiation
result : numpy array
in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
+ """
+ assert result.__class__ is np.ndarray
+ assert tab.__class__ is np.ndarray
+ assert work is not result
+ assert work.size == result.size
+ wk = work.reshape(result.shape)
+ np.subtract(tab[self._a1[cdir]], tab[self._m1[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], 8.,
+ wk[self.output_indices])
+ np.add(wk[self.output_indices], tab[self._m2[cdir]],
+ wk[self.output_indices])
+ np.subtract(wk[self.output_indices], tab[self._a2[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+ return result
+
+ def compute_and_add_abs(self, tab, cdir, result, work):
+ """Apply FD scheme and add np.abs of the derivative inplace.
+
+ Parameters
+ ----------
+ tab : numpy array
+ input field
+ cdir : int
+ direction of differentiation
+ result : numpy array
+ in/out, derivative of tab + result
+ work : numpy array
+ internal workspace. Must be 1D and of the same
+ size as result.
+
"""
assert result.__class__ is np.ndarray
assert tab.__class__ is np.ndarray
- result[self.output_indices] += self._coeff[cdir] * (
- 8 * (tab[self._a1[cdir]] - tab[self._m1[cdir]]) +
- tab[self._m2[cdir]] - tab[self._a2[cdir]])
+ assert work is not result
+ assert work.size == result.size
+ wk = work.reshape(result.shape)
+ np.subtract(tab[self._a1[cdir]], tab[self._m1[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], 8.,
+ wk[self.output_indices])
+ np.add(wk[self.output_indices], tab[self._m2[cdir]],
+ wk[self.output_indices])
+ np.subtract(wk[self.output_indices], tab[self._a2[cdir]],
+ wk[self.output_indices])
+ np.multiply(wk[self.output_indices], self._coeff[cdir],
+ wk[self.output_indices])
+ np.abs(wk[self.output_indices], wk[self.output_indices])
+ np.add(wk[self.output_indices], result[self.output_indices],
+ result[self.output_indices])
+ return result
diff --git a/hysop/numerics/tests/test_finite_differences.py b/hysop/numerics/tests/test_finite_differences.py
new file mode 100644
index 0000000000000000000000000000000000000000..b098eb7b03b12cb43fd210fe543a9e8795af8a60
--- /dev/null
+++ b/hysop/numerics/tests/test_finite_differences.py
@@ -0,0 +1,199 @@
+"""FD schemes tests
+"""
+from hysop.numerics.finite_differences import FDC2, FDC4, FD2C2
+import hysop.tools.numpywrappers as npw
+from hysop import Box, Discretization, Field
+import numpy as np
+from hysop.domain.subsets import SubBox
+
+
+Nx = Ny = Nz = 128
+g = 2
+d2d = Discretization([Nx, Ny], [g, g])
+d3d = Discretization([Nx, Ny, Nz], [g, g, g])
+
+
+def f2d(res, x, y, t):
+ res[0][...] = np.cos(x) * y ** 2 + x * y ** 3
+ return res
+
+
+def ref_f2d(res, x, y, t):
+ res[0][...] = -np.sin(x) * y ** 2 + y ** 3
+ return res
+
+
+def ref2_f2d(res, x, y, t):
+ res[0][...] = 2 * np.cos(x) + 6 * x * y
+ return res
+
+
+def f3d(res, x, y, z, t):
+ res[0][...] = np.cos(x) * y ** 2 + x * y ** 3
+ return res
+
+
+def ref_f3d(res, x, y, z, t):
+ res[0][...] = -np.sin(x) * y ** 2 + y ** 3
+ return res
+
+
+def ref2_f3d(res, x, y, z, t):
+ res[0][...] = 2 * np.cos(x) + 6 * x * y
+ return res
+
+
+lengthes = [1., 1., 1.]
+
+
+def init(discr, formulas):
+ dimension = len(discr.resolution)
+ box = Box(length=lengthes[:dimension])
+ f1 = Field(name='input', domain=box, formula=formulas[0])
+ ref = Field(name='ref', domain=box, formula=formulas[1])
+ topo = box.create_topology(discr)
+ f1d = f1.discretize(topo).data[0]
+ refd = ref.discretize(topo).data[0]
+ f1.initialize(topo=topo)
+ ref.initialize(topo=topo)
+ ref2 = Field(name='ref', domain=box, formula=formulas[2])
+ ref2d = ref2.discretize(topo).data[0]
+ ref2.initialize(topo=topo)
+ return topo, f1d, refd, ref2d
+
+
+def run_schemes(schemes, f1d, result, ind, hh, iout=None):
+
+ work = npw.zeros(result.size)
+ if iout is None:
+ iout = ind
+ for sc in schemes:
+ ref = schemes[sc][0]
+ order = schemes[sc][2]
+ direction = schemes[sc][1]
+ result = sc(f1d, direction, result)
+ assert np.allclose(result[iout], ref[ind], atol=hh[0] ** order)
+ result[...] = 1.
+ result = sc.compute_and_add(f1d, direction, result, work)
+ assert np.allclose(result[iout], 1. + ref[ind], atol=hh[0] ** order)
+ result[...] = 1.
+ result = sc.compute_and_add_abs(f1d, direction, result, work)
+ assert np.allclose(result[iout], 1. + np.abs(ref[ind]),
+ atol=hh[0] ** order)
+
+
+def check_fd_schemes(discr, formulas):
+ topo, f1d, refd, ref2d = init(discr, formulas)
+ ind = topo.mesh.compute_index
+ hh = topo.mesh.space_step
+ sc = {
+ FDC2(hh, ind): [refd, 0, 2],
+ FDC4(hh, ind): [refd, 0, 4],
+ FD2C2(hh, ind): [ref2d, 1, 2]}
+ shape = topo.mesh.resolution
+ result = npw.zeros(shape)
+ run_schemes(sc, f1d, result, ind, hh)
+
+
+def check_fd_schemes_reduce_input(discr, formulas):
+ topo, f1d, refd, ref2d = init(discr, formulas)
+ sl = topo.domain.length * 0.5
+ orig = topo.domain.origin + 0.1 * topo.domain.length
+ subbox = SubBox(parent=topo.domain, origin=orig, length=sl)
+ ind = subbox.discretize(topo)[0]
+ hh = topo.mesh.space_step
+ sc = {
+ FDC2(hh, ind): [refd, 0, 2],
+ FDC4(hh, ind): [refd, 0, 4],
+ FD2C2(hh, ind): [ref2d, 1, 2]}
+ shape = topo.mesh.resolution
+ result = npw.zeros(shape)
+ run_schemes(sc, f1d, result, ind, hh)
+
+
+def check_fd_schemes_reduce_output(discr, formulas):
+ topo, f1d, refd, ref2d = init(discr, formulas)
+ ind = topo.mesh.compute_index
+ hh = topo.mesh.space_step
+ sc = {
+ FDC2(hh, ind, indices_out=True): [refd, 0, 2],
+ FDC4(hh, ind, indices_out=True): [refd, 0, 4],
+ FD2C2(hh, ind, indices_out=True): [ref2d, 1, 2]
+ }
+ shape = np.asarray(topo.mesh.resolution).copy()
+ shape -= 2 * g
+ shape = tuple(shape)
+ result = npw.zeros(shape)
+ iout = [slice(0, shape[i]) for i in xrange(len(shape))]
+ run_schemes(sc, f1d, result, ind, hh, iout)
+
+
+def check_fd_schemes_reduce_all(discr, formulas):
+ topo, f1d, refd, ref2d = init(discr, formulas)
+ hh = topo.mesh.space_step
+ sl = topo.domain.length * 0.5
+ orig = topo.domain.origin + 0.1 * topo.domain.length
+ subbox = SubBox(parent=topo.domain, origin=orig, length=sl)
+ ind = subbox.discretize(topo)[0]
+ sc = {
+ FDC2(hh, ind, indices_out=True): [refd, 0, 2],
+ FDC4(hh, ind, indices_out=True): [refd, 0, 4],
+ FD2C2(hh, ind, indices_out=True): [ref2d, 1, 2]
+ }
+ shape = subbox.mesh[topo].resolution
+ result = npw.zeros(shape)
+ iout = [slice(0, shape[i]) for i in xrange(len(shape))]
+ run_schemes(sc, f1d, result, ind, hh, iout)
+
+
+def check_fd_schemes_setiout(discr, formulas):
+ topo, f1d, refd, ref2d = init(discr, formulas)
+ hh = topo.mesh.space_step
+ sl = topo.domain.length * 0.5
+ orig = topo.domain.origin + 0.1 * topo.domain.length
+ subbox = SubBox(parent=topo.domain, origin=orig, length=sl)
+ ind = subbox.discretize(topo)[0]
+ shape = np.asarray(subbox.mesh[topo].resolution).copy()
+ shape += 6
+ shape = tuple(shape)
+ result = npw.zeros(shape)
+ iout = [slice(3, shape[i] - 3) for i in xrange(len(shape))]
+ sc = {
+ FDC2(hh, ind, indices_out=iout): [refd, 0, 2],
+ FDC4(hh, ind, indices_out=iout): [refd, 0, 4],
+ FD2C2(hh, ind, indices_out=iout): [ref2d, 1, 2]
+ }
+ run_schemes(sc, f1d, result, ind, hh, iout)
+
+
+def test_fd_2d():
+ check_fd_schemes(d2d, [f2d, ref_f2d, ref2_f2d])
+
+
+def test_fd_3d():
+ check_fd_schemes(d3d, [f3d, ref_f3d, ref2_f3d])
+
+
+def test_fd_3d_reduce_input():
+ check_fd_schemes_reduce_input(d3d, [f3d, ref_f3d, ref2_f3d])
+
+
+def test_fd_3d_reduce_output():
+ check_fd_schemes_reduce_output(d3d, [f3d, ref_f3d, ref2_f3d])
+
+
+def test_fd_3d_reduce_all():
+ check_fd_schemes_reduce_all(d3d, [f3d, ref_f3d, ref2_f3d])
+
+
+def test_fd_3d_setiout():
+ check_fd_schemes_setiout(d3d, [f3d, ref_f3d, ref2_f3d])
+
+
+if __name__ == '__main__':
+ test_fd_2d()
+ test_fd_3d()
+ test_fd_3d_reduce_input()
+ test_fd_3d_reduce_output()
+ test_fd_3d_reduce_all()
+ test_fd_3d_setiout()
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 93ba310d003dc1be383e7cd5b9e37ec4e527d807..b298ea7200e214064203ce4cd26b060907e17924 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -8,7 +8,7 @@ Definition of the adaptative time step according to the flow fields.
from hysop.constants import debug
from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation,\
dtCrit
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.operator.discrete.adapt_timestep import AdaptTimeStep_D
from hysop.operator.continuous import opsetup
from hysop.operator.computational import Computational
@@ -102,7 +102,7 @@ class AdaptTimeStep(Computational):
return res
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
nbGhosts = 2
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/baroclinic.py b/hysop/operator/baroclinic.py
index 7609042784860623d860ddb7d85d1a0208bd9312..7937b55a558377a6420c858b106630bfb1aed6c5 100644
--- a/hysop/operator/baroclinic.py
+++ b/hysop/operator/baroclinic.py
@@ -7,7 +7,7 @@ MultiPhase Rot Grad P
from hysop.operator.computational import Computational
from hysop.operator.discrete.baroclinic import Baroclinic as BD
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.constants import debug
import hysop.default_methods as default
from hysop.operator.continuous import opsetup
@@ -51,7 +51,7 @@ class Baroclinic(Computational):
assert SpaceDiscretisation in self.method.keys()
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
nbGhosts = 2
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/baroclinic_from_rhs.py b/hysop/operator/baroclinic_from_rhs.py
index 2678d8aca0b9250faa449de187e5b520dbf91420..a352c4e645ec6579d1ce390fa257a46a3544f066 100644
--- a/hysop/operator/baroclinic_from_rhs.py
+++ b/hysop/operator/baroclinic_from_rhs.py
@@ -7,7 +7,7 @@ MultiPhase baroclinic term
from hysop.operator.computational import Computational
from hysop.operator.discrete.baroclinic_from_rhs import BaroclinicFromRHS as BD
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.constants import debug
import hysop.default_methods as default
from hysop.operator.continuous import opsetup
diff --git a/hysop/operator/discrete/drag_and_lift.py b/hysop/operator/discrete/drag_and_lift.py
index 875b9ee3b38cab61ead102183ac928eb406fd507..73f3fd550a06b5face864dc76309f271423e7137 100644
--- a/hysop/operator/discrete/drag_and_lift.py
+++ b/hysop/operator/discrete/drag_and_lift.py
@@ -10,7 +10,7 @@ from hysop.constants import XDIR, YDIR, ZDIR
from hysop.domain.control_box import ControlBox
from hysop.domain.subsets import Subset
from hysop.numerics.differential_operations import Laplacian
-from hysop.numerics.finite_differences import FD_C_2
+from hysop.numerics.finite_differences import FDC2
import numpy as np
@@ -462,8 +462,8 @@ class NocaForces(Forces):
res[i_n] += self._coeff * self.nu * normal * npw.real_sum(buff)
# 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.
- fd_scheme = FD_C_2(self._topology.mesh.space_step, ind,
+ # Default = FDC2. Todo : set this as an input method value.
+ fd_scheme = FDC2(self._topology.mesh.space_step, ind,
reduce_output_shape=True)
fd_scheme.compute(vd[i_n], i_n, buff)
res[i_n] += 2.0 * normal * self.nu * npw.real_sum(buff)
diff --git a/hysop/operator/discrete/reprojection.py b/hysop/operator/discrete/reprojection.py
index a883d293c968ff713e784013010e3adb405d247e..86506809712df2b94b312f4d3e7aa85d334a86a7 100644
--- a/hysop/operator/discrete/reprojection.py
+++ b/hysop/operator/discrete/reprojection.py
@@ -7,7 +7,7 @@ import numpy as np
from hysop.constants import debug, HYSOP_MPI_REAL
from hysop.methods_keys import SpaceDiscretisation
from hysop.operator.discrete.discrete import DiscreteOperator
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.numerics.differential_operations import GradV
import hysop.tools.numpywrappers as npw
from hysop.numerics.update_ghosts import UpdateGhosts
@@ -29,7 +29,7 @@ class Reprojection(DiscreteOperator):
@param frequency : set frequency of execution of the reprojection
"""
if 'method' in kwds and kwds['method'] is None:
- kwds['method'] = {SpaceDiscretisation: FD_C_4}
+ kwds['method'] = {SpaceDiscretisation: FDC4}
## vorticity field
self.vorticity = vorticity
diff --git a/hysop/operator/drag_and_lift.py b/hysop/operator/drag_and_lift.py
index b20ff77544129586c3a039220ec9271186a00121..f0580dd042885e1873d223325edb6b99d1dcf1f1 100755
--- a/hysop/operator/drag_and_lift.py
+++ b/hysop/operator/drag_and_lift.py
@@ -201,11 +201,11 @@ class NocaForces(Forces):
import hysop.default_methods as default
self.method = default.FORCES
from hysop.methods_keys import SpaceDiscretisation
- from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+ from hysop.numerics.finite_differences import FDC4, FDC2
assert SpaceDiscretisation in self.method.keys()
- if SpaceDiscretisation is FD_C_2:
+ if SpaceDiscretisation is FDC2:
self._min_ghosts = 1
- elif SpaceDiscretisation is FD_C_4:
+ elif SpaceDiscretisation is FDC4:
self._min_ghosts = 2
if surfdir is None:
diff --git a/hysop/operator/forcing.py b/hysop/operator/forcing.py
index bfe03b4e5bf98f0ef60b389dc6608df88f858651..d5cf8314209699b25400d22fb62258d8711dafbc 100644
--- a/hysop/operator/forcing.py
+++ b/hysop/operator/forcing.py
@@ -14,8 +14,8 @@ from hysop.constants import debug
from hysop.operator.continuous import opsetup
import hysop.default_methods as default
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4,\
- FD_C_2
+from hysop.numerics.finite_differences import FDC4,\
+ FDC2
from hysop.operator.differential import Curl
from hysop.fields.continuous import Field
@@ -113,11 +113,11 @@ class ForcingConserv(Forcing):
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
# Finite differences method
# Minimal number of ghost points
nb_ghosts = 2
- elif self.method[SpaceDiscretisation] is FD_C_2:
+ elif self.method[SpaceDiscretisation] is FDC2:
nb_ghosts = 1
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/low_pass_filt.py b/hysop/operator/low_pass_filt.py
index 91f03f28be47654913a9af35b81b612e455cf9ed..a86381ce800cef3ef4aab3a5ecd0dfe2d85b4ecb 100644
--- a/hysop/operator/low_pass_filt.py
+++ b/hysop/operator/low_pass_filt.py
@@ -13,8 +13,8 @@ from hysop.constants import debug
from hysop.operator.continuous import opsetup
import hysop.default_methods as default
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4,\
- FD_C_2
+from hysop.numerics.finite_differences import FDC4,\
+ FDC2
from hysop.operator.differential import Curl
from hysop.fields.continuous import Field
@@ -110,11 +110,11 @@ class LowPassFiltConserv(LowPassFilt):
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
# Finite differences method
# Minimal number of ghost points
nb_ghosts = 2
- elif self.method[SpaceDiscretisation] is FD_C_2:
+ elif self.method[SpaceDiscretisation] is FDC2:
nb_ghosts = 1
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/multiphase_gradp.py b/hysop/operator/multiphase_gradp.py
index f3ad480e5fa1c43fb265cef3e340bf366e6c1091..a9411d3469a67a446ff14ed0c889263d4a5a39d2 100644
--- a/hysop/operator/multiphase_gradp.py
+++ b/hysop/operator/multiphase_gradp.py
@@ -9,7 +9,7 @@ Computation of the pressure gradient in a multiphasic flow:
"""
from hysop.operator.computational import Computational
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.constants import debug, np
import hysop.default_methods as default
from hysop.operator.continuous import opsetup
@@ -51,7 +51,7 @@ class MultiphaseGradP(Computational):
assert SpaceDiscretisation in self.method.keys()
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
nbGhosts = 2
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/penalization.py b/hysop/operator/penalization.py
index 2cabff25dfc17f6885c3aa148dd4b6bce23052f4..bb5e73c66fb772f122e70807f1841cb2aa2eff52 100644
--- a/hysop/operator/penalization.py
+++ b/hysop/operator/penalization.py
@@ -18,8 +18,8 @@ from hysop.operator.continuous import opsetup
from hysop.domain.subsets import Subset
import hysop.default_methods as default
from hysop.methods_keys import SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4,\
- FD_C_2
+from hysop.numerics.finite_differences import FDC4,\
+ FDC2
from hysop.operator.differential import Curl
from hysop.fields.continuous import Field
@@ -139,11 +139,11 @@ class PenalizeVorticity(Penalization):
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
# Finite differences method
# Minimal number of ghost points
nb_ghosts = 2
- elif self.method[SpaceDiscretisation] is FD_C_2:
+ elif self.method[SpaceDiscretisation] is FDC2:
nb_ghosts = 1
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/stretching.py b/hysop/operator/stretching.py
index 5320c192e7912fcf210cbf13eb0e44d2bc06692a..ad6371f311f2d86a9ef0fa8f5985a609ac341948 100755
--- a/hysop/operator/stretching.py
+++ b/hysop/operator/stretching.py
@@ -11,7 +11,7 @@ See also
from hysop.constants import debug
from hysop.methods_keys import TimeIntegrator, Formulation, \
SpaceDiscretisation
-from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.finite_differences import FDC4
from hysop.operator.computational import Computational
from hysop.operator.continuous import opsetup
from hysop.operator.discrete.stretching import Conservative, GradUW
@@ -113,7 +113,7 @@ class Stretching(Computational):
@debug
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
nbghosts = 2
else:
raise ValueError("Unknown method for space discretization of the\
@@ -175,7 +175,7 @@ class StretchingLinearized(Stretching):
@debug
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
+ if self.method[SpaceDiscretisation] is FDC4:
nbghosts = 2
else:
raise ValueError("Unknown method for space discretization of the\
diff --git a/hysop/operator/tests/test_Stretching.py b/hysop/operator/tests/test_Stretching.py
index 159d506352946f3fabbb1be902459aa9a73ebf7c..ad76dca56e977ee3d6bf9f4a2d226d7213b3e0be 100755
--- a/hysop/operator/tests/test_Stretching.py
+++ b/hysop/operator/tests/test_Stretching.py
@@ -7,7 +7,7 @@ from hysop.operator.stretching import Stretching, \
from hysop.problem.simulation import Simulation
from hysop.methods_keys import TimeIntegrator, Formulation,\
SpaceDiscretisation
-from hysop.methods import Euler, RK3, FD_C_4, Conservative
+from hysop.methods import Euler, RK3, FDC4, Conservative
from hysop.tools.parameters import Discretization
import hysop.tools.numpywrappers as npw
pi = np.pi
@@ -116,7 +116,7 @@ def test_stretching():
# Operators
#method = {TimeIntegrator: RK3, Formulation: Conservative,
- # SpaceDiscretisation: FD_C_4}
+ # SpaceDiscretisation: FDC4}
stretch = Stretching(velo, vorti, discretization=nbElem)
stretch.discretize()
topo = stretch.discreteFields[velo].topology
@@ -199,7 +199,7 @@ def test_stretching_external_work():
# Operators
#method = {TimeIntegrator: RK3, Formulation: Conservative,
- # SpaceDiscretisation: FD_C_4}
+ # SpaceDiscretisation: FDC4}
stretch = Stretching(velo, vorti, discretization=nbElem)
stretch.discretize()
wk_p = stretch.get_work_properties()
diff --git a/hysop/operator/tests/test_differential.py b/hysop/operator/tests/test_differential.py
index 2ab50a59af60b251ff00bf1df7ff428f866b7aab..9eaa3895b537fd371434e01bcad3d3fd734697d3 100755
--- a/hysop/operator/tests/test_differential.py
+++ b/hysop/operator/tests/test_differential.py
@@ -6,7 +6,7 @@ from hysop.domain.box import Box
from hysop.fields.continuous import Field
import hysop.tools.numpywrappers as npw
from hysop.methods_keys import SpaceDiscretisation
-from hysop.methods import FD_C_4, FD_C_2
+from hysop.methods import FDC4, FDC2
from hysop.operator.differential import Curl, Grad, DivAdvection
from hysop.tools.parameters import Discretization
@@ -154,28 +154,28 @@ def call_op_fft(class_name, ref_formula, dom, discr,
def test_curl_fd_1():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Curl, vorticity_f, topo1_3d, method=method)
def test_curl_fd_2():
- method = {SpaceDiscretisation: FD_C_2}
+ method = {SpaceDiscretisation: FDC2}
call_op(Curl, vorticity_f, topo1_3d, method=method, order=2)
def test_curl_fd_1_2d():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Curl, vorticity_f2d, topo1_2d, method=method,
op_dim=1, vform=velocity_f2d)
def test_curl_fd_3():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Curl, vorticity_f, topo3_3d, method=method)
def test_curl_fd_3_2d():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Curl, vorticity_f2d, topo3_2d, op_dim=1,
method=method, vform=velocity_f2d)
@@ -214,28 +214,28 @@ def test_curl_fft_1():
def test_grad_1():
- method = {SpaceDiscretisation: FD_C_2}
+ method = {SpaceDiscretisation: FDC2}
call_op(Grad, grad_velo, topo1_3d, op_dim=9, method=method, order=2)
def test_grad_2():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Grad, grad_velo, topo1_3d, op_dim=9, method=method)
def test_grad_3():
- method = {SpaceDiscretisation: FD_C_2}
+ method = {SpaceDiscretisation: FDC2}
call_op(Grad, grad_velo, topo3_3d, op_dim=9, method=method, order=2)
def test_grad_3_2d():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Grad, grad_velo_2d, topo3_2d, op_dim=4, method=method,
vform=velocity_f2d)
def test_curl_fd_work():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(Curl, vorticity_f, topo3_3d, use_work=True, method=method)
@@ -246,5 +246,5 @@ def divadvection_func(res, x, y, z, t):
def test_div_advection():
- method = {SpaceDiscretisation: FD_C_4}
+ method = {SpaceDiscretisation: FDC4}
call_op(DivAdvection, divadvection_func, topo3_3d, op_dim=1, method=method)
diff --git a/hysop/operator/tests/test_multiphase_gradp.py b/hysop/operator/tests/test_multiphase_gradp.py
index 1f8fc2e1f22d1a6a8d5adafc63f3422c415f7c00..d3e583c3ca297d6dc7680e3664f3c7a21b1a5de9 100755
--- a/hysop/operator/tests/test_multiphase_gradp.py
+++ b/hysop/operator/tests/test_multiphase_gradp.py
@@ -4,7 +4,7 @@ from hysop.tools.parameters import Discretization
from hysop.domain.box import Box
from hysop.fields.continuous import Field
from hysop.operator.multiphase_gradp import MultiphaseGradP
-from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+from hysop.numerics.finite_differences import FDC4, FDC2
import hysop.tools.numpywrappers as npw
import numpy as np
from hysop.methods_keys import Support, SpaceDiscretisation, ExtraArgs