diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index 6d70a253c33e1d6b10affc93e8e4fa5f8c6f55d7..1217d0ab50b5615923667e5b0e31cbd26271197c 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -22,18 +22,18 @@ from parmepy.operator.diffusion import Diffusion
from parmepy.operator.penalization import Penalization
from parmepy.operator.differential import Curl
from parmepy.operator.redistribute import Redistribute
-from parmepy.operator.adapt_timestep import AdaptTimeStep
-from parmepy.problem.navier_stokes import NSProblem
from parmepy.operator.monitors.printer import Printer
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
from parmepy.operator.monitors.compute_forces import Forces
from parmepy.problem.simulation import Simulation
-from parmepy.constants import VTK
+from parmepy.constants import VTK, HDF5
from parmepy.domain.obstacle.planes import PlaneBoundaries
from dataNS_bb import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
OUTPUT_FREQ, FILENAME, PROJ, LCFL, CFL, CURL_METHOD, \
TIMESTEP_METHOD, OBST_Ox, OBST_Oy, OBST_Oz, RADIUS
-import os
+from parmepy.operator.monitors.compute_forces import DragAndLift
+from parmepy.domain.obstacle.controlBox import ControlBox
+
## ----------- A 3d problem -----------
print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
@@ -44,33 +44,17 @@ cos = np.cos
sin = np.sin
## Domain
-#box = pp.Box(dim, length=[1., 1., 1.], origin=[0., 0., 0.])
-box = pp.Box(dim, length=[2.0 * pi, pi, pi], origin=[0., 0., 0.])
-#box = pp.Box(dim, length=[12., 10., 10.], origin=[-4., -5., -5.])
-#box = pp.Box(dim, length=[2.0 * pi, 2.0 * pi, 2.0 * pi])
+boxlength = [2.0 * pi, pi, pi]
+boxorigin = [0., 0., 0.]
+box = pp.Box(dim, length=boxlength, origin=boxorigin)
## Global resolution
-#nbElem = [nb] * dim
nbElem = [65, 33, 33]
# Upstream flow velocity
uinf = 1.0
-
# Function to compute potential flow past a sphere
-def computeVel0(res, x, y, z, t):
- module = np.sqrt((x - OBST_Ox) * (x - OBST_Ox) +
- (y - OBST_Oy) * (y - OBST_Oy) +
- (z - OBST_Oz) * (z - OBST_Oz))
- ind = np.where(np.abs(module) < 1e-12)
- module[ind] = RADIUS
- res[0][...] = - ((uinf * x) / module) * (1 - (RADIUS ** 3 / module ** 3))
- res[1][...] = ((uinf * y) / module) * (1 + (RADIUS ** 3
- / (2. * module ** 2)))
- res[2][...] = 0.
- return res
-
-
def computeVel(res, x, y, z, t):
res[0][...] = uinf
res[1][...] = 0.0
@@ -85,21 +69,7 @@ def computeVort(res, x, y, z, t):
res[2][...] = 0.
return res
-## Function to compute TG velocity
-#def computeVel(res, x, y, z, t):
-# res[0][...] = sin(x) * cos(y) * cos(z)
-# res[1][...] = - cos(x) * sin(y) * cos(z)
-# res[2][...] = 0.
-# return res
-
-## Function to compute reference vorticity
-#def computeVort(res, x, y, z, t):
-# res[0][...] = - cos(x) * sin(y) * sin(z)
-# res[1][...] = - sin(x) * cos(y) * sin(z)
-# res[2][...] = 2. * sin(x) * sin(y) * cos(z)
-# return res
-
-## Fields
+## Vector Fields
velo = Field(domain=box, formula=computeVel,
name='Velocity', isVector=True)
vorti = Field(domain=box, formula=computeVort,
@@ -117,82 +87,74 @@ ghosts = np.ones((box.dimension)) * NBGHOSTS
topo = Cartesian(box, box.dimension, nbElem,
ghosts=ghosts)
-## Obstacles (sphere + up and down plates)
-sphere = Sphere(box, position=[OBST_Ox, OBST_Oy, OBST_Oz],
- radius=RADIUS)
+## === Obstacles (sphere + up and down plates) ===
+sphere_pos = (boxlength - boxorigin) * 0.5
+sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
bc = PlaneBoundaries(box, 2, thickness=0.1)
-## Operators
-advec = Advection(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- method=ADVECTION_METHOD
- )
-
-stretch = Stretching(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- topo=topo
- )
-
-diffusion = Diffusion(vorti,
- resolution=nbElem,
- viscosity=VISCOSITY)
-
-poisson = Poisson(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- projection=PROJ)
-
-curl = Curl(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- method=CURL_METHOD)
-# topo=topo)
-
-curl.discretize()
-advec.discretize()
-stretch.discretize()
-diffusion.discretize()
-poisson.discretize()
-
-
-## Topology without ghost points
-topofft = poisson.discreteFields[poisson.vorticity].topology
-## Add monitors and setting solver to Problem
-vorti.setTopoInit(topofft)
-velo.setTopoInit(topofft)
-velo.initialize()
-ind = sphere.discretize(topofft)
-vd = velo.discreteFields[topofft].data
-penal = Penalization(velo, [sphere, bc],
- coeff=[1e8],
- topo=topofft,
- resolutions={velo: nbElem})
-penal.discretize()
-correction = VelocityCorrection(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- uinf=uinf,
- topo=topofft)
-
-correction.discretize()
-
-## Bridges between the different topologies in order to
-## redistribute data.
+## === Computational Operators ===
+op = {}
+op['advection'] = Advection(velo, vorti,
+ resolutions={velo: nbElem, vorti: nbElem},
+ method=ADVECTION_METHOD)
+
+op['stretching'] = Stretching(velo, vorti,
+ resolutions={velo: nbElem, vorti: nbElem},
+ topo=topo)
+
+op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
+
+op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
+ projection=PROJ)
+
+op['curl'] = Curl(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
+ method=CURL_METHOD)
+
+## Discretisation of computational operators
+for ope in op:
+ ope.discretize()
+
+# Get fft topology and use it for penalization and correction operators
+opfft = op['poisson']
+topofft = opfft.discreteFields[opfft.vorticity].topology
+
+op['penal'] = Penalization(velo, [sphere, bc], coeff=[1e8], topo=topofft,
+ resolutions={velo: nbElem})
+
+op['correc'] = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem, vorti: nbElem},
+ uinf=uinf, topo=topofft)
+
+op['penal'].discretize()
+op['correc'].discretize()
+
+## === Bridges between the different topologies ===
+distr = {}
# 1 - Curl_fft to stretching (velocity only)
-distrCurlStr_velo = Redistribute([velo], curl, stretch)
-distrAdvStr_velo = Redistribute([velo], advec, stretch)
+distr['Curl2Str'] = Redistribute([velo], op['curl'], op['stretching'])
# 2 - Advection to stretching (vorticity only)
-distrAdvStr_vorti = Redistribute([vorti], advec, stretch)
+distr['Advec2Str_v'] = Redistribute([velo], op['advection'], op['stretching'])
+distr['Advec2Str_w'] = Redistribute([vorti], op['advection'], op['stretching'])
# 3 - Stretching to Diffusion
-distrStrDiff = Redistribute([vorti], stretch, diffusion)
+distr['stretch2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
# + - Brigdes required if Curl op. is solved using a FD method
distrPoissCurl = Redistribute([vorti, velo], poisson, curl)
distrCurlAdv = Redistribute([vorti, velo], curl, advec)
-## Diagnostics/monitors related to the problem
+distrAdvStr_velo.discretize()
+distrAdvStr_vorti.discretize()
+distrStrDiff.discretize()
+distrCurlStr_velo.discretize()
+distrCurlAdv.discretize()
+
+## === Fields initialization ===
+velo.initialize(topo=topofft)
+
+ind = sphere.discretize(topofft)
+vdfft = velo.discreteFields[topofft].data
+
+## === Monitoring operators ===
outputdir = 'res_' + str(topofft.size) + 'procs'
pref = outputdir + '/vwfft'
@@ -201,10 +163,11 @@ printerFFT = Printer(variables=[velo, vorti],
topo=topofft,
frequency=1,
prefix=pref,
- formattype=VTK)
+ formattype=HDF5)
+
printerFFT.setUp()
-prefE = outputdir + '/ener'
+prefE = outputdir + '/ener'
energy = Energy_enstrophy(velo, vorti,
topo=topofft,
viscosity=VISCOSITY,
@@ -212,16 +175,23 @@ energy = Energy_enstrophy(velo, vorti,
frequency=1,
prefix=prefE)
+# 3D Control box
+ref_step = topofft.mesh.space_step
+cbpos = boxorigin + boxlength + 10 * ref_step
+cblength = boxlength - 20 * ref_step
+cb = ControlBox(box, cbpos, cblength)
+nu = 0.1
+coeffForce = 1.
+# Monitor for forces
+forces = DragAndLift(velo, vorti, nu, coeffForce, topofft, cb,
+ obstacle=[sphere],
+ filename=outputdir + 'forces.dat')
+
## Simulation with fixed time step
simu = Simulation(tinit=0.0,
tend=5., timeStep=0.005,
iterMax=1000000)
-distrAdvStr_velo.discretize()
-distrAdvStr_vorti.discretize()
-distrStrDiff.discretize()
-distrCurlStr_velo.discretize()
-distrCurlAdv.discretize()
curl.setUp()
advec.setUp()
diff --git a/Examples/NS_bluff_bodies.py b/Examples/NS_bluff_bodies.py
index c9af9a64b050a21144123bc77d86c6842d41429f..2fc9b0e491e2d1d7f89d4e513811e04ed99afb1d 100755
--- a/Examples/NS_bluff_bodies.py
+++ b/Examples/NS_bluff_bodies.py
@@ -29,7 +29,7 @@ from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
from parmepy.operator.monitors.compute_forces import Forces
from parmepy.problem.simulation import Simulation
from parmepy.constants import VTK
-from dataNS_bb import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
+from dataNS_bb import dim, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
OUTPUT_FREQ, FILENAME, PROJ, LCFL, CFL, CURL_METHOD, \
TIMESTEP_METHOD, OBST_Ox, OBST_Oy, OBST_Oz, RADIUS
from parmepy.domain.obstacle.planes import PlaneBoundaries
@@ -57,9 +57,10 @@ uinf = 1.0
# Function to compute potential flow past a sphere
def computeVel(res, x, y, z, t):
- module = np.sqrt((x - OBST_Ox) * (x - OBST_Ox) + \
- (y - OBST_Oy) * (y - OBST_Oy) + \
- (z - OBST_Oz) * (z - OBST_Oz))
+ module = np.sqrt((x - OBST_Ox) * (x - OBST_Ox) +
+ (y - OBST_Oy) * (y - OBST_Oy) +
+ (z - OBST_Oz) * (z - OBST_Oz))
+
# print 'module', module, np.shape(module)
for i in range (np.shape(module)[0]):
for j in range (np.shape(module)[1]):
diff --git a/Examples/demo_mpi.py b/Examples/demo_mpi.py
index b85ac15e9adf4598a38662dfa4fa5d31851fa5f2..38d3f8c66527ce1002db8d54b9f883941a1d447d 100644
--- a/Examples/demo_mpi.py
+++ b/Examples/demo_mpi.py
@@ -6,7 +6,7 @@ from parmepy.mpi.topology import Cartesian
from parmepy.operator.analytic import Analytic
from parmepy.operator.redistribute import Redistribute
-assert main_size == 4, 'Use 4 process for this test'
+#assert main_size == 4, 'Use 4 process for this test'
def func(res, x, y, z, t=0.):
@@ -20,18 +20,27 @@ b = pp.Box()
f = Field(domain=b, name="Test_Vec", isVector=True, formula=func)
nb_elem = [17, ] * 3
-topo_q = Cartesian.withResolutionFixed(
- b,
- shape=npw.integerarray([2, 1, 2]),
- globalMeshResolution=npw.integerarray(nb_elem),
- )
-topo_p = Cartesian.withResolutionFixed(
- b,
- shape=npw.integerarray([1, 4, 1]),
- globalMeshResolution=npw.integerarray(nb_elem),
- )
-topo_q.setUp()
-topo_p.setUp()
+## topo_q = Cartesian.withResolutionFixed(
+## b,
+## shape=npw.integerarray([2, 1, 2]),
+## globalMeshResolution=[17, 17, 17]
+## )
+## topo_p = Cartesian.withResolutionFixed(
+## b,
+## shape=npw.integerarray([1, 4, 1]),
+## globalMeshResolution=[17, 17, 17]
+## )
+topo_q = Cartesian(b, dim=2, globalMeshResolution=[17, 17, 17])#, cutdir=[False, True, True])
+
+topo_p = Cartesian(b, dim=1, globalMeshResolution=[17, 17, 17])
+if not topo_p.isNew:
+ topo_p = topo_q
+
+#print topo_q
+
+#print topo_p
+
+
op_p = Analytic([f], {f: nb_elem}, topo=topo_p)
op_q = Analytic([f], {f: nb_elem}, topo=topo_q)
op_p.discretize()
@@ -49,8 +58,7 @@ op_p2q.setUp()
print main_rank, 'topo_p', topo_p
print main_rank, 'topo_q', topo_q
-f.setTopoInit(topo_p)
-f.initialize()
+f.initialize(topo=topo_p)
op_p2q.apply()
op_p2q.wait()
diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py
index 5d14340a780cbb3185963afb9fcbedeb72138777..132c22ca439f8d10f1381757c72b8df818b33d0b 100644
--- a/HySoP/hysop/domain/obstacle/sphere.py
+++ b/HySoP/hysop/domain/obstacle/sphere.py
@@ -13,7 +13,7 @@ class Sphere(Obstacle):
"""
def __init__(self, domain, position, radius=1.0,
- vd=0.0, porousLayers=[]):
+ vd=0.0, porousLayers=None):
"""
Description of a sphere in a domain.
@param domain : the physical domain that contains the sphere.
@@ -41,6 +41,8 @@ class Sphere(Obstacle):
self.chi = [dist]
## List of thicknesses for porous layers
+ if porousLayers is None:
+ porousLayers = []
self.layers = porousLayers
def discretize(self, topo):
@@ -81,7 +83,7 @@ class HemiSphere(Sphere):
x < xs for xs == x position of the center of the sphere.
"""
def __init__(self, domain, position, radius=1.0,
- vd=0.0, porousLayers=[]):
+ vd=0.0, porousLayers=None):
"""
Description of a sphere in a domain.
@param domain : the physical domain that contains the sphere.
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 46e33e45efcd83e16f69c2f2aba0162aef8067d9..a60a3d6611f784786e4ef5d9fbc93065a3811dc5 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -133,6 +133,8 @@ class Field(object):
"""
if topo is None:
topo = self.topoInit
+ else:
+ self.topoInit = topo
df = self.discretization(topo)
df.initialize(self._formula, self.doVectorize, currentTime,
*self.extraParameters)
diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index 82be9a6a5be833b8a937e4a14d4d39f09866ba3f..e3b3fd95b5874a0f9f44a17102585a9aba05fac4 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -207,10 +207,14 @@ class Cartesian(object):
# If it already exist (in the sense of the comparison operator
# of the present class) then isNew is false and
# and the present instance should be destroyed.
- self.domain.register(self) # isNew = ...
- # if isNew >= 0:
- # print "A similar topology already exist in the domain's list.\
- # You'd rather destroy this one and use topo of id", isNew
+ newId = self.domain.register(self)
+ self.isNew = True
+ if newId >= 0:
+ s = "A similar topology already exist in the domain's list."
+ s += "You'd rather destroy this one and use topo of id"
+ s += str(newId)
+ print s
+ self.isNew = False
def getParent(self):
"""
diff --git a/HySoP/hysop/operator/redistribute.py b/HySoP/hysop/operator/redistribute.py
index 1722ccc4e87bbf5c14a6f0ac16f22c0fae84251d..c989275b4955ac8b951dab8522f18dd1b44ae477 100644
--- a/HySoP/hysop/operator/redistribute.py
+++ b/HySoP/hysop/operator/redistribute.py
@@ -176,6 +176,7 @@ class Redistribute(Operator):
self._isUpToDate = True
def _apply_toHost_bridges_toDevice(self, simulation=None):
+
if __VERBOSE__:
print "{"+str(main_rank)+"}", "_apply_toHost_bridges_toDevice"
self._toHost()