cleaned scalar advection example

examples/scalar_advection/scalar_advection.py

@@ -3,25 +3,17 @@ import numpy as np
 
 def compute(args):
 
-    from hysop import IO, IOParams, Field, Box, \
-                  Simulation, Problem, ScalarParameter, \
-                  MPIParams
-
-    from hysop.constants import AutotunerFlags, Implementation
-    from hysop.operator.directional.advection_dir import DirectionalAdvection
-
-    from hysop.methods import StrangOrder, Remesh, TimeIntegrator, \
-                              OpenClKernelConfig, OpenClKernelAutotunerConfig
-
-    from hysop.numerics.splitting.strang import StrangSplitting
-    from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
+    from hysop import Field, Box, Simulation, Problem, \
+                      ScalarParameter, MPIParams
+    from hysop.constants import Implementation
+    from hysop.operators import DirectionalAdvection, StrangSplitting
+    from hysop.methods import Remesh, TimeIntegrator, ComputeGranularity, \
+                              Interpolation
 
     ## Function to compute initial velocity values
     def init_velocity(data, coords):
         for i,d in enumerate(data):
             d[...] = args.velocity[::-1][i]
-        # data[0][...] = np.arange(data[0].shape[0])[:,None]
-        # data[1][...] = -np.arange(data[0].shape[0])[:,None]
 
     ## Function to compute initial scalar values
     def init_scalar(data, coords):
@@ -30,8 +22,9 @@ def compute(args):
             data[0][...] *= np.cos(x)
 
     # Define domain
+    dim  = args.ndim
     npts = args.npts
-    box  = Box(origin=args.box_origin, length=args.box_length, dim=args.ndim)
+    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,
@@ -42,12 +35,18 @@ def compute(args):
     velo   = Field(domain=box, name='V',  is_vector=True,  dtype=args.dtype)
     scalar = Field(domain=box, name='S0', nb_components=1, dtype=args.dtype)
     
+    # Setup operator method dictionnary
     # Advection-Remesh operator discretization parameters
-    method = {TimeIntegrator: Euler, Remesh: Remesh.L2_1}
+    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 = {}
+    extra_op_kwds = { 'mpi_params': mpi_params }
     if (impl is Implementation.PYTHON):
         pass
     elif (impl is Implementation.OPENCL_CODEGEN):
@@ -59,9 +58,10 @@ def compute(args):
         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)
+        # Configure OpenCL kernel generation and tuning method
+        # (already done by HysopArgParser for simplicity)
         from hysop.methods import OpenClKernelConfig
-        method = {OpenClKernelConfig: args.opencl_kernel_config}
+        method[OpenClKernelConfig] = args.opencl_kernel_config
         
         # Setup opencl specific extra operator keyword arguments
         extra_op_kwds['cl_env'] = cl_env
@@ -78,14 +78,12 @@ def compute(args):
             advected_fields = (scalar,),
             dt = dt,
             variables = {velo: npts, scalar: npts},
-            method = method,
-            mpi_params = mpi_params,
             **extra_op_kwds
         )
 
     # Build the directional splitting operator graph
-    splitting = StrangSplitting(splitting_dim=box.dim, 
-                    order=StrangOrder.STRANG_SECOND_ORDER)
+    splitting = StrangSplitting(splitting_dim=dim, 
+                    order=args.strang_order)
     splitting.push_operators(advec)
     
     # Create the problem we want to solve and insert our 
@@ -96,23 +94,33 @@ def compute(args):
     problem.dump_inputs(fields=scalar, filename='S0', frequency=args.dump_freq)
     problem.build()
 
-    # If a visu_rank was provided, display the graph on the given process rank.
-    problem.display(args.visu_rank)
+    # 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
     dfields = problem.input_discrete_fields
     dfields[velo].initialize(formula=init_velocity)
     dfields[scalar].initialize(formula=init_scalar)
-    # import sys
-    # sys.exit(1)
+    
+    # Determine a timestep using the supplied CFL
+    # (velocity is constant for the whole simulation)
+    dx   = dfields[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)
     
     # Create a simulation and solve the problem 
     # (do not forget to specify the dt parameter here)
-    dx  = dfields[scalar].space_step.min()
-    dt0 = 0.99*args.cfl*dx
     simu = Simulation(start=args.tstart, end=args.tend, 
                         dt0=dt0, dt=dt)
-    problem.solve(simu)
+    
+    # Finally solve the problem 
+    problem.solve(simu, dry_run=args.dry_run)
+    
+    # Finalize
     problem.finalize()
 
 
@@ -136,9 +144,6 @@ if __name__=='__main__':
         
         def _add_main_args(self):
             args = super(ScalarAdvectionArgParser, self)._add_main_args()
-            args.add_argument('-cfl', '--velocity-cfl', type=float,
-                                dest='cfl',
-                                help='CFL constant used to determine timestep.')
             args.add_argument('-vel', '--velocity', type=str,
                                 action=self.split, container=tuple, append=False, convert=float,
                                 dest='velocity',
@@ -147,7 +152,6 @@ if __name__=='__main__':
         
         def _check_main_args(self, args):
             super(ScalarAdvectionArgParser, self)._check_main_args(args)
-            self._check_default(args, 'cfl', float, allow_none=False)
             self._check_default(args, 'velocity', tuple, allow_none=False)
     
         def _setup_parameters(self, args):

examples/shear_layer/shear_layer.py

@@ -9,7 +9,6 @@ import numpy as np
 def compute(args):
     from hysop import Box, Simulation, Problem, MPIParams,\
                       ScalarParameter
-    from hysop.tools.string_utils import vprint_banner, vprint
     from hysop.defaults import VelocityField, VorticityField, \
                                TimeParameters
     from hysop.constants import Implementation, AdvectionCriteria
@@ -34,7 +33,7 @@ def compute(args):
     
     # Setup usual implementation specific variables
     impl = args.impl
-    extra_op_kwds = {'mpi_params': mpi_params}
+    extra_op_kwds = { 'mpi_params': mpi_params }
     if (impl is Implementation.PYTHON):
         method = {}
     elif (impl is Implementation.OPENCL_CODEGEN):
@@ -161,11 +160,7 @@ def compute(args):
     dfields[vorti].initialize(formula=init_vorticity)
 
     # Finally solve the problem 
-    if args.dry_run:
-        vprint()
-        vprint_banner('** Dry-run requested, skipping simulation. **')
-    else:
-        problem.solve(simu)
+    problem.solve(simu, dry_run=args.dry_run)
     
     # Finalize
     problem.finalize()

examples/taylor_green/taylor_green.py

@@ -25,7 +25,6 @@ def init_vorticity(data, coords):
 
 def compute(args):
     from hysop import Box, Simulation, Problem, MPIParams
-    from hysop.tools.string_utils import vprint_banner, vprint
     from hysop.defaults import VelocityField, VorticityField, \
                                EnstrophyParameter, TimeParameters
     from hysop.constants import Implementation, AdvectionCriteria
@@ -166,12 +165,9 @@ def compute(args):
     dfields = problem.input_discrete_fields
     dfields[vorti].initialize(formula=init_vorticity)
     
-    if args.dry_run:
-        vprint()
-        vprint_banner('** Dry-run requested, skipping simulation. **')
-    else:
-        problem.solve(simu)
-
+    # Finally solve the problem 
+    problem.solve(simu, dry_run=args.dry_run)
+    
     # Finalize
     problem.finalize()
 

hysop/problem.py

 
-from hysop import vprint
 from hysop.tools.decorators import debug
 from hysop.core.graph.computational_graph import ComputationalGraph
+from hysop.tools.string_utils import vprint_banner, vprint
 
 class Problem(ComputationalGraph):
 
@@ -26,7 +26,12 @@ class Problem(ComputationalGraph):
         self.setup(work)
 
     @debug
-    def solve(self, simu, max_iter=None, report_freq=10, **kargs):
+    def solve(self, simu, dry_run=False,
+                max_iter=None, report_freq=10, **kargs):
+        if dry_run:
+            vprint()
+            vprint_banner('** Dry-run requested, skipping simulation. **')
+            return
         vprint('\nSolving problem...')
         simu.initialize()
         while not simu.is_over: