Everytging concerning Noca forces computation (still some bugs to fix)

Examples/NS_bluff_bodies.py

+#!/usr/bin/python
+
+"""
+Navier Stokes 3D : flow past bluff bodies (using penalization).
+
+All parameters are set and defined in python module dataNS_bb.
+
+"""
+
+import parmepy as pp
+from parmepy.f2py import fftw2py
+import numpy as np
+import math as m
+from parmepy.fields.continuous import Field
+from parmepy.mpi.topology import Cartesian
+from parmepy.domain.obstacle.sphere import Sphere
+from parmepy.domain.obstacle.plates import Plates
+from parmepy.operator.advection import Advection
+from parmepy.operator.stretching import Stretching
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.operator.penalization import Penalization
+from parmepy.operator.curl import Curl
+from parmepy.operator.redistribute import Redistribute
+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 dataNS_bb import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
+    WITH_PROJ, OUTPUT_FREQ, FILENAME, simu, OBST_Ox, OBST_Oy, OBST_Oz, RADIUS
+from parmepy.constants import CONSERVATIVE
+
+
+## ----------- A 3d problem -----------
+print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
+
+## Domain
+box = pp.Box(dim, length=[1., 1., 1.], origin=[0., 0., 0.])
+
+## Global resolution
+nbElem = [nb] * dim
+
+# Upstream flow velocity
+uinf = 1.0
+
+## Function to compute potential flow past a sphere
+def computeVel(x, y, z):
+    vx = 0.
+    module = (x - OBST_Ox) * (x - OBST_Ox) + \
+             (y - OBST_Oy) * (y - OBST_Oy) + \
+             (z - OBST_Oz) * (z - OBST_Oz)
+    if (module >= RADIUS * RADIUS):
+        vx = uinf * (1 - (RADIUS * RADIUS/ module))
+    vy = 0.
+    vz = 0.
+    return vx, vy, vz
+
+## Function to compute initial vorticity
+def computeVort(x, y, z):
+    wx = 0.
+    wy = 0.
+    wz = 0.
+#    module = (x - OBST_Ox) * (x - OBST_Ox) + \
+#             (y - OBST_Oy) * (y - OBST_Oy) + \
+#             (z - OBST_Oz) * (z - OBST_Oz)
+#    if (module >= RADIUS * RADIUS):
+#        wy = (2.0 * uinf * RADIUS * RADIUS * (z - OBST_Oz)) / (module * module)
+#        wz = -(2.0 * uinf * RADIUS * RADIUS * (y - OBST_Oy)) / (module * module)
+    return wx, wy, wz
+
+## Fields
+velo = Field(domain=box, formula=computeVel,
+             name='Velocity', isVector=True)
+vorti = Field(domain=box, formula=computeVort,
+              name='Vorticity', isVector=True)
+
+## Usual Cartesian topology definition
+# At the moment we use two (or three?) topologies :
+# - "topo" for Stretching and all operators based on finite differences.
+#    --> ghost layer = 2
+# - topo from Advection operator for all the other operators.
+#    --> no ghost layer
+# - topo from fftw for Poisson and Diffusion.
+# Todo : check compat between scales and fft operators topologies.
+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)
+
+plates = Plates(box, normal_dir=2, epsilon=0.01)
+
+## Operators
+advec = Advection(velo, vorti,
+                  resolutions={velo: nbElem,
+                               vorti: nbElem},
+                  method=ADVECTION_METHOD,
+                  topoStab=topo
+                  )
+
+stretch = Stretching(velo, vorti,
+                     resolutions={velo: nbElem,
+                                  vorti: nbElem},
+                     formulation=CONSERVATIVE,
+                     topo=topo
+                     )
+
+diffusion = Diffusion(vorti,
+                      resolution=nbElem,
+                      viscosity=VISCOSITY
+                     )
+
+poisson = Poisson(velo, vorti,
+                  resolutions={velo: nbElem,
+                               vorti: nbElem},
+                  projection=WITH_PROJ)
+
+penal = Penalization(velo, [sphere, plates], 
+                     coeff=[1e8],
+                     resolutions={velo: nbElem})
+
+curl = Curl(velo, vorti,
+            resolutions={velo: nbElem,
+                         vorti: nbElem},
+            topo=topo)
+
+## Diagnostics related to the problem
+
+forces = Forces(velo, vorti, topo, ghostUpdate=True, 
+                obstacle=sphere, boxMin= [0.1, 0.1, 0.1], 
+                boxMax= [0.8, 0.8, 0.8], Reynolds=500.,
+                method='', frequency=1, prefix='./res/Noca_sphere2.dat')
+
+printer = Printer(variables=[velo, vorti],
+                  topo=topo,
+                  frequency=1000,
+                  prefix='./res/NS_sphere',
+                  ext='.vtk')
+
+# Bridges between the different topologies in order to
+# redistribute data.
+# 1 -Advection to stretching
+distrAdvStr = Redistribute([vorti, velo], advec, stretch)
+# 2 - Stretching to Poisson/Diffusion
+distrStrPoiss = Redistribute([vorti, velo], stretch, poisson)
+# 3 - Poisson to Curl
+distrPoissCurl = Redistribute([vorti, velo], poisson, curl)
+# 4 - Curl to Advection
+distrCurlAdv = Redistribute([vorti, velo], curl, advec)
+
+## Define the problem to solve
+pb = NSProblem(operators=[distrCurlAdv, advec, distrAdvStr, stretch, distrStrPoiss,
+                          diffusion, poisson, penal, 
+                          distrPoissCurl, curl],
+               simulation=simu, monitors=[forces], dumpFreq=-1)
+
+## Setting solver to Problem
+pb.setUp()
+
+penal.apply(simu)
+distrPoissCurl.apply(simu)
+curl.apply(simu)
+
+for topo in box.topologies.values():
+    print topo
+
+
+## Solve problem
+pb.solve()
+
+## end of time loop ##
+
+## Clean memory buffers
+fftw2py.clean_fftw_solver(box.dimension)

Examples/dataNS_bb.py

+
+
+# problem dimension
+dim = 3
+# resolution
+nb = 65
+# number of ghosts in usual cartesian topo
+NBGHOSTS = 2
+# Advection method
+ADVECTION_METHOD = 'scales_p_M6'
+#
+VISCOSITY = 1. / 500.
+# Vorticity projection?
+WITH_PROJ = False
+
+# Obstacle description
+OBST_Ox = 0.4
+OBST_Oy = 0.5
+OBST_Oz = 0.5
+RADIUS = 0.15
+
+# post treatment ...
+OUTPUT_FREQ = 1
+FILENAME = './res/Noca_sphere.dat'
+
+# simulation parameters
+from parmepy.problem.simulation import Simulation
+simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.01,
+                  iterMax=1000000)

HySoP/hysop/numerics/differential_operations.py

@@ -457,6 +457,89 @@ class DivT(DifferentialOperation):
 
         return result
 
+class DivStressTensor3D(DifferentialOperation):
+    """
+    Computes the 3D stress tensor T defined by 
+    \f$ T=div(\nabla.(U) + \nabla^{T}.(U)) \f$,
+    \f$U\f$ being a velocity vector field.
+
+    Methods :
+    1 - FD_C_4 (default) : 4th order, centered finite differences, with
+    local workspace and based on fd.compute_and_add.
+    init : func = GradS(topo, method)
+    call : result = func(var, result)
+
+    var stands for the vector field U.
+    result must be a vector field (same number of components as var),
+    each component with the same size as var[i].
+    work must be a three component vector field
+    each component with the same size as var[i].
+
+    """
+    def __init__(self, topo, method=None):
+        if method is not None and method.find('FD_order2') >= 0:
+                raise ValueError("2nd order scheme Not yet implemented")
+        else :
+            # - 4th ordered FD
+            # - with work
+            # - fd.compute_and_add method
+
+            # connect call function
+            self.fcall = self.FDCentral4_with_work
+            self.fd_scheme = FD_C_4((topo.mesh.space_step))
+
+    def __call__(self, var, result, work, ind):
+        assert work[0].shape == var[0].shape
+        return self.fcall(var, result, work, ind)
+
+    def FDCentral4_with_work(self, var, result, work, ind):
+        """
+        """
+        # compute indices according to "ind" list
+        self._indices = ind
+        self.fd_scheme.computeIndices(self._indices)
+
+        # 1st component of the stress tensor
+        #--------grad(var[XDIR])-----------
+        self.fd_scheme.compute(var[XDIR], XDIR, work[XDIR])
+        self.fd_scheme.compute(var[XDIR], YDIR, work[YDIR])
+        self.fd_scheme.compute(var[XDIR], ZDIR, work[ZDIR])
+
+        work[XDIR][...] *= 2.
+        self.fd_scheme.compute_and_add(var[YDIR], XDIR, work[YDIR])
+        self.fd_scheme.compute_and_add(var[ZDIR], XDIR, work[ZDIR])
+
+        for cdir in xrange(1, len(var)):
+            self.fd_scheme.compute_and_add(work[cdir], cdir, result[XDIR])
+
+        # 2nd component of the stress tensor
+        #--------grad(var[YDIR])-----------
+        self.fd_scheme.compute(var[YDIR], XDIR, work[XDIR])
+        self.fd_scheme.compute(var[YDIR], YDIR, work[YDIR])
+        self.fd_scheme.compute(var[YDIR], ZDIR, work[ZDIR])
+
+        self.fd_scheme.compute_and_add(var[XDIR], YDIR, work[XDIR])
+        work[YDIR][...] *= 2.
+        self.fd_scheme.compute_and_add(var[ZDIR], YDIR, work[ZDIR])
+
+        for cdir in xrange(1, len(var)):
+            self.fd_scheme.compute_and_add(work[cdir], cdir, result[YDIR])
+
+        # 3rd component of the stress tensor
+        #--------grad(var[ZDIR])-----------
+        self.fd_scheme.compute(var[ZDIR], XDIR, work[XDIR])
+        self.fd_scheme.compute(var[ZDIR], YDIR, work[YDIR])
+        self.fd_scheme.compute(var[ZDIR], ZDIR, work[ZDIR])
+
+        self.fd_scheme.compute_and_add(var[XDIR], ZDIR, work[XDIR])
+        self.fd_scheme.compute_and_add(var[YDIR], ZDIR, work[YDIR])
+        work[ZDIR][...] *= 2.
+
+        for cdir in xrange(1, len(var)):
+            self.fd_scheme.compute_and_add(work[cdir], cdir, result[ZDIR])
+
+        return result
+
 
 class GradS(DifferentialOperation):
     """
@@ -545,19 +628,21 @@ class GradVxW(DifferentialOperation):
 
             if diagnostics is not None:
                 # maxima of partial derivatives : needed for advection stability condition
-                diagnostics[0] = max(diagnostics[1],
-                                     np.max(abs(data[direc])))
+                diagnostics[0] = max(diagnostics[0],
+                                     np.max(abs(data[direc + direc])))
 
                 # maxima of partial derivatives : needed for stretching stability condition
-                diagnostics[1] = max(diagnostics[0],
+                diagnostics[1] = max(diagnostics[1],
                                      np.max(sum([abs(data[i])
                                           for i in xrange(direc,
                                                           self.nbComponents
                                                           + direc)])))
 
-            # compute grad(Vx).var2
-            data[direc] = self.prod(data[direc:self.nbComponents + direc],
-                                    var2)
+            if var2 is not None:
+                # compute grad(Vx).var2
+                data[direc] = self.prod(data[direc:self.nbComponents + direc],
+                                        var2)
+
         return data, diagnostics
 
 

HySoP/hysop/numerics/tests/test_diffOp.py

@@ -5,7 +5,8 @@ from parmepy.constants import PARMES_REAL, ORDER
 from parmepy.domain.box import Box
 from parmepy.fields.continuous import Field
 from parmepy.mpi.topology import Cartesian
-from parmepy.numerics.differential_operations import DivT, GradVxW, GradS, Curl
+from parmepy.numerics.differential_operations import DivT, GradVxW, GradS, Curl, DivStressTensor3D
+import parmepy.tools.numpywrappers as npw
 import numpy.testing as npt
 import math as m
 
@@ -47,8 +48,8 @@ def testDivWV():
 
     # Div operator
     FD_method = 'FD_order4'
-    StretchOp = DivT(topoG, FD_method)
-    result = StretchOp(velo_d.data, vorti_d.data)
+    StretchOp = DivT(topoG, FD_method, memshape=velo_d.data[0].shape)
+    result = StretchOp(velo_d.data, vorti_d.data, result)
 
     # Numerical VS analytical
     ind = topoG.mesh.iCompute
@@ -161,6 +162,55 @@ def testCurl():
     print np.allclose(analyt_d[1], resCurl[1][ind], rtol=errY)
     print np.allclose(analyt_d[2], resCurl[2][ind], rtol=errZ)
 
+def testDivStressTensor():
+    """Basic test for divergence of stress tensor T 
+    (comparison with periodic analytical solution based on Taylor Green benchmark)
+    """
+    # Domain and topologies definitions
+    nb = 64
+    Lx = Ly = Lz = 2. * np.pi
+    box = Box(dimension=3, length=[Lx, Ly, Lz], 
+              origin=[0., 0., 0.])
+    resol = [nb, nb, nb]
+    nbghosts = 2
+    ghosts = np.ones((box.dimension)) * nbghosts
+    topoG = Cartesian(box, box.dimension, resol,
+                     ghosts=ghosts)
+    topoNoG = Cartesian(box, box.dimension, resol)
+
+    # Velocity, work and result fields
+    velo = Field(domain=box, formula=vitesse, 
+                 name='Velocity', isVector=True)
+    velo_d = velo.discretize(topoG)
+    velo.initialize()
+
+    result = npw.zeros_like(velo_d.data)
+    work = npw.zeros_like(velo_d.data)
+
+    # Analytical field
+    analyt = Field(domain=box, formula=analyticalDivStressTensor, 
+                   name='Analytical', isVector=True)
+    analyt_d = analyt.discretize(topoNoG)
+    analyt.initialize()
+
+    # Grad operator
+    FD_method = 'FD_order4'
+    ind = topoG.mesh.iCompute
+    divT = DivStressTensor3D(topoG, FD_method)
+    result = divT(velo_d.data, result, work, ind)
+
+    # Numerical VS analytical
+    errX = errY = errZ = (Lx / (nb - 1)) ** 4
+    print 'err = O(h**4) =', errX
+    print '=============='
+    print 'div(T) test'
+    print '=============='
+    print 'Is the numerical solution equal to the analytical one at the order 4 ? :'
+    print np.allclose(analyt_d[0], result[0][ind], rtol=errX)
+    print np.allclose(analyt_d[1], result[1][ind], rtol=errY)
+    print np.allclose(analyt_d[2], result[2][ind], rtol=errZ)
+
+
 def vitesse(x, y, z):
     vx = m.sin(x) * m.cos(y) * m.cos(z)
     vy = - m.cos(x) * m.sin(y) * m.cos(z)
@@ -187,7 +237,14 @@ def analyticalGradVxW(x, y, z):
     az = 0.
     return ax, ay, az
 
+def analyticalDivStressTensor(x, y, z):
+    ax = - 3. * m.sin(x) * m.cos(y) * m.cos(z)
+    ay = 2. * m.sin(x) * m.cos(y) * m.cos(z) - m.sin(x) * m.sin(y) * m.sin(z)
+    az = 0.
+    return ax, ay, az
+
 if __name__ == "__main__":
-    testDivWV()
-    testGradVxW()
-    testCurl()
+#    testDivWV()
+#    testGradVxW()
+#    testCurl()
+    testDivStressTensor()

HySoP/hysop/operator/tests/test_forces.py

+# -*- coding: utf-8 -*-
+import parmepy as pp
+import numpy as np
+from parmepy.domain.obstacle.sphere import Sphere
+from parmepy.mpi.topology import Cartesian
+from parmepy.operator.monitors.compute_forces import Forces
+from parmepy.fields.continuous import Field
+from parmepy.problem.simulation import Simulation
+
+
+def computeVel(x, y, z):
+    vx = 1.
+    vy = 1.
+    vz = 1.
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    wx = 0.
+    wy = 0.
+    wz = 0.
+    return wx, wy, wz
+
+
+def testForces3D():
+
+    nb = 33
+    Lx = Ly = Lz = 2
+    dom = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[-1., -1., -1.])
+    resol = [nb, nb, nb]
+    velo = Field(domain=dom, name='Velo', isVector=True)
+    vorti = Field(domain=dom, name='Vorti', isVector=True)
+    sphere = Sphere(dom, position=[0., 0., 0.], radius=0.2)
+
+    ghosts = np.ones((dom.dimension)) * 2.0
+    topo = Cartesian(dom, dom.dimension, resol, ghosts=ghosts)
+
+    forces = Forces(velo, vorti, topo, ghostUpdate=True, 
+                    obstacle=sphere, boxMin= [-0.6, -0.6, -0.6], 
+                    boxMax= [0.6, 0.6, 0.6], Reynolds=200.,
+                    method='', frequency=1, prefix='./NS_sphere')
+
+    forces.setUp()
+    velo.initialize(topo)
+    vorti.initialize(topo)
+    forces.apply(Simulation())
+
+if __name__ == "__main__":
+    testForces3D()