From 7bb9bea8373c5a999b39972d6a3965aeac98c8d1 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Tue, 24 Jul 2018 11:24:21 +0200
Subject: [PATCH] Add TaylorGreen Fortran(FFTW+Scales)+Python example
---
ci/scripts/test.sh | 1 +
.../taylor_green/taylor_green_cpuFortran.py | 410 ++++++++++++------
2 files changed, 280 insertions(+), 131 deletions(-)
diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index 1d20dbba8..4a24f0291 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -89,6 +89,7 @@ scalar_advection/scalar_advection.py
shear_layer/shear_layer.py
taylor_green/taylor_green_cpu.py
taylor_green/taylor_green.py
+taylor_green/taylor_green_cpuFortran.py
bubble/periodic_bubble.py
bubble/periodic_bubble_levelset.py
bubble/periodic_bubble_levelset_penalization.py
diff --git a/examples/taylor_green/taylor_green_cpuFortran.py b/examples/taylor_green/taylor_green_cpuFortran.py
index 4f28f64f3..0bb386d80 100644
--- a/examples/taylor_green/taylor_green_cpuFortran.py
+++ b/examples/taylor_green/taylor_green_cpuFortran.py
@@ -23,136 +23,284 @@ def init_vorticity(data, coords):
# initial volume averaged enstrophy: 6*pi^3 / (2*pi)^3 = 0.75
-from hysop import Box, Simulation, Problem, MPIParams
-from hysop.defaults import VelocityField, VorticityField, \
- EnstrophyParameter, TimeParameters
-from hysop.constants import Implementation, AdvectionCriteria, HYSOP_REAL, \
- StretchingFormulation
-from hysop.operators import Advection, DirectionalStretching, Diffusion, \
- PoissonRotational, AdaptiveTimeStep, \
- Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
- ParameterPlotter
-from hysop.numerics.odesolvers.runge_kutta import RK2
-from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
- ComputeGranularity, Interpolation, StrangOrder
-# Define the domain
-dim = 3
-npts = (65, )*dim
-box = Box(origin=(0., )*dim, length=(2*pi,)*dim, dim=dim)
-cfl = 0.5
-lcfl = 0.125
-Re = 1600.
-viscosity = (1.0/Re)
-
-# 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 = Implementation.OPENCL
-extra_op_kwds = {'mpi_params': mpi_params}
-method = {}
-
-# Define parameters and field (time, timestep, velocity, vorticity, enstrophy)
-t, dt = TimeParameters(dtype=HYSOP_REAL)
-velo = VelocityField(domain=box, dtype=HYSOP_REAL)
-vorti = VorticityField(domain=box, dtype=HYSOP_REAL)
-enstrophy = EnstrophyParameter(dtype=HYSOP_REAL)
-
-### Build the directional operators
-#> Directional advection
-advec = Advection(implementation=Implementation.FORTRAN,
- name='advec',
- velocity = velo,
- advected_fields = (vorti,),
- #velocity_cfl = cfl,
- variables = {velo: npts, vorti: npts},
- dt=dt, **extra_op_kwds)
-#> Directional stretching + diffusion
-stretch = DirectionalStretching(implementation=impl,
- name='stretch',
- formulation =StretchingFormulation.CONSERVATIVE,
- velocity = velo,
- vorticity = vorti,
- variables = {velo: npts, vorti: npts},
- dt=dt, **extra_op_kwds)
-#> Directional splitting operator subgraph
-splitting = StrangSplitting(splitting_dim=dim, order=StrangOrder.STRANG_SECOND_ORDER)
-splitting.push_operators(stretch)
-
-diffuse = Diffusion(implementation=Implementation.FORTRAN,
- name='diffuse',
- viscosity = viscosity,
- input_field = vorti,
- variables = {vorti: npts},
- dt=dt, **extra_op_kwds)
-### 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=None,
- implementation=Implementation.FORTRAN, **extra_op_kwds)
-#> We ask to dump 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: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 compute the enstrophy
-enstrophy_op = Enstrophy(name='enstrophy', vorticity=vorti, enstrophy=enstrophy,
- variables={vorti:npts}, implementation=impl, **extra_op_kwds)
-
-### Adaptive timestep operator
-adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True)
-dt_cfl = adapt_dt.push_cfl_criteria(cfl=cfl, Finf=min_max_U.Finf,
- equivalent_CFL=True)
-dt_advec = adapt_dt.push_advection_criteria(lcfl=lcfl, Finf=min_max_W.Finf,
+def compute(args):
+ from hysop import Box, Simulation, Problem, MPIParams, Field
+ from hysop.defaults import VelocityField, VorticityField, \
+ EnstrophyParameter, TimeParameters
+ from hysop.constants import Implementation, AdvectionCriteria, HYSOP_REAL, \
+ StretchingFormulation
+ from hysop.operators import Advection, StaticDirectionalStretching, Diffusion, \
+ PoissonRotational, AdaptiveTimeStep, \
+ Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
+ ParameterPlotter
+ from hysop.numerics.odesolvers.runge_kutta import RK2
+ from hysop.topology.cartesian_topology import CartesianTopology
+ from hysop.tools.parameters import Discretization
+ from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
+ ComputeGranularity, Interpolation, StrangOrder
+ # 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}
+ method = {}
+
+ # Define parameters and field (time, timestep, velocity, vorticity, enstrophy)
+ t, dt = TimeParameters(dtype=HYSOP_REAL)
+ velo = VelocityField(domain=box, dtype=HYSOP_REAL)
+ vorti = VorticityField(domain=box, dtype=HYSOP_REAL)
+ enstrophy = EnstrophyParameter(dtype=HYSOP_REAL)
+ wdotw = Field(domain=box, dtype=HYSOP_REAL, is_vector=False, name="WdotW")
+
+ # Topologies
+ topo_nogh = CartesianTopology(domain=box,
+ discretization=Discretization(npts),
+ mpi_params=mpi_params,
+ cutdirs=[False, False, True])
+
+
+ ### Build the directional operators
+ #> Directional advection
+ advec = Advection(implementation=Implementation.FORTRAN,
+ name='advec',
+ velocity = velo,
+ advected_fields = (vorti,),
+ variables = {velo: npts, vorti: npts},
+ dt=dt, **extra_op_kwds)
+ #> Directional stretching
+ stretch = StaticDirectionalStretching(implementation=impl,
+ name='stretch',
+ formulation = args.stretching_formulation,
+ velocity = velo,
+ vorticity = vorti,
+ variables = {velo: npts, vorti: npts},
+ dt=dt, **extra_op_kwds)
+ #> Directional splitting operator subgraph
+ splitting = StrangSplitting(splitting_dim=dim, order=args.strang_order)
+ splitting.push_operators(stretch)
+ #> Diffusion
+ diffuse = Diffusion(implementation=Implementation.FORTRAN,
+ name='diffuse',
+ viscosity = (1.0/args.Re),
+ input_field = vorti,
+ variables = {vorti: topo_nogh},
+ dt=dt, **extra_op_kwds)
+ ### Build standard operators
+ #> Poisson operator to recover the velocity from the vorticity
+ poisson = PoissonRotational(name='poisson', velocity=velo, vorticity=vorti,
+ variables={velo:topo_nogh, vorti: topo_nogh},
+ projection=args.reprojection_frequency,
+ implementation=Implementation.FORTRAN, **extra_op_kwds)
+ #> We ask to dump 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: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 compute the enstrophy
+ enstrophy_op = Enstrophy(name='enstrophy', vorticity=vorti, enstrophy=enstrophy,
+ variables={vorti:topo_nogh,wdotw:topo_nogh}, implementation=impl, **extra_op_kwds)
+
+ ### Adaptive timestep operator
+ adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True)
+ 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: 0,
- SpaceDiscretization: 4,
- TimeIntegrator: RK2,
- Remesh: Remesh.L4_2,
- Interpolation: Interpolation.LINEAR})
-problem = Problem(method=method)
-problem.insert(poisson, advec, splitting, diffuse, enstrophy_op,
- min_max_U, min_max_W, adapt_dt)
-problem.build()
-
-# Create a simulation
-# (do not forget to specify the t and dt parameters here)
-simu = Simulation(start=0., end=10.01,
- #nb_iter=999999,
- #max_iter=args.max_iter,
- dt0=1e-5, #times_of_interest=args.dump_times,
- t=t, dt=dt)
-simu.write_parameters(t, dt_cfl, dt_advec, dt, enstrophy,
- min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
- filename='parameters.txt', precision=8)
-
-# Attach a field debug dumper if requested
-debug_dumper = None
-# from hysop.tools.debug_dumper import DebugDumper
-# debug_dumper = DebugDumper(
-# path=args.debug_dump_dir,
-# name=args.debug_dump_target,
-# force_overwrite=True, enable_on_op_apply=True)
-
-# Initialize only the vorticity
-problem.initialize_field(vorti, formula=init_vorticity)
-
-# Finally solve the problem
-problem.solve(simu, debug_dumper=debug_dumper)
-
-# Finalize
-problem.finalize()
+ #> Custom operator to plot enstrophy
+ if args.plot_enstrophy:
+ class EnstrophyPlotter(ParameterPlotter):
+ """Custom plotting operator for enstrophy."""
+ def __init__(self, **kwds):
+ import matplotlib.pyplot as plt
+ if all(n==npts[0] for n in npts):
+ snpts='${}^3$'.format(npts[0]-1)
+ else:
+ snpts='x'.join(str(n-1) for n in npts)
+ tag='hysop-{}'.format(snpts)
+ fig = plt.figure(figsize=(30,18))
+ axe0 = plt.subplot2grid((3,2), (0,0), rowspan=3, colspan=1)
+ axe1 = plt.subplot2grid((3,2), (0,1), rowspan=2, colspan=1)
+ axe2 = plt.subplot2grid((3,2), (2,1), rowspan=1, colspan=1)
+ axes = (axe0, axe1, axe2)
+ parameters={axe0:{tag:enstrophy},
+ axe1:{dt_advec.name: dt_advec,
+ dt_cfl.name: dt_cfl,
+ dt.name: dt},
+ axe2:{'CFL*': adapt_dt.equivalent_CFL }}
+ super(EnstrophyPlotter, self).__init__(name='enstrophy_dt',
+ parameters=parameters, fig=fig, axes=axes, **kwds)
+ config='{} {} FD{} PROJECTION_{} {}'.format(args.time_integrator, args.remesh_kernel,
+ args.fd_order, args.reprojection_frequency, args.strang_order)
+ fig = fig.suptitle('HySoP Taylor-Green Example {}\n{}'.format(
+ snpts, config), fontweight='bold')
+ axe0.set_title('Integrated Enstrophy')
+ axe0.set_xlabel('Non-dimensional time', fontweight='bold')
+ axe0.set_ylabel('$\zeta$',
+ rotation=0, fontweight='bold')
+ axe0.set_xlim(args.tstart, args.tend)
+ axe0.set_ylim(0, 26)
+ datadir = os.path.realpath(
+ os.path.join(os.getcwd(), os.path.dirname(__file__)))
+ datadir+='/data'
+ for d in (64,128,256,512):
+ reference=os.path.join(datadir, 'reference_{d}_{d}_{d}.txt'.format(d=d))
+ data = np.loadtxt(reference, usecols=(0,2), dtype=np.float32)
+ axe0.plot(data[:,0], data[:,1]*(1+(d==512)), '--',
+ linewidth=1.0,
+ label='hysop-origin ${}^3$'.format(d) if (d<512) else 'reference-$512^3$')
+ axe0.legend()
+ axe1.set_title('Timesteps (CFL={}, LCFL={})'.format(args.cfl, args.lcfl))
+ axe1.set_xlabel('Non-dimensional time', fontweight='bold')
+ axe1.set_ylabel('Non-dimensional time steps', fontweight='bold')
+ axe1.set_xlim(args.tstart, args.tend)
+ axe1.set_ylim(1e-5, 1e0)
+ axe1.set_yscale('log')
+ axe1.legend()
+ axe2.set_title('Equivalent CFL')
+ axe2.set_xlabel('Non-dimensional time', fontweight='bold')
+ axe2.set_ylabel('CFL*', fontweight='bold')
+ axe2.set_xlim(args.tstart, args.tend)
+ axe2.axhline(y=args.cfl, color='r', linestyle='--')
+ axe2.set_ylim(0., 1.1*args.cfl)
+ plot = EnstrophyPlotter(update_frequency=args.plot_freq,
+ visu_rank=args.visu_rank)
+ else:
+ plot = None
+
+ ## 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, advec, splitting, diffuse, enstrophy_op,
+ min_max_U, min_max_W, adapt_dt, plot)
+ 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,
+ min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
+ filename='parameters.txt', precision=8)
+
+ # Attach a field debug dumper if requested
+ from hysop.tools.debug_dumper import DebugDumper
+ if args.debug_dump_target:
+ debug_dumper = DebugDumper(
+ path=args.debug_dump_dir,
+ name=args.debug_dump_target,
+ force_overwrite=True, enable_on_op_apply=True)
+ else:
+ debug_dumper = None
+
+ # Initialize only the vorticity
+ problem.initialize_field(vorti, formula=init_vorticity)
+
+ # Finally solve the problem
+ problem.solve(simu, debug_dumper=debug_dumper)
+
+ # Finalize
+ problem.finalize()
+
+
+if __name__=='__main__':
+ from examples.example_utils import HysopArgParser, colors
+
+ class TaylorGreenArgParser(HysopArgParser):
+ def __init__(self):
+ prog_name = 'taylor_green'
+ default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(),
+ prog_name)
+
+ description=colors.color('HySoP Taylor-Green Example: ', fg='blue', style='bold')
+ description+=colors.color('[Van Rees 2011] (first part)', fg='yellow', style='bold')
+ description+=colors.color('\nA comparison of vortex and pseudo-spectral methods '
+ +'for the simulation of periodic vortical flows at high Reynolds numbers.',
+ fg='yellow')
+ description+='\n'
+ description+='\nThis example focuses on a validation study for the '
+ description+='hybrid particle-mesh vortex method at Reynolds 1600 for '
+ description+='the 3D Taylor-Green vortex.'
+ description+='\n'
+ description+='\nSee the original paper at '
+ description+='http://vanreeslab.com/wp-content/papercite-data/pdf/rees-2011.pdf.'
+
+ super(TaylorGreenArgParser, self).__init__(
+ prog_name=prog_name,
+ description=description,
+ default_dump_dir=default_dump_dir)
+
+ def _add_main_args(self):
+ args = super(TaylorGreenArgParser, self)._add_main_args()
+ args.add_argument('-Re', '--reynolds-number', type=float,
+ dest='Re',
+ help='Set the simulation Reynolds number.')
+ return args
+
+ def _check_main_args(self, args):
+ super(TaylorGreenArgParser, self)._check_main_args(args)
+ self._check_default(args, 'Re', float, allow_none=False)
+ self._check_positive(args, 'Re', strict=True, allow_none=False)
+
+ def _add_graphical_io_args(self):
+ graphical_io = super(TaylorGreenArgParser, self)._add_graphical_io_args()
+ graphical_io.add_argument('-pe', '--plot-enstrophy', action='store_true',
+ dest='plot_enstrophy',
+ help=('Plot the enstrophy component 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=10,
+ dest='plot_freq',
+ help='Plotting update frequency in terms of iterations.')
+
+ def _check_file_io_args(self, args):
+ super(TaylorGreenArgParser, self)._check_file_io_args(args)
+ self._check_default(args, 'plot_enstrophy', 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 != 3):
+ msg='This example only works for 3D domains.'
+ self.error(msg)
+
+ parser = TaylorGreenArgParser()
+
+ parser.set_defaults(impl='PYTHON', ndim=3, npts=(65,),
+ box_origin=(0.0,), box_length=(2*pi,),
+ tstart=0.0, tend=10.01,
+ dt=1e-5,
+ cfl=0.5, lcfl=0.125,
+ dump_freq=100, dump_times=(),
+ Re=1600.0)
+
+ parser.run(compute)
--
GitLab