diff --git a/hysop/backend/host/fortran/operator/poisson.py b/hysop/backend/host/fortran/operator/poisson.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c7a767efbf758bcca00492c77e0b01bb6ca6e8b
--- /dev/null
+++ b/hysop/backend/host/fortran/operator/poisson.py
@@ -0,0 +1,100 @@
+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
diff --git a/hysop/backend/host/fortran/operator/poisson_rotational.py b/hysop/backend/host/fortran/operator/poisson_rotational.py
index 5620efbb0fade2da9b1ee2831a14fff9c44fe081..a718a487b5482ca9f8a611019dcad833b7cd7e9a 100644
--- a/hysop/backend/host/fortran/operator/poisson_rotational.py
+++ b/hysop/backend/host/fortran/operator/poisson_rotational.py
@@ -1,4 +1,3 @@
-
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()
-
diff --git a/hysop/core/graph/computational_graph.py b/hysop/core/graph/computational_graph.py
index 8c9d2ab0d6c5c65e6920040b156e1c4b9ca2c7ec..417207d514135fe2ad98d62f57bb575a289dda76 100644
--- a/hysop/core/graph/computational_graph.py
+++ b/hysop/core/graph/computational_graph.py
@@ -19,29 +19,29 @@ class ComputationalGraph(ComputationalGraphNode):
"""
Interface of an abstract graph of continuous operators (ie. a computational graph).
"""
-
+
__metaclass__ = ABCMeta
__FORCE_REPORTS__ = False
-
+
@debug
def __init__(self, **kwds):
"""
Parameters
----------
kwds: arguments for base classes.
-
+
Notes
-----
- The following base class variables cannot be specified during graph construction:
+ The following base class variables cannot be specified during graph construction:
variables, input_variables, output_variables
- Order of operation is: add_node, initialize, discretize,
+ Order of operation is: add_node, initialize, discretize,
get_work_properties, setup, apply, finalize.
Nodes can also be added during pre initialization step.
- Graph building is done at the end of the initialization step, after all internal
+ Graph building is done at the end of the initialization step, after all internal
nodes have been initialized.
"""
-
+
if ('input_fields' in kwds.keys()) or ('output_fields' in kwds.keys()):
msg='input_fields or output_fields parameters should not be used in {}, they are \
deduced during graph construction (building step).'.format(cls)
@@ -50,27 +50,27 @@ class ComputationalGraph(ComputationalGraphNode):
msg='input_params or output_params parameters should not be used in {}, they are \
deduced during graph construction (building step).'.format(cls)
raise ValueError(msg)
-
+
super(ComputationalGraph,self).__init__(input_fields=None, output_fields=None,
**kwds)
-
+
self.nodes = []
self.graph = None
self.graph_built = False
self.graph_is_rendering = False
-
+
@graph_built
def _get_profiling_info(self):
"""Collect profiling informations of profiled attributes."""
for node in self.nodes:
self._profiler += node._profiler
-
+
def field_requirements_report(self, requirements):
inputs, outputs = {}, {}
sinputs, soutputs = {}, {}
def sorted_reqs(reqs):
return sorted(reqs, key=lambda x: \
- u'{}::{}'.format(x.operator.pretty_name.decode('utf-8'),
+ u'{}::{}'.format(x.operator.pretty_name.decode('utf-8'),
x.field.pretty_name.decode('utf-8')))
for field, mreqs in requirements.input_field_requirements.iteritems():
for td, reqs in mreqs.requirements.iteritems():
@@ -87,7 +87,7 @@ class ComputationalGraph(ComputationalGraphNode):
ghosts = u'{}<=ghosts<{}'.format(min_ghosts, max_ghosts)
can_split=req.can_split.view(npw.int8)
basis=u'{}'.format(u','.join(str(basis)[:3] for basis in req.basis)) \
- if req.basis else 'ANY'
+ if req.basis else 'ANY'
tstates=u'{}'.format(u','.join(str(ts) for ts in req.tstates)) \
if req.tstates else 'ANY'
sin.append( (opname, fname, ghosts, basis, tstates) )
@@ -106,11 +106,11 @@ class ComputationalGraph(ComputationalGraphNode):
ghosts = u'{}<=ghosts<{}'.format(min_ghosts, max_ghosts)
can_split=req.can_split.view(npw.int8)
basis=u'{}'.format(u','.join(str(basis)[:3] for basis in req.basis)) \
- if req.basis else u'ANY'
+ if req.basis else u'ANY'
tstates=u'{}'.format(u','.join(str(ts) for ts in req.tstates)) \
if req.tstates else u'ANY'
sin.append( (opname, fname, ghosts, basis, tstates) )
-
+
titles = [[(u'OPERATOR', u'FIELD', u'GHOSTS', u'BASIS', u'TSTATES')]]
name_size = max(len(s[0]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
field_size = max(len(s[1]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
@@ -119,7 +119,7 @@ class ComputationalGraph(ComputationalGraphNode):
tstates_size = max(len(s[4]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
template = u'\n {:<{name_size}} {:^{field_size}} {:^{ghosts_size}} {:^{basis_size}} {:^{tstates_size}}'
-
+
ss= u'>INPUTS:'
if sinputs:
for (td, sreqs) in sinputs.iteritems():
@@ -128,12 +128,12 @@ class ComputationalGraph(ComputationalGraphNode):
else:
ss+=u'\n {}'.format(td)
ss+= template.format(*titles[0][0],
- name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
+ name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
basis_size=basis_size, tstates_size=tstates_size)
for (opname, fname, ghosts, basis, tstates) in sreqs:
ss+=template.format(
opname, fname, ghosts, basis, tstates,
- name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
+ name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
basis_size=basis_size, tstates_size=tstates_size)
else:
ss+=u' None'
@@ -145,12 +145,12 @@ class ComputationalGraph(ComputationalGraphNode):
else:
ss+=u'\n {}'.format(td)
ss+= template.format(*titles[0][0],
- name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
+ name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
basis_size=basis_size, tstates_size=tstates_size)
for (opname, fname, ghosts, basis, tstates) in sreqs:
ss+=template.format(
opname, fname, ghosts, basis, tstates,
- name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
+ name_size=name_size, field_size=field_size, ghosts_size=ghosts_size,
basis_size=basis_size, tstates_size=tstates_size)
else:
ss+=u' None'
@@ -173,10 +173,10 @@ class ComputationalGraph(ComputationalGraphNode):
outfields = u'[{}]'.format(foutputs) if foutputs else u''
inparams = u'[{}]'.format(pinputs) if pinputs else u''
outparams = u'[{}]'.format(poutputs) if poutputs else u''
-
+
inputs = u'{}{}{}'.format(infields, u'x' if infields and inparams else u'', inparams)
outputs = u'{}{}{}'.format(outfields, u'x' if outfields and outparams else u'', outparams)
-
+
if inputs == u'':
inputs=u'no inputs'
if outputs == u'':
@@ -198,7 +198,7 @@ class ComputationalGraph(ComputationalGraphNode):
if outputs == '':
outputs=u'no outputs'
ops.setdefault(None, []).append( (op.pretty_name.decode('utf-8'), inputs, outputs, type(op).__name__) )
-
+
name_size = max(strlen(s[0]) for ss in ops.values() for s in ss)
in_size = max(strlen(s[1]) for ss in ops.values() for s in ss)
out_size = max(strlen(s[2]) for ss in ops.values() for s in ss)
@@ -211,7 +211,7 @@ class ComputationalGraph(ComputationalGraphNode):
ss += u'\n>{}'.format(domain.short_description())
ss += u'\n {:<{name_size}} {:<{in_size}} {:<{out_size}} {:<{type_size}}'.format(
'OPERATOR', 'INPUTS', 'OUTPUTS', 'OPERATOR TYPE',
- name_size=name_size, in_size=in_size,
+ name_size=name_size, in_size=in_size,
out_size=out_size, type_size=type_size)
for (opname, inputs, outputs, optype) in dops:
ss += u'\n {:<{name_size}} {:<{in_size}} -> {:<{out_size}} {:<{type_size}}'.format(
@@ -222,17 +222,17 @@ class ComputationalGraph(ComputationalGraphNode):
ss += u'\n>Domainless operators:'
ss += u'\n {:<{name_size}} {:<{in_size}} {:<{out_size}} {:<{type_size}}'.format(
'OPERATOR', 'INPUTS', 'OUTPUTS', 'OPERATOR TYPE',
- name_size=name_size, in_size=in_size,
+ name_size=name_size, in_size=in_size,
out_size=out_size, type_size=type_size)
for (opname, inputs, outputs, optype) in ops[None]:
ss += u'\n {:<{name_size}} {:<{in_size}} -> {:<{out_size}} {:<{type_size}}'.format(
opname, inputs, outputs, optype,
name_size=name_size, in_size=in_size,
out_size=out_size, type_size=type_size)
-
+
title=u' ComputationalGraph {} domain and operator report '.format(self.pretty_name.decode('utf-8'))
return u'\n{}\n'.format(framed_str(title=title, msg=ss[1:])).encode('utf-8')
-
+
def topology_report(self):
ss=''
@@ -247,7 +247,7 @@ class ComputationalGraph(ComputationalGraphNode):
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
fields = self.fields
-
+
topologies = {}
for field in self.fields:
field_topologies = {}
@@ -261,7 +261,7 @@ class ComputationalGraph(ComputationalGraphNode):
topo = op.output_fields[field]
field_topologies.setdefault(topo, []).append(op)
for topo in field_topologies.keys():
- pnames = set(op.pretty_name.decode('utf-8') for op in field_topologies[topo])
+ pnames = set(op.pretty_name.decode('utf-8') for op in field_topologies[topo])
sops=u', '.join(sorted(pnames))
entries = (str(topo.backend.kind).lower(), topo.tag, sops)
topologies.setdefault(field, []).append(entries)
@@ -272,7 +272,7 @@ class ComputationalGraph(ComputationalGraphNode):
template = u'\n {:<{backend_size}} {:<{topo_size}} {}'
sizes = {'backend_size': backend_size,
'topo_size': topo_size}
-
+
ss = u''
for field in sorted(self.fields, key=lambda x: x.tag):
ss += u'\n>FIELD {}::{}'.format(field.name, field.pretty_name.decode('utf-8'))
@@ -284,7 +284,7 @@ class ComputationalGraph(ComputationalGraphNode):
title = u' ComputationalGraph {} fields report '.format(self.pretty_name.decode('utf-8'))
ss = u'\n{}\n'.format(framed_str(title=title, msg=ss[1:]))
return ss.encode('utf-8')
-
+
def operator_report(self):
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
@@ -299,16 +299,16 @@ class ComputationalGraph(ComputationalGraphNode):
foutputs = u','.join( sorted([u'{}.{}'.format(f.pretty_name.decode('utf-8'),
t.pretty_tag.decode('utf-8'))
for (f,t) in op.output_fields.iteritems()]))
- pinputs = u','.join( sorted([p.pretty_name.decode('utf-8')
+ pinputs = u','.join( sorted([p.pretty_name.decode('utf-8')
for p in op.input_params.values()]))
- poutputs = u','.join( sorted([p.pretty_name.decode('utf-8')
+ poutputs = u','.join( sorted([p.pretty_name.decode('utf-8')
for p in op.output_params.values()]))
infields = u'[{}]'.format(finputs) if finputs else u''
outfields = u'[{}]'.format(foutputs) if foutputs else u''
inparams = u'[{}]'.format(pinputs) if pinputs else u''
outparams = u'[{}]'.format(poutputs) if poutputs else u''
-
+
inputs = u'{}{}{}'.format(infields, u'x' if infields and inparams else u'', inparams)
outputs = u'{}{}{}'.format(outfields, u'x' if outfields and outparams else u'', outparams)
if inputs == '':
@@ -317,7 +317,7 @@ class ComputationalGraph(ComputationalGraphNode):
outputs=u'no outputs'
ops.append( (op.pretty_name.decode('utf-8'), inputs, outputs, type(op).__name__) )
-
+
isize = len(self.sorted_nodes)
isize = int(npw.ceil(npw.log10(isize)))
name_size = max(len(s[0]) for s in ops)
@@ -328,7 +328,7 @@ class ComputationalGraph(ComputationalGraphNode):
ss = u' {:<{isize}} {:<{name_size}} {:<{in_size}} {:<{out_size}} {:<{type_size}}'.format(
'ID', 'OPERATOR', 'INPUTS', 'OUTPUTS', 'OPERATOR TYPE',
isize=isize,
- name_size=name_size, in_size=in_size,
+ name_size=name_size, in_size=in_size,
out_size=out_size, type_size=type_size)
for i, (opname, inputs, outputs, optype) in enumerate(ops):
ss += u'\n {:>{isize}} {:<{name_size}} {:<{in_size}} -> {:<{out_size}} {:<{type_size}}'.format(
@@ -336,7 +336,7 @@ class ComputationalGraph(ComputationalGraphNode):
isize=isize,
name_size=name_size, in_size=in_size,
out_size=out_size, type_size=type_size)
-
+
title = u' ComputationalGraph {} discrete operator report '.format(self.pretty_name.decode('utf-8'))
return u'\n{}\n'.format(framed_str(title=title, msg=ss)).encode('utf-8')
@@ -346,7 +346,7 @@ class ComputationalGraph(ComputationalGraphNode):
for (domain, ops) in node.get_domains().iteritems():
domains.setdefault(domain, set()).update(ops)
return domains
-
+
@graph_built
def get_topologies(self):
topologies = {}
@@ -358,7 +358,7 @@ class ComputationalGraph(ComputationalGraphNode):
topologies[backend] = set()
topologies[backend].update(topos)
return topologies
-
+
@debug
@not_initialized
def push_nodes(self, *args):
@@ -374,7 +374,7 @@ class ComputationalGraph(ComputationalGraphNode):
else:
msg = 'Given node is not an instance of ComputationalGraphNode (got a {}).'
raise ValueError(msg.format(node.__class__))
-
+
def available_methods(self):
avail_methods = {}
if not self.nodes:
@@ -388,33 +388,33 @@ class ComputationalGraph(ComputationalGraphNode):
else:
avail_methods[k] = v
return avail_methods
-
+
def default_method(self):
return {}
@debug
- def initialize(self,
+ def initialize(self,
is_root=True,
- topgraph_method=None,
+ topgraph_method=None,
outputs_are_inputs=False,
**kwds):
if self.initialized:
return
self.is_root = is_root
-
+
if is_root:
self.pre_initialize(**kwds)
msg=u'ComputationalGraph {} is empty.'
assert len(self.nodes) > 0, msg.format(self.pretty_name.decode('utf-8')).encode('utf-8')
-
+
for node in self.nodes:
node.pre_initialize(**kwds)
-
+
method = self._setup_method(topgraph_method)
self.handle_method(method)
-
+
for node in self.nodes:
node.initialize(is_root=False,
outputs_are_inputs=False,
@@ -430,7 +430,7 @@ class ComputationalGraph(ComputationalGraphNode):
if is_root:
field_requirements = self.get_and_set_field_requirements()
field_requirements.build_topologies()
-
+
self._build_graph(outputs_are_inputs=outputs_are_inputs, current_level=0)
# fix for auto generated nodes
@@ -439,15 +439,15 @@ class ComputationalGraph(ComputationalGraphNode):
node.get_and_set_field_requirements()
assert (node._field_requirements is not None)
- # from now on, all nodes contained in self.nodes are ComputationalGraphOperator
+ # from now on, all nodes contained in self.nodes are ComputationalGraphOperator
# (ordered for sequential execution)
self.handle_topologies(self.input_topology_states, self.output_topology_states)
self.check()
-
+
self.initialized=True
if is_root:
self.post_initialize(**kwds)
-
+
@debug
def check(self):
super(ComputationalGraph, self).check()
@@ -467,7 +467,7 @@ class ComputationalGraph(ComputationalGraphNode):
if (self.is_root and __VERBOSE__ or __DEBUG__ or self.__FORCE_REPORTS__):
print self.field_requirements_report(requirements)
return requirements
-
+
@debug
def handle_topologies(self, input_topology_states, output_topology_states):
# do not call super method
@@ -489,17 +489,17 @@ class ComputationalGraph(ComputationalGraphNode):
from hysop.core.graph.graph_builder import GraphBuilder
builder = GraphBuilder(node=self)
- builder.configure(current_level=current_level,
+ builder.configure(current_level=current_level,
outputs_are_inputs=outputs_are_inputs, **kwds)
-
+
builder.build_graph()
-
+
# keep variables
- self._init_base(input_fields=builder.input_fields,
+ self._init_base(input_fields=builder.input_fields,
output_fields=builder.output_fields,
input_params=builder.input_params,
output_params=builder.output_params)
-
+
self.graph = builder.graph
self.reduced_graph = builder.reduced_graph
self.sorted_nodes = builder.sorted_nodes
@@ -507,12 +507,12 @@ class ComputationalGraph(ComputationalGraphNode):
self.input_topology_states = builder.op_input_topology_states
self.output_topology_states = builder.op_output_topology_states
-
+
self.initial_input_topology_states = builder.input_topology_states
self.final_output_topology_states = builder.output_topology_states
-
+
self.level = current_level
-
+
self.initialized = True
self.graph_built = True
@@ -541,13 +541,13 @@ class ComputationalGraph(ComputationalGraphNode):
else:
command_queues = None
active_ops = None
-
+
def draw():
import time
from gi.repository import Gtk, GObject
from graph_tool.draw import GraphWindow, sfdp_layout
-
- pos_layout = sfdp_layout(graph)
+
+ pos_layout = sfdp_layout(graph)
win = GraphWindow(graph, pos_layout, geometry=(800,600),
vertex_text = vertex_text,
@@ -565,19 +565,19 @@ class ComputationalGraph(ComputationalGraphNode):
win.graph.queue_draw()
time.sleep(0.01)
return True
-
+
GObject.idle_add(update_window)
win.connect("delete_event", Gtk.main_quit)
win.show_all()
Gtk.main()
self.graph_is_rendering = False
-
+
self.graph_is_rendering = True
from threading import Thread
display_thread = Thread(target=draw)
display_thread.start()
-
+
@debug
@graph_built
def discretize(self):
@@ -588,7 +588,7 @@ class ComputationalGraph(ComputationalGraphNode):
op = operators[vertex]
if not op.discretized:
op.discretize()
-
+
if self.is_root:
input_discrete_fields = {}
for (field,topo) in self.input_fields.iteritems():
@@ -596,7 +596,7 @@ class ComputationalGraph(ComputationalGraphNode):
istate = istate.copy(is_read_only=False) # problem inputs are writeable for initialization
dfield = field.discretize(topo, istate)
input_discrete_fields[field] = dfield
-
+
output_discrete_fields = {}
for field,topo in self.output_fields.iteritems():
ostate = self.final_output_topology_states[field][1]
@@ -608,15 +608,15 @@ class ComputationalGraph(ComputationalGraphNode):
self.input_discrete_fields = input_discrete_fields
self.output_discrete_fields = output_discrete_fields
-
+
self.discretized = True
-
-
+
+
@debug
@discretized
def get_work_properties(self):
requests = MultipleOperatorMemoryRequests()
-
+
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
for vid in self.sorted_nodes:
@@ -630,8 +630,8 @@ class ComputationalGraph(ComputationalGraphNode):
ss = (srequests if (srequests != u'') else u' *no extra work requested*')
title= u' ComputationalGraph {} work properties report '.format(self.pretty_name.decode('utf-8'))
vprint(u'\n{}\n'.format(framed_str(title=title, msg=ss)).encode('utf-8'))
- return requests
-
+ return requests
+
@debug
@discretized
def setup(self, work=None):
@@ -649,9 +649,9 @@ class ComputationalGraph(ComputationalGraphNode):
op.setup(work=work)
self.ready = True
- def build(self, outputs_are_inputs=True, method=None, allow_subbuffers=False):
+ def build(self, outputs_are_inputs=True, method=None, allow_subbuffers=False):
"""
- Shortcut for initialize(), discretize(), get_work_properties(), setup()
+ Shortcut for initialize(), discretize(), get_work_properties(), setup()
for quick graph initialization.
"""
self.initialize(outputs_are_inputs=outputs_are_inputs,
@@ -669,7 +669,7 @@ class ComputationalGraph(ComputationalGraphNode):
drawing = self.graph_is_rendering
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
-
+
if drawing:
active_ops = reduced_graph.vertex_properties['active_ops']
old_color = None
@@ -690,14 +690,14 @@ class ComputationalGraph(ComputationalGraphNode):
@debug
@ready
- def finalize(self):
+ def finalize(self, **kwds):
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
for vid in self.sorted_nodes:
vertex = reduced_graph.vertex(vid)
op = operators[vertex]
if op.ready:
- op.finalize()
+ op.finalize(**kwds)
self.ready = False
@classmethod
diff --git a/hysop/numerics/fftw_f/fft2d.f90 b/hysop/numerics/fftw_f/fft2d.f90
index 4a0a708d6c6989cbab7142aafebec699694ea541..98f79f55ee20467eef16980b1ea6848bbda1a9ce 100755
--- a/hysop/numerics/fftw_f/fft2d.f90
+++ b/hysop/numerics/fftw_f/fft2d.f90
@@ -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
diff --git a/hysop/numerics/fftw_f/fft3d.f90 b/hysop/numerics/fftw_f/fft3d.f90
index c839f1752a9a2d510bc00faca0738283be900c83..c87435d7f02e4a842cdd3da3290d14a0f587bc23 100755
--- a/hysop/numerics/fftw_f/fft3d.f90
+++ b/hysop/numerics/fftw_f/fft3d.f90
@@ -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()
diff --git a/hysop/numerics/fftw_f/fftw2py.f90 b/hysop/numerics/fftw_f/fftw2py.f90
index 09517a4ca3034a9800a76e99c94ff9656963e2a9..230d6455b160e938f0ccbf0988bda6a1060a4929 100755
--- a/hysop/numerics/fftw_f/fftw2py.f90
+++ b/hysop/numerics/fftw_f/fftw2py.f90
@@ -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.
diff --git a/hysop/numerics/fftw_f/fftw2py.pyf b/hysop/numerics/fftw_f/fftw2py.pyf
index 74ac41ee1c12458b18ea26880d940ab6f0a52d9f..e7338d4c38ac4f67d61712dfc4d4b211e4fc59ad 100644
--- a/hysop/numerics/fftw_f/fftw2py.pyf
+++ b/hysop/numerics/fftw_f/fftw2py.pyf
@@ -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
diff --git a/hysop/operator/poisson.py b/hysop/operator/poisson.py
index a07d8cfa58e2ad4a07989f8cba82e46a5a27f450..fe575f60d1e643a49f8f6445ec9a414df9b23115 100644
--- a/hysop/operator/poisson.py
+++ b/hysop/operator/poisson.py
@@ -1,4 +1,3 @@
-
"""
@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)
diff --git a/hysop/operator/poisson_rotational.py b/hysop/operator/poisson_rotational.py
index 4c30890232104945f8fd98b2ab08b15b0e5f0330..599c82d4cb4271dbbf6ebd8b4ad45f2c4ee4075a 100644
--- a/hysop/operator/poisson_rotational.py
+++ b/hysop/operator/poisson_rotational.py
@@ -1,4 +1,3 @@
-
"""
@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)
diff --git a/hysop/operator/tests/test_poisson.py b/hysop/operator/tests/test_poisson.py
index d5cf80c30852a67d46033ec6a42b98bfa48b210f..c345b8c4514507df6675dd65a11fe9777300c4f7 100644
--- a/hysop/operator/tests/test_poisson.py
+++ b/hysop/operator/tests/test_poisson.py
@@ -1,4 +1,3 @@
-
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()
diff --git a/hysop/operator/tests/test_poisson_rotational.py b/hysop/operator/tests/test_poisson_rotational.py
index 6223d9b5613f5d8b1b02945563010fb49696a9b3..971eaba7971c1a7e5a055ce322254bb1f753138a 100644
--- a/hysop/operator/tests/test_poisson_rotational.py
+++ b/hysop/operator/tests/test_poisson_rotational.py
@@ -1,3 +1,4 @@
+# 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()