bubble example

examples/bubble/periodic_bubble.py

+
+## HySoP Example: Bubble nD
+## Osher1995 (first part): 
+##  A Level Set Formulation of Eulerian Interface Capturing Methods for 
+##  Incompressible Fluid flows.
+
+import os
+import numpy as np
+
+Bc = ((0.5,0.5,0.5),)
+Br = (0.25,)
+
+def bubble_mask(coords, Bc, Br):
+    assert len(Bc)==len(Br)>=1
+    dim = len(coords)
+    shape = tuple(c.size for c in coords[::-1])
+    mask = np.zeros(shape=shape, dtype=np.int8)
+    for (center, radius) in zip(Bc, Br):
+        center = center[:dim]
+        D = np.zeros(shape=shape, dtype=np.float64)
+        for (Xi,Ci) in zip(coords, center):
+            D += (Xi-Ci)**2
+        mask |= (D<=radius**2)
+    return mask
+
+def init_rho(data, coords, Bc, Br, rho1, rho2):
+    mask = bubble_mask(coords, Bc, Br)
+    data[0][...] = mask*rho1 + (1-mask)*rho2
+    
+def init_mu(data, coords, Bc, Br, mu1, mu2):
+    mask = bubble_mask(coords, Bc, Br)
+    data[0][...] = mask*mu1 + (1-mask)*mu2
+
+def init_vorticity(data, coords):
+    for d in data:
+        d[...] = 0.0
+
+def compute(args):
+    from hysop import Box, Simulation, Problem, MPIParams, IOParams
+    from hysop.defaults import VelocityField, VorticityField, \
+                               DensityField, ViscosityField, \
+                               EnstrophyParameter, TimeParameters, \
+                               VolumicIntegrationParameter
+    from hysop.constants import Implementation, AdvectionCriteria
+
+    from hysop.operators import DirectionalAdvection, DirectionalDiffusion,        \
+                                DirectionalStretchingDiffusion,                    \
+                                PoissonRotational, AdaptiveTimeStep,               \
+                                Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
+                                ParameterPlotter, Integrate, HDF_Writer
+
+    from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
+                              ComputeGranularity, Interpolation
+
+    # Define the domain
+    dim  = args.ndim
+    npts = args.npts
+    box  = Box(origin=args.box_origin, length=args.box_length, dim=dim)
+    
+    # Get default MPI Parameters from domain (even for serial jobs)
+    mpi_params = MPIParams(comm=box.task_comm,
+                           task_id=box.current_task())
+    
+    # Setup usual implementation specific variables
+    impl = args.impl
+    extra_op_kwds = {'mpi_params': mpi_params}
+    if (impl is Implementation.PYTHON):
+        method = {}
+    elif (impl is Implementation.OPENCL):
+        # For the OpenCL implementation we need to setup the compute device
+        # and configure how the code is generated and compiled at runtime.
+                
+        # Create an explicit OpenCL context from user parameters
+        from hysop.backend.device.opencl.opencl_tools import get_or_create_opencl_env
+        cl_env = get_or_create_opencl_env(mpi_params=mpi_params, 
+                                          platform_id=args.cl_platform_id, 
+                                          device_id=args.cl_device_id)
+        
+        # Configure OpenCL kernel generation and tuning (already done by HysopArgParser)
+        from hysop.methods import OpenClKernelConfig
+        method = { OpenClKernelConfig: args.opencl_kernel_config }
+        
+        # Setup opencl specific extra operator keyword arguments
+        extra_op_kwds['cl_env'] = cl_env
+    else:
+        msg='Unknown implementation \'{}\'.'.format(impl)
+        raise ValueError(msg)
+    
+    # Define parameters and field (time, timestep, velocity, vorticity, enstrophy)
+    t, dt = TimeParameters(dtype=args.dtype)
+    velo  = VelocityField(domain=box, dtype=args.dtype)
+    vorti = VorticityField(domain=box, dtype=args.dtype)
+    rho   = DensityField(domain=box, dtype=args.dtype) 
+    mu    = ViscosityField(domain=box, dtype=args.dtype, mu=True)
+
+    enstrophy = EnstrophyParameter(dtype=args.dtype)
+    rhov = VolumicIntegrationParameter(field=rho)
+    muv  = VolumicIntegrationParameter(field=mu)
+    
+    ### Build the directional operators
+    #> Directional advection 
+    advec = DirectionalAdvection(implementation=impl,
+            name='advec',
+            velocity = velo,       
+            advected_fields = (vorti, rho, mu),
+            velocity_cfl = args.cfl,
+            variables = {velo: npts, vorti: npts, rho: npts, mu: npts},
+            dt=dt, **extra_op_kwds)
+    #> Directional stretching + diffusion
+    if (dim==3):
+        stretch_diffuse = DirectionalStretchingDiffusion(implementation=impl,
+                 name='stretch_diffuse',
+                 formulation = args.stretching_formulation,
+                 viscosity = mu,
+                 velocity  = velo,       
+                 vorticity = vorti,
+                 variables = {velo: npts, vorti: npts, mu: npts},
+                 dt=dt, **extra_op_kwds)
+    elif (dim==2):
+        stretch_diffuse = DirectionalDiffusion(implementation=impl,
+                 name='stretch_diffuse',
+                 coeffs = mu,
+                 fields  = vorti,
+                 variables = {vorti: npts, mu: npts},
+                 dt=dt, **extra_op_kwds)
+    else:
+        msg='Unsupported dimension {}.'.format(dim)
+        raise RuntimeError(msg)
+
+    #> Directional splitting operator subgraph
+    splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
+    splitting.push_operators(advec, stretch_diffuse)
+    
+    ### Build standard operators
+    #> Poisson operator to recover the velocity from the vorticity
+    poisson = PoissonRotational(name='poisson', velocity=velo, vorticity=vorti, 
+                            variables={velo:npts, vorti: npts}, 
+                            projection=args.reprojection_frequency,
+                            implementation=impl, **extra_op_kwds)
+
+    #> Operators to dump rho and mu
+    io_params = IOParams(filename='fields', frequency=args.dump_freq)
+    dump_fields = HDF_Writer(name='hdf_writer',
+                             io_params=io_params,
+                             variables={velo: npts, 
+                                        vorti: npts,    
+                                        rho: npts, 
+                                        mu: npts})
+    
+    #> Operator to compute the infinite norm of the velocity
+    min_max_U = MinMaxFieldStatistics(name='min_max_U', field=velo,
+            Finf=True, implementation=impl, variables={velo:npts},
+            **extra_op_kwds)
+    #> Operator to compute the infinite norm of the vorticity
+    min_max_W = MinMaxFieldStatistics(name='min_max_W', field=vorti,
+            Finf=True, implementation=impl, variables={vorti:npts},
+            **extra_op_kwds)
+    #> Operator to track min and max density and viscosity
+    min_max_rho = MinMaxFieldStatistics(name='min_max_rho', field=rho,
+            Fmin=True, Fmax=True,
+            implementation=impl, variables={rho:npts},
+            **extra_op_kwds)
+    min_max_mu = MinMaxFieldStatistics(name='min_max_mu', field=mu,
+            Fmin=True, Fmax=True,
+            implementation=impl, variables={mu:npts},
+            **extra_op_kwds)
+    #> Operators to compute the integrated enstrophy, density and viscosity
+    integrate_enstrophy = Enstrophy(name='enstrophy', vorticity=vorti, enstrophy=enstrophy,
+            variables={vorti:npts}, implementation=impl, **extra_op_kwds)
+    integrate_rho = Integrate(name='integrate_rho', field=rho, variables={rho: npts},
+                                    parameter=rhov, scaling='normalize',
+                                    implementation=impl, **extra_op_kwds)
+    integrate_mu = Integrate(name='integrate_mu', field=mu, variables={mu: npts},
+                                    parameter=muv, scaling='normalize',
+                                    implementation=impl, **extra_op_kwds)
+    
+
+    ### Adaptive timestep operator
+    adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=0.1)
+    dt_cfl   = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
+                                            equivalent_CFL=True)
+    dt_advec = adapt_dt.push_advection_criteria(lcfl=args.lcfl, Finf=min_max_W.Finf,
+                                                    criteria=AdvectionCriteria.W_INF)
+    
+    ## Create the problem we want to solve and insert our 
+    # directional splitting subgraph and the standard operators.
+    # The method dictionnary passed to this graph will be dispatched
+    # accross all operators contained in the graph.
+    method.update(
+            { 
+               ComputeGranularity:    args.compute_granularity,
+               SpaceDiscretization:   args.fd_order,
+               TimeIntegrator:        args.time_integrator,
+               Remesh:                args.remesh_kernel,
+               Interpolation:         args.interpolation
+            }
+    )
+    problem = Problem(method=method)
+    problem.insert(poisson, 
+                   dump_fields,
+                   splitting, 
+                   # integrate_enstrophy, integrate_rho, integrate_mu,
+                   min_max_rho, min_max_mu, min_max_U, min_max_W, 
+                   adapt_dt)
+    problem.build()
+
+    # If a visu_rank was provided, and show_graph was set,
+    # display the graph on the given process rank.
+    if args.display_graph:
+        problem.display(args.visu_rank)
+    
+    # Create a simulation
+    # (do not forget to specify the t and dt parameters here)
+    simu = Simulation(start=args.tstart, end=args.tend, 
+                      nb_iter=args.nb_iter,
+                      max_iter=args.max_iter,
+                      dt0=args.dt, times_of_interest=args.dump_times,
+                      t=t, dt=dt)
+    simu.write_parameters(t, dt_cfl, dt_advec, dt, 
+            enstrophy, rhov, muv,
+            min_max_U.Finf, min_max_W.Finf, 
+            min_max_rho.Fmin, min_max_rho.Fmax,
+            min_max_mu.Fmin, min_max_mu.Fmax,
+            adapt_dt.equivalent_CFL,
+            filename='parameters.txt', precision=8)
+    
+    # Initialize only the vorticity
+    dfields = problem.input_discrete_fields
+    dfields[vorti].initialize(formula=init_vorticity)
+    dfields[rho].initialize(formula=init_rho, rho1=args.rho1, rho2=args.rho2, Bc=Bc, Br=Br)
+    dfields[mu].initialize(formula=init_mu, mu1=args.mu1, mu2=args.mu2, Bc=Bc, Br=Br)
+
+    # Finally solve the problem 
+    problem.solve(simu, dry_run=args.dry_run)
+    
+    # Finalize
+    problem.finalize()
+
+
+if __name__=='__main__':
+    from examples.example_utils import HysopArgParser, colors
+
+    class PeriodicBubbleArgParser(HysopArgParser):
+        def __init__(self):
+            prog_name = 'periodic_bubble'
+            default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(), 
+                    prog_name)
+
+            description=colors.color('HySoP Periodic Bubble Example: ', fg='blue', style='bold')
+            description+=colors.color('[Osher 1995] (first part)', fg='yellow', style='bold')
+            description+=colors.color('\nA Level Set Formulation of Eulerian Interface Capturing Methods for '
+                    +'Incompressible Fluid flows.',
+                    fg='yellow')
+            description+='\n'
+            description+='\nThis example focuses on a validation study for the '
+            description+='hybrid particle-mesh vortex method for varying densities without using a levelset function.'
+    
+            super(PeriodicBubbleArgParser, self).__init__(
+                 prog_name=prog_name,
+                 description=description,
+                 default_dump_dir=default_dump_dir)
+
+        def _add_main_args(self):
+            args = super(PeriodicBubbleArgParser, self)._add_main_args()
+            args.add_argument('-rho1', '--bubble-density', type=float,
+                                dest='rho1',
+                                help='Set the density of the bubbles.')
+            args.add_argument('-rho2', '--fluid-density', type=float,
+                                dest='rho2',
+                                help='Set the density of the fluid surrounding the bubbles.')
+            args.add_argument('-mu1', '--bubble-viscosity', type=float,
+                                dest='mu1',
+                                help='Set the viscosity of the bubbles.')
+            args.add_argument('-mu2', '--fluid-viscosity', type=float,
+                                dest='mu2',
+                                help='Set the viscosity of the fluid surrounding the bubbles.')
+            return args
+
+        def _check_main_args(self, args):
+            super(PeriodicBubbleArgParser, self)._check_main_args(args)
+            vars_ = ('rho1', 'rho2', 'mu1', 'mu2')
+            self._check_default(args, vars_, float, allow_none=False)
+            self._check_positive(args, vars_, strict=False, allow_none=False)
+            
+        def _add_graphical_io_args(self):
+            graphical_io = super(PeriodicBubbleArgParser, self)._add_graphical_io_args()
+            graphical_io.add_argument('-pp', '--plot-parameters', action='store_true',
+                    dest='plot_parameters',
+                    help=('Plot the density and viscosity integrals during simulation. '+
+                         'Simulation will stop at each time of interest and '+
+                         'the plot will be updated every specified freq iterations.'))
+            graphical_io.add_argument('-pf', '--plot-freq', type=int, default=100, 
+                    dest='plot_freq',
+                    help='Plotting update frequency in terms of iterations.')
+        
+        def _check_file_io_args(self, args):
+            super(PeriodicBubbleArgParser, self)._check_file_io_args(args)
+            self._check_default(args, 'plot_parameters', bool, allow_none=False)
+            self._check_default(args, 'plot_freq', int, allow_none=False)
+            self._check_positive(args, 'plot_freq', strict=True, allow_none=False)
+            
+        def _setup_parameters(self, args):
+            if (args.ndim not in (2,3)):
+                msg='This example only works for 2D and 3D domains.'
+                self.error(msg)
+
+    parser = PeriodicBubbleArgParser()
+
+    parser.set_defaults(impl='cl', ndim=2, npts=(257,),
+                        box_origin=(0.0,), box_length=(1.0,), 
+                        tstart=0.0, tend=0.50, 
+                        dt=1e-5, cfl=0.5, lcfl=0.125,
+                        dump_freq=0, 
+                        dump_times=(0.0, 0.1, 0.20, 0.30, 0.325, 0.4, 0.45),
+                        rho1=1.0, rho2=10.0, mu1=0.00025, mu2=0.00050)
+
+    parser.run(compute)

hysop/backend/device/opencl/operator/integrate.py

+
+from hysop.tools.decorators import debug
+from hysop.backend.device.opencl.opencl_operator import OpenClOperator, op_apply
+from hysop.operator.base.integrate import IntegrateBase
+
+class OpenClIntegrate(IntegrateBase, OpenClOperator):
+
+    @debug
+    def __init__(self, **kwds):
+        super(OpenClIntegrate, self).__init__(**kwds)
+    
+    @debug
+    def get_field_requirements(self):
+        # force 0 ghosts for the reduction (pyopencl reduction kernel)
+        requirements = super(OpenClIntegrate, self).get_field_requirements()
+        topo, req = requirements.get_input_requirement(self.field)
+        req.max_ghosts = (0,)*self.field.dim
+        return requirements
+    
+    @debug
+    def setup(self, work):
+        self.sum_kernels = tuple(self.dF.backend.sum(a=self.dF[i], build_kernel_launcher=True, async=True)
+                                    for i in xrange(self.dF.nb_components))
+    
+    @op_apply
+    def apply(self, **kwds):
+        queue = self.cl_env.default_queue
+        evts = ()
+        for knl in self.sum_kernels:
+            evt = knl(queue=queue)
+            evts += (evt,)
+        for (i,evt) in enumerate(evts):
+            evt.wait()
+            Pi = self.sum_kernel[i].out.get()[0]
+            if (self.scaling_coeff[i] is None):
+                self.scaling_coeff[i] = 1.0 / Pi
+            self.parameter[i] = self.scaling_coeff[i] * Pi

hysop/core/memory/memory_request.py

@@ -452,7 +452,7 @@ class MultipleOperatorMemoryRequests(object):
                     precision = ''
                 ss+= '\n {}{}:'.format(backend.full_tag, precision)
                 ss+= template.format(*titles, **sizes)
-                for op in sorted(backend_requests.keys(), key=lambda op: op.name):
+                for op in sorted(backend_srequests.keys(), key=lambda op: op.name):
                     sop_reqs = backend_srequests[op]
                     for sreq in sop_reqs:
                         ss+= template.format(*sreq, **sizes)

hysop/defaults.py

 
-from hysop.fields.default_fields import VelocityField, VorticityField
+from hysop.fields.default_fields import VelocityField, VorticityField, \
+                                        DensityField, ViscosityField
 from hysop.parameters.default_parameters import TimeParameters, \
-        EnstrophyParameter, KineticEnergyParameter
+        EnstrophyParameter, KineticEnergyParameter, VolumicIntegrationParameter

hysop/fields/default_fields.py

@@ -26,3 +26,18 @@ def VorticityField(domain, name=None, pretty_name=None,
     return Field(domain=domain, name=name, pretty_name=pretty_name, 
                     nb_components=nb_components, **kwds)
 
+
+def DensityField(domain, name=None, pretty_name=None, **kwds):
+    name        = first_not_None(name, 'rho')
+    pretty_name = first_not_None(pretty_name, greak[16])
+    return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)
+
+
+def ViscosityField(domain, name=None, pretty_name=None, mu=False, **kwds):
+    if mu:
+        name = first_not_None(name, 'mu')
+        pretty_name = first_not_None(pretty_name, greak[11])
+    else:
+        name = first_not_None(name, 'nu')
+        pretty_name = first_not_None(pretty_name, greak[12])
+    return Field(domain=domain, name=name, pretty_name=pretty_name, **kwds)

hysop/operator/base/integrate.py

+
+
+from abc import ABCMeta
+
+from hysop.tools.types       import check_instance, to_tuple, first_not_None
+from hysop.tools.decorators  import debug
+from hysop.fields.continuous_field import Field
+
+from hysop.core.memory.memory_request import MemoryRequest
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.parameters.scalar_parameter import ScalarParameter, TensorParameter
+
+class IntegrateBase(object):
+    """
+    Common implementation interface for field integration.
+    """
+
+    __metaclass__ = ABCMeta
+
+    @debug
+    def __init__(self, field, variables, 
+                    parameter=None, scaling=None, **kwds):
+        """
+        Initialize a Integrate operator base.
+
+        Integrate a field on it compute domain and put the result in a parameter.
+
+        in:  field
+            Possibly as multi-component field that should be integrated.
+        out: parameter
+             P = scaling * integral_V(field)
+             where V is the field domain volume
+             and scaling depends on specified scaling method.
+        
+        parameter
+        ----------
+        field: Field
+            Input continuous field to be integrated.
+        variables: dict
+            dictionary of fields as keys and topologies as values.
+        parameter: ScalarParameter or TensorParameter
+            The output parameter that will contain the integral.
+            Should match field.nb_components. 
+            A default parameter will be created if not specified.
+        scaling: None, float, str or array-like of str, optional
+            Scaling method used after integration.
+            'volumic':   scale by domain size (product of mesh space steps)
+            'normalize': scale by first integration (first value will be 1.0) 
+            Can also be a custom float value of tuple of float values.
+            Defaults to volumic integration.
+        kwds:
+            Extra keywords arguments that will be passed towards implementation 
+            enstrophy operator __init__.
+        """
+        
+        check_instance(field, Field)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+
+        scaling = first_not_None(scaling, 'volumic')
+        if isinstance(scaling, float):
+            scaling = (scaling,)*field.nb_components
+        
+        # Generate parameter if not supplied.
+        if (parameter is None):
+            parameter = VolumicIntegrationParameter(field=field)
+        if (parameter.size != field.nb_components):
+            msg='Expected a parameter of size {} but got a parameter of size {}.'
+            msg=msg.format(field.nb_components, parameter.size)
+            raise RuntimeError(msg)
+        
+        check_instance(parameter, (ScalarParameter, TensorParameter))
+        check_instance(scaling, (str,tuple))
+        if isinstance(scaling, tuple):
+            check_instance(tuple, values=float, size=field.nb_components)
+        
+        input_fields  = { field: variables[field] }
+        output_params = { parameter.name: parameter }
+
+        self.field     = field
+        self.parameter = parameter
+        self.scaling   = scaling
+        self.scaling_coeff = None
+        
+        super(IntegrateBase, self).__init__(input_fields=input_fields, output_params=output_params, **kwds)
+    
+    @debug
+    def discretize(self):
+        if self.discretized:
+            return
+        super(IntegrateBase, self).discretize()
+        dF = self.input_discrete_fields[self.field]
+
+        scaling = self.scaling
+        if (scaling == 'volumic'):
+            scaling_coeff = npw.prod(dF.space_step) / npw.prod(dF.domain.length)
+            scaling_coeff = (scaling_coeff,)*dF.nb_components
+        elif (scaling == 'normalize'):
+            scaling_coeff = [None,]*dF.nb_components
+        elif isinstance(scaling, tuple):
+            scaling_coeff = tuple(scaling)
+        else:
+            msg='Unknown scaling method {}'.format(self.scaling)
+            raise ValueError(msg)
+
+        assert len(scaling_coeff)==dF.nb_components
+
+        self.dF = dF
+        self.scaling_coeff = scaling_coeff
+

hysop/operator/base/min_max.py

@@ -71,7 +71,8 @@ class MinMaxFieldStatisticsBase(object):
                 param = make_param(k, quiet=True)
             else:
                 param = None
-            assert npw.prod(param.shape) == nb_components
+            if (param is not None):
+                assert npw.prod(param.shape) == nb_components
             parameters[k] = param
         return parameters 
 

hysop/operator/directional/diffusion_dir.py

@@ -105,8 +105,10 @@ class DirectionalDiffusion(DirectionalSymbolic):
 
         if isinstance(coeffs, npw.ndarray):
             coeffs = (coeffs,)
-        if isinstance(coeffs, (Field, Parameter)):
+        elif isinstance(coeffs, Parameter):
             coeffs = (coeffs.s,)
+        elif isinstance(coeffs, Field):
+            coeffs = (coeffs.s(*coeffs.domain.frame.vars),)
 
         fields     = to_tuple(fields)
         coeffs     = to_tuple(coeffs)

hysop/operator/integrate.py

+
+
+"""
+@file enstrophy.py
+Enstrophy solver frontend.
+"""
+from hysop.constants         import Implementation
+from hysop.tools.types       import check_instance
+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.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
+from hysop.parameters.scalar_parameter import ScalarParameter, TensorParameter
+from hysop.parameters.default_parameters import VolumicIntegrationParameter
+
+
+class Integrate(ComputationalGraphNodeFrontend):
+    """
+    Interface for integrating fields on their domain
+    Available implementations are: 
+        *OPENCL (gpu based implementation)
+    """
+    
+    @classmethod
+    def implementations(cls):
+        from hysop.backend.device.opencl.operator.integrate import OpenClIntegrate
+        implementations = {
+                Implementation.OPENCL: OpenClIntegrate
+        }
+        return implementations
+    
+    @classmethod
+    def default_implementation(cls):
+        return Implementation.OPENCL
+    
+    @debug
+    def __init__(self, field, variables,
+                parameter=None, scaling=None,
+                implementation=None, base_kwds=None, **kwds):
+        """
+        Initialize a Integrate operator frontend.
+
+        Integrate a field on it compute domain and put the result in a parameter.
+
+        in:  field
+            Possibly as multi-component field that should be integrated.
+        out: parameter
+             P = scaling * integral_V(field)
+             where V is the field domain volume
+             and scaling depends on specified scaling method.
+        
+        parameter
+        ----------
+        field: Field
+            Input continuous field to be integrated.
+        variables: dict
+            dictionary of fields as keys and topologies as values.
+        parameter: ScalarParameter or TensorParameter
+            The output parameter that will contain the integral.
+            Should match field.nb_components. 
+            A default parameter will be created if not specified.
+        scaling: None, float, str or array-like of str, optional
+            Scaling method used after integration.
+            'volumic':   scale by domain size (product of mesh space steps)
+            'normalize': scale by first integration (first value will be 1.0) 
+            Defaults to volumic integration.
+        implementation: Implementation, optional, defaults to None
+            target implementation, should be contained in available_implementations().
+            If None, implementation will be set to default_implementation().
+        base_kwds: dict, optional, defaults to None
+            Base class keywords arguments.
+            If None, an empty dict will be passed.
+        kwds:
+            Extra keywords arguments that will be passed towards implementation 
+            enstrophy operator __init__.
+
+        Notes
+        -----
+        An Integrate operator implementation should at least support 
+        the hysop.operator.base.integrate.IntegrateBase interface.
+        """
+        base_kwds = base_kwds or dict()
+        
+        check_instance(field, Field)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+        check_instance(parameter, (ScalarParameter, TensorParameter), allow_none=True)
+        check_instance(scaling, str, allow_none=True)
+        check_instance(base_kwds, dict, keys=str)
+        
+        # Pregenerate parameter so that we can directly store it in self.
+        if (parameter is None):
+            parameter = VolumicIntegrationParameter(field=field)
+        if (parameter.size != field.nb_components):
+            msg='Expected a parameter of size {} but got a parameter of size {}.'
+            msg=msg.format(field.nb_components, parameter.size)
+            raise RuntimeError(msg)
+        
+        super(Integrate, self).__init__(field=field, variables=variables,
+                parameter=parameter, scaling=scaling, base_kwds=base_kwds, **kwds)

hysop/operators.py

@@ -16,6 +16,7 @@ from hysop.operator.adapt_timestep    import AdaptiveTimeStep
 from hysop.operator.hdf_io            import HDF_Writer, HDF_Reader
 from hysop.operator.custom_symbolic   import CustomSymbolicOperator
 from hysop.operator.parameter_plotter import ParameterPlotter
+from hysop.operator.integrate         import Integrate
 
 from hysop.operator.derivative import SpaceDerivative,       \
                                       MultiSpaceDerivatives, \

hysop/parameters/default_parameters.py

+from hysop.constants import HYSOP_REAL
 from hysop.tools.types import first_not_None
 from hysop.tools.sympy_utils import greak, Greak, subscripts
+from hysop.parameters.tensor_parameter import TensorParameter
 from hysop.parameters.scalar_parameter import ScalarParameter
 
 def TimeParameters(dtype, **kwds):
@@ -16,4 +18,19 @@ def KineticEnergyParameter(name=None, pretty_name=None, **kwds):
     name        = first_not_None(name, 'kinetic_energy')
     pretty_name = first_not_None(pretty_name, u'E\u2096')
     return ScalarParameter(name=name, pretty_name=pretty_name, **kwds)
-    
+
+def VolumicIntegrationParameter(name=None, pretty_name=None, field=None, dtype=None, **kwds):
+    if (field is not None):
+        pretty_name = first_not_None(pretty_name, name, field.pretty_name)
+        name        = first_not_None(name, field.name)
+        dtype       = first_not_None(dtype, field.dtype)
+    dtype = first_not_None(dtype, HYSOP_REAL)
+    name += '_v'
+    pretty_name += u'\u1d65'.encode('utf-8')
+    if (field.nb_components>1):
+        return TensorParameter(name=name, pretty_name=pretty_name, 
+                               shape=(field.nb_components,), 
+                               dtype=dtype, **kwds)
+    else:
+        return ScalarParameter(name=name, pretty_name=pretty_name, 
+                                dtype=dtype, **kwds)

hysop/parameters/tensor_parameter.py

@@ -165,6 +165,9 @@ class TensorParameter(Parameter):
     def _get_shape(self):
         """Get parameter shape."""
         return self._value.shape
+    def _get_size(self):
+        """Get parameter size."""
+        return self._value.size
     def _get_dtype(self):
         """Get parameter dtype."""
         return self._value.dtype
@@ -253,6 +256,7 @@ TensorParameter[name={}, pname={}]
         return ss
 
     shape = property(_get_shape)
+    size = property(_get_size)
     dtype = property(_get_dtype)
     ctype = property(_get_ctype)
     min_value = property(_get_min_value)