diff --git a/examples/multiresolution/scalar_advection.py b/examples/multiresolution/scalar_advection.py
index bf87e882ccc267e996bcd2f000ef4fe7af849d7a..26d500cf59e5c033c78616b5d967c2ffa2497e20 100644
--- a/examples/multiresolution/scalar_advection.py
+++ b/examples/multiresolution/scalar_advection.py
@@ -122,6 +122,8 @@ def compute(args):
lpfilter = LowpassFilter(input_variables={scalar: snpts},
output_variables={scalar: npts},
filtering_method=FilteringMethod.SPECTRAL,
+ plot_input_energy=IOParams(filename='E_fine_{fname}_{it}', frequency=args.dump_freq),
+ plot_output_energy=IOParams(filename='E_coarse_{fname}_{it}', frequency=args.dump_freq),
implementation=impl,
**extra_op_kwds)
@@ -212,7 +214,7 @@ if __name__=='__main__':
parser.set_defaults(box_origin=(0.0,), box_length=(2*np.pi,),
tstart=0.0, tend=2*np.pi, npts=(16,16,16,),
- dump_freq=1, cfl=0.5, velocity=(1.0,),
+ dump_freq=10, cfl=0.5, velocity=(1.0,),
ndim=3, compute_precision='fp64')
parser.run(compute)
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index ca812e4a68c6ff6f20ed4afd879b66db27606e82..425230bc0b901790a9e7ede20a5d65fc3688386b 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -5,6 +5,7 @@ LowpassFilter operator generator.
from hysop.constants import Implementation
from hysop.methods import Remesh
from hysop.numerics.remesh.remesh import RemeshKernel
+from hysop.tools.io_utils import IOParams
from hysop.tools.types import check_instance, to_list, first_not_None, InstanceOf
from hysop.tools.numpywrappers import npw
from hysop.tools.decorators import debug
@@ -68,7 +69,28 @@ class SpectralLowpassFilterBase(LowpassFilterBase, SpectralOperatorBase):
using the spectral method.
"""
@debug
- def __init__(self, **kwds):
+ def __init__(self, plot_input_energy=None,
+ plot_output_energy=None,
+ **kwds):
+ """
+ Initialize a SpectralLowpassFilterBase.
+
+ Parameters
+ ----------
+ plot_input_energy: IOParams, optional, defaults to None
+ Plot input field energy in a custom file.
+ plot_output_energy: IOParams, optional, defaults to None
+ Plot output field energy in a custom file.
+
+ Notes
+ -----
+ IOParams filename is formatted before being used:
+ {fname} is replaced with field name
+ {it} is replaced with simulation iteration id
+ If None is passed, no plots are generated.
+ """
+ check_instance(plot_input_energy, IOParams, allow_none=True)
+ check_instance(plot_output_energy, IOParams, allow_none=True)
super(SpectralLowpassFilterBase, self).__init__(**kwds)
@@ -84,8 +106,9 @@ class SpectralLowpassFilterBase(LowpassFilterBase, SpectralOperatorBase):
# build spectral transforms
tg_fine = self.new_transform_group()
tg_coarse = self.new_transform_group()
- Ft = tg_fine.require_forward_transform(Fin, custom_output_buffer='auto')
- Bt = tg_coarse.require_backward_transform(Fout, custom_input_buffer='B0')
+
+ Ft = tg_fine.require_forward_transform(Fin, custom_output_buffer='auto', plot_energy=plot_input_energy)
+ Bt = tg_coarse.require_backward_transform(Fout, custom_input_buffer='B0', plot_energy=plot_output_energy)
self.tg_fine = tg_fine
self.tg_coarse = tg_coarse
diff --git a/hysop/operator/base/spectral_operator.py b/hysop/operator/base/spectral_operator.py
index fa13dca8cdcb39dba36783b1d0280c077e061217..3e3a433c8bcab0aaefb07608129275494e043933 100644
--- a/hysop/operator/base/spectral_operator.py
+++ b/hysop/operator/base/spectral_operator.py
@@ -12,6 +12,7 @@ from hysop.tools.decorators import debug
from hysop.tools.units import bytes2str
from hysop.tools.numerics import is_fp, is_complex, complex_to_float_dtype, \
float_to_complex_dtype, determine_fp_types
+from hysop.tools.io_utils import IOParams
from hysop.tools.spectral_utils import SpectralTransformUtils as STU
from hysop.core.arrays.array_backend import ArrayBackend
from hysop.core.arrays.array import Array
@@ -550,7 +551,8 @@ class SpectralTransformGroup(object):
@_not_initialized
def require_forward_transform(self, field, axes=None, transform_tag=None,
- custom_output_buffer=None, action=None):
+ custom_output_buffer=None, action=None,
+ plot_energy=None):
"""
Tells this SpectralTransformGroup to build a forward SpectralTransform
on given field. Only specified axes are transformed.
@@ -594,6 +596,12 @@ class SpectralTransformGroup(object):
compute slices of the output buffer.
SpectralTransformAction.ACCUMULATE will sum the current content of the buffer
with the result of the forward transform.
+ plot_energy: IOParams, optional, defaults to None
+ Plot field energy after each call to the forward transform to a custom file.
+ IOParams filename is formatted before being used:
+ {fname} is replaced with field name
+ {it} is replaced with simulation iteration id
+ If None is passed, no plots are generated (default behaviour).
"""
transform_tag = first_not_None(transform_tag, 'default')
action = first_not_None(action, SpectralTransformAction.OVERWRITE)
@@ -614,7 +622,7 @@ class SpectralTransformGroup(object):
tag=self.tag + '_' + transform_tag + '_' + f.name,
symbolic_transform=transforms[idx],
custom_output_buffer=custom_output_buffer,
- action=action)
+ action=action, plot_energy=plot_energy)
self._forward_transforms[(f,axes,transform_tag)] = planned_transforms[idx]
else:
assert (field,axes,transform_tag) not in self._forward_transforms, msg.format(field.name, axes, transform_tag)
@@ -625,7 +633,7 @@ class SpectralTransformGroup(object):
tag=self.tag + '_' + transform_tag + '_' + field.name,
symbolic_transform=transforms,
custom_output_buffer=custom_output_buffer,
- action=action)
+ action=action, plot_energy=plot_energy)
self._forward_transforms[(field,axes,transform_tag)] = planned_transforms
return planned_transforms
@@ -633,7 +641,7 @@ class SpectralTransformGroup(object):
def require_backward_transform(self, field, axes=None, transform_tag=None,
custom_input_buffer=None,
matching_forward_transform=None,
- action=None):
+ action=None, plot_energy=None):
"""
Same as require_forward_transform but for backward transforms.
This corresponds to the following backward transform mappings:
@@ -677,8 +685,6 @@ class SpectralTransformGroup(object):
Extra tag to register the backward transform (a single scalar field can be
transformed multiple times). Default tag is 'default'.
custom_input_buffer: None or str or F, optional
-
- # self.Ft.input_buffer is just a pypencl Array.
Force this transform to take as input one of the two common transform group buffers.
Default None value will force the user to supply an input buffer.
Specifying 'B0' or 'B1' will tell the planner to take as transform input
@@ -692,7 +698,12 @@ class SpectralTransformGroup(object):
compute slices of the given output field.
SpectralTransformAction.ACCUMULATE will sum the current content of the field
with the result of the backward transform.
-
+ plot_energy: IOParams, optional, defaults to None
+ Plot field energy before each call to the backward transform to a custom file.
+ IOParams filename is formatted before being used:
+ {fname} is replaced with field name
+ {it} is replaced with simulation iteration id
+ If None is passed, no plots are generated (default behaviour).
"""
transform_tag = first_not_None(transform_tag, 'default')
action = first_not_None(action, SpectralTransformAction.OVERWRITE)
@@ -788,13 +799,14 @@ class PlannedSpectralTransform(object):
def __init__(self, transform_group, tag, symbolic_transform, action,
custom_input_buffer=None, custom_output_buffer=None,
- matching_forward_transform=None):
+ matching_forward_transform=None, plot_energy=None):
check_instance(transform_group, SpectralTransformGroup)
check_instance(transform_group.op, SpectralOperatorBase)
check_instance(tag, str)
check_instance(symbolic_transform, AppliedSpectralTransform)
check_instance(action, SpectralTransformAction)
+ check_instance(plot_energy, IOParams, allow_none=True)
assert custom_input_buffer in (None, 'B0', 'B1', 'auto'), custom_input_buffer
assert custom_output_buffer in (None, 'B0', 'B1', 'auto'), custom_output_buffer
@@ -806,6 +818,10 @@ class PlannedSpectralTransform(object):
self._custom_output_buffer = custom_output_buffer
self._matching_forward_transform = matching_forward_transform
self._action = action
+
+ self._plot_energy_ioparams = plot_energy
+ self._do_plot_energy = (plot_energy is not None)
+ del plot_energy
field = self.s.field
is_forward = self.s.is_forward
@@ -1128,7 +1144,7 @@ class PlannedSpectralTransform(object):
self._zero_fill_output_slices = zero_fill_output_slices
self._backend = dfield.backend
-
+
if self.DEBUG:
def axis_format(info):
prefix='\n'+' '*4
@@ -1171,6 +1187,8 @@ class PlannedSpectralTransform(object):
print 'transpose_info: {}'.format(axis_format(transpose_info))
print ':ZERO FILL:'
print 'zero_fill_output_slices: {}'.format(slc_format(self._zero_fill_output_slices))
+
+
def get_mapped_input_buffer(self):
return self.get_mapped_full_input_buffer()[self.input_slices]
@@ -1487,14 +1505,27 @@ class PlannedSpectralTransform(object):
mem_tag = self.transform_group.mem_tag
kind = backend.kind
- B0_tag = '{}_{}_B0'.format(mem_tag, kind)
- B1_tag = '{}_{}_B1'.format(mem_tag, kind)
- TMP_tag = '{}_{}_TMP'.format(mem_tag, kind)
- self.B0_tag, self.B1_tag, self.TMP_tag = B0_tag, B1_tag, TMP_tag
-
- return {B0_tag: nbytes,
- B1_tag: nbytes,
- TMP_tag: tmp_nbytes}
+ B0_tag = '{}_{}_B0'.format(mem_tag, kind)
+ B1_tag = '{}_{}_B1'.format(mem_tag, kind)
+ TMP_tag = '{}_{}_TMP'.format(mem_tag, kind)
+ PLOT_tag = '{}_{}_PLOT'.format(mem_tag, kind)
+ self.B0_tag, self.B1_tag, self.TMP_tag, self.PLOT_tag = B0_tag, B1_tag, TMP_tag, PLOT_tag
+
+ requests = {B0_tag: nbytes,
+ B1_tag: nbytes,
+ TMP_tag: tmp_nbytes}
+
+ # determine arrays to compute energy for plotting
+ if self._do_plot_energy:
+ transforms = self.transforms
+ print transforms
+ import sys
+ sys.exit(1)
+ plot_nbytes = compute_nbytes(maxk, self.dfield.dtype)
+ requests[PLOT_tag] = plot_nbytes
+
+ return requests
+
@_discretized
def setup(self, work):
@@ -1804,8 +1835,13 @@ SPECTRAL TRANSFORM SETUP
self._queue = queue
self._ready = True
- def __call__(self):
+ def __call__(self, simu=None):
assert (self._ready)
assert (self._queue is not None)
- return self._queue.execute()
+ if (self.is_backward and self._do_plot_energy):
+ evt = self.plot_energy(simu=simu)
+ evt = self._queue.execute()
+ if (self.is_forward and self._do_plot_energy):
+ evt = self.plot_energy(simu=simu)
+ return evt