Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • particle_methods/hysop
1 result
Show changes
Commits on Source (7)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 872 additions and 371 deletions
ci/scripts/test.sh
View file @ a5549c54
......@@ -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
......
hysop/__init__.py.in
View file @ a5549c54
......@@ -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()
hysop/backend/host/fortran/operator/diffusion.py
View file @ a5549c54
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
hysop/backend/host/fortran/operator/fortran_fftw.py
View file @ a5549c54
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
......
hysop/backend/host/fortran/operator/poisson.py 0 → 100644
View file @ a5549c54
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.constants import FieldProjection
from hysop.backend.host.fortran.operator.fortran_fftw import fftw2py, FortranFFTWOperator
from hysop.core.graph.graph import op_apply
import numpy as np
class PoissonFFTW(FortranFFTWOperator):
def __init__(self, Fin, Fout, variables,
extra_input_kwds=None, **kwds):
"""Operator to solve Poisson equation using FFTW in Fortran.
Solves:
Laplacian(Fout) = Fin
Parameters
----------
Fout: Field
Input continuous field (rhs).
Fin: Field
Output continuous field (lhs), possibly inplace.
variables: dict
Dictionary of fields as keys and topology descriptors as values.
kwds: dict, optional
Base class arguments.
"""
check_instance(Fin, Field)
check_instance(Fout, Field)
check_instance(variables, dict, keys=Field,
values=CartesianTopologyDescriptors)
assert Fin.nb_components == Fout.nb_components
assert Fin.nb_components == 1, \
"Poisson operator is implemented only for scalar fields"
assert Fin.domain is Fout.domain, \
'only one domain is supported'
assert variables[Fin] is variables[Fout], \
'only one topology is supported'
input_fields = {Fin: variables[Fin]}
output_fields = {Fout: variables[Fout]}
super(PoissonFFTW, self).__init__(input_fields=input_fields,
output_fields=output_fields,
**kwds)
self.Fin = Fin
self.Fout = Fout
def initialize(self, **kwds):
super(PoissonFFTW, self).initialize(**kwds)
dim = self.dim
if (dim == 2):
self._solve = fftw2py.solve_laplace_2d
elif (dim == 3):
self._solve = fftw2py.solve_laplace_3d
@debug
def discretize(self):
if self.discretized:
return
super(PoissonFFTW, self).discretize()
self.dFin = self.Fin.discrete_fields[self.topology]
self.dFout = self.Fout.discrete_fields[self.topology]
assert (self.dFin.ghosts == self.dFout.ghosts).all(), \
"Input and output fields must have the same ghosts."
@op_apply
def apply(self, **kargs):
super(PoissonFFTW, self).apply(**kargs)
gh = self.input_fields[self.Fin].ghosts
changeLayout, di, do = self._initialize_mem_layout(self.dFin, self.dFout)
# Vectors are given in ZYX layout to Fortran
do = self._solve(di[0], gh)
if changeLayout:
self.dFout.data[0][...] = -np.ascontiguousarray(do)
self.dFout.exchange_ghosts()
def _initialize_mem_layout(self, fi, fo):
"""Convert C-contiguous numpy arrays to Fortran-ordered data if needed."""
changeLayout = False
if all([__.flags['C_CONTIGUOUS']
for _ in (fi, fo) for __ in _.data ]):
changeLayout = True
print "CHange layout"
new_i = tuple(np.asfortranarray(_) for _ in fi.buffers)
new_o = tuple(np.asfortranarray(_) for _ in fo.buffers)
elif all([__.flags['F_CONTIGUOUS']
for _ in (fi, fo) for __ in _.data ]):
new_i, new_o = fi, fo
else:
raise RuntimeError("Memory layout comatibility issue")
return changeLayout, new_i, new_o
hysop/backend/host/fortran/operator/poisson_rotational.py
View file @ a5549c54
from hysop.tools.types import check_instance, InstanceOf
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
......@@ -7,13 +6,14 @@ from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.constants import FieldProjection
from hysop.backend.host.fortran.operator.fortran_fftw import fftw2py, FortranFFTWOperator
from hysop.core.graph.graph import op_apply
import numpy as np
class PoissonRotationalFFTW(FortranFFTWOperator):
def __init__(self, velocity, vorticity, variables,
projection=FieldProjection.NONE, **kwds):
def __init__(self, velocity, vorticity, variables,
projection=FieldProjection.NONE, **kwds):
"""PoissonRotational operator to solve incompressible flows using FFTW in Fortran.
Parameters
----------
velocity : :class:`~hysop.fields.continuous_field.Field
......@@ -24,13 +24,13 @@ class PoissonRotationalFFTW(FortranFFTWOperator):
dictionary of fields as keys and topologies as values.
projection: hysop.constants.FieldProjection or positive integer, optional
Project vorticity such that resolved velocity is divergence free (for 3D fields).
When active, projection is done prior to every solve, unless projection is
When active, projection is done prior to every solve, unless projection is
an integer in which case it is done every n applies.
This parameter is ignored for 2D fields and defaults to no projection.
kwds :
kwds :
base class parameters.
"""
check_instance(velocity, Field)
check_instance(vorticity, Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
......@@ -48,14 +48,14 @@ class PoissonRotationalFFTW(FortranFFTWOperator):
else:
assert (velocity.dim==3) and (vorticity.dim==3), 'Projection only available in 3D.'
output_fields[vorticity] = wtopology
super(PoissonRotationalFFTW,self).__init__(input_fields=input_fields, output_fields=output_fields,
**kwds)
self.velocity = velocity
self.vorticity = vorticity
self.projection = projection
def initialize(self, **kwds):
super(PoissonRotationalFFTW,self).initialize(**kwds)
dim = self.dim
......@@ -63,22 +63,22 @@ class PoissonRotationalFFTW(FortranFFTWOperator):
self._solve = self._solve_2d
elif (dim==3):
self._solve = self._solve_3d
if (dim!=3) and (self.projection!=FieldProjection.NONE):
raise ValueError('Velocity reprojection only available in 3D.')
projection = self.projection
if (projection == FieldProjection.NONE):
self._do_project = lambda simu: False
elif (projection == FieldProjection.EVERY_STEP):
self._do_project = lambda simu: True
else: # projection is an integer representing frenquency
else: # projection is an integer representing frenquency
freq = projection
assert freq>=1
self._do_project = lambda simu: ((simu.current_iteration % freq)==0)
self.projection = projection
@debug
def discretize(self):
......@@ -92,24 +92,51 @@ class PoissonRotationalFFTW(FortranFFTWOperator):
def apply(self, simulation=None, **kargs):
super(PoissonRotationalFFTW, self).apply(simulation=simulation, **kargs)
if self._do_project(simulation):
self._project(simulation)
self._solve(simulation)
def _solve_2d(self, simu=None):
self._project()
self._solve()
def _initialize_mem_layout(self):
"""Convert C-contiguous numpy arrays to Fortran-ordered data if needed."""
changeLayout = False
if all([__.flags['C_CONTIGUOUS']
for _ in (self.dvorticity,self.dvelocity) for __ in _.data ]):
changeLayout = True
print "CHange layout"
dw = tuple(np.asfortranarray(_) for _ in self.dvorticity.buffers)
dv = tuple(np.asfortranarray(_) for _ in self.dvelocity.buffers)
elif all([__.flags['F_CONTIGUOUS']
for _ in (self.dvorticity,self.dvelocity) for __ in _.data ]):
dw, dv = self.dvorticity, self.dvelocity
else:
raise RuntimeError("Memory layout comatibility issue")
return changeLayout, dv, dw
def _finalize_mem_layout(self, dv):
"""Get back with C-contiguous data.
We enforce a sign change for data because of the change
of direction ordering ZYX."""
for dvi, dvi_f in zip(self.dvelocity.buffers, dv[::-1]):
dvi[...] = -np.ascontiguousarray(dvi_f)
def _solve_2d(self):
"""
Solve 2D poisson problem, no projection, no correction.
"""
assert self.dim==2
ghosts_w = self.input_fields[self.vorticity].ghosts
ghosts_v = self.output_fields[self.velocity].ghosts
self.dvelocity_data =\
fftw2py.solve_poisson_2d(self.dvorticity.data[0],
self.dvelocity.data[0],
self.dvelocity.data[1],
ghosts_w, ghosts_v)
changeLayout, dv, dw = self._initialize_mem_layout()
# Vectors are given in ZYX layout to Fortran
dv = fftw2py.solve_poisson_2d(dw[0],
dv[1], dv[0],
ghosts_w, ghosts_v)
if changeLayout:
self._finalize_mem_layout(dv)
self.dvelocity.exchange_ghosts()
def _solve_3d(self,simu=None):
def _solve_3d(self):
"""
Solve 3D poisson problem, no projection, no correction
"""
......@@ -117,26 +144,27 @@ class PoissonRotationalFFTW(FortranFFTWOperator):
assert self.dim==3
ghosts_w = self.input_fields[self.vorticity].ghosts
ghosts_v = self.output_fields[self.velocity].ghosts
self.dvelocity.data =\
fftw2py.solve_poisson_3d(self.dvorticity.data[0],
self.dvorticity.data[1],
self.dvorticity.data[2],
self.dvelocity.data[0],
self.dvelocity.data[1],
self.dvelocity.data[2],
ghosts_w, ghosts_v)
changeLayout, dv, dw = self._initialize_mem_layout()
# Vectors are given in ZYX layout to Fortran
dv = fftw2py.solve_poisson_3d(dw[2], dw[1], dw[0],
dv[2], dv[1], dv[0],
ghosts_w, ghosts_v)
if changeLayout:
self._finalize_mem_layout(dv)
self.dvelocity.exchange_ghosts()
def _project(self, simulation):
def _project(self):
"""
Apply projection on vorticity such that the
resolved velocity will have div(U)=0.
"""
assert self.dim==3
ghosts_w = self.output_fields[self.vorticity].ghosts
changeLayout, dv, dw = self._initialize_mem_layout()
# Vectors are given in ZYX layout to Fortran
self.dvorticity.data =\
fftw2py.projection_om_3d(self.dvorticity.data[0],
self.dvorticity.data[1],
self.dvorticity.data[2], ghosts_w)
fftw2py.projection_om_3d(dw[2], dw[1], dw[0], ghosts_w)
if changeLayout:
self._finalize_mem_layout(dv)
self.dvorticity.exchange_ghosts()
hysop/backend/host/host_array.py
View file @ a5549c54
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)
hysop/core/graph/computational_graph.py
View file @ a5549c54
This diff is collapsed.
hysop/core/graph/computational_node_frontend.py
View file @ a5549c54
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.
......
hysop/numerics/fftw_f/fft2d.f90
View file @ a5549c54
......@@ -29,7 +29,7 @@ module fft2d
public :: init_c2c_2d,init_r2c_2d, r2c_scalar_2d, c2c_2d,c2r_2d,c2r_scalar_2d,&
r2c_2d,cleanFFTW_2d, filter_poisson_2d, filter_curl_2d, getParamatersTopologyFFTW2d, &
filter_diffusion_2d, init_r2c_2dBIS
filter_diffusion_2d, init_r2c_2dBIS, filter_laplace_2d
!> plan for fftw "c2c" forward or r2c transform
......@@ -289,7 +289,7 @@ contains
do i = 1, fft_resolution(c_X)
ig = i + ghosts(c_X)
rdatain1(i, j) = input_x(ig,jg)
rdatain2(i, j) = input_y(ig,jg)
rdatain2(i, j) = input_y(ig,jg)
end do
end do
......@@ -484,6 +484,20 @@ contains
end subroutine filter_poisson_2d
subroutine filter_laplace_2d()
integer(C_INTPTR_T) :: i, j
complex(C_DOUBLE_COMPLEX) :: coeff
do i = 1,local_resolution(c_X)
do j = 1, fft_resolution(c_Y)
coeff = ONE/(kx(i)**2+ky(j)**2)
dataout1(j,i) = coeff*dataout1(j,i)
end do
end do
dataout1(1,1) = 0
end subroutine filter_laplace_2d
subroutine filter_diffusion_2d(nudt)
real(C_DOUBLE), intent(in) :: nudt
......
hysop/numerics/fftw_f/fft3d.f90
View file @ a5549c54
......@@ -29,7 +29,8 @@ module fft3d
getParamatersTopologyFFTW3d,filter_poisson_3d,filter_curl_diffusion_3d, &
init_r2c_3d_many, r2c_3d_many, c2r_3d_many, filter_diffusion_3d_many,&
filter_poisson_3d_many, filter_diffusion_3d, filter_curl_3d, filter_projection_om_3d,&
filter_multires_om_3d, filter_pressure_3d, r2c_3d_scal, filter_spectrum_3d
filter_multires_om_3d, filter_pressure_3d, r2c_3d_scal, filter_spectrum_3d, &
filter_laplace_3d
!> plan for fftw "c2c" forward or r2c transform
type(C_PTR) :: plan_forward1, plan_forward2, plan_forward3
......@@ -315,7 +316,7 @@ contains
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
plan_backward1 = fftw_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
call computeKy(lengths(c_Y))
call computeKz(lengths(c_Z))
......@@ -326,7 +327,7 @@ contains
is3DUpToDate = .true.
end subroutine init_r2c_scalar_3d
!> forward transform - The result is saved in local buffers
!! @param[in] omega_x 3d scalar field, x-component of the input vector field
!! @param[in] omega_y 3d scalar field, y-component of the input vector field
......@@ -749,6 +750,21 @@ contains
end subroutine filter_diffusion_3d
subroutine filter_laplace_3d()
integer(C_INTPTR_T) :: i, j, k
complex(C_DOUBLE_COMPLEX) :: coeff
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
do i = 1, local_resolution(c_X)
coeff = one/((kx(i)**2+ky(j)**2+kz(k)**2))
dataout1(i,k,j) = coeff*dataout1(i,k,j)
end do
end do
end do
dataout1(1,1,1) = 0
end subroutine filter_laplace_3d
!> Solve curl problem in the Fourier space :
!! \f{eqnarray*} \omega &=& \nabla \times v
subroutine filter_curl_3d()
......
hysop/numerics/fftw_f/fftw2py.f90
View file @ a5549c54
......@@ -15,7 +15,7 @@ module fftw2py
implicit none
contains
!> Initialisation of fftw context : create plans and memory buffers
!! @param[in] resolution global resolution of the discrete domain
!! @param[in] lengths width of each side of the domain
......@@ -89,18 +89,19 @@ contains
! Duplicate comm into client_data::main_comm (used later in fft2d and fft3d)
call mpi_comm_dup(comm, main_comm, ierr)
!print*, "Init fftw/poisson solver for a 3d problem"
call init_r2c_scalar_3d(resolution,lengths)
call getParamatersTopologyFFTW3d(datashape,offset)
end subroutine init_fftw_solver_scalar
!> Free memory allocated for fftw-related objects (plans and buffers)
subroutine clean_fftw_solver(dim)
integer, intent(in) :: dim
print *, "Clean FFTW fortran"
if(dim == 2) then
call cleanFFTW_2d()
else
......@@ -116,7 +117,7 @@ contains
real(wp),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
integer(kind=ip), dimension(2), intent(in) :: ghosts_vort, ghosts_velo
!f2py intent(in,out) :: velocity_x,velocity_y
call r2c_scalar_2d(omega, ghosts_vort)
call filter_poisson_2d()
......@@ -126,6 +127,34 @@ contains
end subroutine solve_poisson_2d
!> Solve
!! \f[ \nabla (u) = f \f]
!! u being a 2D scalar field and f a 2D scalar field.
subroutine solve_laplace_2d(f, u, ghosts)
real(wp),dimension(:,:),intent(in):: f
real(wp),dimension(size(f,1),size(f,2)),intent(out) :: u
integer(kind=ip), dimension(2), intent(in) :: ghosts
call r2c_scalar_2d(f, ghosts)
call filter_laplace_2d()
call c2r_scalar_2d(u, ghosts)
end subroutine solve_laplace_2d
!> Solve
!! \f[ \nabla (u) = f \f]
!! u being a 3D scalar field and f a 3D scalar field.
subroutine solve_laplace_3d(f, u, ghosts)
real(wp),dimension(:,:,:),intent(in):: f
real(wp),dimension(size(f,1),size(f,2),size(f,3)),intent(out) :: u
integer(kind=ip), dimension(3), intent(in) :: ghosts
call r2c_scalar_3d(f, ghosts)
call filter_laplace_3d()
call c2r_scalar_3d(u, ghosts)
end subroutine solve_laplace_3d
!> Solve
!! \f{eqnarray*} \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
!! omega being a 2D scalar field.
......
hysop/numerics/fftw_f/fftw2py.pyf
View file @ a5549c54
......@@ -35,6 +35,16 @@ module fftw2py ! in fftw2py.f90
integer(kind=ip) dimension(2),intent(in) :: ghosts_vort
integer(kind=ip) dimension(2),intent(in) :: ghosts_velo
end subroutine solve_poisson_2d
subroutine solve_laplace_2d(f, u, ghosts)
real(kind=wp),dimension(:,:),intent(in):: f
real(kind=wp),dimension(size(f,1),size(f,2)),intent(out) :: u
integer(kind=ip), dimension(2), intent(in) :: ghosts
end subroutine solve_laplace_2d
subroutine solve_laplace_3d(f, u, ghosts)
real(kind=wp),dimension(:,:,:),intent(in):: f
real(kind=wp),dimension(size(f,1),size(f,2),size(f,3)),intent(out) :: u
integer(kind=ip), dimension(3), intent(in) :: ghosts
end subroutine solve_laplace_3d
subroutine solve_diffusion_2d(nudt,omega,ghosts_vort) ! in fftw2py.f90:fftw2py
real(kind=wp) intent(in) :: nudt
real(kind=wp) dimension(:,:),intent(in,out) :: omega
......
hysop/operator/diffusion.py
View file @ a5549c54
"""
@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)
hysop/operator/poisson.py
View file @ a5549c54
"""
@file poisson.py
Poisson solver frontend.
......@@ -13,37 +12,41 @@ from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeF
from hysop.backend.host.python.operator.poisson import PythonPoisson
from hysop.backend.device.opencl.operator.poisson import OpenClPoisson
from hysop.backend.host.fortran.operator.poisson import PoissonFFTW
class Poisson(ComputationalGraphNodeFrontend):
"""
Interface the poisson solver.
Available implementations are:
Available implementations are:
*PYTHON (numpy local solver)
*OPENCL
*FORTRAN (FFTW solver)
"""
__implementations = {
Implementation.PYTHON: PythonPoisson,
Implementation.OPENCL: OpenClPoisson
Implementation.PYTHON: PythonPoisson,
Implementation.OPENCL: OpenClPoisson,
Implementation.FORTRAN: PoissonFFTW
}
@classmethod
def implementations(cls):
return cls.__implementations
@classmethod
def default_implementation(cls):
return Implementation.PYTHON
@debug
def __init__(self, Fin, Fout, variables,
def __init__(self, Fin, Fout, variables,
implementation=None, base_kwds=None, **kwds):
"""
Initialize a n-dimensional Poisson operator frontend.
Solves:
Laplacian(Fout) = Fin
Parameters
----------
Fout: Field
......@@ -59,7 +62,7 @@ class Poisson(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__.
"""
base_kwds = base_kwds or dict()
......@@ -68,6 +71,6 @@ class Poisson(ComputationalGraphNodeFrontend):
check_instance(Fin, Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str)
super(Poisson, self).__init__(Fin=Fin, Fout=Fout,
super(Poisson, self).__init__(Fin=Fin, Fout=Fout,
variables=variables, base_kwds=base_kwds, implementation=implementation, **kwds)
hysop/operator/poisson_rotational.py
View file @ a5549c54
"""
@file poisson.py
PoissonRotational solver frontend.
......@@ -11,40 +10,42 @@ from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
# from hysop.backend.host.fortran.operator.poisson_rotational import PoissonRotationalFFTW
from hysop.backend.host.fortran.operator.poisson_rotational import PoissonRotationalFFTW
from hysop.backend.host.python.operator.poisson_rotational import PythonPoissonRotational
from hysop.backend.device.opencl.operator.poisson_rotational import OpenClPoissonRotational
class PoissonRotational(ComputationalGraphNodeFrontend):
"""
Interface the poisson solver.
Available implementations are:
Available implementations are:
*FORTRAN (fftw based solver)
*PYTHON (gpyfft based solver)
*OPENCL (clfft based solver)
"""
__implementations = {
Implementation.PYTHON: PythonPoissonRotational,
Implementation.OPENCL: OpenClPoissonRotational,
#Implementation.FORTRAN: PoissonRotationalFFTW
Implementation.PYTHON: PythonPoissonRotational,
Implementation.OPENCL: OpenClPoissonRotational,
Implementation.FORTRAN: PoissonRotationalFFTW
}
@classmethod
def implementations(cls):
return cls.__implementations
@classmethod
def default_implementation(cls):
return Implementation.PYTHON
@debug
def __init__(self, velocity, vorticity, variables,
def __init__(self, velocity, vorticity, variables,
implementation=None, base_kwds=None, **kwds):
"""
Initialize a PoissonRotational operator frontend for 2D or 3D streamfunction-vorticity formulations.
in = W (vorticity)
out = U (velocity)
About dimensions:
- if velocity is a 2D vector field, vorticity should have only one component Wx.
- if velocity is a 3D vector field, vortcitiy should have three components (Wx,Wy,Wz).
......@@ -53,11 +54,11 @@ class PoissonRotational(ComputationalGraphNodeFrontend):
laplacian(psi) = -W
Ux = +dpsi/dy
Uy = -dpsi/dx
In 3D:
laplacian(psi) = -W
U = rot(psi)
Parameters
----------
velocity: Field
......@@ -73,12 +74,12 @@ class PoissonRotational(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 PoissonRotational operator implementation should at least support
A PoissonRotational operator implementation should at least support
the following __init__ attributes: velocity, vorticity, variables
"""
base_kwds = base_kwds or dict()
......@@ -87,11 +88,11 @@ class PoissonRotational(ComputationalGraphNodeFrontend):
check_instance(vorticity, Field)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str)
dim = velocity.domain.dim
vcomp = velocity.nb_components
wcomp = vorticity.nb_components
if (velocity.domain != vorticity.domain):
msg = 'Velocity and vorticity do not share the same domain.'
raise ValueError(msg)
......@@ -109,5 +110,5 @@ class PoissonRotational(ComputationalGraphNodeFrontend):
msg='Vorticity component mistmach, got {} components but expected 3.'.format(wcomp)
raise RuntimeError(msg)
super(PoissonRotational, self).__init__(velocity=velocity, vorticity=vorticity,
super(PoissonRotational, self).__init__(velocity=velocity, vorticity=vorticity,
variables=variables, base_kwds=base_kwds, implementation=implementation, **kwds)
hysop/operator/tests/test_diffusion.py 0 → 100644
View file @ a5549c54
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()
hysop/operator/tests/test_directional_diffusion.py
View file @ a5549c54
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()
......
hysop/operator/tests/test_poisson.py
View file @ a5549c54
import random, primefac
from hysop.deps import it, sm, random
from hysop.constants import HYSOP_REAL
......@@ -16,19 +15,19 @@ from hysop import Field, Box
class TestPoissonOperator(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_poisson')
if enable_debug_mode:
cls.size_min = 15
cls.size_max = 16
else:
cls.size_min = 23
cls.size_max = 87
cls.enable_extra_tests = enable_extra_tests
cls.enable_debug_mode = enable_debug_mode
......@@ -41,9 +40,9 @@ class TestPoissonOperator(object):
from hysop.symbolic.frame import SymbolicFrame
from hysop.symbolic.field import curl, laplacian
enable_pretty_printing()
#at this moment we just test a periodic solver
analytic_solutions = {}
analytic_functions = {}
for dim in (1,2,3,4):
......@@ -51,7 +50,7 @@ class TestPoissonOperator(object):
def gen_psi():
kis = tuple(random.randint(1,5) for _ in xrange(dim))
qis = tuple(npw.random.rand(dim).round(decimals=3).tolist())
basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi))
basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi))
for (ki,qi,xi) in zip(kis, qis, frame.coords) )
psi = sm.Integer(1)
for i in xrange(dim):
......@@ -73,14 +72,14 @@ class TestPoissonOperator(object):
def teardown_class(cls):
pass
def _test(self, dim, dtype,
size_min=None, size_max=None):
enable_extra_tests = self.enable_extra_tests
size_min = first_not_None(size_min, self.size_min)
size_max = first_not_None(size_max, self.size_max)
valid_factors = {2,3,5,7,11,13}
factors = {1}
while (factors-valid_factors):
......@@ -88,14 +87,14 @@ class TestPoissonOperator(object):
shape = tuple(npw.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
for Si in shape:
factors.update( set(primefac.primefac(int(Si-1))) )
domain = Box(length=(2*npw.pi,)*dim)
Psi = Field(domain=domain, name='Psi', dtype=dtype,
Psi = Field(domain=domain, name='Psi', dtype=dtype,
nb_components=1, register_object=False)
W = Field(domain=domain, name='W', dtype=dtype,
nb_components=1, register_object=False)
self._test_one(shape=shape, dim=dim, dtype=dtype,
self._test_one(shape=shape, dim=dim, dtype=dtype,
domain=domain, Psi=Psi, W=W)
print
......@@ -114,8 +113,8 @@ class TestPoissonOperator(object):
else:
msg = 'Unknown dtype {}.'.format(dtype)
raise NotImplementedError(msg)
def _test_one(self, shape, dim, dtype,
def _test_one(self, shape, dim, dtype,
domain, Psi, W):
print '\nTesting periodic {}D Poisson: dtype={} shape={}'.format(
......@@ -129,16 +128,20 @@ class TestPoissonOperator(object):
print ' >Testing all implementations:'
implementations = Poisson.implementations()
# Compute reference solution
variables = { Psi:shape, W:shape }
def iter_impl(impl):
base_kwds = dict(Fout=Psi, Fin=W, variables=variables,
implementation=impl,
implementation=impl,
name='poisson_{}'.format(str(impl).lower()))
if impl is Implementation.PYTHON:
if impl is Implementation.FORTRAN:
msg=' *Fortran FFTW: '
print msg,
yield Poisson(**base_kwds)
elif impl is Implementation.PYTHON:
msg=' *Python FFTW: '
print msg,
yield Poisson(**base_kwds)
......@@ -146,14 +149,14 @@ class TestPoissonOperator(object):
msg=' *OpenCl CLFFT: '
print msg
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())
print msg,
yield Poisson(cl_env=cl_env, **base_kwds)
else:
msg='Unknown implementation to test {}.'.format(impl)
raise NotImplementedError(msg)
# Compare to analytic solution
Psiref = None
Wref = None
......@@ -161,27 +164,33 @@ class TestPoissonOperator(object):
if (dim>3) and (impl is Implementation.OPENCL):
print ' *OpenCl CLFFT: NO SUPPORT'
continue
if (impl is Implementation.FORTRAN):
if ((dim<2) or (dim>3) or (not dtype is HYSOP_REAL)):
print ' *Fortran FFTW: NO SUPPORT'
continue
for op in iter_impl(impl):
op = op.build()
dw = op.input_discrete_fields[W]
dpsi = op.output_discrete_fields[Psi]
dw.initialize(self.__analytic_init, dtype=dtype,
dw.initialize(self.__analytic_init, dtype=dtype,
fns=self.analytic_functions[dim]['Ws'])
if (Psiref is None):
dpsi.initialize(self.__analytic_init, dtype=dtype,
dpsi.initialize(self.__analytic_init, dtype=dtype,
fns=self.analytic_functions[dim]['Psis'])
Wref = tuple( data.get().handle.copy() for data in dw.data )
Psiref = tuple( data.get().handle.copy() for data in dpsi.data )
dpsi.initialize(self.__random_init, dtype=dtype)
op.apply()
Wout = tuple( data.get().handle.copy() for data in dw.data )
Psiout = tuple( data.get().handle.copy() for data in dpsi.data )
self._check_output(impl, op, Wref, Psiref, Wout, Psiout)
if (impl is Implementation.FORTRAN):
op.finalize(clean_fftw_solver=True)
print
@classmethod
def _check_output(cls, impl, op, Wref, Psiref, Wout, Psiout):
check_instance(Wref, tuple, values=npw.ndarray)
......@@ -202,14 +211,14 @@ class TestPoissonOperator(object):
print
msg = msg0.format(iname)
raise ValueError(msg)
for (out_buffers, ref_buffers, name) in zip((Wout, Psiout),
for (out_buffers, ref_buffers, name) in zip((Wout, Psiout),
(Wref, Psiref), ('W', 'Psi')):
for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
iname = '{}{}'.format(name,i)
assert fout.dtype == fref.dtype, iname
assert fout.shape == fref.shape, iname
eps = npw.finfo(fout.dtype).eps
dist = npw.abs(fout-fref)
dinf = npw.max(dist)
......@@ -219,7 +228,7 @@ class TestPoissonOperator(object):
continue
has_nan = npw.any(npw.isnan(fout))
has_inf = npw.any(npw.isinf(fout))
print
print
print 'Test output comparisson for {} failed for component {}:'.format(name, i)
......@@ -253,10 +262,10 @@ class TestPoissonOperator(object):
print w
print
print
msg = 'Test failed for {} on component {} for implementation {}.'
msg=msg.format(name, i, impl)
raise RuntimeError(msg)
raise RuntimeError(msg)
def test_1d_float32(self):
......@@ -267,7 +276,7 @@ class TestPoissonOperator(object):
self._test(dim=3, dtype=npw.float32)
def test_4d_float32(self):
self._test(dim=4, dtype=npw.float32)
def test_1d_float64(self):
self._test(dim=1, dtype=npw.float64)
def test_2d_float64(self):
......@@ -288,13 +297,13 @@ class TestPoissonOperator(object):
self.test_2d_float64()
self.test_3d_float64()
self.test_4d_float64()
if __name__ == '__main__':
TestPoissonOperator.setup_class(enable_extra_tests=False,
TestPoissonOperator.setup_class(enable_extra_tests=False,
enable_debug_mode=False)
test = TestPoissonOperator()
with test_context():
test.perform_tests()
......
hysop/operator/tests/test_poisson_rotational.py
View file @ a5549c54
# coding: utf-8
import random, primefac
from hysop.deps import it, sm, random
......@@ -16,19 +17,19 @@ from hysop import Field, Box
class TestPoissonRotationalOperator(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_poisson_rotational')
if enable_debug_mode:
cls.size_min = 5
cls.size_max = 5
else:
cls.size_min = 53
cls.size_max = 87
cls.enable_extra_tests = enable_extra_tests
cls.enable_debug_mode = enable_debug_mode
......@@ -41,7 +42,7 @@ class TestPoissonRotationalOperator(object):
from hysop.symbolic.frame import SymbolicFrame
from hysop.symbolic.field import curl, laplacian
enable_pretty_printing()
#at this moment we just test a periodic solver
frame = SymbolicFrame(dim=3)
def gen_psi():
......@@ -53,11 +54,11 @@ class TestPoissonRotationalOperator(object):
for i in xrange(dim):
psi *= random.choice(basis[i])
return psi
analytic_solutions = {}
analytic_functions = {}
for dim in (2,3):
psis = npw.asarray([gen_psi() for _ in xrange(3)],
psis = npw.asarray([gen_psi() for _ in xrange(3)],
dtype=object).view(TensorBase)
if (dim==2):
psis[:2] = 0
......@@ -84,7 +85,7 @@ class TestPoissonRotationalOperator(object):
def teardown_class(cls):
pass
def _test(self, dim, dtype,
size_min=None, size_max=None):
enable_extra_tests = self.enable_extra_tests
......@@ -99,15 +100,15 @@ class TestPoissonRotationalOperator(object):
shape = tuple(npw.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
for Si in shape:
factors.update( set(primefac.primefac(int(Si-1))) )
domain = Box(length=(2*npw.pi,)*dim)
U = Field(domain=domain, name='U', dtype=dtype,
U = Field(domain=domain, name='U', dtype=dtype,
nb_components=dim, register_object=False)
W = Field(domain=domain, name='W', dtype=dtype,
nb_components=(3 if (dim==3) else 1), register_object=False)
self._test_one(shape=shape, dim=dim, dtype=dtype,
domain=domain, U=U, W=W)
self._test_one(shape=shape, dim=dim, dtype=dtype,
domain=domain, U=U, W=W)
print
@classmethod
......@@ -126,7 +127,7 @@ class TestPoissonRotationalOperator(object):
msg = 'Unknown dtype {}.'.format(dtype)
raise NotImplementedError(msg)
def _test_one(self, shape, dim, dtype,
def _test_one(self, shape, dim, dtype,
domain, U, W):
print '\nTesting periodic {}D PoissonRotational: dtype={} shape={}'.format(
......@@ -140,14 +141,14 @@ class TestPoissonRotationalOperator(object):
print ' >Testing all implementations:'
implementations = PoissonRotational.implementations()
# Compute reference solution
variables = { U:shape, W:shape }
def iter_impl(impl):
base_kwds = dict(velocity=U, vorticity=W, variables=variables,
implementation=impl,
implementation=impl,
name='poisson_{}'.format(str(impl).lower()))
if impl is Implementation.FORTRAN:
msg=' *Fortran FFTW: '
......@@ -161,18 +162,21 @@ class TestPoissonRotationalOperator(object):
msg=' *OpenCl CLFFT: '
print msg
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())
print msg,
yield PoissonRotational(cl_env=cl_env, projection=0, **base_kwds)
else:
msg='Unknown implementation to test {}.'.format(impl)
raise NotImplementedError(msg)
# Compare to analytic solution
Uref = None
Wref = None
for impl in implementations:
if (impl is Implementation.FORTRAN) and not dtype is HYSOP_REAL:
print ' *FORTRAN: NO SUPPORT for ' + str(dtype)
continue
for (i,op) in enumerate(iter_impl(impl)):
from hysop.tools.debug_dumper import DebugDumper
name='{}_{}'.format(impl, i)
......@@ -181,24 +185,24 @@ class TestPoissonRotationalOperator(object):
op = op.build()
dw = op.input_discrete_fields[W]
du = op.output_discrete_fields[U]
dw.initialize(self.__analytic_init, dtype=dtype,
dw.initialize(self.__analytic_init, dtype=dtype,
fns=self.analytic_functions[dim]['Ws'])
if (Uref is None):
du.initialize(self.__analytic_init, dtype=dtype,
du.initialize(self.__analytic_init, dtype=dtype,
fns=self.analytic_functions[dim]['Us'])
Wref = tuple( data.get().handle.copy() for data in dw.data )
Uref = tuple( data.get().handle.copy() for data in du.data )
du.initialize(self.__random_init, dtype=dtype)
op.apply(simulation=None, debug_dumper=dbg)
Wout = tuple( data.get().handle.copy() for data in dw.data )
Uout = tuple( data.get().handle.copy() for data in du.data )
# dbg(1, 1.0, 'U', Uout)
self._check_output(impl, op, Wref, Uref, Wout, Uout)
print
@classmethod
def _check_output(cls, impl, op, Wref, Uref, Wout, Uout):
check_instance(Wref, tuple, values=npw.ndarray)
......@@ -219,13 +223,13 @@ class TestPoissonRotationalOperator(object):
print
msg = msg0.format(iname)
raise ValueError(msg)
for (out_buffers, ref_buffers, name) in zip((Wout, Uout), (Wref, Uref), ('W', 'U')):
for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
iname = '{}{}'.format(name,i)
assert fout.dtype == fref.dtype, iname
assert fout.shape == fref.shape, iname
eps = npw.finfo(fout.dtype).eps
dist = npw.abs(fout-fref)
dinf = npw.max(dist)
......@@ -235,7 +239,7 @@ class TestPoissonRotationalOperator(object):
continue
has_nan = npw.any(npw.isnan(fout))
has_inf = npw.any(npw.isinf(fout))
print
print
print 'Test output comparisson for {} failed for component {}:'.format(name, i)
......@@ -269,10 +273,10 @@ class TestPoissonRotationalOperator(object):
print w
print
print
msg = 'Test failed for {} on component {} for implementation {}.'
msg = msg.format(name, i, impl)
raise RuntimeError(msg)
raise RuntimeError(msg)
def test_2d_float32(self):
......@@ -291,15 +295,15 @@ class TestPoissonRotationalOperator(object):
self.test_2d_float64()
self.test_3d_float64()
if __name__ == '__main__':
TestPoissonRotationalOperator.setup_class(enable_extra_tests=False,
TestPoissonRotationalOperator.setup_class(enable_extra_tests=False,
enable_debug_mode=False)
test = TestPoissonRotationalOperator()
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()
......