From 2313f125751ee19478126a177b6948827768c82c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Chlo=C3=A9=20Mimeau?= <Chloe.Mimeau@imag.fr>
Date: Tue, 27 Oct 2015 15:19:01 +0100
Subject: [PATCH] add SFD stuff (ie: method to get the filtered solution (=the
base flow) needed to perform global linear stability analysis)
---
Examples/FlowAroundSphere_linearized.py | 498 ++++++++++++++++++
HySoP/hysop/operator/discrete/forcing.py | 124 +++++
.../hysop/operator/discrete/low_pass_filt.py | 123 +++++
.../operator/discrete/monitoringPoints.py | 134 +++++
HySoP/hysop/operator/discrete/residual.py | 117 ++++
HySoP/hysop/operator/forcing.py | 155 ++++++
HySoP/hysop/operator/low_pass_filt.py | 153 ++++++
HySoP/hysop/operator/monitoringPoints.py | 68 +++
HySoP/hysop/operator/residual.py | 53 ++
9 files changed, 1425 insertions(+)
create mode 100644 Examples/FlowAroundSphere_linearized.py
create mode 100644 HySoP/hysop/operator/discrete/forcing.py
create mode 100644 HySoP/hysop/operator/discrete/low_pass_filt.py
create mode 100644 HySoP/hysop/operator/discrete/monitoringPoints.py
create mode 100644 HySoP/hysop/operator/discrete/residual.py
create mode 100644 HySoP/hysop/operator/forcing.py
create mode 100644 HySoP/hysop/operator/low_pass_filt.py
create mode 100644 HySoP/hysop/operator/monitoringPoints.py
create mode 100644 HySoP/hysop/operator/residual.py
diff --git a/Examples/FlowAroundSphere_linearized.py b/Examples/FlowAroundSphere_linearized.py
new file mode 100644
index 000000000..d41cf782e
--- /dev/null
+++ b/Examples/FlowAroundSphere_linearized.py
@@ -0,0 +1,498 @@
+#!/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
+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.monitoringPoints import MonitoringPoints
+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
+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, IOParams
+
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+
+# ====== pi constant and trigonometric functions ======
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+# ====== Flow constants =======
+VISCOSITY = 1. / 300.
+
+# ======= Domain =======
+dim = 3
+Nx = 129
+Ny = Nz = 65
+#Nx = 257
+#Ny = Nz = 129
+#Nx = 513
+#Ny = Nz = 257
+#Nx = 1025
+#Ny = Nz = 513
+g = 2
+boxlength = npw.asrealarray([10.24, 5.12, 5.12])
+boxorigin = npw.asrealarray([-2.0, -2.56, -2.56])
+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 set initial velocity BF =======
+def setZeroVel(res, x, y, z, t):
+ res[0][...] = 0.
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+
+# ======= Function to set initial vorticity BF =======
+def setZeroVort(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 =======
+# Base flow
+veloBF = Field(domain=box, formula=setZeroVel,
+ name='Velocity', is_vector=True)
+vortiBF = Field(domain=box, formula=setZeroVort,
+ name='Vorticity', is_vector=True)
+# Small perturbation whose stability is analyzed
+velo = Field(domain=box, formula=setZeroVel,
+ name='Velocity_fluc', is_vector=True)
+vorti = Field(domain=box, formula=setZeroVort,
+ name='Vorticity_fluc', is_vector=True)
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=2000.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")
+# 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['advection1'] = Advection(velo, vortiBF,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
+
+op['advection2'] = Advection(veloBF, vorti,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
+
+op['stretching1'] = Stretching(veloBF, vorti,
+ discretization=topo_with_ghosts)
+
+op['stretching2'] = Stretching(velo, vortiBF,
+ 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['advection2'].discreteFields[vorti].topology
+
+# ===== Smooth vorticity absorption at the outlet =====
+op['vort_absorption'] = AbsorptionBC(velo, vorti, discretization=topofft,
+ req_flowrate=rate,
+ x_coords_absorp=[7.24, 8.24])
+# 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['advec22str1'] = RedistributeIntra(source=op['advection2'],
+ target=op['stretching1'],
+ variables=[veloBF])
+distr['advec12str2'] = RedistributeIntra(source=op['advection1'],
+ target=op['stretching2'],
+ variables=[vortiBF])
+
+distr['fft2str1'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretching1'],
+ variables=[vorti])
+distr['fft2str2'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretching2'],
+ variables=[velo])
+distr['str12fft'] = RedistributeIntra(source=op['stretching1'],
+ target=op['poisson'],
+ variables=[vorti])
+distr['str22fft'] = RedistributeIntra(source=op['stretching2'],
+ target=op['poisson'],
+ variables=[velo])
+
+distr['fft2advec1'] = RedistributeIntra(source=op['poisson'],
+ target=op['advection1'],
+ variables=[velo])
+distr['fft2advec2'] = RedistributeIntra(source=op['poisson'],
+ target=op['advection2'],
+ variables=[vorti])
+distr['advec12fft'] = RedistributeIntra(source=op['advection1'],
+ target=op['poisson'],
+ variables=[velo])
+distr['advec22fft'] = RedistributeIntra(source=op['advection2'],
+ target=op['poisson'],
+ variables=[vorti])
+
+distr['fft2penal'] = RedistributeIntra(source=op['poisson'],
+ target=op['penalVort'],
+ variables=[velo, vorti])
+
+# ========= Monitoring operators =========
+monitors = {}
+#iop = IOParams('fields', frequency=100)
+#monitors['writer'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+# io_params=iop)
+
+io_ener = IOParams('energy_enstrophy')
+monitors['energy'] = EnergyEnstrophy(velo, vorti, discretization=topofft,
+ io_params=io_ener, is_normalized=False)
+
+rk = 0
+if (0.0 in topofft.mesh.coords[2]):
+ rk = main_rank
+io_prof = IOParams('profile_X_axis', frequency=8400, io_leader=rk)
+monitors['profile'] = Profiles(velo, vorti, discretization=topofft,
+ io_params=io_prof, prof_coords=[0.0, 0.0],
+ direction=0, beginMeanComput=25.0)
+
+coordsMonit = [4.0, 0.0, 0.0]
+rk = 0
+if (coordsMonit[2] in topofft.mesh.coords[2]):
+ rk = main_rank
+io_monit = IOParams('monit', frequency=1, io_leader=rk)
+monitors['monit'] = MonitoringPoints(velo, vorti, discretization=topofft,
+ io_params=io_monit,
+ monitPt_coords=coordsMonit)
+
+io_resi = IOParams('residual')
+monitors['residual'] = Residual(vorti, discretization=topofft,
+ io_params=io_resi)
+
+from hysop.domain.control_box import ControlBox
+from hysop.operator.drag_and_lift import MomentumForces, NocaForces
+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(parent=box, origin=cbpos, length=cblength)
+#coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
+
+#io_forces=IOParams('drag_and_lift_NocaII')
+#monitors['forcesNoca'] = NocaForces(velo, vorti,
+# discretization=topo_with_ghosts,
+# nu=VISCOSITY,
+# volume_of_control=cb,
+# normalization=coeffForce,
+# obstacles=[sphere],
+# io_params=io_forces)
+
+#io_forcesPenal=IOParams('drag_and_lift_Mom')
+#monitors['forcesMom'] = MomentumForces(velocity=velo,
+# discretization=topo_with_ghosts,
+# normalization=coeffForce,
+# obstacles=[sphere],
+# penalisation_coeff=[1e8],
+# io_params=io_forcesPenal)
+
+#io_forcesPenal=IOParams('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 = IOParams('sliceXY', frequency=2000)
+thickSliceXY = ControlBox(parent=box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ length=[10.24- step_dir, 5.12- step_dir, 4.0 * step_dir])
+#thickSliceXY = ControlBox(parent=box, origin=[-2.56, -2.56, -2.0 * step_dir],
+# length=[5.12 - step_dir, 5.12 - step_dir, 4.0 * step_dir])
+monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXY, subset=thickSliceXY,
+ xmfalways=True)
+
+io_sliceXZ = IOParams('sliceXZ', frequency=2000)
+thickSliceXZ = ControlBox(parent=box, origin=[-2.0, -2.0 * step_dir, -2.56],
+ length=[10.24- step_dir, 4.0 * step_dir, 5.12- step_dir])
+monitors['writerSliceXZ'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXZ, subset=thickSliceXZ,
+ xmfalways=True)
+
+io_subBox = IOParams('subBox', frequency=5000)
+subBox = ControlBox(parent=box, origin=[-0.7, -2.0, -2.0], length=[8.0, 4.0, 4.0])
+monitors['writerSubBox'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_subBox, subset=subBox,
+ xmfalways=True)
+
+# ========= 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.discretize()
+for monit in monitors.values():
+ monit.setup()
+
+print '[', main_rank, '] total time for setup:', MPI.Wtime() - time_setup
+
+# ========= Fields initialization =========
+# - initialize (veloBF, vortiBF) on topofft
+# - initialize (veloBF, vortiBF) and (velo, vort) on topostr
+# - penalize vorti from velo on topostr
+# - redistribute topostr --> topofft
+
+time_init = MPI.Wtime()
+ind = sphere.discretize(topofft)
+def initFields():
+ veloBF.initialize(topo=topofft)
+ vortiBF.initialize(topo=topofft)
+ veloBF.initialize(topo=topo_with_ghosts)
+ vortiBF.initialize(topo=topo_with_ghosts)
+ velo.initialize(topo=topo_with_ghosts)
+ vorti.initialize(topo=topo_with_ghosts)
+ op['penalVort'].apply(simu)
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].wait()
+
+# /!\ CAREFULL /!\ : INITIALIZATION COMMENTED OUT !!!!!!!!!!!!!!!!!
+#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
+ # monitors['forcesMom'].apply(simu) # Forces Heloise
+ distr['fft2penal'].apply(simu)
+ distr['fft2penal'].wait()
+ op['penalVort'].apply(simu) # Vorticity penalization
+ op['stretching1'].apply(simu) # Stretching 1
+ op['stretching2'].apply(simu) # Stretching 2
+ distr['str12fft'].apply(simu)
+ distr['str12fft'].wait()
+ distr['str22fft'].apply(simu)
+ distr['str22fft'].wait()
+ op['diffusion'].apply(simu) # Diffusion
+ distr['fft2advec1'].apply(simu)
+ distr['fft2advec1'].wait()
+ distr['fft2advec2'].apply(simu)
+ distr['fft2advec2'].wait()
+ op['advection1'].apply(simu) # Advection 1 (scales)
+ op['advection2'].apply(simu) # Advection 2 (scales)
+ distr['advec12fft'].apply(simu)
+ distr['advec12fft'].wait()
+ distr['advec22fft'].apply(simu)
+ distr['advec22fft'].wait()
+ monitors['writerSliceXY'].apply(simu)
+ monitors['writerSliceXZ'].apply(simu)
+ monitors['writerSubBox'].apply(simu)
+ monitors['energy'].apply(simu) # Energy/enstrophy
+ monitors['profile'].apply(simu) # Profile
+ monitors['monit'].apply(simu) # Monitoring points in the wake
+ monitors['residual'].apply(simu) # Vorticity residual
+ distr['fft2penal'].apply(simu)
+ distr['fft2penal'].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 = True
+doRestart = True
+dumpFreq = 8000
+io_default=IOParams('restart')
+dump_filename = io.Writer(io_params=io_default).filename
+#===== Restart (if needed) =====
+if doRestart:
+ # Load data from dumped files
+ simu = restart(dump_filename)
+ iop_vel = IOParams('velo_00000.h5')
+ veloBF.hdf_load(topofft, io_params=iop_vel)
+ iop_vort = IOParams('vorti_00000.h5')
+ vortiBF.hdf_load(topofft, io_params=iop_vort)
+ # Initialize velo and vorti to a small perturbation
+ # equal to baseFlow * 1e-8
+ # (cf Florian Guiho's thesis p.30 eq. (2.13))
+ velo.initialize(topo=topofft)
+ vorti.initialize(topo=topofft)
+ velo.discreteFields[topofft].data[0][...] = \
+ veloBF.discreteFields[topofft].data[0][...] * 1e-8
+ vorti.discreteFields[topofft].data[0][...] = \
+ vortiBF.discreteFields[topofft].data[0][...] * 1e-8
+ # Set up for monitors and redistribute
+ for ope in distr.values():
+ ope.setup()
+ for monit in monitors.values():
+ monit.setup()
+ # Redistribute veloBF and vortBF on topoGhosts
+ distr['advec22str1'].apply(simu)
+ distr['advec22str1'].wait()
+ distr['advec12str2'].apply(simu)
+ distr['advec12str2'].wait()
+
+# ======= 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:
+ print 'dump ...'
+ dump(dump_filename)
+ iop_vel = IOParams('veloFluc')
+ velo.hdf_dump(topofft, io_params=iop_vel)
+ iop_vort = IOParams('vortiFluc')
+ vorti.hdf_dump(topofft, io_params=iop_vort)
+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()
+for monit in monitors.values():
+ monit.finalize()
diff --git a/HySoP/hysop/operator/discrete/forcing.py b/HySoP/hysop/operator/discrete/forcing.py
new file mode 100644
index 000000000..55a00fa9e
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/forcing.py
@@ -0,0 +1,124 @@
+# -*- coding: utf-8 -*-
+"""DOperator implementing the forcing term in the NS,
+ depending on the filtered field
+ (--> computation of base flow).
+
+.. currentmodule:: hysop.operator.discrete.forcing
+
+"""
+from hysop.constants import debug
+from hysop.operator.discrete.discrete import DiscreteOperator
+from hysop.tools.profiler import profile
+
+class Forcing(DiscreteOperator):
+ """Discretized forcing operator.
+ i.e. solve (with an implicit Euler scheme)
+ \f{eqnarray*}
+ \frac{\partial \omega}{\partial t} &=& -\chi(\omega - \bar{\omega}))
+ \Leftrightarrow \omega^{n+1} &=&
+ \frac{\omega^n + \Delta t \chi \bar{\omega^{n+1}}}{1+\Delta t \chi}
+ \f}
+ """
+
+ @debug
+ def __init__(self, strength=None, **kwds):
+ """
+ Parameters
+ ----------
+ strength : strength of the forcing, chosen in the order
+ of the amplification rate related to the unstable flow.
+ **kwds : extra parameters for parent class.
+
+ """
+ super(Forcing, self).__init__(**kwds)
+ topo = self.variables[0].topology
+ for v in self.variables:
+ msg = 'Multiresolution not implemented for penalization.'
+ assert v.topology == topo, msg
+
+ ## variable
+ self.var = self.variables[0]
+ ## forced variable
+ self.varFiltered = self.variables[1]
+ ## strength of the forcing
+ self.strength = strength
+
+ def _apply(self, dt):
+ nbc = self.var.nb_components
+ for d in xrange(nbc):
+ self.var[d][...] += self.varFiltered[d][...] * \
+ (dt * self.strength)
+ self.var[d][...] *= 1.0 / (1.0 + dt * self.strength)
+ print 'forcing: non conservative formulation'
+
+ def apply(self, simulation=None):
+ assert simulation is not None, \
+ "Simulation parameter is required."
+ dt = simulation.timeStep
+ self._apply(dt)
+
+
+class ForcingConserv(Forcing):
+ """Discretized forcing operator.
+ i.e. solve (with an implicit Euler scheme and a CONSERVATIVE formulation)
+ \f{eqnarray*}
+ \frac{\partial \omega}{\partial t} &=& -\chi(\omega - \bar{\omega}))
+ \Leftrightarrow \omega^{n+1} &=&
+ \frac{\omega^n + \Delta t \chi \bar{\omega^{n+1}}}{1+\Delta t \chi}
+ \f}
+ """
+
+ @debug
+ def __init__(self, vorticity, velocity, velocityFilt, curl, **kwds):
+ """
+ Parameters
+ ----------
+ velocity, vorticity, velocityFilt: :class:`~hysop.fields.continuous.Field`
+ curl : :class:`~hysop.operator.differential`
+ internal operator to compute the curl of the forced velocity
+ **kwds : extra parameters for parent class.
+
+ Notes
+ -----
+ velocity and velocityFilt are not modified by this operator.
+ vorticity is an in-out parameter.
+ input and ouput variables of the curl are some local buffers.
+ """
+
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(ForcingConserv, self).__init__(variables=[vorticity,
+ velocity,
+ velocityFilt],
+ **kwds)
+ # Input vector fields
+ self.velocity = velocity
+ self.vorticity = vorticity
+ self.velocityFilt = velocityFilt
+ # warning : a buffer is added for invar variable in curl
+ topo = self.velocity.topology
+ msg = 'Multiresolution not implemented for vort forcing.'
+ assert self.vorticity.topology == topo, msg
+ self._curl = curl
+
+ def _apply(self, dt):
+ # Vorticity forcing
+ # warning : the buff0 array ensures "invar" to be 0
+ # outside the obstacle for the curl evaluation
+ invar = self._curl.invar
+ nbc = invar.nb_components
+ for d in xrange(nbc):
+ invar.data[d][...] = 0.0
+ coeff = (dt * self.strength) / (1.0 + dt * self.strength)
+ for d in xrange(nbc):
+ invar.data[d][...] = \
+ (self.velocityFilt[d][...] -
+ self.velocity[d][...]) * coeff
+ self._curl.apply()
+ for d in xrange(self.vorticity.nb_components):
+ self.vorticity[d][...] += self._curl.outvar[d][...]
+ print 'forcing: conservative formulation'
+
+
+
+
+
diff --git a/HySoP/hysop/operator/discrete/low_pass_filt.py b/HySoP/hysop/operator/discrete/low_pass_filt.py
new file mode 100644
index 000000000..45104de99
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/low_pass_filt.py
@@ -0,0 +1,123 @@
+# -*- coding: utf-8 -*-
+"""Discrete operator for vorticity low-pass filtering
+ (--> computation of base flow).
+.. currentmodule:: hysop.operator.discrete.low_pass_filt_vort
+
+"""
+from hysop.constants import debug
+from hysop.operator.discrete.discrete import DiscreteOperator
+from hysop.tools.profiler import profile
+
+class LowPassFilt(DiscreteOperator):
+ """Discretized vorticity filtering operator.
+ i.e. solve (with an explicit Euler scheme)
+ \f{eqnarray*}
+ \frac{\partial \bar{\omega}}{\partial t} &=& \Omega_c(\omega - \bar{\omega}))
+ \Leftrightarrow \bar{\omega^{n+1}} &=& \bar{\omega^{n}} +
+ \Delta t \Omega_c (\omega^n - \bar{\omega^n})
+ \Leftrightarrow \bar{\omega^{n+1}} &=&
+ \bar{\omega^{n}} (1 - \Delta t \Omega_c) +
+ \Delta t \Omega_c \omega^n
+ \f}
+ """
+
+ @debug
+ def __init__(self, cutFreq=None, **kwds):
+ """
+ Parameters
+ ----------
+ cutFreq : cutting circular frequency corresponding to the half of
+ the eigenfrequency of the flow instability
+ **kwds : extra parameters for parent class.
+
+ """
+ super(LowPassFilt, self).__init__(**kwds)
+ topo = self.variables[0].topology
+ for v in self.variables:
+ msg = 'Multiresolution not implemented for penalization.'
+ assert v.topology == topo, msg
+
+ ## variable
+ self.var = self.variables[0]
+ ## filtered variable
+ self.varFiltered = self.variables[1]
+ ## cutting circular frequency
+ self.cutFreq = cutFreq
+
+ def _apply(self, dt):
+ nbc = self.varFiltered.nb_components
+ coeff = dt * self.cutFreq
+ for d in xrange(nbc):
+ self.varFiltered[d][...] *= (1.0 - coeff)
+ self.varFiltered[d][...] += self.var[d][...] * coeff
+ print 'filtering: non conservative formulation'
+
+ def apply(self, simulation=None):
+ assert simulation is not None, \
+ "Simulation parameter is required."
+ dt = simulation.timeStep
+ self._apply(dt)
+
+
+class LowPassFiltConserv(LowPassFilt):
+ """Discretized vorticity filtering operator.
+ i.e. solve (with an explicit Euler scheme and a CONSERVATIVE formulation)
+ \f{eqnarray*}
+ \frac{\partial \bar{\omega}}{\partial t} &=& \Omega_c(\omega - \bar{\omega}))
+ \Leftrightarrow \bar{\omega^{n+1}} &=& \bar{\omega^{n}} +
+ \nabla \times (\Delta t \Omega_c (u^n - \bar{u^n}))
+ \f}
+ """
+
+ @debug
+ def __init__(self, vorticityFilt, velocity, velocityFilt, curl, **kwds):
+ """
+ Parameters
+ ----------
+ vorticityFilt, vorticity, velocityFilt: :class:`~hysop.fields.continuous.Field`
+ curl : :class:`~hysop.operator.differential`
+ internal operator to compute the curl of the forced velocity
+ **kwds : extra parameters for parent class.
+
+ Notes
+ -----
+ velocity and velocityFilt are not modified by this operator.
+ vorticityFilt is an in-out parameter.
+ input and ouput variables of the curl are some local buffers.
+ """
+
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(LowPassFiltConserv, self).__init__(variables=[vorticityFilt,
+ velocity,
+ velocityFilt],
+ **kwds)
+ # Input vector fields
+ self.vorticityFilt = vorticityFilt
+ self.velocity = velocity
+ self.velocityFilt = velocityFilt
+ # warning : a buffer is added for invar variable in curl
+ topo = self.velocity.topology
+ msg = 'Multiresolution not implemented for vort forcing.'
+ assert self.vorticityFilt.topology == topo, msg
+ assert self.velocityFilt.topology == topo, msg
+ self._curl = curl
+
+ def _apply(self, dt):
+ # Vorticity filtering
+ # warning : the buff0 array ensures "invar" to be 0
+ # outside the obstacle for the curl evaluation
+ invar = self._curl.invar
+ nbc = invar.nb_components
+ for d in xrange(nbc):
+ invar.data[d][...] = 0.0
+ coeff = dt * self.cutFreq
+ for d in xrange(nbc):
+ invar.data[d][...] = \
+ (self.velocity[d][...] -
+ self.velocityFilt[d][...]) * coeff
+ self._curl.apply()
+ for d in xrange(self.vorticityFilt.nb_components):
+ self.vorticityFilt[d][...] += self._curl.outvar[d][...]
+ print 'filtering: conservative formulation'
+
+
diff --git a/HySoP/hysop/operator/discrete/monitoringPoints.py b/HySoP/hysop/operator/discrete/monitoringPoints.py
new file mode 100644
index 000000000..eac5f9154
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/monitoringPoints.py
@@ -0,0 +1,134 @@
+# -*- coding: utf-8 -*-
+"""
+@file monitoringPoints.py
+Print time evolution of flow variables (velo, vorti)
+at a particular monitoring point in the wake
+"""
+from hysop.constants import debug
+from hysop.tools.profiler import profile
+import hysop.tools.numpywrappers as npw
+import scitools.filetable as ft
+import numpy as np
+from hysop.operator.discrete.discrete import DiscreteOperator
+
+
+class MonitoringPoints(DiscreteOperator):
+ """
+ Print time evolution of flow variables at a given position in the wake.
+ """
+ def __init__(self, velocity, vorticity, monitPt_coords, **kwds):
+ """
+ Constructor.
+ @param velocity field
+ @param vorticity field
+ @param monitPt_coords : list of coordinates corresponding
+ to the space location of the different monitoring points in the wake
+ """
+ ## velocity field
+ self.velocity = velocity
+ ## vorticity field
+ self.vorticity = vorticity
+ ## Monitoring point coordinates
+ self.monitPt_coords = monitPt_coords
+
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(MonitoringPoints, self).__init__(variables=[velocity, vorticity],
+ **kwds)
+ topo_v = self.velocity.topology
+ self.shape_v = self.velocity.data[0][topo_v.mesh.iCompute].shape
+ self.space_step = topo_v.mesh.space_step
+ self.length = topo_v.domain.length
+ self.origin = topo_v.domain.origin
+ self.coords = topo_v.mesh.coords
+ self.nbIter = 0
+ ## Normalized flow values at monitoring position (velNorm, vortNorm)
+ self.velNorm = 0.0
+ self.vortNorm = 0.0
+
+ # Is current processor working ? (Is monitPt_coords(z) in z-coords ?)
+ self.is_rk_computing = False
+ s = self._dim - 1
+ if (self.monitPt_coords[s] in self.coords[s]):
+ self.is_rk_computing = True
+
+ def _set_work_arrays(self, rwork=None, iwork=None):
+
+ v_ind = self.velocity.topology.mesh.iCompute
+ w_ind = self.vorticity.topology.mesh.iCompute
+ shape_v = self.velocity.data[0][v_ind].shape
+ shape_w = self.velocity.data[0][w_ind].shape
+ # setup for rwork, iwork is useless.
+ if rwork is None:
+ # --- Local allocation ---
+ if shape_v == shape_w:
+ self._rwork = [npw.zeros(shape_v)]
+ else:
+ self._rwork = [npw.zeros(shape_v), npw.zeros(shape_w)]
+ else:
+ assert isinstance(rwork, list), 'rwork must be a list.'
+ # --- External rwork ---
+ self._rwork = rwork
+ if shape_v == shape_w:
+ assert len(self._rwork) == 1
+ assert self._rwork[0].shape == shape_v
+ else:
+ assert len(self._rwork) == 2
+ assert self._rwork[0].shape == shape_v
+ assert self._rwork[1].shape == shape_w
+
+ def get_work_properties(self):
+
+ v_ind = self.velocity.topology.mesh.iCompute
+ w_ind = self.vorticity.topology.mesh.iCompute
+ shape_v = self.velocity.data[0][v_ind].shape
+ shape_w = self.velocity.data[0][w_ind].shape
+ if shape_v == shape_w:
+ return {'rwork': [shape_v], 'iwork': None}
+ else:
+ return {'rwork': [shape_v, shape_w], 'iwork': None}
+
+ @debug
+ @profile
+ def apply(self, simulation=None):
+ assert simulation is not None, \
+ "Missing simulation value for computation."
+
+ time = simulation.time
+ ite = simulation.currentIteration
+ filename = self._writer.filename #+ '_ite' + format(ite)
+
+ if self.is_rk_computing :
+ self.nbIter += 1
+ vd = self.velocity.data
+ vortd = self.vorticity.data
+ nbc = self.velocity.nb_components
+ tab = [self.monitPt_coords[0], self.monitPt_coords[1],
+ self.monitPt_coords[2]]
+
+ ind = []
+ for d in xrange(nbc):
+ cond = np.where(abs(self.coords[d] - tab[d])
+ < (self.space_step[d] * 0.5))
+ if cond[0].size > 0:
+ ind.append(cond[d][0])
+ else:
+ raise ValueError("Wrong set of coordinates.")
+
+
+ for i in xrange (self.shape_v[0]):
+ self.velNorm = np.sqrt(vd[0][ind[0],ind[1],ind[2]] ** 2 +
+ vd[1][ind[0],ind[1],ind[2]] ** 2 +
+ vd[2][ind[0],ind[1],ind[2]] ** 2)
+ self.vortNorm = np.sqrt(vortd[0][ind[0],ind[1],ind[2]] ** 2 +
+ vortd[1][ind[0],ind[1],ind[2]] ** 2 +
+ vortd[2][ind[0],ind[1],ind[2]] ** 2)
+
+
+ if self._writer is not None and self._writer.do_write(ite) :
+ self._writer.buffer[0, 0] = time
+ self._writer.buffer[0, 1] = self.velNorm
+ self._writer.buffer[0, 2] = self.vortNorm
+ self._writer.write()
+
+
+
diff --git a/HySoP/hysop/operator/discrete/residual.py b/HySoP/hysop/operator/discrete/residual.py
new file mode 100644
index 000000000..498924fd0
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/residual.py
@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+"""
+@file residual.py
+Compute and print the time evolution of the residual
+"""
+from hysop.constants import debug
+from hysop.tools.profiler import profile
+import hysop.tools.numpywrappers as npw
+import scitools.filetable as ft
+import numpy as np
+from hysop.operator.discrete.discrete import DiscreteOperator
+
+
+class Residual(DiscreteOperator):
+ """
+ Compute and print the residual as a function of time
+ """
+ def __init__(self, vorticity, **kwds):
+ """
+ Constructor.
+ @param vorticity field
+ """
+ ## vorticity field
+ self.vorticity = vorticity
+
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(Residual, self).__init__(variables=[vorticity], **kwds)
+ topo_w = self.vorticity.topology
+ self.shape_w = self.vorticity.data[0][topo_w.mesh.iCompute].shape
+ self.space_step = topo_w.mesh.space_step
+ self.length = topo_w.domain.length
+ self.origin = topo_w.domain.origin
+ self.coords = topo_w.mesh.coords
+ self.nbIter = 0
+ ## Global residual
+ self.residual = 0.0
+ # Time stem of the previous iteration
+ self._old_dt = None
+ # Define array to store vorticity field at previous iteration
+ self._vortPrev = [npw.zeros_like(d) for d in self.vorticity.data]
+
+ def _set_work_arrays(self, rwork=None, iwork=None):
+
+ w_ind = self.vorticity.topology.mesh.iCompute
+ shape_w = self.vorticity.data[0][w_ind].shape
+ # setup for rwork, iwork is useless.
+ if rwork is None:
+ # --- Local allocation ---
+ self._rwork = [npw.zeros(shape_w)]
+ else:
+ assert isinstance(rwork, list), 'rwork must be a list.'
+ # --- External rwork ---
+ self._rwork = rwork
+ assert len(self._rwork) == 1
+ assert self._rwork[0].shape == shape_w
+
+ def get_work_properties(self):
+
+ w_ind = self.vorticity.topology.mesh.iCompute
+ shape_w = self.vorticity.data[0][w_ind].shape
+ return {'rwork': [shape_w], 'iwork': None}
+
+ def initialize_vortPrev(self):
+
+ w_ind = self.vorticity.topology.mesh.iCompute
+ for d in xrange(self.vorticity.dimension):
+ self._vortPrev[d][w_ind] = self.vorticity[d][w_ind]
+
+ @debug
+ @profile
+ def apply(self, simulation=None):
+ assert simulation is not None, \
+ "Missing simulation value for computation."
+
+ time = simulation.time
+ ite = simulation.currentIteration
+ dt = simulation.timeStep
+ if self._old_dt is None:
+ self._old_dt = dt
+
+ filename = self._writer.filename
+
+ # Compute on local proc (w^(n+1)-w^n) ** 2
+ local_res = 0.
+ # get the list of computation points (no ghosts)
+ wd = self.vorticity
+ nbc = wd.nb_components
+ w_ind = self.vorticity.topology.mesh.iCompute
+ for i in xrange(nbc):
+ self._rwork[0][...] = (wd[i][w_ind] -
+ self._vortPrev[i][w_ind]) ** 2
+ local_res += npw.real_sum(self._rwork[0])
+
+ # --- Reduce local_res over all proc ---
+ sendbuff = npw.zeros((1))
+ recvbuff = npw.zeros((1))
+ sendbuff[:] = [local_res]
+ #
+ self.vorticity.topology.comm.Allreduce(sendbuff, recvbuff)
+
+ # Update global residual
+ self.residual = np.sqrt(recvbuff[0])
+
+ # Print results, if required
+ if self._writer is not None and self._writer.do_write(ite) :
+ self._writer.buffer[0, 0] = time
+ self._writer.buffer[0, 1] = ite
+ self._writer.buffer[0, 2] = self.residual
+ self._writer.write()
+
+ # update vort(n-1) for next iteration
+ for i in xrange(nbc):
+ self._vortPrev[i][w_ind] = self.vorticity.data[i][w_ind]
+ self._old_dt = dt
+
+
+
diff --git a/HySoP/hysop/operator/forcing.py b/HySoP/hysop/operator/forcing.py
new file mode 100644
index 000000000..bfe03b4e5
--- /dev/null
+++ b/HySoP/hysop/operator/forcing.py
@@ -0,0 +1,155 @@
+# -*- coding: utf-8 -*-
+"""Operator implementing the forcing term in the NS,
+ depending on the filtered field
+ (--> computation of base flow).
+
+.. currentmodule:: hysop.operator.forcing
+
+"""
+
+from hysop.operator.computational import Computational
+from hysop.operator.discrete.forcing import Forcing as DForcing
+from hysop.operator.discrete.forcing import ForcingConserv as DForcingCsv
+from hysop.constants import debug
+from hysop.operator.continuous import opsetup
+import hysop.default_methods as default
+from hysop.methods_keys import SpaceDiscretisation
+from hysop.numerics.finite_differences import FD_C_4,\
+ FD_C_2
+from hysop.operator.differential import Curl
+from hysop.fields.continuous import Field
+
+
+class Forcing(Computational):
+ """
+ Integrate a forcing term in the right hand side
+ of the NS equations, depending on the filtered field.
+ i.e. solve
+ \f{eqnarray*}
+ \frac{\partial \omega}{\partial t} &=& -\chi(\omega - \bar{\omega}))
+ \f}
+ The strength of the forcing is chosen in the order
+ of the amplification rate related to the unstable flow.
+
+ """
+
+ @debug
+ def __init__(self, strength=None, **kwds):
+ """
+ Parameters
+ ----------
+ @param strength : strength of the filter
+ **kwds : extra parameters for parent class
+
+ """
+ super(Forcing, self).__init__(**kwds)
+
+ ## strength of the filter
+ self.strength = strength
+
+ self.input = self.output = self.variables
+# self.output = [self.variables[0]]
+
+ def discretize(self):
+ super(Forcing, self)._standard_discretize()
+ # all variables must have the same resolution
+ assert self._single_topo, 'multi-resolution case not allowed.'
+
+ @debug
+ def setup(self, rwork=None, iwork=None):
+ self.discrete_op = DForcing(
+ variables=self.discreteFields.values(),
+ strength=self.strength,
+ rwork=rwork, iwork=iwork)
+
+ self._is_uptodate = True
+
+
+class ForcingConserv(Forcing):
+ """
+ Integrate a forcing term in the right hand side
+ of the NS equations, depending on the filtered vorticity.
+ i.e. solve
+ \f{eqnarray*}
+ \frac{\partial \omega}{\partial t} &=& -\chi(\omega - \bar{\omega}))
+ \f}
+ The equation is solved using a CONSERVATIVE formulation.
+ The strength of the forcing is chosen in the order
+ of the amplification rate related to the unstable flow.
+
+ """
+ @debug
+ def __init__(self, velocity, vorticity, velocityFilt, **kwds):
+ """
+ Parameters
+ ----------
+ @param[in] velocity, vorticity, velocityFilt fields
+ @param[in, out] forced vorticity field
+ @param strength : strength of the filter
+ **kwds : extra parameters for parent class
+
+ Notes
+ -----
+ velocity and velocityFilt are not modified by this operator.
+ vorticity is an in-out parameter.
+ """
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(ForcingConserv, self).__init__(
+ variables=[velocity, vorticity, velocityFilt], **kwds)
+ ## velocity variable
+ self.velocity = velocity
+ ## vorticity variable
+ self.vorticity = vorticity
+ ## filtered velocity variable
+ self.velocityFilt = velocityFilt
+ # A method is required to set how the curl will be computed.
+ if self.method is None:
+ self.method = default.DIFFERENTIAL
+ # operator to compute buffer = curl(penalised velocity)
+ self._curl = None
+
+ self.input = [self.velocity, self.vorticity, self.velocityFilt]
+ self.output = [self.vorticity]
+
+ def discretize(self):
+
+ if self.method[SpaceDiscretisation] is FD_C_4:
+ # Finite differences method
+ # Minimal number of ghost points
+ nb_ghosts = 2
+ elif self.method[SpaceDiscretisation] is FD_C_2:
+ nb_ghosts = 1
+ else:
+ raise ValueError("Unknown method for space discretization of the\
+ differential operator in penalization.")
+ super(ForcingConserv, self)._standard_discretize(nb_ghosts)
+ # all variables must have the same resolution
+ assert self._single_topo, 'multi-resolution case not allowed.'
+ topo = self.variables[self.velocity]
+ invar = Field(domain=self.velocity.domain,
+ name='curl_in', is_vector=True)
+ dimension = self.domain.dimension
+ outvar = Field(domain=self.velocity.domain,
+ name='curl_out',
+ is_vector=dimension == 3)
+ self._curl = Curl(invar=invar, outvar=outvar,
+ discretization=topo, method=self.method)
+ self._curl.discretize()
+
+ def get_work_properties(self):
+ return self._curl.get_work_properties()
+
+ @debug
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ self._curl.setup(rwork, iwork)
+ self.discrete_op = DForcingCsv(
+ vorticity=self.discreteFields[self.vorticity],
+ velocity=self.discreteFields[self.velocity],
+ velocityFilt=self.discreteFields[self.velocityFilt],
+ curl=self._curl.discrete_op,
+ strength=self.strength,
+ rwork=rwork, iwork=iwork)
+ self._is_uptodate = True
+
+
diff --git a/HySoP/hysop/operator/low_pass_filt.py b/HySoP/hysop/operator/low_pass_filt.py
new file mode 100644
index 000000000..91f03f28b
--- /dev/null
+++ b/HySoP/hysop/operator/low_pass_filt.py
@@ -0,0 +1,153 @@
+# -*- coding: utf-8 -*-
+"""Operator for any field low-pass filtering
+ (--> computation of base flow).
+
+.. currentmodule:: hysop.operator.low_pass_filt
+
+"""
+
+from hysop.operator.computational import Computational
+from hysop.operator.discrete.low_pass_filt import LowPassFilt as DFilt
+from hysop.operator.discrete.low_pass_filt import LowPassFiltConserv as DFiltCsv
+from hysop.constants import debug
+from hysop.operator.continuous import opsetup
+import hysop.default_methods as default
+from hysop.methods_keys import SpaceDiscretisation
+from hysop.numerics.finite_differences import FD_C_4,\
+ FD_C_2
+from hysop.operator.differential import Curl
+from hysop.fields.continuous import Field
+
+
+class LowPassFilt(Computational):
+ """
+ Provides a filtered field by low-pass filtering
+ the flow with the frequency of the instability divided by 2.
+ i.e. solve
+ \f{eqnarray*}
+ \frac{\partial \bar{\omega}}{\partial t} &=& \Omega_c(\omega - \bar{\omega}))
+ \f}
+
+ """
+
+ @debug
+ def __init__(self, cutFreq=None, **kwds):
+ """
+ Parameters
+ ----------
+ @param cutFreq : cutting circular frequency corresponding to the half of
+ the eigenfrequency of the flow instability
+ **kwds : extra parameters for parent class
+
+ """
+ super(LowPassFilt, self).__init__(**kwds)
+
+ ## cutting circular frequency
+ self.cutFreq = cutFreq
+
+ self.input = self.output = self.variables
+
+ def discretize(self):
+ super(LowPassFilt, self)._standard_discretize()
+ # all variables must have the same resolution
+ assert self._single_topo, 'multi-resolution case not allowed.'
+
+ @debug
+ def setup(self, rwork=None, iwork=None):
+ self.discrete_op = DFilt(
+ variables=self.discreteFields.values(),
+ cutFreq=self.cutFreq,
+ rwork=rwork, iwork=iwork)
+
+ self._is_uptodate = True
+
+
+class LowPassFiltConserv(LowPassFilt):
+ """
+ Provides a filtered field by low-pass filtering
+ the flow with the frequency of the instability divided by 2.
+ i.e. solve
+ \f{eqnarray*}
+ \frac{\partial \bar{\omega}}{\partial t} &=& \Omega_c(\omega - \bar{\omega}))
+ \f}
+ The equation is solved using a CONSERVATIVE formulation.
+
+ """
+
+ @debug
+ def __init__(self, velocity, vorticityFilt, velocityFilt, **kwds):
+ """
+ Parameters
+ ----------
+ @param[in] velocity, vorticityFilt, velocityFilt fields
+ @param[in, out] vorticityFilt field (i.e. filtered vorticity field)
+ @param strength : strength of the filter
+ **kwds : extra parameters for parent class
+
+ Notes
+ -----
+ velocity and velocityFilt are not modified by this operator.
+ vorticityFilt is an in-out parameter.
+ """
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(LowPassFiltConserv, self).__init__(
+ variables=[velocity, vorticityFilt, velocityFilt], **kwds)
+
+ ## velocity variable
+ self.velocity = velocity
+ ## filtered vorticity variable
+ self.vorticityFilt = vorticityFilt
+ ## filtered velocity variable
+ self.velocityFilt = velocityFilt
+ # A method is required to set how the curl will be computed.
+ if self.method is None:
+ self.method = default.DIFFERENTIAL
+ # operator to compute buffer = curl(penalised velocity)
+ self._curl = None
+
+ self.input = [self.velocity, self.vorticityFilt, self.velocityFilt]
+ self.output = [self.vorticityFilt]
+
+ def discretize(self):
+
+ if self.method[SpaceDiscretisation] is FD_C_4:
+ # Finite differences method
+ # Minimal number of ghost points
+ nb_ghosts = 2
+ elif self.method[SpaceDiscretisation] is FD_C_2:
+ nb_ghosts = 1
+ else:
+ raise ValueError("Unknown method for space discretization of the\
+ differential operator in penalization.")
+ super(LowPassFiltConserv, self)._standard_discretize(nb_ghosts)
+ # all variables must have the same resolution
+ assert self._single_topo, 'multi-resolution case not allowed.'
+ topo = self.variables[self.velocity]
+ invar = Field(domain=self.velocity.domain,
+ name='curl_in', is_vector=True)
+ dimension = self.domain.dimension
+ outvar = Field(domain=self.velocity.domain,
+ name='curl_out',
+ is_vector=dimension == 3)
+ self._curl = Curl(invar=invar, outvar=outvar,
+ discretization=topo, method=self.method)
+ self._curl.discretize()
+
+ def get_work_properties(self):
+ return self._curl.get_work_properties()
+
+ @debug
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ self._curl.setup(rwork, iwork)
+ self.discrete_op = DFiltCsv(
+ vorticityFilt=self.discreteFields[self.vorticityFilt],
+ velocity=self.discreteFields[self.velocity],
+ velocityFilt=self.discreteFields[self.velocityFilt],
+ curl=self._curl.discrete_op,
+ cutFreq=self.cutFreq,
+ rwork=rwork, iwork=iwork)
+ self._is_uptodate = True
+
+
+
diff --git a/HySoP/hysop/operator/monitoringPoints.py b/HySoP/hysop/operator/monitoringPoints.py
new file mode 100644
index 000000000..3066fda28
--- /dev/null
+++ b/HySoP/hysop/operator/monitoringPoints.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+@file monitoringPoints.py
+Print time evolution of flow variables (velo, vorti)
+at a particular monitoring point in the wake
+"""
+from hysop.operator.discrete.monitoringPoints import MonitoringPoints as MonitD
+from hysop.operator.computational import Computational
+from hysop.operator.continuous import opsetup
+
+
+class MonitoringPoints(Computational):
+ """
+ Compute and print velo/vorti profiles
+ """
+
+ def __init__(self, velocity, vorticity, monitPt_coords, **kwds):
+ """
+ Constructor.
+ @param velocity field
+ @param vorticity field
+ @param monitPts_coords : coordinates corresponding
+ to the space location of the monitoring point in the wake
+
+ Default file name = 'monit.dat'
+ See hysop.tools.io_utils.Writer for details
+ """
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(MonitoringPoints, self).__init__(variables=
+ [velocity, vorticity],
+ **kwds)
+ ## velocity field
+ self.velocity = velocity
+ ## vorticity field
+ self.vorticity = vorticity
+ ## coordinates of the monitoring point
+ self.monitPt_coords = monitPt_coords
+ self.input = [velocity, vorticity]
+ self.output = []
+
+ def get_work_properties(self):
+ if not self._is_discretized:
+ msg = 'The operator must be discretized '
+ msg += 'before any call to this function.'
+ raise RuntimeError(msg)
+ vd = self.discreteFields[self.velocity]
+ wd = self.discreteFields[self.vorticity]
+ v_ind = vd.topology.mesh.iCompute
+ w_ind = wd.topology.mesh.iCompute
+ shape_v = vd[0][v_ind].shape
+ shape_w = wd[0][w_ind].shape
+ if shape_v == shape_w:
+ return {'rwork': [shape_v], 'iwork': None}
+ else:
+ return {'rwork': [shape_v, shape_w], 'iwork': None}
+
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ if not self._is_uptodate:
+
+ self.discrete_op = MonitD(self.discreteFields[self.velocity],
+ self.discreteFields[self.vorticity],
+ self.monitPt_coords, rwork=rwork)
+ # Output setup
+ self._set_io('monit', (1, 3))
+ self.discrete_op.setWriter(self._writer)
+ self._is_uptodate = True
+
diff --git a/HySoP/hysop/operator/residual.py b/HySoP/hysop/operator/residual.py
new file mode 100644
index 000000000..1a5de67ae
--- /dev/null
+++ b/HySoP/hysop/operator/residual.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+"""
+@file residual.py
+Compute and print the time evolution of the residual
+"""
+from hysop.operator.discrete.residual import Residual as ResD
+from hysop.operator.computational import Computational
+from hysop.operator.continuous import opsetup
+
+
+class Residual(Computational):
+ """
+ Compute and print the residual time evolution
+ """
+
+ def __init__(self, vorticity, **kwds):
+ """
+ Constructor.
+ @param vorticity field
+
+ Default file name = 'residual.dat'
+ See hysop.tools.io_utils.Writer for details
+ """
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(Residual, self).__init__(variables=[vorticity], **kwds)
+ ## vorticity field
+ self.vorticity = vorticity
+ self.input = [vorticity]
+ self.output = []
+
+ def get_work_properties(self):
+ if not self._is_discretized:
+ msg = 'The operator must be discretized '
+ msg += 'before any call to this function.'
+ raise RuntimeError(msg)
+ wd = self.discreteFields[self.vorticity]
+ w_ind = wd.topology.mesh.iCompute
+ shape_w = wd[0][w_ind].shape
+ return {'rwork': [shape_w], 'iwork': None}
+
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ if not self._is_uptodate:
+
+ self.discrete_op = ResD(self.discreteFields[self.vorticity],
+ rwork=rwork)
+ # Initialization of w^(n-1) vorticity value
+ self.discrete_op.initialize_vortPrev()
+ # Output setup
+ self._set_io('residual', (1, 3))
+ self.discrete_op.setWriter(self._writer)
+ self._is_uptodate = True
+
--
GitLab