From a1ff716e91e2892de18a719271f267edf7c0f310 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Tue, 13 Mar 2018 17:07:32 +0100
Subject: [PATCH] shear layer example rework
---
examples/example_utils.py | 19 +-
examples/shear_layer/shear_layer.py | 295 +++++++++++++++++++---------
2 files changed, 212 insertions(+), 102 deletions(-)
diff --git a/examples/example_utils.py b/examples/example_utils.py
index c9a2cfc7f..e654bbee5 100644
--- a/examples/example_utils.py
+++ b/examples/example_utils.py
@@ -240,7 +240,7 @@ class HysopArgParser(argparse.ArgumentParser):
help='Backend implementation (either python or opencl).')
args.add_argument('-cp', '--compute-precision', type=str, default='fp32',
dest='compute_precision',
- help='Floating-point precision used to compute the Field.')
+ help='Floating-point precision used to discretize the parameters and fields.')
return args
def _check_main_args(self, args):
@@ -528,7 +528,7 @@ class HysopArgParser(argparse.ArgumentParser):
term_io.add_argument('-stdout', '--std-out',
type=str, default='{dpath}/p{size}/{rank}.out',
dest='stdout',
- help='Redirect stdout to this file.' + msg0)
+ help='Redirect stdout to this file. ' + msg0)
term_io.add_argument('-stderr', '--std-err',
type=str, default='{dpath}/p{size}/{rank}.out',
dest='stderr',
@@ -866,7 +866,7 @@ class HysopHelpFormatter(ColorHelpFormatter):
def _format_usage(self, usage, actions, groups, prefix):
def predicate(action):
- blacklist = ('-stdout', '-stderr', '-V', '-D', '-P', '-T', '-h',
+ blacklist = ('-tee', '-stdout', '-stderr', '-V', '-D', '-K', '-P', '-T', '-h',
'--version', '--hardware-info', '--hardware-statistics')
p = not action.option_strings[0].startswith('--opencl')
p &= not action.option_strings[0].startswith('--autotuner')
@@ -877,18 +877,19 @@ class HysopHelpFormatter(ColorHelpFormatter):
usage = super(HysopHelpFormatter, self)._format_usage(usage=usage, actions=actions,
groups=groups, prefix=prefix)
usage = usage.rstrip()
-
+ s=' '*(8+len(self._prog))
opencl_parameters = colors.color('OPENCL_PARAMETERS', fg='green', style='bold')
autotuner_parameters = colors.color('AUTOTUNER_PARAMETERS', fg='green', style='bold')
trace = colors.color('TRACE', fg='green', style='bold')
+ tee = colors.color('RANKS', fg='green', style='bold')
stdout = colors.color('STDOUT', fg='green', style='bold')
stderr = colors.color('STDERR', fg='green', style='bold')
- usage +='\n\t\t[--opencl-{{...}} {}]\n\t\t[--autotuner-{{...}} {}]'.format(
- opencl_parameters, autotuner_parameters)
- usage +='\n\t\t[-stdout {}] [-stderr {}]'.format(stdout, stderr)
- usage +='\n\t\t[-V] [-D] [-P] [-T {}]'.format(trace)
- usage +='\n\t\t[--help] [--version] [--hardware-info] [--hardware-summary]'
+ usage +='\n{s}[--opencl-{{...}} {}]\n{s}[--autotuner-{{...}} {}]'.format(
+ opencl_parameters, autotuner_parameters, s=s)
+ usage +='\n{s}[-tee {}] [-stdout {}] [-stderr {}]'.format(tee, stdout, stderr, s=s)
+ usage +='\n{s}[-V] [-D] [-K] [-P] [-T {}]'.format(trace, s=s)
+ usage +='\n{s}[--help] [--version] [--hardware-info] [--hardware-summary]'.format(s=s)
usage +='\n\n'
return usage
diff --git a/examples/shear_layer/shear_layer.py b/examples/shear_layer/shear_layer.py
index 5042c6e78..3fd1bcba6 100644
--- a/examples/shear_layer/shear_layer.py
+++ b/examples/shear_layer/shear_layer.py
@@ -1,83 +1,95 @@
-from hysop import IO, Domain, Field, Box, Discretization, \
- Simulation, Problem, IOParams
-from hysop.deps import np, sm, sys
-from hysop.symbolic import space_symbols
-from hysop.symbolic.field import curl
-from hysop.tools.sympy_utils import greak
-from hysop.parameters.scalar_parameter import ScalarParameter
-from hysop.parameters.default_parameters import EnstrophyParameter, KineticEnergyParameter
-from hysop.constants import Backend, AutotunerFlags, Implementation, AdvectionCriteria
-from hysop.tools.contexts import printoptions
-
-from hysop.operators import DirectionalAdvection, DirectionalDiffusion, \
- PoissonRotational, Diffusion, AdaptiveTimeStep, Enstrophy, KineticEnergy, \
- MinMaxFieldStatistics
-
-from hysop.methods import StrangOrder, SpaceDiscretization, Remesh, \
- TimeIntegrator, ComputeGranularity, \
- OpenClKernelConfig, OpenClKernelAutotunerConfig
-
-from hysop.numerics.splitting.strang import StrangSplitting
-from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
-from hysop.operators import SpaceDerivative
-
-IO.set_default_path('/tmp/hysop_examples/shear_layer')
-
-
+## HySoP Example: Shear Layer
## See Brown 1995:
-## Performance of under-resolved two dimensional incrompressible flow simulations
-delta = 0.05
-rho = 100 # 30 for the thick layer (case 1), 100 for the thin layer (case 1 and 2).
-viscosity = 0.25*1e-4 # viscosity = 1e-4 for case 1, 0.5*1e-4 for case 2, 0.25*1e-4 for case 3.
-
-## Compute initial vorticity symbolically
-(x,y) = space_symbols[:2]
-u = sm.tanh(rho*(0.25-sm.sqrt((y-0.5)*(y-0.5))))
-v = delta*sm.sin(2*sm.pi*x)
-w = v.diff(x) - u.diff(y)
-W0 = sm.utilities.lambdify((x,y), w)
-
-## Function to compute initial vorticity
-def init_vorticity(data, coords):
- data[0][...] = W0(*coords)
- data[0][np.isnan(data[0])] = 0.0
-
-
-def run(dim=2, npts=512+1, lcfl=0.125, cfl=0.7):
- d3d=(npts,)*dim
- impl = Implementation.OPENCL_CODEGEN
-
- method = { ComputeGranularity: 0,
- SpaceDiscretization: 4,
- TimeIntegrator: RK4 }
-
- if (impl is Implementation.OPENCL_CODEGEN):
- autotuner_config = OpenClKernelAutotunerConfig(
- max_candidates=1,
- autotuner_flag=AutotunerFlags.PATIENT,
-
- prune_threshold=1.2, override_cache=False, verbose=0)
- kernel_config = OpenClKernelConfig(autotuner_config=autotuner_config)
- method.update({OpenClKernelConfig : kernel_config})
-
- box = Box(length=(1,)*dim)
-
- velo = Field(domain=box, name='V', is_vector=True, dtype=np.float64)
- vorti = Field(domain=box, name='W', nb_components=(1 if dim==2 else dim), dtype=np.float64)
+## Performance of under-resolved two dimensional incompressible flow simulations
+
+def compute(args):
+ from hysop import IO, Domain, Field, Box, Discretization, \
+ Simulation, Problem, IOParams
+ from hysop.deps import np, sm, sys
+ from hysop.symbolic import space_symbols
+ from hysop.symbolic.field import curl
+ from hysop.tools.sympy_utils import greak
+ from hysop.parameters.scalar_parameter import ScalarParameter
+ from hysop.constants import Backend, AutotunerFlags, Implementation, AdvectionCriteria
+ from hysop.tools.contexts import printoptions
+
+ from hysop.operators import DirectionalAdvection, DirectionalDiffusion, \
+ PoissonRotational, AdaptiveTimeStep, \
+ MinMaxFieldStatistics
+
+ from hysop.methods import StrangOrder, SpaceDiscretization, Remesh, \
+ TimeIntegrator, ComputeGranularity, \
+ OpenClKernelConfig, OpenClKernelAutotunerConfig
+
+ from hysop.numerics.splitting.strang import StrangSplitting
+ from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
+ from hysop.operators import SpaceDerivative
- kinetic_energy = KineticEnergyParameter(dtype=velo.dtype)
- enstrophy = EnstrophyParameter(dtype=vorti.dtype)
+ # Check that we are in 2D.
+ if (args.dim != 2):
+ msg='This example only works for 2D domains.'
+ raise ValueError(msg)
- simu = Simulation(start=0.0, end=1.20, dt0=1e-4)
- t = simu.t
- dt = simu.dt
-
- diffusion = DirectionalDiffusion(implementation=impl,
- name='diffusion',
- fields=vorti, coeffs=viscosity,
- variables={vorti: d3d}, dt=dt)
-
+ # Define the domain
+ npts = args.npts
+ box = Box(origin=args.box_origin, length=args.box_length, dim=args.ndim)
+
+ # 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 = {}
+ if (impl is Implementation.PYTHON):
+ pass
+ elif (impl is Implementation.OPENCL_CODEGEN):
+ # 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)
+
+ # Get back user paramaters
+ # rho: thickness of the shear layers
+ # delta: the strength of the initial perturbation
+ rho = args.rho
+ delta = args.delta
+
+ # Compute initial vorticity fomula symbolically
+ # and define the function to compute initial vorticity.
+ (x,y) = space_symbols[:2]
+ u = sm.tanh(rho*(0.25-sm.sqrt((y-0.5)*(y-0.5))))
+ v = delta*sm.sin(2*sm.pi*x)
+ w = v.diff(x) - u.diff(y)
+ W0 = sm.utilities.lambdify((x,y), w)
+
+ def init_vorticity(data, coords):
+ data[0][...] = W0(*coords)
+ data[0][np.isnan(data[0])] = 0.0
+
+ # Define parameters and field (time, timestep, viscosity, velocity, vorticity)
+ t = ScalarParameter('t', dtype=args.dtype)
+ dt = ScalarParameter('dt', dtype=args.dtype)
+ nu = ScalarParameter('nu', initial_value=args.nu, const=True, dtype=args.dtype)
+ velo = Field(domain=box, name='V', is_vector=True, dtype=args.dtype)
+ vorti = Field(domain=box, name='W', nb_components=(1 if dim==2 else dim), dtype=args.dtype)
+
+ ### Build the directional operators
+ #> Directional advection
advec = DirectionalAdvection(implementation=impl,
name='advec',
velocity = velo,
@@ -86,42 +98,139 @@ def run(dim=2, npts=512+1, lcfl=0.125, cfl=0.7):
variables = {velo: d3d, vorti: d3d},
method = {Remesh: Remesh.L4_2},
dt=dt)
-
- min_max_U = MinMaxFieldStatistics(name='min_max_U', field=velo,
- Finf=True, implementation=impl, variables={velo:d3d})
- min_max_W = MinMaxFieldStatistics(name='min_max_W', field=vorti,
- Finf=True, implementation=impl, variables={vorti:d3d})
-
+ #> Directional diffusion
+ diffusion = DirectionalDiffusion(implementation=impl,
+ name='diffusion',
+ fields=vorti, coeffs=nu,
+ variables={vorti: d3d}, dt=dt)
+ #> Directional splitting operator graph
splitting = StrangSplitting(splitting_dim=dim,
order=StrangOrder.STRANG_SECOND_ORDER)
splitting.push_operators(advec, diffusion)
-
+
+
+ ### Build standard operators
+ #> Poisson operator to recover the velocity from the vorticity
poisson = PoissonRotational(name='poisson', velocity=velo, vorticity=vorti,
- variables={velo:d3d, vorti: d3d}, projection=None)
- poisson.dump_outputs(fields=(vorti,), frequency=10)
- poisson.dump_outputs(fields=(velo,), frequency=10)
+ variables={velo:d3d, vorti: d3d}, projection=None,
+ implementation=impl)
+ #> We ask to dump the inputs and the outputs of this operator
+ poisson.dump_outputs(fields=(vorti,), frequency=args.dump_freq)
+ poisson.dump_outputs(fields=(velo,), frequency=args.dump_freq)
+ #> 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:d3d})
+ #> 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:d3d})
+ ### Adaptive timestep operator
adapt_dt = AdaptiveTimeStep(dt)
- adapt_dt.push_cfl_criteria(cfl=cfl, Finf=min_max_U.Finf)
- adapt_dt.push_advection_criteria(lcfl=lcfl, Finf=min_max_W.Finf,
+ adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf)
+ 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
+ # accrossed all operators contained in the graph.
+ method = { ComputeGranularity: 0,
+ SpaceDiscretization: 4,
+ TimeIntegrator: RK4 }
problem = Problem(method=method)
problem.insert(poisson, splitting, min_max_U, min_max_W, adapt_dt)
problem.build()
- problem.display()
+ # If a visu_rank was provided, display the graph on the given process rank.
+ 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,
+ dt0=args.dt0, nb_iter=args.nb_iter,
+ t=t, dt=dt)
+ simu.write_parameters(t, dt, filename='parameters.txt', precision=4)
+
+ # Initialize only the vorticity
dfields = problem.input_discrete_fields
dfields[vorti].initialize(formula=init_vorticity)
-
- simu.write_parameters(t, dt, filename='parameters.txt', precision=4)
+ # Finally solve the problem
with printoptions(threshold=10000, linewidth=240,
nanstr='nan', infstr='inf',
formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
problem.solve(simu)
problem.finalize()
-if __name__=='__main__':
- run()
+if __name__=='__main__':
+ from examples.example_utils import HysopArgParser, colors
+
+ class ShearLayerArgParser(HysopArgParser):
+ def __init__(self):
+ prog_name = 'shear_layer'
+ default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(),
+ prog_name)
+
+ description=colors.color('HySoP Shear Layer Example: ', fg='blue', style='bold')
+ description+=colors.color('[Brown 1995]', fg='yellow', style='bold')
+ description+=colors.color('\nPerformance of under-resolved two dimensional '
+ +'incrompressible flow simulation (Navier-Stokes 2D).', fg='yellow')
+ description+='\n'
+ description+='\nThis example focuses on a doubly periodic shear layer to '
+ description+='which a sinusoidal perturbation perpendicular to the '
+ description+='orientation of the shear layers is made.'
+ description+='\n'
+ description+='\nEach of the shear layers rolls up in a single vortex '
+ description+='as the flow evolves.'
+ description+='\n'
+ description+='\nDefault example parameters depends on the chosen case:'
+ description+='\n CASE 0 1 2'
+ description+='\n delta 0.5 0.5 0.5'
+ description+='\n rho 30 100 100'
+ description+='\n visco 1.0e-4 0.5e-4 0.25e-4'
+ description+='\n comment thick thin thin'
+ description+='\n'
+ description+='\nSee the original paper at '
+ description+='http://crd.lbl.gov/assets/pubs_presos/underIIJCP.pdf.'
+
+ super(ShearLayerArgParser, self).__init__(
+ prog_name=prog_name,
+ description=description,
+ default_dump_dir=default_dump_dir)
+
+ def _add_main_args(self):
+ args = super(ShearLayerArgParser, self)._add_main_args()
+ args.add_argument('-c', '--case', type=int,
+ dest='case',
+ help='Set (rho,delta,nu) to case-specific default values.')
+ args.add_argument('-rho', '--layer-thickness', type=float,
+ dest='rho',
+ help='Initial thickness of the layers.')
+ args.add_argument('-delta', '--initial-perturbation-stength', type=float,
+ dest='delta',
+ help='Initial perturbation strengh.')
+ args.add_argument('-nu', '--viscosity', type=float,
+ dest='nu',
+ help='Viscosity coefficient.')
+ return args
+ def _check_main_args(self, args):
+ super(ShearLayerArgParser, self)._check_main_args(args)
+ self._check_default(args, 'c', int, allow_none=False)
+ self._check_default(args, 'rho', float, allow_none=False)
+ self._check_default(args, 'delta', float, allow_none=False)
+ self._check_default(args, 'nu', float, allow_none=False)
+
+
+ parser = ShearLayerArgParser()
+
+ parser.set_defaults(impl='python', dim=2, npts=256,
+ box_origin=(0.0,), box_length=(1.0,),
+ tstart=0.0, tend=1.20,
+ nb_iter=None, nb_max_iter=None, dt=None,
+ dt0=1e-4, cfl=0.5, lcfl=0.125,
+ case=0, nu=None, rho=None, delta=None,
+ dump_freq=5)
+
+ parser.run(compute)
--
GitLab