From a8223323e4cc9e51d5c7be426302461f73a09627 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Tue, 19 Jun 2018 11:07:28 +0200
Subject: [PATCH] Fix hostarray nanmin and nanmax. Fix default timestep
parameter.Fix Fortran FFTW diffusion operator. Add fftw diffusion test.
---
ci/scripts/test.sh | 1 +
hysop/__init__.py.in | 6 +-
.../host/fortran/operator/diffusion.py | 50 ++--
.../host/fortran/operator/fortran_fftw.py | 32 +-
hysop/backend/host/host_array.py | 43 +--
.../core/graph/computational_node_frontend.py | 39 ++-
hysop/operator/diffusion.py | 26 +-
hysop/operator/tests/test_diffusion.py | 277 ++++++++++++++++++
.../tests/test_directional_diffusion.py | 91 +++---
hysop/operators.py | 1 +
hysop/parameters/default_parameters.py | 9 +-
11 files changed, 432 insertions(+), 143 deletions(-)
create mode 100644 hysop/operator/tests/test_diffusion.py
diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index be5da252b..40cf839ca 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -61,6 +61,7 @@ python "$HYSOP_DIR/operator/tests/test_solenoidal_projection.py"
python "$HYSOP_DIR/operator/tests/test_custom_symbolic.py"
python "$HYSOP_DIR/operator/tests/test_directional_advection.py"
python "$HYSOP_DIR/operator/tests/test_directional_diffusion.py"
+python "$HYSOP_DIR/operator/tests/test_diffusion.py"
python "$HYSOP_DIR/operator/tests/test_directional_stretching.py"
if [ "$HAS_CACHE_DIR" = true ]; then
diff --git a/hysop/__init__.py.in b/hysop/__init__.py.in
index f63bf48fa..a12fafcf2 100644
--- a/hysop/__init__.py.in
+++ b/hysop/__init__.py.in
@@ -32,7 +32,6 @@ __OPTIMIZE__ = not __debug__
__VERBOSE__ = get_env('VERBOSE', ("@VERBOSE@" is "ON"))
__DEBUG__ = get_env('DEBUG', ("@DEBUG@" is "ON"))
__PROFILE__ = get_env('PROFILE', ("@PROFILE@" is "ON"))
-__KERNEL_DEBUG__ = get_env('KERNEL_DEBUG', False)
__TRACE_CALLS__ = get_env('TRACE_CALLS', False)
__TRACE_WARNINGS__ = get_env('TRACE_WARNINGS', False)
@@ -80,7 +79,7 @@ def sprint(*args, **kargs):
"""prints only shared memory master rank 0"""
if (shm_rank==0):
vprint(*args, **kargs)
-
+
def reset():
"""Reset hysop state (registered objects)."""
from hysop.tools.handle import RegisteredObject
@@ -114,7 +113,7 @@ from hysop.problem import Problem
from hysop.tools.io_utils import IO, IOParams
__all__ = ['Box', 'Field', 'DiscreteField', 'ScalarParameter', 'TensorParameter',
'Domain', 'Discretization', 'Simulation', 'MPIParams',
- 'Problem', 'IO', 'IOParams',
+ 'Problem', 'IO', 'IOParams',
'Topology', 'CartesianTopology', 'TopologyDescriptor']
if __MPI_ENABLED__:
__all__ += ['MPI', 'main_rank', 'main_size']
@@ -138,4 +137,3 @@ msg_io = '\n*Default path for all i/o is \'{}\'.'.format(default_path)
msg_io += '\n*Default path for caching is \'{}\'.'.format(cache_path)
mprint(msg_io)
mprint()
-
diff --git a/hysop/backend/host/fortran/operator/diffusion.py b/hysop/backend/host/fortran/operator/diffusion.py
index 2dea7de14..7d7e8b474 100644
--- a/hysop/backend/host/fortran/operator/diffusion.py
+++ b/hysop/backend/host/fortran/operator/diffusion.py
@@ -1,20 +1,20 @@
-
from hysop.backend.host.fortran.operator.fortran_fftw import FortranFFTWOperator, fftw2py
from hysop.tools.types import check_instance, InstanceOf
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.core.graph.graph import op_apply
class DiffusionFFTW(FortranFFTWOperator):
__default_method = {}
__available_methods = {}
-
- def __init__(self, input_field, output_field, viscosity, variables, **kargs):
+
+ def __init__(self, input_field, output_field, viscosity, variables, dt, **kargs):
"""Diffusion operator using FFTW in Fortran.
-
+
Parameters
----------
input_field : :class:`~hysop.fields.continuous_field.Field
@@ -22,9 +22,11 @@ class DiffusionFFTW(FortranFFTWOperator):
variables: dictionary of fields:topology
viscosity : double
constant viscosity value
- kargs :
+ dt: ScalarParameter
+ Timestep parameter that will be used for time integration.
+ kargs :
base class parameters.
-
+
Notes:
*Equations:
dW/dt = K*Laplacian(W)
@@ -35,7 +37,8 @@ class DiffusionFFTW(FortranFFTWOperator):
check_instance(input_field, Field)
check_instance(output_field, Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
-
+ check_instance(dt, ScalarParameter)
+
assert variables[input_field] == variables[output_field], \
'input and output topology mismatch'
assert input_field.domain is output_field.domain,\
@@ -43,13 +46,16 @@ class DiffusionFFTW(FortranFFTWOperator):
input_fields = { input_field: variables[input_field] }
output_fields = { output_field: variables[output_field] }
+ input_params = {dt.name: dt}
super(DiffusionFFTW,self).__init__(input_fields=input_fields, output_fields=output_fields,
- **kargs)
-
+ input_params=input_params,
+ **kargs)
+
self.input_field = input_field
self.output_field = output_field
self.viscosity = viscosity
-
+ self.dt = dt
+
def initialize(self, **kwds):
super(DiffusionFFTW,self).initialize(**kwds)
@@ -59,9 +65,9 @@ class DiffusionFFTW(FortranFFTWOperator):
elif (dim==3):
self._solve = self._solve_3d
else:
- raise AttributeError(dim + "D case not yet implemented.")
+ raise AttributeError(str(dim) + "D case not yet implemented.")
+
-
@debug
def discretize(self):
if self.discretized:
@@ -71,22 +77,20 @@ class DiffusionFFTW(FortranFFTWOperator):
self.dout = self.output_field.discrete_fields[self.topology]
@op_apply
- def apply(self, simulation=None, **kargs):
- super(DiffusionFFTW,self).apply(simulation=simulation,**kargs)
- self._solve(simulation)
-
- def _solve_2d(self, simulation):
+ def apply(self, **kargs):
+ super(DiffusionFFTW,self).apply(**kargs)
+ self._solve(self.dt())
+
+ def _solve_2d(self, dt):
""" Solve 2D diffusion problem
"""
- dt = simulation.time_step
ghosts = self.topology.ghosts
self.dout.data =\
fftw2py.solve_diffusion_2d(self.viscosity * dt, self.din.data[0], ghosts)
-
- def _solve_3d(self,simulation):
+
+ def _solve_3d(self, dt):
""" Solve 3D diffusion problem
"""
- dt = simulation.time_step
ghosts = self.topology.ghosts
self.dout.data = \
fftw2py.solve_diffusion_3d(self.viscosity * dt,
@@ -94,12 +98,10 @@ class DiffusionFFTW(FortranFFTWOperator):
self.din.data[1],
self.din.data[2],
ghosts)
-
+
def available_methods(self):
return self.__available_methods
def default_method(self):
return self.__default_method
def handle_method(self, method):
pass
-
-
diff --git a/hysop/backend/host/fortran/operator/fortran_fftw.py b/hysop/backend/host/fortran/operator/fortran_fftw.py
index a551c5879..f1a730cac 100644
--- a/hysop/backend/host/fortran/operator/fortran_fftw.py
+++ b/hysop/backend/host/fortran/operator/fortran_fftw.py
@@ -1,4 +1,3 @@
-
import hysop
try:
@@ -21,9 +20,9 @@ class FortranFFTWOperator(ComputationalGraphOperator):
@debug
def __init__(self, input_fields, output_fields, **kwds):
- super(FortranFFTWOperator, self).__init__(input_fields=input_fields,
+ super(FortranFFTWOperator, self).__init__(input_fields=input_fields,
output_fields=output_fields, **kwds)
-
+
check_instance(input_fields, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(output_fields, dict, keys=Field, values=CartesianTopologyDescriptors)
@@ -34,14 +33,14 @@ class FortranFFTWOperator(ComputationalGraphOperator):
msg+='\nHYSOP_REAL is {} but field {} has dtype {}.'
msg=msg.format(HYSOP_REAL.__name__, f.name, f.dtype)
raise RuntimeError(msg)
-
+
domain = self.input_fields.keys()[0].domain
self.dim = domain.dim
self.domain = domain
-
+
def handle_topologies(self, input_topology_states, output_topology_states):
super(FortranFFTWOperator,self).handle_topologies(input_topology_states, output_topology_states)
-
+
topology = self.input_fields.values()[0].topology
domain = self.domain
for (field,topoview) in self.input_fields.iteritems():
@@ -51,11 +50,12 @@ class FortranFFTWOperator(ComputationalGraphOperator):
assert topoview.topology is topology, 'topology mismatch'
assert field.domain is domain, 'domain mismatch'
self.topology = topology
-
+
@debug
def get_field_requirements(self):
requirements = super(FortranFFTWOperator, self).get_field_requirements()
-
+ dim = self.domain.dim
+
# set can_split to True in all directions except the contiguous one
# for inputs and outputs
for is_input, (field, td, req) in requirements.iter_requirements():
@@ -63,30 +63,32 @@ class FortranFFTWOperator(ComputationalGraphOperator):
can_split[:-1] = False
can_split[-1] = True
req.can_split = can_split
+ req.min_ghosts = (0, ) * dim
+ req.max_ghosts = (0, ) * dim
return requirements
-
+
@debug
def discretize(self):
if self.discretized:
return
super(FortranFFTWOperator,self).discretize()
self._fftw_discretize()
-
+
@debug
def get_work_properties(self):
return super(FortranFFTWOperator,self).get_work_properties()
-
+
@debug
def setup(self, work=None):
super(FortranFFTWOperator,self).setup(work=work)
-
+
@debug
def finalize(self, clean_fftw_solver=False, **kwds):
super(FortranFFTWOperator,self).finalize(**kwds)
if clean_fftw_solver:
fftw2py.clean_fftw_solver(self.dim)
-
+
def _fftw_discretize(self):
"""
fftw specific way to discretize variables for a given
@@ -94,12 +96,12 @@ class FortranFFTWOperator(ComputationalGraphOperator):
It is assumed that in fft case, only one topology must be used
for all variables of all fortran fftw operators.
"""
-
+
self._set_domain_and_tasks()
topo = self.topology
comm = self.mpi_params.comm
size = self.mpi_params.size
-
+
msg = 'input topology is not compliant with fftw.'
assert topo.cart_dim == 1, msg
diff --git a/hysop/backend/host/host_array.py b/hysop/backend/host/host_array.py
index d333f137e..bd2b2370b 100644
--- a/hysop/backend/host/host_array.py
+++ b/hysop/backend/host/host_array.py
@@ -1,4 +1,3 @@
-
from hysop.deps import np
from hysop.core.arrays import MemoryType, MemoryOrdering
from hysop.core.arrays import default_order
@@ -8,29 +7,29 @@ from hysop.backend.host.host_array_backend import HostArrayBackend
class HostArray(Array):
"""
Host memory array wrapper.
- An HostArray is a numpy.ndarray work-alike that stores its data and performs
+ An HostArray is a numpy.ndarray work-alike that stores its data and performs
its computations on CPU with numpy.
"""
def __init__(self, handle, backend, **kwds):
"""
Build an HostArray instance.
-
+
Parameters
----------
handle: numpy.ndarray
implementation of this array
backend: HostArrayBackend
backend used to build this handle
- kwds:
+ kwds:
arguments for base classes.
Notes
-----
This should never be called directly by the user.
- Arrays should be constructed using array backend facilities, like zeros or empty.
+ Arrays should be constructed using array backend facilities, like zeros or empty.
The parameters given here refer to a low-level method for instantiating an array.
"""
-
+
if not isinstance(handle, np.ndarray):
msg='Handle should be a np.ndarray but got a {}.'
msg=msg.format(handle.__class__.__name__)
@@ -39,14 +38,14 @@ class HostArray(Array):
msg='Backend should be a HostArrayBackend but got a {}.'
msg=msg.format(handle.__class__.__name__)
raise ValueError(msg)
-
+
# if handle.dtype in [np.bool]:
- # msg='{} unsupported yet, use HYSOP_BOOL={} instead.'.format(handle.dtype,
+ # msg='{} unsupported yet, use HYSOP_BOOL={} instead.'.format(handle.dtype,
# HYSOP_BOOL)
# raise TypeError(msg)
-
+
super(HostArray,self).__init__(handle=handle, backend=backend, **kwds)
-
+
# array interface
self.__array_interface__ = handle.__array_interface__
@@ -56,14 +55,14 @@ class HostArray(Array):
"""
from hysop.symbolic.array import HostSymbolicArray
return HostSymbolicArray(memory_object=self, name=name, **kwds)
-
+
def as_symbolic_buffer(self, name, **kwds):
"""
Return a symbolic buffer variable that contain a reference to this array.
"""
from hysop.symbolic.array import HostSymbolicBuffer
return HostSymbolicBuffer(memory_object=self, name=name, **kwds)
-
+
def get_ndim(self):
return self.handle.ndim
def get_shape(self):
@@ -98,8 +97,8 @@ class HostArray(Array):
return self.handle.nbytes
def get_int_ptr(self):
return int(self.handle.ctypes.data)
-
- # array properties
+
+ # array properties
ndim = property(get_ndim)
shape = property(get_shape)
offset = property(get_offset)
@@ -116,7 +115,7 @@ class HostArray(Array):
itemsize = property(get_itemsize)
nbytes = property(get_nbytes)
int_ptr = property(get_int_ptr)
-
+
def get(self, handle=False):
"""
Returns equivalent array on host device, ie. self.
@@ -139,12 +138,12 @@ class HostArray(Array):
"""
self.handle.flags.writeable = True
return self
-
+
def tofile(self, fid, sep='', format='%s'):
"""
Write array to a file as text or binary (default).
- This is a convenience function for quick storage of array data.
+ This is a convenience function for quick storage of array data.
Information on endianness and precision is lost.
"""
self.handle.tofile(fid=fid,sep=sep,format=format)
@@ -171,12 +170,12 @@ class HostArray(Array):
"""
handle = self.handle.take(indices=indices, axis=axis, out=out, mode=mode)
return self.wrap(handle)
- def astype(self, dtype, order=MemoryOrdering.SAME_ORDER,
+ def astype(self, dtype, order=MemoryOrdering.SAME_ORDER,
casting='unsafe', subok=True, copy=True):
"""
Copy of the array, cast to a specified type.
"""
- handle = self.handle.astype(dtype=dtype, order=order, casting=casting,
+ handle = self.handle.astype(dtype=dtype, order=order, casting=casting,
subok=subok, copy=copy)
return self.wrap(handle)
def view(self, dtype):
@@ -185,3 +184,9 @@ class HostArray(Array):
"""
handle = self.handle.view(dtype=dtype)
return self.wrap(handle)
+
+ def nanmin(self, axis=None, out=None):
+ return self.backend.nanmin(a=self, axis=None, out=None)
+
+ def nanmax(self, axis=None, out=None):
+ return self.backend.nanmax(a=self, axis=None, out=None)
diff --git a/hysop/core/graph/computational_node_frontend.py b/hysop/core/graph/computational_node_frontend.py
index 1b8df9cac..bfb7d5de9 100644
--- a/hysop/core/graph/computational_node_frontend.py
+++ b/hysop/core/graph/computational_node_frontend.py
@@ -1,4 +1,3 @@
-
from hysop.constants import Implementation, Backend, implementation_to_backend
from hysop.tools.decorators import debug
from hysop.tools.types import check_instance, first_not_None
@@ -23,8 +22,8 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
base_kwds: dict, optional, defaults to None
Base class keywords arguments.
If None, an empty dict will be passed.
- impl_kwds:
- Keywords arguments that will be passed towards implementation implemention
+ impl_kwds:
+ Keywords arguments that will be passed towards implementation implemention
__init__.
Attributes
@@ -35,19 +34,19 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
the backend corresponding to implementation
impl: ComputationalGraphNodeGenerator or ComputationalGraphNode
the implementation class
- impl_kwds:
- Keywords arguments that will be passed towards implementation implemention
+ impl_kwds:
+ Keywords arguments that will be passed towards implementation implemention
impl.__init__ during a call to _generate.
"""
-
+
base_kwds = first_not_None(base_kwds, {})
super(ComputationalGraphNodeFrontend,self).__init__(**base_kwds)
check_instance(implementation, Implementation, allow_none=True)
-
+
default_implementation = self.__class__.default_implementation()
- available_implementations = self.__class__.available_implementations()
-
+ available_implementations = self.__class__.available_implementations()
+
if not isinstance(default_implementation, Implementation):
msg='default_implementation is not a instance of hysop.backend.Implementation.'
raise TypeError(msg)
@@ -60,7 +59,7 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
if default_implementation not in available_implementations:
msg='default_implementation is not contained in available_implementations.'
raise ValueError(msg)
-
+
if implementation is None:
implementation = default_implementation
elif implementation not in available_implementations:
@@ -72,7 +71,7 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
msg = 'Specified implementation \'{}\' is not an available implementation, '
msg+= 'available implementations are:\n {}'
- msg=msg.format(implementation, '\n '.join(implementations))
+ msg=msg.format(implementation, '\n '.join(simplementations))
raise ValueError(msg)
self.implementation = implementation
@@ -82,13 +81,13 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
self._generated = False
self._input_fields_to_dump = []
self._output_fields_to_dump = []
-
+
@debug
def _generate(self):
try:
op = self.impl(**self.impl_kwds)
except:
- sargs = ['*{} = {}'.format(k,v.__class__)
+ sargs = ['*{} = {}'.format(k,v.__class__)
for (k,v) in self.impl_kwds.iteritems()]
msg = 'FATAL ERROR during {}.generate():\n'
msg += ' => failed to call {}.__init__()\n with the following keywords:'
@@ -96,7 +95,7 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
msg = msg.format(self.__class__, self.impl)
print '\n{}\n'.format(msg)
raise
-
+
for kwds in self._input_fields_to_dump:
op.dump_inputs(**kwds)
for kwds in self._output_fields_to_dump:
@@ -104,13 +103,13 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
self._generated = True
return (op,)
-
+
@classmethod
@not_implemented
def implementations(cls):
"""
Should return all implementations as a dictionnary.
- Keys are Implementation instances and values are either ComputationalGraphNode
+ Keys are Implementation instances and values are either ComputationalGraphNode
or ComputationalGraphNodeGenerator.
"""
pass
@@ -122,7 +121,7 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
Return the default Implementation, should be compatible with available_implementations.
"""
pass
-
+
@classmethod
def available_implementations(cls):
"""
@@ -134,7 +133,7 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
"""
Tell the generated operator to dump some of its inputs before
apply is called.
-
+
Target folder, file, dump frequency and other io pameters
are passed trough io_params or as keywords.
@@ -144,12 +143,12 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
msg=msg.format(self.name)
assert (self._generated is False), msg
self._input_fields_to_dump.append(kwds)
-
+
def dump_outputs(self, **kwds):
"""e
Tell the generated operator to dump some of its outputs after
apply is called.
-
+
Target folder, file, dump frequency and other io pameters
are passed trough instance io_params of this parameter or
as keywords.
diff --git a/hysop/operator/diffusion.py b/hysop/operator/diffusion.py
index c594ed67d..dabce77bc 100644
--- a/hysop/operator/diffusion.py
+++ b/hysop/operator/diffusion.py
@@ -1,4 +1,3 @@
-
"""
@file advection_dir.py
Directional advection operator generator.
@@ -9,6 +8,7 @@ from hysop.tools.enum import EnumFactory
from hysop.tools.decorators import debug
from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
from hysop.backend.host.fortran.operator.diffusion import DiffusionFFTW
@@ -18,7 +18,7 @@ class Diffusion(ComputationalGraphNodeFrontend):
Interface the poisson solver.
Available implementations are: FFTW
"""
-
+
__implementations = {
Implementation.FORTRAN: DiffusionFFTW
}
@@ -26,18 +26,18 @@ class Diffusion(ComputationalGraphNodeFrontend):
@classmethod
def implementations(cls):
return cls.__implementations
-
+
@classmethod
def default_implementation(cls):
return Implementation.FORTRAN
-
+
@debug
- def __init__(self, input_field, variables, viscosity,
+ def __init__(self, input_field, variables, viscosity, dt,
output_field=None, implementation=None, base_kwds=None, **kwds):
"""
Initialize a Poisson operator frontend.
Solves dF/dt = viscosity * Laplacian(F)
-
+
Parameters
----------
input_field: field
@@ -49,6 +49,8 @@ class Diffusion(ComputationalGraphNodeFrontend):
dictionary of fields as keys and topologies as values.
viscosity: scalar or Field
Warning: some implementations may only offer scalar viscosity.
+ dt: ScalarParameter
+ Timestep parameter that will be used for time integration.
implementation: Implementation, optional, defaults to None
target implementation, should be contained in available_implementations().
If None, implementation will be set to default_implementation().
@@ -56,13 +58,13 @@ class Diffusion(ComputationalGraphNodeFrontend):
Base class keywords arguments.
If None, an empty dict will be passed.
kwds:
- Keywords arguments that will be passed towards implementation
+ Keywords arguments that will be passed towards implementation
poisson operator __init__.
Notes
-----
- A diffusion operator implementation should at least support
- the following __init__ attributes:
+ A diffusion operator implementation should at least support
+ the following __init__ attributes:
input_field, output_field, variables, viscosity
"""
output_field = output_field or input_field
@@ -72,7 +74,9 @@ class Diffusion(ComputationalGraphNodeFrontend):
check_instance(output_field, Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str)
+ check_instance(dt, ScalarParameter)
super(Diffusion, self).__init__(input_field=input_field, output_field=output_field,
- variables=variables, viscosity=viscosity, implementation=implementation, base_kwds=base_kwds,
- **kwds)
+ variables=variables, viscosity=viscosity, dt=dt,
+ implementation=implementation, base_kwds=base_kwds,
+ **kwds)
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
new file mode 100644
index 000000000..50166e67b
--- /dev/null
+++ b/hysop/operator/tests/test_diffusion.py
@@ -0,0 +1,277 @@
+from hysop.deps import sys
+from hysop.tools.numpywrappers import npw
+import numpy as np
+from hysop.testsenv import __ENABLE_LONG_TESTS__
+from hysop.testsenv import iter_clenv
+from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.tools.contexts import printoptions
+from hysop.tools.types import first_not_None, to_tuple
+from hysop.tools.numerics import is_integer, is_fp
+from hysop.tools.io_utils import IO
+from hysop.tools.sympy_utils import enable_pretty_printing
+from hysop.operator.directional.diffusion_dir import DirectionalDiffusion
+from hysop.operator.diffusion import Diffusion
+
+from hysop import Field, Box, Simulation
+from hysop.methods import Remesh, TimeIntegrator
+from hysop.constants import Implementation, DirectionLabels, ComputeGranularity, SpaceDiscretization
+from hysop.numerics.splitting.strang import StrangSplitting, StrangOrder
+from hysop.numerics.remesh.remesh import RemeshKernel
+from hysop.numerics.stencil.stencil_generator import StencilGenerator, CenteredStencilGenerator, MPQ
+from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator, Euler, RK2, RK3, RK4, RK4_38
+
+class TestDirectionalDiffusionOperator(object):
+
+ @classmethod
+ def setup_class(cls,
+ enable_extra_tests=__ENABLE_LONG_TESTS__,
+ enable_debug_mode=False):
+
+ IO.set_default_path('/tmp/hysop_tests/test_directional_diffusion')
+
+ if enable_debug_mode:
+ cls.size_min = 8
+ cls.size_max = 9
+ else:
+ cls.size_min = 11
+ cls.size_max = 23
+
+ cls.enable_extra_tests = enable_extra_tests
+ cls.enable_debug_mode = enable_debug_mode
+
+ @classmethod
+ def teardown_class(cls):
+ pass
+
+ def _test(self, dim, is_inplace,
+ size_min=None, size_max=None):
+ assert dim > 0
+
+ # periodic boundaries removes one computational point
+ # so we add one here.
+ size_min = first_not_None(size_min, self.size_min) + 1
+ size_max = first_not_None(size_max, self.size_max) + 1
+
+ shape = tuple(npw.random.randint(low=size_min, high=size_max, size=dim).tolist())
+ shape = (33,33,33)
+ shape = (65,65,65)
+
+ if self.enable_extra_tests:
+ flt_types = (npw.float32, npw.float64)
+ time_integrators = (Euler, RK2,)
+ orders = (2, 4, 6)
+ else:
+ flt_types = (npw.float64,)
+ time_integrators = (RK2, )
+ orders = (4,)
+
+ domain = Box(length=(2*npw.pi,)*dim)
+ for dtype in flt_types:
+ Fin = Field(domain=domain, name='Fin', dtype=dtype,
+ nb_components=3, register_object=False)
+ Fout = Field(domain=domain, name='Fout', dtype=dtype,
+ nb_components=3, register_object=False)
+ coeffs = npw.random.rand(Fin.nb_components, dim).astype(dtype)
+ for order in orders:
+ for time_integrator in time_integrators:
+ print
+ self._test_one(time_integrator=time_integrator, order=order,
+ shape=shape, dim=dim, dtype=dtype,
+ is_inplace=is_inplace, domain=domain,
+ Fin=Fin, Fout=Fout, coeffs=coeffs)
+
+ @staticmethod
+ def __random_init(data, coords, dtype):
+ if is_integer(dtype):
+ for d in data:
+ d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+ elif is_fp(dtype):
+ for d in data:
+ d[...] = npw.random.random(size=d.shape)
+ else:
+ msg = 'Unknown dtype {}.'.format(dtype)
+ raise NotImplementedError(msg)
+
+ @staticmethod
+ def __scalar_init(data, coords, dtype):
+ if is_fp(dtype):
+ for (i,d) in enumerate(data):
+ d[...] = 1
+ for xi in coords:
+ d[...] *= npw.cos(xi*(i+1))
+ else:
+ msg = 'Unknown dtype {}.'.format(dtype)
+ raise NotImplementedError(msg)
+
+ def _test_one(self, time_integrator, order, shape, dim,
+ dtype, is_inplace, domain, Fin, Fout, coeffs):
+
+ print '\nTesting {}D DirectionalDiffusion_{}_FCD{}: inplace={} dtype={} shape={}'.format(
+ dim, time_integrator.name(), order,
+ is_inplace, dtype.__name__, shape),
+
+ dt = ScalarParameter('dt0', initial_value=0.1, const=True,
+ dtype=dtype)
+
+ fin = Fin
+ if is_inplace:
+ fout = Fin
+ variables = { Fin : shape }
+ else:
+ fout = Fout
+ variables = { Fin : shape, Fout: shape }
+ msg='Out of place diffusion has not been implemented yet.'
+ raise NotImplementedError(msg)
+
+ implementations = (Implementation.OPENCL, Implementation.FORTRAN)
+ ref_impl = Implementation.OPENCL
+
+ method = {TimeIntegrator: time_integrator, SpaceDiscretization: order}
+
+ def iter_impl(impl):
+ base_kwds = dict(fields=fin, coeffs=coeffs, dt=dt,
+ variables=variables, implementation=impl, laplacian_formulation=False,
+ method=method, name='test_diffusion_{}'.format(str(impl).lower()))
+ if (impl is Implementation.OPENCL):
+ for cl_env in iter_clenv():
+ msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
+ cl_env.device.name.strip())
+ da = DirectionalDiffusion(**base_kwds)
+ split = StrangSplitting(splitting_dim=dim,
+ extra_kwds=dict(cl_env=cl_env),
+ order=StrangOrder.STRANG_SECOND_ORDER)
+ split.push_operators(da)
+ yield msg, split
+ elif (impl is Implementation.FORTRAN):
+ diff = Diffusion(input_field=fin, viscosity=coeffs, dt=dt,
+ variables=variables, implementation=impl,
+ method=method, name='test_diffusion_{}'.format(str(impl).lower()))
+ yield "", diff.to_graph()
+ else:
+ msg='Unknown implementation to test {}.'.format(impl)
+ raise NotImplementedError(msg)
+
+ # Compare to other implementations
+ reference_fields = {}
+ outputs = {}
+
+ for impl in implementations:
+ print '\n >Implementation {}:'.format(impl),
+ for sop,op in iter_impl(impl):
+ print '\n *{}: '.format(sop),
+ if (not is_inplace):
+ op.push_copy(dst=fin, src=fout)
+ op.build()
+
+ dfin = op.input_discrete_fields[fin]
+ dfout = op.output_discrete_fields[fout]
+
+ view = dfin.compute_slices
+
+ sys.stdout.write('.')
+ sys.stdout.flush()
+ try:
+ dfout.initialize(self.__random_init, dtype=dfout.dtype)
+ dfin.initialize(self.__scalar_init, dtype=dfin.dtype)
+ S0 = dfin.integrate()
+
+ if not npw.all(npw.isfinite(S0)):
+ msg='Integral is not finite. Got {}.'.format(S0)
+ raise RuntimeError(msg)
+
+ _input = tuple(dfin[i].get().handle.copy()
+ for i in xrange(dfin.nb_components))
+
+ op.apply()
+ if impl is ref_impl:
+ for _ in range(100):
+ op.apply()
+
+ outputs[impl] = tuple(dfout[i].get().handle.copy()[view]
+ for i in xrange(dfout.nb_components))
+ if ref_impl in outputs.keys() and not impl is ref_impl:
+ reference = outputs[ref_impl]
+ output = outputs[impl]
+ for i in xrange(dfout.nb_components):
+ if npw.may_share_memory(reference[i], output[i]):
+ msg='Output and reference arrays may share the same memory.'
+ raise RuntimeError(msg)
+ if (reference[i].dtype != output[i].dtype):
+ msg='Output dtype {} differs from the reference one {}.'
+ msg=msg.format(output[i].dtype, reference[i].dtype)
+ raise RuntimeError(msg)
+ if (reference[i].shape != output[i].shape):
+ msg='Output shape {} differs from the reference one {}.'
+ msg=msg.format(output[i].shape, reference[i].shape)
+ raise RuntimeError(msg)
+ mask = npw.isfinite(output[i])
+ if not mask.all():
+ msg='\nFATAL ERROR: Output is not finite on axis {}.\n'.format(i)
+ print msg
+ npw.fancy_print(output[i], replace_values={(lambda a: npw.isfinite(a)): '.'})
+ raise RuntimeError(msg)
+ di = npw.abs(reference[i] - output[i])
+ max_di = npw.max(di)
+ neps = 1000
+ max_tol = neps*npw.finfo(dfout.dtype).eps
+ if (max_di>max_tol):
+ print
+ print 'SCALAR INPUT'
+ print _input[i]
+ print 'SCALAR REFERENCE'
+ print reference[i]
+ print 'SCALAR OUTPUT'
+ print output[i]
+ print 'ABS(REF - OUT)'
+ npw.fancy_print(di, replace_values={(lambda a: a<max_tol): '.'})
+ print
+ print reference[i].shape
+ msg='Output did not match reference output for component {}'.format(i)
+ msg+='\n > max computed distance was {}.'.format(max_di)
+ msg+='\n > max tolerence was set to {} ({} eps).'.format(max_tol, neps)
+ raise RuntimeError(msg)
+
+ S1 = dfout.integrate()
+ if not npw.all(npw.isfinite(S1)):
+ msg='Integral is not finite. Got {}.'.format(S1)
+ raise RuntimeError(msg)
+ if (npw.abs(S0-S1)>max_tol).any():
+ msg='Scalar was not conserved expected {} but got {}.'
+ msg=msg.format(to_tuple(S0, cast=float),
+ to_tuple(S1, cast=float))
+ raise RuntimeError(msg)
+ except:
+ sys.stdout.write('\bx\n\n')
+ sys.stdout.flush()
+ raise
+
+
+ def test_diffusion_1D_inplace(self):
+ self._test(dim=1, is_inplace=True)
+ def test_diffusion_2D_inplace(self):
+ self._test(dim=2, is_inplace=True)
+ def test_diffusion_3D_inplace(self):
+ self._test(dim=3, is_inplace=True)
+
+ def perform_tests(self):
+ #self.test_diffusion_1D_inplace()
+ #self.test_diffusion_2D_inplace()
+ self.test_diffusion_3D_inplace()
+
+ print
+
+
+if __name__ == '__main__':
+ TestDirectionalDiffusionOperator.setup_class(enable_extra_tests=False,
+ enable_debug_mode=True)
+
+ test = TestDirectionalDiffusionOperator()
+
+ enable_pretty_printing()
+
+ with printoptions(threshold=10000, linewidth=240,
+ nanstr='nan', infstr='inf',
+ formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
+ test.perform_tests()
+
+ TestDirectionalDiffusionOperator.teardown_class()
diff --git a/hysop/operator/tests/test_directional_diffusion.py b/hysop/operator/tests/test_directional_diffusion.py
index 04c40f0dd..02b0ff34a 100644
--- a/hysop/operator/tests/test_directional_diffusion.py
+++ b/hysop/operator/tests/test_directional_diffusion.py
@@ -1,4 +1,3 @@
-
from hysop.deps import sys
from hysop.tools.numpywrappers import npw
from hysop.testsenv import __ENABLE_LONG_TESTS__
@@ -22,19 +21,19 @@ from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator, Euler, RK2, RK
class TestDirectionalDiffusionOperator(object):
@classmethod
- def setup_class(cls,
+ def setup_class(cls,
enable_extra_tests=__ENABLE_LONG_TESTS__,
enable_debug_mode=False):
IO.set_default_path('/tmp/hysop_tests/test_directional_diffusion')
-
+
if enable_debug_mode:
cls.size_min = 8
cls.size_max = 9
else:
cls.size_min = 11
cls.size_max = 23
-
+
cls.enable_extra_tests = enable_extra_tests
cls.enable_debug_mode = enable_debug_mode
@@ -45,12 +44,12 @@ class TestDirectionalDiffusionOperator(object):
def _test(self, dim, is_inplace,
size_min=None, size_max=None):
assert dim > 0
-
+
# periodic boundaries removes one computational point
# so we add one here.
size_min = first_not_None(size_min, self.size_min) + 1
size_max = first_not_None(size_max, self.size_max) + 1
-
+
shape = tuple(npw.random.randint(low=size_min, high=size_max, size=dim).tolist())
if self.enable_extra_tests:
@@ -61,20 +60,20 @@ class TestDirectionalDiffusionOperator(object):
flt_types = (npw.float32,)
time_integrators = (Euler, RK2)
orders = (4,)
-
+
domain = Box(length=(2*npw.pi,)*dim)
for dtype in flt_types:
- Fin = Field(domain=domain, name='Fin', dtype=dtype,
+ Fin = Field(domain=domain, name='Fin', dtype=dtype,
nb_components=3, register_object=False)
- Fout = Field(domain=domain, name='Fout', dtype=dtype,
+ Fout = Field(domain=domain, name='Fout', dtype=dtype,
nb_components=3, register_object=False)
coeffs = npw.random.rand(Fin.nb_components, dim).astype(dtype)
for order in orders:
for time_integrator in time_integrators:
print
- self._test_one(time_integrator=time_integrator, order=order,
+ self._test_one(time_integrator=time_integrator, order=order,
shape=shape, dim=dim, dtype=dtype,
- is_inplace=is_inplace, domain=domain,
+ is_inplace=is_inplace, domain=domain,
Fin=Fin, Fout=Fout, coeffs=coeffs)
@staticmethod
@@ -82,14 +81,14 @@ class TestDirectionalDiffusionOperator(object):
shape = data[0].shape
if is_integer(dtype):
for d in data:
- d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+ d[...] = npw.random.random_integers(low=0, high=255, size=shape)
elif is_fp(dtype):
for d in data:
- d[...] = npw.random.random(size=d.shape)
+ d[...] = npw.random.random(size=d.shape)
else:
msg = 'Unknown dtype {}.'.format(dtype)
raise NotImplementedError(msg)
-
+
@staticmethod
def __scalar_init(data, coords, dtype):
if is_fp(dtype):
@@ -120,11 +119,11 @@ class TestDirectionalDiffusionOperator(object):
variables = { Fin : shape, Fout: shape }
msg='Out of place diffusion has not been implemented yet.'
raise NotImplementedError(msg)
-
+
implementations = DirectionalDiffusion.implementations().keys()
method = {TimeIntegrator: time_integrator, SpaceDiscretization: order}
-
+
def iter_impl(impl):
base_kwds = dict(fields=fin, coeffs=coeffs, dt=dt,
variables=variables, implementation=impl, laplacian_formulation=False,
@@ -132,16 +131,16 @@ class TestDirectionalDiffusionOperator(object):
da = DirectionalDiffusion(**base_kwds)
if (impl is Implementation.OPENCL):
for cl_env in iter_clenv():
- msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
+ msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
cl_env.device.name.strip())
- split = StrangSplitting(splitting_dim=dim,
+ split = StrangSplitting(splitting_dim=dim,
extra_kwds=dict(cl_env=cl_env),
order=StrangOrder.STRANG_SECOND_ORDER)
yield msg, da, split
else:
msg='Unknown implementation to test {}.'.format(impl)
raise NotImplementedError(msg)
-
+
# Compare to other implementations
reference_fields = {}
@@ -153,33 +152,33 @@ class TestDirectionalDiffusionOperator(object):
if (not is_inplace):
split.push_copy(dst=fin, src=fout)
split.build()
-
+
dfin = split.input_discrete_fields[fin]
dfout = split.output_discrete_fields[fout]
-
+
view = dfin.compute_slices
reference = None
-
+
sys.stdout.write('.')
sys.stdout.flush()
try:
dfout.initialize(self.__random_init, dtype=dfout.dtype)
dfin.initialize(self.__scalar_init, dtype=dfin.dtype)
S0 = dfin.integrate()
-
+
if not npw.all(npw.isfinite(S0)):
msg='Integral is not finite. Got {}.'.format(S0)
raise RuntimeError(msg)
-
- _input = tuple(dfin[i].get().handle.copy()
+
+ _input = tuple(dfin[i].get().handle.copy()
for i in xrange(dfin.nb_components))
if (reference is None):
reference = self._compute_reference(fin=dfin, dfin=_input, dt=dt(), dim=dim,
time_integrator=time_integrator, order=order, coeffs=coeffs, view=view)
split.apply()
-
+
output = tuple(dfout[i].get().handle.copy()[view]
for i in xrange(dfout.nb_components))
@@ -220,7 +219,7 @@ class TestDirectionalDiffusionOperator(object):
msg+='\n > max computed distance was {}.'.format(max_di)
msg+='\n > max tolerence was set to {} ({} eps).'.format(max_tol, neps)
raise RuntimeError(msg)
-
+
S1 = dfout.integrate()
if not npw.all(npw.isfinite(S1)):
msg='Integral is not finite. Got {}.'.format(S1)
@@ -234,11 +233,11 @@ class TestDirectionalDiffusionOperator(object):
sys.stdout.write('\bx\n\n')
sys.stdout.flush()
raise
-
+
def _compute_reference(self, fin, dfin, dt,
- time_integrator, order, coeffs, dim,
+ time_integrator, order, coeffs, dim,
view):
-
+
dfin = tuple(f.copy() for f in dfin)
csg = CenteredStencilGenerator()
@@ -249,10 +248,10 @@ class TestDirectionalDiffusionOperator(object):
dt_coeffs = (0.5,)*(2*dim)
Xin = { 'F{}'.format(i): dfin[i] for i in xrange(len(dfin)) }
views = { 'F{}'.format(i): view for i in xrange(len(dfin)) }
-
+
ghost_exchangers = tuple(fin.build_ghost_exchanger(directions=i, data=dfin)
for i in xrange(dim))
-
+
for (j, dt_coeff) in zip(directions, dt_coeffs):
ndir = (dim-j-1)
def rhs(out, X, **kwds):
@@ -261,53 +260,53 @@ class TestDirectionalDiffusionOperator(object):
fi = 'F{}'.format(i)
Fi_in = X[fi]
d2Fi_dxj = out[fi]
- d2Fi_dxj[...] = D2(a=Fi_in, out=None, axis=ndir,
+ d2Fi_dxj[...] = D2(a=Fi_in, out=None, axis=ndir,
symbols={D2.dx:fin.space_step[ndir]}) * coeff[i]
out = time_integrator(Xin=Xin, RHS=rhs, dt=(dt_coeff*dt))
for i in xrange(len(dfin)):
fi = 'F{}'.format(i)
dfin[i][...] = out[fi]
- ghost_exchangers[ndir].exchange_ghosts()
+ ghost_exchangers[ndir].exchange_ghosts()
out = tuple(f[view] for f in dfin)
return out
-
+
def test_diffusion_1D_out_of_place(self):
self._test(dim=1, is_inplace=False)
def test_diffusion_2D_out_of_place(self):
self._test(dim=2, is_inplace=False)
def test_diffusion_3D_out_of_place(self):
self._test(dim=3, is_inplace=False)
-
+
def test_diffusion_1D_inplace(self):
self._test(dim=1, is_inplace=True)
def test_diffusion_2D_inplace(self):
self._test(dim=2, is_inplace=True)
def test_diffusion_3D_inplace(self):
self._test(dim=3, is_inplace=True)
-
+
def perform_tests(self):
self.test_diffusion_1D_inplace()
- # self.test_diffusion_2D_inplace()
+ self.test_diffusion_2D_inplace()
self.test_diffusion_3D_inplace()
-
+
# self.test_diffusion_1D_out_of_place()
# self.test_diffusion_2D_out_of_place()
# self.test_diffusion_3D_out_of_place()
print
-
-
+
+
if __name__ == '__main__':
- TestDirectionalDiffusionOperator.setup_class(enable_extra_tests=False,
+ TestDirectionalDiffusionOperator.setup_class(enable_extra_tests=False,
enable_debug_mode=True)
-
+
test = TestDirectionalDiffusionOperator()
enable_pretty_printing()
-
- with printoptions(threshold=10000, linewidth=240,
- nanstr='nan', infstr='inf',
+
+ with printoptions(threshold=10000, linewidth=240,
+ nanstr='nan', infstr='inf',
formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
test.perform_tests()
diff --git a/hysop/operators.py b/hysop/operators.py
index 27efe79d4..6ace01759 100644
--- a/hysop/operators.py
+++ b/hysop/operators.py
@@ -7,6 +7,7 @@ from hysop.operators import DirectionalAdvection
from hysop.operator.poisson import Poisson
from hysop.operator.poisson_rotational import PoissonRotational
+from hysop.operator.diffusion import Diffusion # FFTW diffusion
from hysop.operator.redistribute import Redistribute
from hysop.operator.analytic import AnalyticField
diff --git a/hysop/parameters/default_parameters.py b/hysop/parameters/default_parameters.py
index 88fc81cf2..b6776362a 100644
--- a/hysop/parameters/default_parameters.py
+++ b/hysop/parameters/default_parameters.py
@@ -4,7 +4,8 @@ from hysop.tools.sympy_utils import greak, Greak, subscripts
from hysop.parameters.tensor_parameter import TensorParameter
from hysop.parameters.scalar_parameter import ScalarParameter
-def TimeParameters(dtype, **kwds):
+def TimeParameters(dtype=None, **kwds):
+ dtype = first_not_None(dtype, HYSOP_REAL)
t = ScalarParameter('t', dtype=dtype, **kwds)
dt = ScalarParameter('dt', dtype=dtype, **kwds)
return (t, dt)
@@ -28,9 +29,9 @@ def VolumicIntegrationParameter(name=None, pretty_name=None, field=None, dtype=N
name += '_v'
pretty_name += u'\u1d65'.encode('utf-8')
if (field.nb_components>1):
- return TensorParameter(name=name, pretty_name=pretty_name,
- shape=(field.nb_components,),
+ return TensorParameter(name=name, pretty_name=pretty_name,
+ shape=(field.nb_components,),
dtype=dtype, **kwds)
else:
- return ScalarParameter(name=name, pretty_name=pretty_name,
+ return ScalarParameter(name=name, pretty_name=pretty_name,
dtype=dtype, **kwds)
--
GitLab