add example for scalability tests

Examples/FlowAroundSphere_DNS.py

+#!/usr/bin/python
+
+"""
+Taylor Green 3D : see paper van Rees 2011.
+
+All parameters are set and defined in python module dataTG.
+
+"""
+
+from hysop import Box, IOParams, IO
+from hysop.f2hysop import fftw2py
+import numpy as np
+import cPickle
+#from scitools.NumPyDB import NumPyDB_cPickle as hysopPickle
+from hysop.fields.continuous import Field
+from hysop.fields.variable_parameter import VariableParameter
+from hysop.mpi.topology import Cartesian
+from hysop.operator.advection import Advection
+from hysop.operator.stretching import Stretching
+from hysop.operator.absorption_BC import AbsorptionBC
+from hysop.operator.poisson import Poisson
+from hysop.operator.diffusion import Diffusion
+from hysop.operator.adapt_timestep import AdaptTimeStep
+from hysop.operator.redistribute_intra import RedistributeIntra
+from hysop.operator.hdf_io import HDF_Writer
+from hysop.operator.energy_enstrophy import EnergyEnstrophy
+from hysop.operator.profiles import Profiles
+from hysop.operator.residual import Residual
+from hysop.problem.simulation import Simulation
+from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\
+    Remesh, Support, Splitting, dtCrit, SpaceDiscretisation
+from hysop.numerics.integrators.runge_kutta2 import RK2 as RK2
+from hysop.numerics.integrators.runge_kutta3 import RK3 as RK3
+from hysop.numerics.integrators.runge_kutta4 import RK4 as RK4
+from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+from hysop.numerics.interpolation import Linear
+from hysop.numerics.remeshing import L6_4 as rmsh
+from hysop.tools.profiler import Profiler, FProfiler
+import hysop.tools.io_utils as io
+import hysop.tools.numpywrappers as npw
+from hysop.mpi import main_rank, MPI
+from hysop.tools.parameters import Discretization
+
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+
+# ====== pi constant and trigonometric functions ======
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+# ====== Flow constants =======
+uinf = 1.0
+VISCOSITY = 1. / 300.
+
+# ======= Domain =======
+dim = 3
+#-------- choose resolution -----
+#Nz = 129
+Nz = 257
+#Nz = 513
+Ny = Nx = 129
+#Ny = Nx = 257
+#--------------------------------
+g = 2
+
+boxorigin = npw.asrealarray([-2.56, -2.56, -2.56])
+#-------- chose domain lenght --------
+#boxlength = npw.asrealarray([5.12, 5.12, 5.12])
+boxlength = npw.asrealarray([5.12, 5.12, 10.24])
+#boxlength = npw.asrealarray([5.12, 5.12, 20.48])
+#-------------------------------------
+box = Box(length=boxlength, origin=boxorigin)
+
+# A global discretization with ghost points
+d3dg = Discretization([Nx, Ny, Nz], [g, g, g])
+# A global discretization, without ghost points
+d3d = Discretization([Nx, Ny, Nz])
+
+# ====== Sphere inside the domain ======
+RADIUS = 0.5
+pos = [0., 0., 0.]
+from hysop.domain.subsets import Sphere, HemiSphere
+sphere = Sphere(origin=pos, radius=RADIUS, parent=box)
+
+
+# ======= Function to compute initial velocity  =======
+def computeVel(res, x, y, z, t):
+    res[0][...] = uinf
+    res[1][...] = 0.
+    res[2][...] = 0.
+    return res
+
+
+# ======= Function to compute initial vorticity =======
+def computeVort(res, x, y, z, t):
+    res[0][...] = 0.
+    res[1][...] = 0.
+    res[2][...] = 0.
+    return res
+
+#  ====== Time-dependant required-flowrate (Variable Parameter) ======
+def computeFlowrate(simu):
+    # === Time-dependant flow rate ===
+    t = simu.tk
+    Tstart = 3.0
+    flowrate = np.zeros(3)
+    flowrate[0] = uinf * box.length[1] * box.length[2]
+    if t >= Tstart and t <= Tstart + 1.0:
+        flowrate[1] = sin(pi * (t - Tstart)) * \
+                      box.length[1] * box.length[2]
+    # === Constant flow rate ===
+#    flowrate = np.zeros(3)
+#    flowrate[0] = uinf * box.length[1] * box.length[2]
+    return flowrate
+
+
+# ======= Fields =======
+velo = Field(domain=box, formula=computeVel,
+             name='Velocity', is_vector=True)
+vorti = Field(domain=box, formula=computeVort,
+              name='Vorticity', is_vector=True)
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=300.0, timeStep=0.0125, iterMax=100)
+
+
+# Adaptative timestep method : dt = min(values(dtCrit))
+# where dtCrit is a list of criterions on which the computation
+# of the adaptative time step is based
+# ex : dtCrit = ['gradU', 'cfl', 'stretch'], means :
+# dt = min (dtAdv, dtCfl, dtStretch), where dtAdv is equal to LCFL / |gradU|
+# For dtAdv, the possible choices are the following:
+# 'vort' (infinite norm of vorticity) : dtAdv = LCFL / |vort|
+# 'gradU' (infinite norm of velocity gradient), dtAdv = LCFL / |gradU|
+# 'deform' (infinite norm of deformation tensor),
+# dtAdv = LCFL / (0.5(gradU + gradU^T))
+op = {}
+iop = IOParams("time_step", fileformat=IO.ASCII)
+# Default topology (i.e. 3D, with ghosts)
+topo_with_ghosts = box.create_topology(d3dg)
+
+
+op['dtAdapt'] = AdaptTimeStep(velo, vorti, simulation=simu,
+                              discretization=topo_with_ghosts,
+                              method={TimeIntegrator: RK3,
+                                      SpaceDiscretisation: FD_C_4,
+                                      dtCrit: ['gradU', 'stretch']},
+                              io_params=iop,
+                              lcfl=0.125,
+                              cfl=0.5)
+
+op['advection'] = Advection(velo, vorti,
+                            discretization=d3d,
+                            method={Scales: 'p_M6',
+                                    Splitting: 'classic'}
+                            )
+
+op['stretching'] = Stretching(velo, vorti,
+                              discretization=topo_with_ghosts)
+
+op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
+                            discretization=d3d)
+
+rate = VariableParameter(formula=computeFlowrate)
+op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)
+
+# ===== Discretization of computational operators ======
+for ope in op.values():
+    ope.discretize()
+
+topofft = op['poisson'].discreteFields[vorti].topology
+topoadvec = op['advection'].discreteFields[vorti].topology
+
+# =====  Smooth vorticity absorption at the outlet =====
+op['vort_absorption'] = AbsorptionBC(velo, vorti, discretization=topofft, 
+                                     req_flowrate=rate,
+                                     x_coords_absorp=[1.56, 2.56])
+op['vort_absorption'].discretize()
+
+# =====  Penalization of the vorticity on a sphere inside the domain =====
+from hysop.operator.penalization import PenalizeVorticity
+op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
+                                    discretization=topo_with_ghosts,
+                                    obstacles=[sphere], coeff=1e8,
+                                    method={SpaceDiscretisation: FD_C_4})
+op['penalVort'].discretize()
+
+# ==== Operators to map data between the different computational operators ===
+# (i.e. between topologies)
+distr = {}
+distr['fft2str'] = RedistributeIntra(source=op['poisson'],
+                                     target=op['stretching'],
+                                     variables=[velo, vorti])
+distr['str2fft'] = RedistributeIntra(source=op['stretching'],
+                                     target=op['poisson'],
+                                     variables=[velo, vorti])
+distr['fft2advec'] = RedistributeIntra(source=op['poisson'],
+                                       target=op['advection'],
+                                       variables=[velo, vorti])
+distr['advec2fft'] = RedistributeIntra(source=op['advection'],
+                                       target=op['poisson'],
+                                       variables=[velo, vorti])
+
+# ========= Setup for all declared operators =========
+time_setup = MPI.Wtime()
+for ope in op.values():
+    ope.setup()
+for ope in distr.values():
+    ope.setup()
+
+
+print '[', main_rank, '] total time for setup:', MPI.Wtime() - time_setup
+
+# ========= Fields initialization =========
+# - initialize velo + vort on topostr
+# - penalize vorticity
+# - redistribute topostr --> topofft
+
+time_init = MPI.Wtime()
+ind = sphere.discretize(topofft)
+def initFields():
+    velo.initialize(topo=topo_with_ghosts)
+    vorti.initialize(topo=topo_with_ghosts)
+    op['penalVort'].apply(simu)
+    distr['str2fft'].apply(simu)
+    distr['str2fft'].wait()
+
+initFields()
+print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init
+
+fullseq = []
+
+def run(sequence):
+    op['vort_absorption'].apply(simu)
+    op['poisson'].apply(simu)               # Poisson + correction
+    distr['fft2str'].apply(simu)
+    distr['fft2str'].wait()
+    op['penalVort'].apply(simu)             # Vorticity penalization
+    op['stretching'].apply(simu)            # Stretching
+    distr['str2fft'].apply(simu)
+    distr['str2fft'].wait()
+    op['diffusion'].apply(simu)             # Diffusion
+    distr['fft2advec'].apply(simu)
+    distr['fft2advec'].wait()
+    op['advection'].apply(simu)             # Advection (scales)
+    distr['advec2fft'].apply(simu)
+    distr['advec2fft'].wait()
+    distr['fft2str'].apply(simu)
+    distr['fft2str'].wait()
+    op['dtAdapt'].apply(simu)               # Update timestep
+    op['dtAdapt'].wait()
+
+# ==== Serialize the simulation data of the problem to a "restart" file ====
+def dump(filename):
+    """
+    Serialize some data of the problem to file
+    (only data required for a proper restart, namely fields in self.input
+    and simulation).
+    @param filename : prefix for output file. Real name = filename_rk_N,
+    N being current process number. If None use default value from problem
+    parameters (self.filename)
+    """
+    if filename is not None:
+        filedump = filename + '_rk_' + str(main_rank)
+    db = open(filedump, 'wb')
+    cPickle.dump(simu, db)
+
+# ====== Load the simulation data of the problem from a "restart" file ======
+def restart(filename):
+    """
+    Load serialized data to restart from a previous state.
+    self.input variables and simulation are loaded.
+    @param  filename : prefix for downloaded file.
+    Real name = filename_rk_N, N being current process number.
+    If None use default value from problem
+    parameters (self.filename)
+    """
+    if filename is not None:
+        filedump = filename + '_rk_' + str(main_rank)
+    db = open(filedump, 'r')
+    simu = cPickle.load(db)
+    simu.start = simu.time - simu.timeStep
+    ite = simu.currentIteration
+    simu.initialize()
+    simu.currentIteration = ite
+    print 'simu', simu
+    print ("load ...", filename)
+    return simu
+
+seq = fullseq
+
+simu.initialize()
+doDump = False
+doRestart = False
+dumpFreq = 8000
+io_default=IOParams('restart', fileformat=IO.ASCII)
+dump_filename = io_default.filename
+print dump_filename
+#===== Restart (if needed) =====
+if doRestart:
+    simu = restart(dump_filename)
+    iop_vel = IOParams('velo_00000.h5')
+    velo.hdf_load(topofft, io_params=iop_vel)
+    iop_vort = IOParams('vorti_00000.h5')
+    vorti.hdf_load(topofft, io_params=iop_vort)
+    # Set up for redistribute
+    for ope in distr.values():
+        ope.setup()
+
+# ======= Time loop =======
+total_time = FProfiler("Total")
+solve_time = FProfiler("SolveTime")
+cttime = MPI.Wtime()
+while not simu.isOver:
+    ctime = MPI.Wtime()
+    #    if topofft.rank == 0:
+    #        simu.printState()
+    run(seq)
+    solve_time += MPI.Wtime() - ctime # Mesure le temps d execution d une iteration complete
+    simu.advance()
+simu.finalize()
+total_time += MPI.Wtime() - cttime # Mesure le temps total de la boucle (ne devrait pas etre tres different de 'solve_time')
+
+prof = Profiler(None, box.comm_task)
+prof += solve_time
+prof += total_time
+for ope in op.values():
+    ope.finalize()
+    ope.get_profiling_info() # permet a l operateur de collecter les infos de son operateur discret
+    prof += ope.profiler # on ajoute les info de op a l objet prof
+for v in (velo, vorti):
+    prof += v.profiler  # On recupere aussi les info des variables
+prof.summarize() # on resume le profile fait des moyennes a travers les procesus
+print str(prof) # On affiche les resultats, on peut faire aussi prof.write()

HySoP/setup.py.in

@@ -79,7 +79,6 @@ def create_fortran_extension(name, pyf_file=None, src_dirs=None, sources=None,
     if pyf_file is not None:
         sources.append(pyf_file)
     sources = sort_f90.sort(sources)
-    print "FORTRAN SOURCES : ", sources
     if debug_mode == 0:
         options.append(('F2PY_REPORT_ON_ARRAY_COPY', '1'))
         if os.uname()[0] == 'Linux':
@@ -97,8 +96,6 @@ def create_fortran_extension(name, pyf_file=None, src_dirs=None, sources=None,
             #inc_dir.append(exti)
     #libs += hysop_link_libraries_names
     #libdir += hysop_link_libraries_dirs
-    print libs
-    print libdir
     # we trust cmake for external libraries and
     # add them to linker, without using libraries option
     extra_link_args = hysop_link_libraries
@@ -156,18 +153,6 @@ def create_swig_extension(name, inc_dirs, src_dirs=None, sources=None):
     extra_compile_args = parseCMakeVar("@CXX_FLAGS@")
     extra_link_args = parseCMakeVar("@CXX_LINKER_FLAGS@")
     define_macros = parseCMakeDefines("@CXX_EXTRA_DEFINES@")
-
-    #print "INCLUDE_DIRS=",include_dirs
-    #print "LIBRARIES=",libraries
-    #print "LIBRARY_DIRS=",library_dirs
-    #print "EXTRA_COMPILE_ARGS=",extra_compile_args
-    #print "EXTRA_LINK_ARGS=",extra_link_args
-    #print "DEFINE_MACROS=",define_macros
-    #  To avoid -I -I in compiler call, which results in a bug:
-    #while inc_dir.count('') > 0:
-    #    inc_dir.remove('')
-    #inc_dir.append('@CMAKE_BINARY_DIR@/Modules')
-    
     swig_ext = Extension(name, sources=sources, language='c++',
                          swig_opts=swig_opts,
                          include_dirs=include_dirs,
@@ -330,7 +315,7 @@ config = Configuration(
     author=authors,
     author_email='hysop-members@lists.forge.imag.fr',
     url='https://forge.imag.fr/projects/hysop/',
-    license='GNU public license',
+    license='GNU General Public License (GPLv3)',
     package_dir={'': '@CMAKE_SOURCE_DIR@'},
     ext_modules=ext_modules,
     packages=packages,