From 1169d73d58f3723f76d9ca8d6345cbd73f02713e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Chlo=C3=A9=20Mimeau?= <Chloe.Mimeau@imag.fr> Date: Thu, 29 Jan 2015 11:09:06 +0100 Subject: [PATCH] renaming --- Examples/FlowAroundSphere.py | 8 +- Examples/{ => TaylorGreen}/TaylorGreen3D.py | 0 .../{ => TaylorGreen}/TaylorGreen3D_GPU.py | 0 .../{ => TaylorGreen}/TaylorGreen3D_debug.py | 10 +- .../TaylorGreen/TaylorGreen3D_debug_filter.py | 299 ++++++++++++++++++ 5 files changed, 310 insertions(+), 7 deletions(-) rename Examples/{ => TaylorGreen}/TaylorGreen3D.py (100%) rename Examples/{ => TaylorGreen}/TaylorGreen3D_GPU.py (100%) rename Examples/{ => TaylorGreen}/TaylorGreen3D_debug.py (95%) create mode 100644 Examples/TaylorGreen/TaylorGreen3D_debug_filter.py diff --git a/Examples/FlowAroundSphere.py b/Examples/FlowAroundSphere.py index 6a5afdc17..4a58fdc61 100644 --- a/Examples/FlowAroundSphere.py +++ b/Examples/FlowAroundSphere.py @@ -49,10 +49,10 @@ VISCOSITY = 1. / 300. # ======= Domain ======= dim = 3 -Nx = 513 -Ny = Nz = 257 -#Nx = 257 -#Ny = Nz = 129 +#Nx = 513 +#Ny = Nz = 257 +Nx = 257 +Ny = Nz = 129 g = 2 boxlength = npw.asrealarray([10.24, 5.12, 5.12]) boxorigin = npw.asrealarray([-2.0, -2.56, -2.56]) diff --git a/Examples/TaylorGreen3D.py b/Examples/TaylorGreen/TaylorGreen3D.py similarity index 100% rename from Examples/TaylorGreen3D.py rename to Examples/TaylorGreen/TaylorGreen3D.py diff --git a/Examples/TaylorGreen3D_GPU.py b/Examples/TaylorGreen/TaylorGreen3D_GPU.py similarity index 100% rename from Examples/TaylorGreen3D_GPU.py rename to Examples/TaylorGreen/TaylorGreen3D_GPU.py diff --git a/Examples/TaylorGreen3D_debug.py b/Examples/TaylorGreen/TaylorGreen3D_debug.py similarity index 95% rename from Examples/TaylorGreen3D_debug.py rename to Examples/TaylorGreen/TaylorGreen3D_debug.py index 27f26dd8d..203a4a12f 100644 --- a/Examples/TaylorGreen3D_debug.py +++ b/Examples/TaylorGreen/TaylorGreen3D_debug.py @@ -37,7 +37,7 @@ sin = np.sin VISCOSITY = 1. / 1600. # ======= Domain ======= dim = 3 -Nx = Ny = Nz = 65 +Nx = Ny = Nz = 257 g = 2 box = Box(length=[2.0 * pi, 2.0 * pi, 2.0 * pi]) @@ -70,7 +70,7 @@ vorti = Field(domain=box, formula=computeVort, # ========= Simulation setup ========= -simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=50) +simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=10000000) # Adaptative timestep method : dt = min(values(dtCrit)) # where dtCrit is a list of criterions on which the computation @@ -135,7 +135,7 @@ distr['str2fft'] = RedistributeIntra(source=op['stretching'], variables=[vorti]) # ## ========= Monitoring operators ========= -iop = IOParams('TG_io', frequency=100) +iop = IOParams('TG_io', frequency=200) writer = HDF_Writer(variables={velo: topofft, vorti: topofft}, io_params=iop) @@ -177,11 +177,15 @@ fullseq = [] def run(sequence): op['advection'].apply(simu) +# print 'enstrophy 1' +# energy.apply(simu) distr['adv2str'].apply(simu) distr['adv2str'].wait() op['stretching'].apply(simu) distr['str2fft'].apply(simu) distr['str2fft'].wait() +# print 'enstrophy 2' +# energy.apply(simu) op['diffusion'].apply(simu) op['poisson'].apply(simu) energy.apply(simu) diff --git a/Examples/TaylorGreen/TaylorGreen3D_debug_filter.py b/Examples/TaylorGreen/TaylorGreen3D_debug_filter.py new file mode 100644 index 000000000..bac8348fa --- /dev/null +++ b/Examples/TaylorGreen/TaylorGreen3D_debug_filter.py @@ -0,0 +1,299 @@ +#!/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 +from hysop.f2py import fftw2py +import numpy as np +from hysop.fields.continuous import Field +from hysop.operator.advection import Advection +from hysop.operator.stretching import Stretching +from hysop.operator.poisson import Poisson +from hysop.operator.filtering import Filtering +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.problem.simulation import Simulation +from hysop.methods_keys import Scales, TimeIntegrator,\ + Splitting, dtCrit, SpaceDiscretisation +from hysop.numerics.integrators.runge_kutta3 import RK3 as RK3 +from hysop.numerics.finite_differences import FD_C_4 +from hysop.mpi import main_rank, MPI +from hysop.tools.parameters import Discretization, IOParams + +print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) =========" + +# ====== pi constant and trigonometric functions ====== +pi = np.pi +cos = np.cos +sin = np.sin + +VISCOSITY = 1. / 1600. +# ======= Domain ======= +dim = 3 +Nx = Ny = Nz = 129 +g = 2 +box = Box(length=[2.0 * pi, 2.0 * pi, 2.0 * pi]) + +# 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]) +# Default topology (i.e. 3D, with ghosts) +topo_with_ghosts = box.create_topology(d3dg) + +# ======= Function to compute TG velocity ======= +def computeVel(res, x, y, z, t): + res[0][...] = sin(x) * cos(y) * cos(z) + res[1][...] = - cos(x) * sin(y) * cos(z) + res[2][...] = 0. + return res + +# ======= Function to compute reference vorticity ======= +def computeVort(res, x, y, z, t): + res[0][...] = - cos(x) * sin(y) * sin(z) + res[1][...] = - sin(x) * cos(y) * sin(z) + res[2][...] = 2. * sin(x) * sin(y) * cos(z) + return res + +# ======= Function to compute reference vorticity ======= +def computeVortFilt(res, x, y, z, t): + res[0][...] = 0. + res[1][...] = 0. + res[2][...] = 0. + return res + +# ======= Fields ======= +velo = Field(domain=box, formula=computeVel, + name='Velocity', is_vector=True) +vorti = Field(domain=box, formula=computeVort, + name='Vorticity', is_vector=True) +vortiFiltered = Field(domain=box, formula=computeVortFilt, + name='VorticityFiltered', is_vector=True) + + +# ========= Simulation setup ========= +simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=10000000) + +# 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") + +op['dtAdapt'] = AdaptTimeStep(velo, vorti, simulation=simu, + discretization=topo_with_ghosts, + method={TimeIntegrator: RK3, + SpaceDiscretisation: FD_C_4, + dtCrit: ['deform', 'cfl', '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) + +op['poisson'] = Poisson(velo, vorti, discretization=d3d, projection=1) + +op['filtering'] = Filtering(vorti, vortiFiltered, discretization=d3d) + +# ===== Discretization of computational operators ====== +for ope in op.values(): + ope.discretize() + +topofft = op['poisson'].discreteFields[vorti].topology +topoadvec = op['advection'].discreteFields[vorti].topology + +# ==== Operators to map data between the different computational operators === +# (i.e. between topologies) +distr = {} +distr['adv2str'] = RedistributeIntra(source=op['advection'], + target=op['stretching'], + variables=[velo, vorti]) +distr['str2adv'] = RedistributeIntra(source=op['stretching'], + target=op['advection'], + variables=[velo, vorti]) +distr['fft2str'] = RedistributeIntra(source=op['poisson'], + target=op['stretching'], + variables=[velo]) +distr['fft2str2'] = RedistributeIntra(source=op['poisson'], + target=op['stretching'], + variables=[velo, vorti]) +distr['str2fft'] = RedistributeIntra(source=op['stretching'], + target=op['diffusion'], + variables=[vorti]) +# ## ========= Monitoring operators ========= + +iop = IOParams('TG_io', frequency=100) +writer = HDF_Writer(variables={velo: topofft, vorti: topofft}, + io_params=iop) + +io_ener=IOParams('energy_enstrophy') +energy = EnergyEnstrophy(velo, vorti, discretization=topofft, + io_params=io_ener) + +io_enerFilt=IOParams('energy_enstrophy_filt') +energyFilt = EnergyEnstrophy(velo, vortiFiltered, discretization=topofft, + io_params=io_enerFilt) + +writer.discretize() +energy.discretize() +energyFilt.discretize() + +# ========= Setup for all declared operators ========= +time_setup = MPI.Wtime() +for ope in op.values(): + ope.setup() +for ope in distr.values(): + ope.setup() +writer.setup() +energy.setup() +energyFilt.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() + + +def initFields_mode1(): + velo.initialize(topo=topoadvec) + vorti.initialize(topo=topoadvec) + vortiFiltered.initialize(topo=topofft) + + +initFields_mode1() +print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init + +fullseq = [] + + +def run(sequence): + op['advection'].apply(simu) +# print 'enstrophy 1' +# energy.apply(simu) + distr['adv2str'].apply(simu) + distr['adv2str'].wait() + op['stretching'].apply(simu) + distr['str2fft'].apply(simu) + distr['str2fft'].wait() +# print 'enstrophy 2' +# energy.apply(simu) + op['diffusion'].apply(simu) + op['poisson'].apply(simu) + op['filtering'].apply(simu) + energy.apply(simu) + energyFilt.apply(simu) + distr['fft2str2'].apply(simu) + distr['fft2str2'].wait() + op['dtAdapt'].apply(simu) + op['dtAdapt'].wait() + distr['str2adv'].apply(simu) + distr['str2adv'].wait() + +# ==== Serialize some 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: +# filename = filename +# filedump = filename + '_rk_' + str(main_rank) +# db = parmesPickle(filedump, mode='store') +# db.dump(simu, 'simulation') +# velo.dump(filename, mode='append') +# vorti.dump(filename, mode='append') + +# ## ====== Load some 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: +# filename = filename +# filedump = filename + '_rk_' + str(main_rank) +# db = parmesPickle(filedump, mode='load') +# simu = db.load('simulation')[0] +# simu.start = simu.time - simu.timeStep +# ite = simu.currentIteration +# simu.initialize() +# simu.currentIteration = ite +# print 'simu', simu +# print ("load ...", filename) +# velo.load(filename) +# vorti.load(filename) +# return simu + +seq = fullseq + +simu.initialize() +#doDump = False +#doRestart = False +#dumpFreq = 5000 +#io_default={"filename":'restart'} +#dump_filename = io.Writer(params=io_default).filename +#===== Restart (if needed) ===== +# if doRestart: +# simu = restart(dump_filename) +# # Set up for monitors and redistribute +# for ope in distr.values(): +# ope.setUp() +# for monit in monitors.values(): +# monit.setUp() + +# ======= Time loop ======= +time_run = MPI.Wtime() +while not simu.isOver: + if topofft.rank == 0: + simu.printState() + run(seq) + simu.advance() +# # testdump = simu.currentIteration % dumpFreq is 0 +# # if doDump and testdump: +# # dump(dump_filename) +print '[', main_rank, '] total time for run :', MPI.Wtime() - time_run + +# ======= Finalize ======= +fftw2py.clean_fftw_solver(box.dimension) +for ope in distr.values(): + ope.finalize() +writer.finalize() +energy.finalize() +energyFilt.finalize() -- GitLab