Skip to content
Snippets Groups Projects
Commit 3ae9d1f3 authored by Chloé Mimeau's avatar Chloé Mimeau
Browse files

general updates

parent b7679ed6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 218 additions and 145 deletions
Examples/FlowAroundSphere_linearized.py
+ 113
− 108
View file @ 3ae9d1f3
......@@ -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()
......
HySoP/hysop/operator/discrete/forcing.py
+ 3
− 3
View file @ 3ae9d1f3
......@@ -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
......
HySoP/hysop/operator/discrete/monitoringPoints.py
+ 8
− 8
View file @ 3ae9d1f3
......@@ -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) :
......
HySoP/hysop/operator/discrete/stretching.py
+ 33
− 10
View file @ 3ae9d1f3
......@@ -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'
HySoP/hysop/operator/stretching.py
+ 2
− 6
View file @ 3ae9d1f3
......@@ -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'
HySoP/hysop/operator/tests/test_Stretching.py
+ 52
− 10
View file @ 3ae9d1f3
......@@ -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():
......
HySoP/hysop/tools/numpywrappers.py
+ 7
− 0
View file @ 3ae9d1f3
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment