From 3ae9d1f397b098fb14f41d501867a82c76422a07 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Chlo=C3=A9=20Mimeau?= <Chloe.Mimeau@imag.fr>
Date: Mon, 29 Feb 2016 11:42:59 +0100
Subject: [PATCH] general updates
---
Examples/FlowAroundSphere_linearized.py | 221 +++++++++---------
HySoP/hysop/operator/discrete/forcing.py | 6 +-
.../operator/discrete/monitoringPoints.py | 16 +-
HySoP/hysop/operator/discrete/stretching.py | 43 +++-
HySoP/hysop/operator/stretching.py | 8 +-
HySoP/hysop/operator/tests/test_Stretching.py | 62 ++++-
HySoP/hysop/tools/numpywrappers.py | 7 +
7 files changed, 218 insertions(+), 145 deletions(-)
diff --git a/Examples/FlowAroundSphere_linearized.py b/Examples/FlowAroundSphere_linearized.py
index d41cf782e..d2bc313bd 100644
--- a/Examples/FlowAroundSphere_linearized.py
+++ b/Examples/FlowAroundSphere_linearized.py
@@ -16,7 +16,9 @@ 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.stretching import Stretching, \
+ StretchingLinearized
+from hysop.operator.dissip_filter import DissipFilter
from hysop.operator.absorption_BC import AbsorptionBC
from hysop.operator.poisson import Poisson
from hysop.operator.diffusion import Diffusion
@@ -33,7 +35,7 @@ from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\
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.finite_differences import FD_C_4, FD_C_2, Filter_C_4
from hysop.numerics.interpolation import Linear
from hysop.numerics.remeshing import L6_4 as rmsh
import hysop.tools.io_utils as io
@@ -53,10 +55,10 @@ VISCOSITY = 1. / 300.
# ======= Domain =======
dim = 3
-Nx = 129
-Ny = Nz = 65
-#Nx = 257
-#Ny = Nz = 129
+#Nx = 129
+#Ny = Nz = 65
+Nx = 257
+Ny = Nz = 129
#Nx = 513
#Ny = Nz = 257
#Nx = 1025
@@ -78,7 +80,7 @@ from hysop.domain.subsets import Sphere, HemiSphere
sphere = Sphere(origin=pos, radius=RADIUS, parent=box)
-# ======= Function to set initial velocity BF =======
+# ======= Function to set initial velocity =======
def setZeroVel(res, x, y, z, t):
res[0][...] = 0.
res[1][...] = 0.
@@ -86,7 +88,7 @@ def setZeroVel(res, x, y, z, t):
return res
-# ======= Function to set initial vorticity BF =======
+# ======= Function to set initial vorticity =======
def setZeroVort(res, x, y, z, t):
res[0][...] = 0.
res[1][...] = 0.
@@ -122,7 +124,7 @@ 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)
+simu = Simulation(tinit=0.0, tend=2000.0, timeStep=0.005, iterMax=10000000)
# Adaptative timestep method : dt = min(values(dtCrit))
@@ -150,36 +152,44 @@ op['dtAdapt'] = AdaptTimeStep(velo, vorti, simulation=simu,
lcfl=0.125,
cfl=0.5)
-op['advection1'] = Advection(velo, vortiBF,
- discretization=d3d,
- method={Scales: 'p_M6',
- Splitting: 'classic'}
- )
+op['advectionBridge'] = Advection(veloBF, vortiBF,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
-op['advection2'] = Advection(veloBF, vorti,
- discretization=d3d,
- method={Scales: 'p_M6',
- Splitting: 'classic'}
- )
+op['advection'] = Advection(veloBF, vorti,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
-op['stretching1'] = Stretching(veloBF, vorti,
- discretization=topo_with_ghosts)
+op['stretchingLin'] = StretchingLinearized(velocity=velo,
+ vorticity=vorti,
+ velocity_BF=veloBF,
+ vorticity_BF=vortiBF,
+ discretization=topo_with_ghosts)
-op['stretching2'] = Stretching(velo, vortiBF,
- discretization=topo_with_ghosts)
+
+op['artifDissip'] = DissipFilter(velo, vorti,
+ discretization=topo_with_ghosts,
+ method={SpaceDiscretisation: Filter_C_4})
op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
discretization=d3d)
rate = VariableParameter(formula=computeFlowrate)
-op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)
+op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)#, projection=1)
+
+op['poissonProj'] = Poisson(velo, vorti, discretization=d3d,
+ flowrate=rate, projection=1)
# ===== Discretization of computational operators ======
for ope in op.values():
ope.discretize()
topofft = op['poisson'].discreteFields[vorti].topology
-topoadvec = op['advection2'].discreteFields[vorti].topology
+topoadvec = op['advection'].discreteFields[vorti].topology
# ===== Smooth vorticity absorption at the outlet =====
op['vort_absorption'] = AbsorptionBC(velo, vorti, discretization=topofft,
@@ -199,42 +209,17 @@ 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'],
+distr['advec2str'] = RedistributeIntra(source=op['advectionBridge'],
+ target=op['stretchingLin'],
+ variables=[veloBF, vortiBF])
+
+distr['fft2str'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretchingLin'],
+ variables=[velo, vorti])
+
+distr['str2fft'] = RedistributeIntra(source=op['stretchingLin'],
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])
+ variables=[velo, vorti])
# ========= Monitoring operators =========
monitors = {}
@@ -254,12 +239,14 @@ 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]
+coordsMonit = [1.0, 0.0, 0.0]
rk = 0
-if (coordsMonit[2] in topofft.mesh.coords[2]):
+if (coordsMonit[2] in topo_with_ghosts.mesh.coords[2]):
rk = main_rank
+print 'rk ... ', rk
io_monit = IOParams('monit', frequency=1, io_leader=rk)
-monitors['monit'] = MonitoringPoints(velo, vorti, discretization=topofft,
+monitors['monit'] = MonitoringPoints(velo, vorti,
+ discretization=topo_with_ghosts,
io_params=io_monit,
monitPt_coords=coordsMonit)
@@ -303,7 +290,7 @@ ref_step = topo_with_ghosts.mesh.space_step
# io_params=io_forcesPenal)
step_dir = ref_step[0]
-io_sliceXY = IOParams('sliceXY', frequency=2000)
+io_sliceXY = IOParams('sliceXY', frequency=1)
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],
@@ -312,7 +299,7 @@ monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft}
io_params=io_sliceXY, subset=thickSliceXY,
xmfalways=True)
-io_sliceXZ = IOParams('sliceXZ', frequency=2000)
+io_sliceXZ = IOParams('sliceXZ', frequency=1)
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},
@@ -365,40 +352,38 @@ print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init
fullseq = []
def run(sequence):
+ print 'norm vort perturb 1 :', vorti.norm(topofft)
op['vort_absorption'].apply(simu)
+ print 'norm vort perturb 2 (apres absorption ):', vorti.norm(topofft)
op['poisson'].apply(simu) # Poisson + correction
+ print 'norm vort perturb 3 (apres Poisson ):', vorti.norm(topofft)
# monitors['forcesMom'].apply(simu) # Forces Heloise
- distr['fft2penal'].apply(simu)
- distr['fft2penal'].wait()
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].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()
+ print 'norm vort perturb 4 (apres penalisation):', vorti.norm(topo_with_ghosts)
+ op['stretchingLin'].apply(simu) # Stretching linearized
+ print 'norm vort perturb 5 (apres stretching ):', vorti.norm(topo_with_ghosts)
+ op['artifDissip'].apply(simu) # Dissipation filter
+ print 'norm vort perturb 6 (apres dissp filter):', vorti.norm(topo_with_ghosts)
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].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)
+ print 'norm vort perturb 7 (apres diffusion ):', vorti.norm(topofft)
+ op['advection'].apply(simu) # Advection
+ print 'norm vort perturb 8 (apres advection ):', vorti.norm(topofft)
+
+# 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()
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ monitors['monit'].apply(simu) # Monitoring points in the wake
+# op['dtAdapt'].apply(simu) # Update timestep
+# op['dtAdapt'].wait()
# ==== Serialize the simulation data of the problem to a "restart" file ====
def dump(filename):
@@ -429,10 +414,12 @@ def restart(filename):
filedump = filename + '_rk_' + str(main_rank)
db = open(filedump, 'r')
simu = cPickle.load(db)
- simu.start = simu.time - simu.timeStep
- ite = simu.currentIteration
+# simu.start = simu.time - simu.timeStep
+# ite = simu.currentIteration
+ simu.start = 0.0
+ simu.timeStep = 0.005
simu.initialize()
- simu.currentIteration = ite
+# simu.currentIteration = ite
print 'simu', simu
print ("load ...", filename)
return simu
@@ -447,31 +434,49 @@ io_default=IOParams('restart')
dump_filename = io.Writer(io_params=io_default).filename
#===== Restart (if needed) =====
if doRestart:
- # Load data from dumped files
+ # Fields initialization
+ veloBF.initialize(topo=topofft)
+ vortiBF.initialize(topo=topofft)
+ veloBF.initialize(topo=topo_with_ghosts)
+ vortiBF.initialize(topo=topo_with_ghosts)
+ velo.initialize(topo=topofft)
+ vorti.initialize(topo=topofft)
+ # Load data from dumped files --> base flow
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()
+ # vortBF projection + Poisson
+ op['poissonProj'].apply(simu)
+ # Initialize velo and vorti to a small perturbation
+ # equal to baseFlow * 1e-8
+ # (cf Florian Guiho's thesis p.30 eq. (2.13))
+ for d in xrange(dim):
+ velo.discreteFields[topofft].data[d][...] = \
+ veloBF.discreteFields[topofft].data[d][...].copy()
+ velo.discreteFields[topofft].data[d][...] *= 1e-8
+ vorti.discreteFields[topofft].data[d][...] = \
+ vortiBF.discreteFields[topofft].data[d][...].copy()
+ vorti.discreteFields[topofft].data[d][...] *= 1e-8
# Redistribute veloBF and vortBF on topoGhosts
- distr['advec22str1'].apply(simu)
- distr['advec22str1'].wait()
- distr['advec12str2'].apply(simu)
- distr['advec12str2'].wait()
+ distr['advec2str'].apply(simu)
+ distr['advec2str'].wait()
+
+ print '======= INITIALIZATION ======='
+ print 'vortiBF', vortiBF.norm(topofft)
+# print 'vortiBF norm Ghosts', vortiBF.norm(topo_with_ghosts)
+ print 'vorti', vorti.norm(topofft)
+
+ print 'veloBF', veloBF.norm(topofft)
+# print 'veloBF norm Ghosts', veloBF.norm(topo_with_ghosts)
+ print 'velo', velo.norm(topofft)
+
# ======= Time loop =======
time_run = MPI.Wtime()
diff --git a/HySoP/hysop/operator/discrete/forcing.py b/HySoP/hysop/operator/discrete/forcing.py
index 55a00fa9e..7d9121db5 100644
--- a/HySoP/hysop/operator/discrete/forcing.py
+++ b/HySoP/hysop/operator/discrete/forcing.py
@@ -64,7 +64,7 @@ class ForcingConserv(Forcing):
\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}
+ \omega^n + \nabla \times (\frac{\Delta t \chi}{1+\Delta t \chi} (\bar{u^{n+1}}-u^n))
\f}
"""
@@ -87,8 +87,8 @@ class ForcingConserv(Forcing):
assert 'variables' not in kwds, 'variables parameter is useless.'
super(ForcingConserv, self).__init__(variables=[vorticity,
- velocity,
- velocityFilt],
+ velocity,
+ velocityFilt],
**kwds)
# Input vector fields
self.velocity = velocity
diff --git a/HySoP/hysop/operator/discrete/monitoringPoints.py b/HySoP/hysop/operator/discrete/monitoringPoints.py
index eac5f9154..74738987b 100644
--- a/HySoP/hysop/operator/discrete/monitoringPoints.py
+++ b/HySoP/hysop/operator/discrete/monitoringPoints.py
@@ -114,14 +114,14 @@ class MonitoringPoints(DiscreteOperator):
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)
+ # Compute velocity and vorticity norm
+ # at the monitoring point
+ 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) :
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 52688d693..675961bae 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -127,7 +127,6 @@ class Stretching(DiscreteOperator):
self._synchronize(self.velocity.data + self.vorticity.data)
# Compute stretching term (rhs) and update vorticity
self._compute(dt, t)
- print 'passage dans Stretch 1'
def apply(self, simulation=None):
"""
@@ -251,8 +250,8 @@ class GradUW(Stretching):
return nb_cycles, subdt
-class StretchingLinearized(Conservative):
- """Conservative formulation of (\om_BF \cdot \nabla)\bu
+class StretchingLinearized(Stretching):
+ """By default: Conservative formulation
"""
def __init__(self, vorticity_BF, usual_op, **kwds):
@@ -265,24 +264,48 @@ class StretchingLinearized(Conservative):
self.vorticity_BF.nb_components)
self.usual_op = usual_op
- super(StretchingLinearized, self).__init__(**kwds)
+ # Right-hand side for time integration
+ def rhs(t, y, result, form):
+ if form == "div(w:u)":
+ result = self.strFunc(y, self.velocity.data, result)
+ else:
+ result = self.strFunc(y, self.vorticity_BF.data, result)
+ return result
+
+ super(StretchingLinearized, self).__init__(formulation=diff_op.DivWV,
+ rhs=rhs, **kwds)
def _integrate(self, dt, t):
- # - Call time integrator -
- # Init workspace with a first evaluation of the
+ # - Call time integrator (1st term over 3) -
+ # Init workspace with a first evaluation of the div(wb:u') term in the
# rhs of the integrator
self._ti_work[:self.nb_components] = \
self.timeIntegrator.f(t, self.vorticity_BF.data,
- self._ti_work[:self.nb_components])
+ self._ti_work[:self.nb_components], "div(w:u)")
+ # perform integration and save result in-place
+ self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
+ result=self.vorticity.data)
+ # - Call time integrator (2nd term over 3) -
+ # Init workspace with a first evaluation of the div(u':wb) term in the
+ # rhs of the integrator
+ self._ti_work[:self.nb_components] = \
+ self.timeIntegrator.f(t, self.velocity.data,
+ self._ti_work[:self.nb_components], "div(u:w)")
# perform integration and save result in-place
- self.vorticity.data = self.timeIntegrator(t, self.vorticity_BF.data, dt,
+ self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
result=self.vorticity.data)
+ def _compute(self, dt, t):
+ # No subcycling for this formulation
+ self._integrate(dt, t)
+
def _apply(self, t, dt):
# Synchronize ghost points
self._synchronize(self.velocity.data + self.vorticity.data)
self._synchronize_vort_BF(self.vorticity_BF.data)
- # Compute stretching term (rhs) and update vorticity
+ # Compute the 2 first "stretching" terms (div(wb:u') and div(u':wb))
+ # and update vorticity for each of them
self._compute(dt, t)
+ # Compute the 3rd stretching term (div(w':ub)) and update vorticity
self.usual_op._apply(t, dt)
- print 'passage dans Stretch 2'
+
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 802533422..5320c192e 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -120,7 +120,6 @@ class Stretching(Computational):
stretching operator.")
super(Stretching, self)._standard_discretize(nbghosts)
- print 'passage dans discretize'
@debug
@opsetup
@@ -130,7 +129,6 @@ class Stretching(Computational):
vorticity=self.discreteFields[self.vorticity],
method=self.method, rwork=rwork, iwork=iwork)
self._is_uptodate = True
- print 'passage dans setup classique'
class StretchingLinearized(Stretching):
@@ -168,7 +166,7 @@ class StretchingLinearized(Stretching):
self.vorticity_BF = vorticity_BF
# Usual stretching operator to compute the
- # first term of the linearized rhs: (w.grad)ub
+ # first term of the linearized rhs: (w'.grad)ub
self.usual_stretch = Stretching(self.velocity_BF, self.vorticity,
discretization=self._discretization)
@@ -184,7 +182,6 @@ class StretchingLinearized(Stretching):
stretching operator.")
self.usual_stretch.discretize()
super(StretchingLinearized, self)._standard_discretize(nbghosts)
- print 'passage dans discretize 2'
@debug
@@ -194,7 +191,7 @@ class StretchingLinearized(Stretching):
self.usual_stretch.setup()
# Setup of a second stretching operator (receiving 3 input variables)
- # to compute the 2nd term of the linearized rhs: (wb.grad)u
+ # to compute the 2nd term of the linearized rhs: (wb.grad)u'
method_lin = self.method.copy()
method_lin[TimeIntegrator] = Euler
self.discrete_op =\
@@ -204,6 +201,5 @@ class StretchingLinearized(Stretching):
usual_op=self.usual_stretch.discrete_op,
method=method_lin, rwork=rwork, iwork=iwork)
self._is_uptodate = True
- print 'passage dans setup bis'
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
index 7a9803a9d..deef2a022 100755
--- a/HySoP/hysop/operator/tests/test_Stretching.py
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -5,9 +5,9 @@ from hysop.fields.continuous import Field
from hysop.operator.stretching import Stretching, \
StretchingLinearized
from hysop.problem.simulation import Simulation
-#from hysop.methods_keys import TimeIntegrator, Formulation,\
-# SpaceDiscretisation
-#from hysop.methods import RK3, FD_C_4, Conservative
+from hysop.methods_keys import TimeIntegrator, Formulation,\
+ SpaceDiscretisation
+from hysop.methods import Euler, RK3, FD_C_4, Conservative
from hysop.tools.parameters import Discretization
import hysop.tools.numpywrappers as npw
pi = np.pi
@@ -54,6 +54,45 @@ def computeVort(res, x, y, z, t):
return res
+def computeVelBF(res, x, y, z, t):
+ amodul = cos(pi * 1. / 3.)
+ pix = pi * x
+ piy = pi * y
+ piz = pi * z
+ pi2x = 2. * pix
+ pi2y = 2. * piy
+ pi2z = 2. * piz
+ res[0][...] = 2. * sin(pix) * sin(pix) \
+ * sin(pi2y) * sin(pi2z) * amodul
+ res[1][...] = - sin(pi2x) * sin(piy) \
+ * sin(piy) * sin(pi2z) * amodul
+ res[2][...] = - sin(pi2x) * sin(piz) \
+ * sin(piz) * sin(pi2y) * amodul
+ return res
+
+
+def computeVortBF(res, x, y, z, t):
+ amodul = cos(pi * 1. / 3.)
+ pix = pi * x
+ piy = pi * y
+ piz = pi * z
+ pi2x = 2. * pix
+ pi2y = 2. * piy
+ pi2z = 2. * piz
+ res[0][...] = 2. * pi * sin(pi2x) * amodul *\
+ (- cos(pi2y) * sin(piz) * sin(piz)
+ + sin(piy) * sin(piy) * cos(pi2z))
+
+ res[1][...] = 2. * pi * sin(pi2y) * amodul *\
+ (2. * cos(pi2z) * sin(pix) * sin(pix)
+ + sin(piz) * sin(piz) * cos(pi2x))
+
+ res[2][...] = -2. * pi * sin(pi2z) * amodul *\
+ (cos(pi2x) * sin(piy) * sin(piy)
+ + sin(pix) * sin(pix) * cos(pi2y))
+
+ return res
+
def test_stretching():
@@ -107,33 +146,36 @@ def test_stretchingLinearized():
domain=box, formula=computeVort,
name='Vorticity', is_vector=True)
veloBF = Field(
- domain=box, formula=computeVel,
+ domain=box, formula=computeVelBF,
name='VelocityBF', is_vector=True)
vortiBF = Field(
- domain=box, formula=computeVort,
+ domain=box, formula=computeVortBF,
name='VorticityBF', is_vector=True)
# Operators
# Usual stretching operator
-# stretch1 = Stretching(velo, vorti, discretization=nbElem)
+ stretch1 = Stretching(velo, vorti, discretization=nbElem)
# Linearized stretching operator
stretch2 = StretchingLinearized(velocity=velo,
vorticity=vorti,
velocity_BF=veloBF,
vorticity_BF=vortiBF,
discretization=nbElem)
-# stretch1.discretize()
+ stretch1.discretize()
stretch2.discretize()
- topo = stretch2.discreteFields[velo].topology
+ topo = stretch1.discreteFields[velo].topology
velo.initialize(topo=topo)
vorti.initialize(topo=topo)
veloBF.initialize(topo=topo)
vortiBF.initialize(topo=topo)
-# stretch1.setup()
+ stretch1.setup()
stretch2.setup()
simulation = Simulation(tinit=0, tend=1., timeStep=timeStep)
-# stretch1.apply(simulation)
+ stretch1.apply(simulation)
+ print 'norm vorti (usual):', vorti.norm(topo)
+
stretch2.apply(simulation)
+ print 'norm vorti (lin):', vorti.norm(topo)
def test_stretching_external_work():
diff --git a/HySoP/hysop/tools/numpywrappers.py b/HySoP/hysop/tools/numpywrappers.py
index f126de7fd..b0868f27e 100644
--- a/HySoP/hysop/tools/numpywrappers.py
+++ b/HySoP/hysop/tools/numpywrappers.py
@@ -38,6 +38,13 @@ def ones_like(tab):
return np.ones_like(tab, dtype=tab.dtype, order=ORDER)
+def reshape(tab, shape):
+ """
+ Wrapper to numpy.reshape, force order to hysop.constants.ORDER
+ """
+ return np.reshape(tab, shape, dtype=tab.dtype, order=ORDER)
+
+
def realempty(tab):
"""
Wrapper to numpy.empty, force order to hysop.constants.ORDER
--
GitLab