Update topo

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()

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]):

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()
 

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.

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)

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):
         """

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()