From 6327a5ac8a30858380435de971180eb9619418f4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Mon, 1 Dec 2014 11:17:14 +0100
Subject: [PATCH] Add flow around sphere example (first draft with new version)
---
Examples/FlowAroundSphere.py | 373 +++++++++++++++++++++++++++++++++++
1 file changed, 373 insertions(+)
create mode 100644 Examples/FlowAroundSphere.py
diff --git a/Examples/FlowAroundSphere.py b/Examples/FlowAroundSphere.py
new file mode 100644
index 000000000..6e43288b5
--- /dev/null
+++ b/Examples/FlowAroundSphere.py
@@ -0,0 +1,373 @@
+#!/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.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.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.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
+import hysop.tools.io_utils as io
+from hysop.mpi import main_rank, MPI
+from hysop.tools.parameters import Discretization, IO_params
+
+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
+Nx = 129
+Ny = Nz = 65
+g = 2
+boxlength = npw.realarray([10.24, 5.12, 5.12])
+boxorigin = npw.realarray([-2.0, -2.56, -2.56])
+box = Box(dim, 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.sphere import Sphere
+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', isVector=True)
+vorti = Field(domain=box, formula=computeVort,
+ name='Vorticity', isVector=True)
+
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=50)
+
+
+# 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 = IO_params("time_step")
+# 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
+
+# ===== Penalization of the vorticity on a sphere inside the domain =====
+from hysop.operator.penalize_vorticity 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['diffusion'],
+ 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])
+# ========= Monitoring operators =========
+
+iop = IO_params('fields', frequency=100)
+monitors['writer'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=iop)
+
+io_ener = IO_params('energy_enstrophy')
+monitors['energy'] = EnergyEnstrophy(velo, vorti, discretization=topofft,
+ io_params=io_ener)
+
+ref_step = topo_with_ghosts.mesh.space_step
+cbpos = npw.zeros(dim)
+cblength = npw.zeros(dim)
+cbpos[...] = boxorigin[...]
+cbpos += 15 * ref_step
+cblength[...] = boxlength[...]
+cblength -= 30 * ref_step
+cb = ControlBox(box, cbpos, cblength)
+coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
+
+io_forces=IO_params('drag_and_lift')
+monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+ discretization=topo_with_ghosts, cb,
+ obstacles=[sphere], io_params=io_forces)
+
+io_forcesPenal=IO_params('drag_and_lift_penal')
+monitors['forcesPenal'] = DragAndLiftPenal(velo, vorti, coeffForce,
+ discretization=topofft,
+ obstacles=[sphere], factor=[1e8],
+ io_params=io_forcesPenal)
+
+step_dir = ref_step[0]
+io_sliceXY = IO_params('sliceXY', frequency=200)
+thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ lengths=[10.24, 5.12, 4.0 * step_dir])
+monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXY, subset=thickSliceXY,
+ xmfalways=True)
+
+io_sliceXZ = IO_params('sliceXZ', frequency=2000)
+thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
+ lengths=[10.24, 4.0 * step_dir, 5.12])
+monitors['writerSliceXZ'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXZ, subset=thickSliceXZ,
+ xmfalways=True)
+
+io_subBox = IO_params('subBox', frequency=2000)
+subBox = ControlBox(box, origin=[-0.7, -2.0, -2.0], lengths=[8.0, 4.0, 4.0])
+monitors['writerSubBox'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_subBox, subset=subBox,
+ xmfalways=True)
+
+## ========= Monitors discretization=========
+for monit in monitors.values():
+ monit.discretize()
+
+# ========= Setup for all declared operators/monitors =========
+time_setup = MPI.Wtime()
+for ope in op.values():
+ ope.setup()
+for ope in distr.values():
+ ope.setup()
+for monit in monitors.values():
+ monit.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['poisson'].apply(simu) # Poisson
+ op['correction'].apply(simu) # Velo correction
+ monitors['forcesPenal'].apply(simu) # Forces Heloise
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ op['penalVort'].apply(simu) # Vorticity penalization
+ op['stretching'].apply(simu) # Stretching
+ monitors['forces'].apply(simu) # Forces Noca
+ 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()
+# monitors['writerSliceXY'].apply(simu)
+# monitors['writerSliceXZ'].apply(simu)
+# monitors['writerSubBox'].apply(simu)
+ monitors['energy'].apply(simu) # Energy/enstrophy
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ op['dtAdapt'].apply(simu) # Update timestep
+ op['dtAdapt'].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()
--
GitLab