Resolve "Fine to coarse grid filters"

Closes #33 (closed)

examples/multiresolution/scalar_advection.py

+
+import numpy as np
+
+def compute(args):
+
+    from hysop import Field, Box, Simulation, Problem, \
+                      ScalarParameter, MPIParams, IOParams
+    from hysop.constants import Implementation
+    from hysop.operators import Advection, DirectionalAdvection, StrangSplitting, \
+                                LowpassFilter, HDF_Writer
+    from hysop.methods import Remesh, TimeIntegrator, ComputeGranularity, \
+                              Interpolation
+
+    ## Function to compute initial velocity values
+    def init_velocity(data, coords, component=None):
+        if (component is None):
+            for i in xrange(len(data)):
+                data[i][...] = args.velocity[::-1][i]
+        else:
+            assert len(data)==1
+            i, = component
+            data[0][...] = args.velocity[::-1][i]
+
+    ## Function to compute initial scalar values
+    def init_scalar(data, coords):
+        data[0][...] = 1.0
+        for x in coords[0]:
+            data[0][...] *= np.cos(x)
+
+        data[1][...] = 1.0
+        for x in coords[0]:
+            data[1][...] *= np.sin(x)
+
+    # Define domain
+    dim  = args.ndim
+    npts  = args.npts
+    snpts = tuple(4*(x-1)+1 for x in 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())
+    
+    # Define parameters and field (time and analytic field)
+    dt     = ScalarParameter('dt', dtype=args.dtype)
+    velo   = Field(domain=box, name='V', is_vector=True,  dtype=args.dtype)
+    scalar = Field(domain=box, name='S', nb_components=2, dtype=args.dtype)
+    
+    # Setup operator method dictionnary
+    # Advection-Remesh operator discretization parameters
+    method = { 
+               ComputeGranularity:  args.compute_granularity,
+               TimeIntegrator:      args.time_integrator,
+               Remesh:              args.remesh_kernel,
+               Interpolation:       args.interpolation
+    }
+    
+    # Setup implementation specific variables
+    impl = args.impl
+    extra_op_kwds = { 'mpi_params': mpi_params }
+    if (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 method
+        # (already done by HysopArgParser for simplicity)
+        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
+    elif (impl in  (Implementation.PYTHON, Implementation.FORTRAN)):
+        pass
+    else:
+        msg='Unknown implementation \'{}\'.'.format(impl)
+        raise ValueError(msg)
+    
+    # Create the problem we want to solve
+    problem = Problem(method=method)
+    
+    if (impl is Implementation.FORTRAN) or (impl is Implementation.PYTHON):
+        # The fortran scales implementation is a special case.
+        # Here directional advection is a black box.
+        advec = Advection(implementation=Implementation.FORTRAN,
+                name='advec',
+                velocity = velo,       
+                advected_fields = (scalar,),
+                variables = {velo: npts, scalar: snpts},
+                dt = dt, **extra_op_kwds)
+        
+        # Finally insert our advection into the problem
+        problem.insert(advec)
+    else:
+        # Build the directional advection operator
+        # here the cfl determines the maximum number of ghosts
+        advec = DirectionalAdvection(implementation=impl,
+                name='advec',
+                velocity = velo,       
+                velocity_cfl = args.cfl,
+                advected_fields = (scalar,),
+                variables = {velo: npts, scalar: snpts},
+                dt=dt, **extra_op_kwds)
+
+        # Build the directional splitting operator graph
+        splitting = StrangSplitting(splitting_dim=dim, 
+                        order=args.strang_order)
+        splitting.push_operators(advec)
+
+        # Finally insert our splitted advection into the problem
+        problem.insert(splitting)
+    
+    #> Operators to dump all fields
+    io_params = IOParams(filename='fine', frequency=args.dump_freq)
+    df0 = HDF_Writer(name='S',
+                             io_params=io_params,
+                             variables={scalar: snpts},
+                             **extra_op_kwds)
+    io_params = IOParams(filename='coarse', frequency=args.dump_freq)
+    df1 = HDF_Writer(name='V',
+                         io_params=io_params,
+                         variables={velo: npts},
+                         **extra_op_kwds)
+    
+    # Add a writer of input field at given frequency.
+    problem.insert(df0, df1)
+    problem.build(args)
+
+    # 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)
+    
+    # Initialize discrete velocity and scalar field
+    problem.initialize_field(velo,   formula=init_velocity)
+    problem.initialize_field(scalar, formula=init_scalar)
+    
+    # Determine a timestep using the supplied CFL
+    # (velocity is constant for the whole simulation)
+    dx   = problem.get_input_discrete_field(scalar).space_step.min()
+    Vinf = max(abs(vi) for vi in args.velocity)
+    dt0  = (args.cfl*dx)/Vinf
+    if (args.dt is not None):
+        dt0 = min(args.dt, dt0)
+    dt0  = 0.99*dt0
+    
+    # Create a simulation and solve the problem 
+    # (do not forget to specify the dt parameter here)
+    simu = Simulation(start=args.tstart, end=args.tend, 
+                      nb_iter=args.nb_iter,
+                      max_iter=args.max_iter,
+                      times_of_interest=args.dump_times,
+                      dt=dt, dt0=dt0)
+    
+    # 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 ScalarAdvectionArgParser(HysopArgParser):
+        def __init__(self):
+            prog_name = 'scalar_advection'
+            default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(), prog_name)
+
+            description=colors.color('HySoP Scalar Advection Example: ', fg='blue', style='bold')
+            description+='Advect a scalar by a given constant velocity. '
+            description+='\n\nThe advection operator is directionally splitted resulting '
+            description+='in the use of one or more advection-remesh operators per direction.'
+    
+            super(ScalarAdvectionArgParser, self).__init__(
+                 prog_name=prog_name,
+                 description=description,
+                 default_dump_dir=default_dump_dir)
+        
+        def _add_main_args(self):
+            args = super(ScalarAdvectionArgParser, self)._add_main_args()
+            args.add_argument('-vel', '--velocity', type=str,
+                                action=self.split, container=tuple, append=False, convert=float,
+                                dest='velocity',
+                                help='Velocity components.')
+            return args
+        
+        def _check_main_args(self, args):
+            super(ScalarAdvectionArgParser, self)._check_main_args(args)
+            self._check_default(args, 'velocity', tuple, allow_none=False)
+    
+        def _setup_parameters(self, args):
+            if len(args.velocity) == 1:
+                args.velocity *= args.ndim
+
+    parser = ScalarAdvectionArgParser()
+
+    parser.set_defaults(box_origin=(0.0,), box_length=(2*np.pi,), 
+                       tstart=0.0, tend=2*np.pi, npts=(16,),
+                       dump_freq=1, cfl=0.5, velocity=(1.0,), 
+                       ndim=3, compute_precision='fp64')
+
+    parser.run(compute)