improve scalar diffusion example with fortran implementation

examples/scalar_diffusion/scalar_diffusion.py

@@ -15,15 +15,16 @@ def compute(args):
 
     from hysop.constants import Implementation, ComputeGranularity
     from hysop.operator.directional.diffusion_dir import DirectionalDiffusion
+    from hysop.operator.diffusion import Diffusion
 
     from hysop.methods import StrangOrder, TimeIntegrator, SpaceDiscretization
     from hysop.numerics.splitting.strang import StrangSplitting
-    
+
     # Define 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())
@@ -32,75 +33,84 @@ def compute(args):
     dt     = ScalarParameter('dt', dtype=args.dtype)
     nu     = ScalarParameter('nu', initial_value=args.nu, const=True, dtype=args.dtype)
     scalar = Field(domain=box, name='S0', nb_components=1, dtype=args.dtype)
-    
+
     # Diffusion operator discretization parameters
-    method = { 
+    method = {
                ComputeGranularity:    args.compute_granularity,
                SpaceDiscretization:   args.fd_order,
                TimeIntegrator:        args.time_integrator,
     }
-    
+
     # Setup implementation specific variables
     impl = args.impl
     extra_op_kwds = { 'mpi_params': mpi_params }
-    if (impl is Implementation.PYTHON):
-        pass
+
+    # Create the problem we want to solve and insert our
+    # directional splitting subgraph.
+    # Add a writer of input field at given frequency.
+    problem = Problem(method=method)
+
+    if (impl is Implementation.FORTRAN):
+        # Build the directional diffusion operator
+        diffusion = Diffusion(implementation=impl,
+                              name='diffusion', dt=dt,
+                              Fin=scalar, nu=nu,
+                              variables = {scalar: npts},
+                              **extra_op_kwds)
+        problem.insert(diffusion)
     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, 
+        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
+
+        # Build the directional diffusion operator
+        diffusion = DirectionalDiffusion(implementation=impl,
+                                         name='diffusion', dt=dt,
+                                         fields=scalar, coeffs=nu,
+                                         variables = {scalar: npts},
+                                         **extra_op_kwds)
+
+        # Build the directional splitting operator graph
+        splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
+        splitting.push_operators(diffusion)
+
+        problem.insert(splitting)
     else:
         msg='Unknown implementation \'{}\'.'.format(impl)
         raise ValueError(msg)
-    
-    # Build the directional diffusion operator
-    diffusion = DirectionalDiffusion(implementation=impl,
-                name='diffusion', dt=dt,
-                fields=scalar, coeffs=nu,
-                variables = {scalar: npts},
-                **extra_op_kwds)
-    
-    # Build the directional splitting operator graph
-    splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
-    splitting.push_operators(diffusion)
-    
-    # Create the problem we want to solve and insert our 
-    # directional splitting subgraph.
-    # Add a writer of input field at given frequency.
-    problem = Problem(method=method)
-    problem.insert(splitting)
+
     problem.dump_inputs(fields=scalar, filename='S0', frequency=args.dump_freq, **extra_op_kwds)
     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 scalar field
     problem.initialize_field(scalar, formula=init_scalar)
-    
-    # Create a simulation and solve the problem 
+
+    # Create a simulation and solve the problem
     # (do not forget to specify the dt parameter here)
-    simu = Simulation(start=args.tstart, end=args.tend, 
+    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=args.dt)
-    
-    # Finally solve the problem 
+
+    # Finally solve the problem
     problem.solve(simu, dry_run=args.dry_run)
-    
+
     # Finalize
     problem.finalize()
 
@@ -117,7 +127,7 @@ if __name__=='__main__':
             description+='Diffuse a scalar with a given diffusion coefficient. '
             description+='\n\nThe diffusion operator is directionally splitted resulting '
             description+='in the use of one or more diffusion operators per direction.'
-    
+
             super(ScalarDiffusionArgParser, self).__init__(
                  prog_name=prog_name,
                  description=description,
@@ -136,8 +146,8 @@ if __name__=='__main__':
 
     parser = ScalarDiffusionArgParser()
 
-    parser.set_defaults(box_origin=(0.0,), box_length=(1.0,), 
-                       tstart=0.0, tend=1.0, nb_iter=500, 
+    parser.set_defaults(box_origin=(0.0,), box_length=(1.0,),
+                       tstart=0.0, tend=1.0, nb_iter=500,
                        dump_freq=5, nu=0.01, impl='cl')
 
     parser.run(compute)