diff --git a/Examples/FlowAroundSphere.py b/Examples/FlowAroundSphere.py
index eee8df58db35c431026ac2951c2ff9043e7a9325..e769394f1f3566c48883faac1f90db10ef77b70c 100644
--- a/Examples/FlowAroundSphere.py
+++ b/Examples/FlowAroundSphere.py
@@ -15,6 +15,7 @@ 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
@@ -48,8 +49,8 @@ VISCOSITY = 1. / 300.
# ======= Domain =======
dim = 3
-Nx = 129
-Ny = Nz = 65
+Nx = 257
+Ny = Nz = 129
g = 2
boxlength = npw.asrealarray([10.24, 5.12, 5.12])
boxorigin = npw.asrealarray([-2.0, -2.56, -2.56])
@@ -106,7 +107,7 @@ vorti = Field(domain=box, formula=computeVort,
# ========= Simulation setup =========
-simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=50)
+simu = Simulation(tinit=0.0, tend=80.0, timeStep=0.0125, iterMax=10000000)
# Adaptative timestep method : dt = min(values(dtCrit))
@@ -146,10 +147,12 @@ op['stretching'] = Stretching(velo, vorti,
op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
discretization=d3d)
-
rate = VariableParameter(formula=computeFlowrate)
+op['vort_absorption'] = AbsorptionBC(velo, vorti, discretization=d3d,
+ req_flowrate=rate,
+ x_coords_absorp=[7.24, 8.24])
+
op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)
-#op['poisson'] = Poisson(velo, vorti, discretization=d3d)
# ===== Discretization of computational operators ======
for ope in op.values():
@@ -159,10 +162,8 @@ topofft = op['poisson'].discreteFields[vorti].topology
topoadvec = op['advection'].discreteFields[vorti].topology
# ===== Penalization of the vorticity on a sphere inside the domain =====
-#from hysop.operator.penalize_vorticity import PenalizeVorticity
-from hysop.operator.penalization_and_curl import PenalizationAndCurl
-#op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
-op['penalVort'] = PenalizationAndCurl(velocity=velo, vorticity=vorti,
+from hysop.operator.penalize_vorticity import PenalizeVorticity
+op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
discretization=topo_with_ghosts,
obstacles=[sphere], coeff=1e8,
method={SpaceDiscretisation: FD_C_4})
@@ -206,7 +207,7 @@ cblength -= 30 * ref_step
cb = ControlBox(parent=box, origin=cbpos, length=cblength)
coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
-io_forces=IO_params('drag_and_lift_Noca')
+io_forces=IO_params('drag_and_lift_NocaII')
#monitors['forcesNoca'] = NocaForces(velo, vorti,
# discretization=topo_with_ghosts,
# nu=VISCOSITY,
@@ -220,6 +221,7 @@ monitors['forcesMom'] = MomentumForces(velocity=velo,
discretization=topo_with_ghosts,
normalization=coeffForce,
obstacles=[sphere],
+ coeff=[1e8],
io_params=io_forcesPenal)
#io_forcesPenal=IO_params('drag_and_lift_penal')
@@ -228,13 +230,13 @@ monitors['forcesMom'] = MomentumForces(velocity=velo,
# obstacles=[sphere], factor=[1e8],
# io_params=io_forcesPenal)
-#step_dir = ref_step[0]
-#io_sliceXY = IO_params('sliceXY', frequency=200)
-#thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir],
-# lengths=[10.24, 5.12, 4.0 * step_dir])
-#monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
-# io_params=io_sliceXY, subset=thickSliceXY,
-# xmfalways=True)
+step_dir = ref_step[0]
+io_sliceXY = IO_params('sliceXY', frequency=10)
+thickSliceXY = ControlBox(parent=box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ length=[10.24, 5.12, 4.0 * step_dir])
+monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXY, subset=thickSliceXY,
+ xmfalways=True)
#io_sliceXZ = IO_params('sliceXZ', frequency=2000)
#thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
@@ -249,10 +251,6 @@ monitors['forcesMom'] = MomentumForces(velocity=velo,
# io_params=io_subBox, subset=subBox,
# xmfalways=True)
-## ========= Monitors discretization=========
-for monit in monitors.values():
- monit.discretize()
-
# ========= Setup for all declared operators/monitors =========
time_setup = MPI.Wtime()
for ope in op.values():
@@ -286,11 +284,17 @@ 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['fft2str'].apply(simu)
distr['fft2str'].wait()
op['penalVort'].apply(simu) # Vorticity penalization
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].wait()
+ op['poisson'].apply(simu)
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
op['stretching'].apply(simu) # Stretching
# monitors['forcesNoca'].apply(simu) # Forces Noca
distr['str2fft'].apply(simu)
@@ -301,7 +305,7 @@ def run(sequence):
op['advection'].apply(simu) # Advection (scales)
distr['advec2fft'].apply(simu)
distr['advec2fft'].wait()
-# monitors['writerSliceXY'].apply(simu)
+ monitors['writerSliceXY'].apply(simu)
# monitors['writerSliceXZ'].apply(simu)
# monitors['writerSubBox'].apply(simu)
monitors['energy'].apply(simu) # Energy/enstrophy
@@ -386,5 +390,5 @@ print '[', main_rank, '] total time for run :', MPI.Wtime() - time_run
fftw2py.clean_fftw_solver(box.dimension)
for ope in distr.values():
ope.finalize()
-writer.finalize()
-energy.finalize()
+for monit in monitors.values():
+ monit.finalize()
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index ebe521df43a3008709ee5cbbdbe1e6dfc7ae942e..9fddce6d0b53237905a55501756ec68dc54b607c 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -20,7 +20,7 @@ class Forces(DiscreteOperator):
"""
__metaclass__ = ABCMeta
- def __init__(self, obstacles, **kwds):
+ def __init__(self, obstacles, normalization=1., **kwds):
"""
@param velocity field
@param vorticity field
@@ -58,6 +58,14 @@ class Forces(DiscreteOperator):
# Ghost points synchronizer
self._synchronize = UpdateGhosts(self._topology, nbc)
+ # Normalizing coefficient for forces
+ # (based on the physics of the flow)
+ # \warning : a normalization is applied for
+ # all the forces (Noca's + momentum).
+ # If the normalization parameter is not specified,
+ # then it is set to 1 (ie : no normalization)
+ self._normalization = normalization
+
# Which formula must be used to compute the forces.
# Must be set in derived class.
self._formula = lambda dt: 0
@@ -94,7 +102,10 @@ class Forces(DiscreteOperator):
Performs MPI reduction (sum result value over all process)
All process get the result of the sum.
"""
- self.force = self._subcomm.allreduce(self.force)
+ # \warning : the MPI communicator has been changed
+ # here for the allreduce operation : it is now the same
+ # as the one used for the energy/enstrophy computation
+ self.force = self._topology.comm.allreduce(self.force)
def _mpi_sum(self, root=0):
"""
@@ -102,7 +113,7 @@ class Forces(DiscreteOperator):
Result send only to 'root' process.
@param root : number of the process which get the result.
"""
- self.force = self._subcomm.reduce(self.force, root=root)
+ self.force = self._topology.comm.reduce(self.force, root=root)
def apply(self, simulation=None):
"""
@@ -124,6 +135,8 @@ class Forces(DiscreteOperator):
self._formula(dt)
# Reduce results over MPI processes
self.mpi_sum()
+ # normalization of the forces --> cD, cL, cZ
+ self.force *= self._normalization
# Print results, if required
ite = simulation.currentIteration
if self._writer is not None and self._writer.do_write(ite):
@@ -139,7 +152,7 @@ class MomentumForces(Forces):
The present class implements formula (52) of Plouhmans2002.
Integral inside the obstacle is not taken into account.
"""
- def __init__(self, velocity, **kwds):
+ def __init__(self, velocity, coeff, **kwds):
"""
@param velocity field
@param vorticity field
@@ -154,11 +167,23 @@ class MomentumForces(Forces):
assert 'variables' not in kwds, 'variables parameter is useless.'
## discrete velocity field
self.velocity = velocity
+ ## penalieation coefficient (lambda)
+ self._coeff = coeff
#
self._obstacle = None
super(MomentumForces, self).__init__(variables=[velocity], **kwds)
self._formula = self._momentum
self._subcomm = self._obstacle.subcomm[self._topology]
+ self._dvol = npw.prod(self._topology.mesh.space_step)
+ self._init_buffers()
+
+ def _init_buffers(self):
+ """
+ Allocate memory for local buffers (used to compute
+ volume integrals of the velocity)
+ """
+ shape = self.velocity.data[0][...].shape
+ self._buff = npw.zeros(shape)
def _init_indices(self, obstacles):
msg = 'obstacles arg must be a list.'
@@ -171,16 +196,21 @@ class MomentumForces(Forces):
return self._obstacle.ind[self._topology]
def _momentum(self, dt):
- # -- Integration over the volume of control --
- nbc = self.velocity.nbComponents
- self._buffer[...] = \
- [self._obstacle.integrate_dfield_on_proc(self.velocity,
- component=d)
- for d in xrange(nbc)]
- self.force[...] = 1. / dt * (self._buffer - self._previous)
- # Update previous for next time step ...
- # Which value must we choose for previous??
- self._previous[...] = self._buffer[...]
+ # -- Integration over the obstacle --
+ # \warning : the force has to be set to 0 before computation
+ self.force[...] = 0.0
+ ind = self._indices
+ for d in xrange(self._dim):
+ i = 0
+ self._buff[...] = 0.0
+ self._buff[ind] = self.velocity.data[d][ind] / \
+ (1.0 + self._coeff[i] * dt)
+ self._buff[ind] *= -1.0
+ self._buff[ind] += self.velocity.data[d][ind]
+ i += 1
+ self.force[d] += npw.real_sum(self._buff[ind])
+ self.force[...] *= self._dvol
+ self.force[...] /= dt
class NocaForces(Forces):
@@ -194,7 +224,7 @@ class NocaForces(Forces):
__metaclass__ = ABCMeta
def __init__(self, velocity, vorticity, nu, volume_of_control,
- normalization=1., surfdir=None, **kwds):
+ surfdir=None, **kwds):
"""
@param vorticity field
@param nu viscosity
@@ -219,9 +249,6 @@ class NocaForces(Forces):
self._coeff = 1. / (self._dim - 1)
## viscosity
self.nu = nu
- # Normalizing coefficient for forces
- # (based on the physics of the flow)
- self._normalization = normalization
#
self._dvol = npw.prod(self._topology.mesh.space_step)
# function to compute the laplacian of a
@@ -338,9 +365,11 @@ class NocaForces(Forces):
j1 = [YDIR, ZDIR, XDIR]
j2 = [ZDIR, XDIR, YDIR]
x0 = coords[i_n].flat[0]
+ print 'buff', buff
for j in i_t:
- buff[...] = vd[j2[j]][ind] * self.vorticity.data[j1[j]][ind]\
- - vd[j1[j]][ind] * self.vorticity.data[j2[j]][ind]
+# buff[...] = vd[j2[j]][ind] * self.vorticity.data[j1[j]][ind]\
+# - vd[j1[j]][ind] * self.vorticity.data[j2[j]][ind]
+ buff[...] = vd[j2[j]][ind]
res[i_n] -= self._coeff * normal * npw.real_sum(coords[j] * buff)
res[j] += self._coeff * normal * x0 * npw.real_sum(buff)
@@ -389,7 +418,7 @@ class NocaForces(Forces):
class NocaI(NocaForces):
def _set_work_arrays(self, rwork=None, iwork=None):
- self._allocate_work_for_controlbox(self, rwork, iwork)
+ self._allocate_work_for_controlbox(rwork, iwork)
def _noca(self, dt):
"""
@@ -431,13 +460,11 @@ class NocaI(NocaForces):
self._compute_gamma_common(s_x, s_id, self._buffer)
self.force += self._buffer * dsurf
- self.force *= self._normalization
-
class NocaII(NocaForces):
def _set_work_arrays(self, rwork=None, iwork=None):
- self._allocate_work_for_controlbox(self, rwork, iwork)
+ self._allocate_work_for_controlbox(rwork, iwork)
def _init_buffers(self):
"""
@@ -505,7 +532,6 @@ class NocaII(NocaForces):
# Then the 'common' part (same in all Noca's formula)
self._compute_gamma_common(s_x, s_id, self._buffer)
self.force += self._buffer * dsurf
- self.force *= self._normalization
self._update_buffers()
def _compute_gamma_momentum(self, dt, surf, s_id, res):
@@ -602,7 +628,6 @@ class NocaIII(NocaForces):
self._compute_gamma_common(s_xn, 1, 2 * d + 1, self._buffer)
self.force += self._buffer * dsurf
- self.force *= self._normalization
self._update_buffers()
def _compute_gamma_flux(self, dt, surf, normal, s_id, res):
diff --git a/HySoP/hysop/operator/discrete/penalize_vorticity.py b/HySoP/hysop/operator/discrete/penalize_vorticity.py
index f0fa15ddf5d60fd51be889e724134909b7b20a03..70c45739727c3727bacba222a40245cdfd7bf22b 100644
--- a/HySoP/hysop/operator/discrete/penalize_vorticity.py
+++ b/HySoP/hysop/operator/discrete/penalize_vorticity.py
@@ -14,7 +14,7 @@ class PenalizeVorticity(Penalization):
"""
@debug
- def __init__(self, vorticity, curl, **kwds):
+ def __init__(self, vorticity, velocity, buff, curl, **kwds):
"""
@param[in,out] vorticity : vorticity field
@param[in] curl : curl operator used to compute
@@ -22,24 +22,33 @@ class PenalizeVorticity(Penalization):
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
super(PenalizeVorticity, self).__init__(variables=[vorticity,
- curl.invar],
+ velocity, buff],
**kwds)
- self.velocity = curl.invar
+ self.velocity = velocity
self.vorticity = vorticity
+ # \warning : a buffer is added for invar variable in curl
+ self.buff0 = buff
topo = self.velocity.topology
msg = 'Multiresolution not implemented for penalization.'
assert self.vorticity.topology == topo, msg
self._curl = curl
def _apply_single_coeff(self, dt):
+ # Vorticity penalization
+ # \warning : the buff0 array ensures "invar" to be 0
+ # outside the obstacle for the curl evaluation
+ for d in xrange(self.buff0.nbComponents):
+ self.buff0.data[d][...] = 0.0
coeff = -dt * self._coeff / (1.0 + dt * self._coeff)
- print coeff
- self._penalize(self.velocity, coeff, self._cond)
+ for d in xrange(self.buff0.nbComponents):
+ self.buff0.data[d][self._cond] = self.velocity[d][self._cond] * coeff
self._curl.apply()
for d in xrange(self.vorticity.nbComponents):
- self.vorticity[d][self._cond] += self._curl.outvar[d][self._cond]
- coeff = -1. / (dt * self._coeff)
- self._penalize(self.velocity, coeff, self._cond)
+ self.vorticity[d][...] += self._curl.outvar[d][...]
+ # Velocity penalization (nor used for the moment since
+ # it implies a non divergence free velocity field)
+# coeff = 1. / (1. + dt * self._coeff)
+# self._penalize(self.velocity, coeff, self._cond)
def _apply_multi_coeff(self, dt):
for i in xrange(len(self._cond)):
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
index cea72bb6159df9264f401174c5561babdd8a4dc2..5652130f89de4c37e61e3e8f926471d123e5a1af 100644
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -16,15 +16,21 @@ class Forces(Computational):
__metaclass__ = ABCMeta
- def __init__(self, obstacles, **kwds):
+ def __init__(self, obstacles, normalization=1., **kwds):
"""
@param obstacles a list of hysop.domain.obstacles inside
the box
+ @param normalization : a normalization coefficient applied to the force
+ (default = 1.)
"""
+ # \warning : the normalization paramater is now an instance of the "Forces" class
super(Forces, self).__init__(**kwds)
self.input = self.variables
## List of hysop.domain.subsets, obstacles to the flow
self.obstacles = obstacles
+ ## Normalizing coefficient for forces
+ ## (based on the physics of the flow)
+ self.normalization = normalization
# Minimal number of ghost points required. Defined in derived class
self._min_ghosts = 0
@@ -67,12 +73,14 @@ class MomentumForces(Forces):
__metaclass__ = ABCMeta
- def __init__(self, velocity, **kwds):
+ def __init__(self, velocity, coeff, **kwds):
"""
@param velocity : velocity continuous vector field
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
self.velocity = velocity
+ # penalization coefficient (lambda)
+ self.coeff = coeff
super(MomentumForces, self).__init__(variables=[velocity], **kwds)
def get_work_properties(self):
@@ -85,7 +93,9 @@ class MomentumForces(Forces):
import MomentumForces as DMF
self.discrete_op = DMF(
velocity=self.discreteFields[self.velocity],
- obstacles=self.obstacles)
+ coeff=self.coeff,
+ obstacles=self.obstacles,
+ normalization=self.normalization)
# output setup
self._set_io('drag_and_lift', (1, 4))
@@ -103,15 +113,13 @@ class NocaForces(Forces):
__metaclass__ = ABCMeta
def __init__(self, velocity, vorticity, nu, volume_of_control=None,
- normalization=1., surfdir=None, **kwds):
+ surfdir=None, **kwds):
"""
@param velocity field
@param vorticity field
@param nu viscosity
@param a volume of control
(hysop.domain.obstacle.controlBox.ControlBox object)
- @param normalization : a normalization coefficient applied to the force
- (default = 1.)
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
## velocity field
@@ -122,9 +130,6 @@ class NocaForces(Forces):
**kwds)
## viscosity
self.nu = nu
- ## Normalizing coefficient for forces
- ## (based on the physics of the flow)
- self.normalization = normalization
# setup for finite differences
if self.method is None:
@@ -147,9 +152,9 @@ class NocaForces(Forces):
if volume_of_control is None:
lr = self.domain.length * 0.9
xr = self.domain.origin + 0.04 * self.domain.length
- from hysop.domain.subsets.control_box import ControlBox
- volume_of_control = ControlBox(parent=self.domain,
- origin=xr, length=lr)
+ from hysop.domain.subsets.control_box import ControlBox
+ volume_of_control = ControlBox(parent=self.domain,
+ origin=xr, length=lr)
self._voc = volume_of_control
def get_work_properties(self):
@@ -173,17 +178,17 @@ class NocaForces(Forces):
def setup(self, rwork=None, iwork=None):
if not self._is_uptodate:
from hysop.operator.discrete.drag_and_lift \
- import NocaForces as DNF
+ import NocaForces, NocaI, NocaII, NocaIII
topo = self.discreteFields[self.velocity].topology
self._voc.discretize(topo)
- self.discrete_op = DNF(
+ self.discrete_op = NocaII(
velocity=self.discreteFields[self.velocity],
vorticity=self.discreteFields[self.vorticity],
nu=self.nu,
volume_of_control=self._voc,
- normalization=self.normalization,
+ surfdir=self._surfdir,
obstacles=self.obstacles,
- surfdir=self._surfdir)
+ normalization=self.normalization)
# output setup
self._set_io('drag_and_lift', (1, 4))
diff --git a/HySoP/hysop/operator/penalize_vorticity.py b/HySoP/hysop/operator/penalize_vorticity.py
index 1cd39ab43f03cd606eba48b3c62ec6130443c8ac..d221a62a72221dfd73a1c30bf833a0f50b8e9284 100644
--- a/HySoP/hysop/operator/penalize_vorticity.py
+++ b/HySoP/hysop/operator/penalize_vorticity.py
@@ -33,13 +33,17 @@ class PenalizeVorticity(Penalization):
@param[in,out] vorticity : vorticity field
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
+ penal_buffer0 = Field(domain=velocity.domain,
+ name='penal_buffer0', isVector=True)
super(PenalizeVorticity, self).__init__(variables=[velocity,
- vorticity],
+ vorticity, penal_buffer0],
**kwds)
- ## velocity field
+ # velocity field
self.velocity = velocity
- ## vorticity field
+ # vorticity field
self.vorticity = vorticity
+ # \warning : a buffer is added for invar variable in curl
+ self.penal_buffer0 = penal_buffer0
# A method is required to set how the curl will be computed.
if self.method is None:
self.method = default.DIFFERENTIAL
@@ -67,7 +71,7 @@ class PenalizeVorticity(Penalization):
obs.discretize(topo)
penal_buffer = Field(domain=self.velocity.domain,
name='penal_buffer', isVector=True)
- self._curl = Curl(invar=self.velocity, outvar=penal_buffer,
+ self._curl = Curl(invar=self.penal_buffer0, outvar=penal_buffer,
discretization=topo, method=self.method)
self._curl.discretize()
@@ -80,6 +84,8 @@ class PenalizeVorticity(Penalization):
self._curl.setup(rwork, iwork)
self.discrete_op = PW(
vorticity=self.discreteFields[self.vorticity],
+ velocity=self.discreteFields[self.velocity],
+ buff=self.discreteFields[self.penal_buffer0],
curl=self._curl.discrete_op,
obstacles=self.obstacles,
coeff=self.coeff, rwork=rwork, iwork=iwork)
diff --git a/HySoP/hysop/tools/numpywrappers.py b/HySoP/hysop/tools/numpywrappers.py
index 67c8c0cf5c50eb8396b083cad8d03eab9cf0b786..f126de7fd463847e91ef4da685b5ed9ffec7e497 100644
--- a/HySoP/hysop/tools/numpywrappers.py
+++ b/HySoP/hysop/tools/numpywrappers.py
@@ -31,6 +31,12 @@ def zeros_like(tab):
"""
return np.zeros_like(tab, dtype=tab.dtype, order=ORDER)
+def ones_like(tab):
+ """
+ Wrapper to numpy.ones_like, force order to hysop.constants.ORDER
+ """
+ return np.ones_like(tab, dtype=tab.dtype, order=ORDER)
+
def realempty(tab):
"""
@@ -202,4 +208,4 @@ def unlock(tab):
"""
set tab as a writeable array
"""
- tabs.flags.writeable = True
+ tab.flags.writeable = True