No commit message

No commit message

Examples/NavierStokes3d_sphere.py

+#!/usr/bin/python
+import parmepy as pp
+import parmepy.f2py
+import numpy as np
+import mpi4py.MPI as MPI
+import math as m
+
+PARMES_REAL = pp.constants.PARMES_REAL
+ORDER = pp.constants.ORDER
+
+#from numpy import linalg as LA
+
+pi = m.pi
+
+## Import scales and fftw solvers
+ppfft = parmepy.f2py.fftw2py
+scales = parmepy.f2py.scales2py
+
+rank = MPI.COMM_WORLD.Get_rank()
+print "Mpi process number ", rank
+
+## ----------- A 3d problem -----------
+print " ========= Start test for Navier-Stokes 3D ========="
+
+## Opening/reading input data file
+inputData = open("input.dat", "r")
+
+Ox = np.float64(inputData.readline().rstrip('\n\r'))
+Oy = np.float64(inputData.readline().rstrip('\n\r'))
+Oz = np.float64(inputData.readline().rstrip('\n\r'))
+Lx = np.float64(inputData.readline().rstrip('\n\r'))
+Ly = np.float64(inputData.readline().rstrip('\n\r'))
+Lz = np.float64(inputData.readline().rstrip('\n\r'))
+resX = np.uint32(inputData.readline().rstrip('\n\r'))
+resY = np.uint32(inputData.readline().rstrip('\n\r'))
+resZ = np.uint32(inputData.readline().rstrip('\n\r'))
+timeStep = np.float64(inputData.readline().rstrip('\n\r'))
+finalTime = np.float64(inputData.readline().rstrip('\n\r'))
+outputModulo = np.uint32(inputData.readline().rstrip('\n\r'))
+outputFilePrefix = inputData.readline().rstrip('\n\r')
+nbGhosts = np.uint32(inputData.readline().rstrip('\n\r'))
+visco = np.float64(inputData.readline().rstrip('\n\r'))
+uinf = np.float64(inputData.readline().rstrip('\n\r'))
+obstName = inputData.readline().rstrip('\n\r')
+lambdFluid = np.float64(inputData.readline().rstrip('\n\r'))
+lambdPorous = np.float64(inputData.readline().rstrip('\n\r'))
+lambdSolid = np.float64(inputData.readline().rstrip('\n\r'))
+obstRadius = np.float64(inputData.readline().rstrip('\n\r'))
+obstOx = np.float64(inputData.readline().rstrip('\n\r'))
+obstOy = np.float64(inputData.readline().rstrip('\n\r'))
+obstOz = np.float64(inputData.readline().rstrip('\n\r'))
+obstOrientation = inputData.readline().rstrip('\n\r')
+thicknPorous = np.float64(inputData.readline().rstrip('\n\r'))
+thicknZbound = np.float64(inputData.readline().rstrip('\n\r'))
+
+inputData.close()
+
+## Physical Domain description
+myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[Ox, Oy, Oz])
+resolution3d = np.asarray((resX, resY, resZ))
+ncells = resolution3d - 1
+hx = Lx / ncells[0]
+hy = Ly / ncells[1]
+hz = Lz / ncells[2]
+
+## Obstacle
+lambd = np.array([lambdFluid, lambdPorous, lambdSolid], 
+                 dtype=PARMES_REAL, order=ORDER)
+sphere = pp.Obstacle(myDomain3d,
+                     name=obstName,
+                     zlayer=thicknZbound,
+                     radius=obstRadius,
+                     center=[obstOx, obstOy, obstOz],
+                     orientation=obstOrientation,
+                     porousLayer=thicknPorous)
+
+### Post
+#outputFilePrefix = './res/NS_'
+#outputModulo = 50
+
+### Simulation parameters
+#timeStep = 1.#1e-1
+#finalTime = 10000.#20.*timeStep
+
+
+## Topologies definitions
+# 1---- FFT topology + diffusion/poisson solvers initialization -----------
+localres, localoffset = ppfft.init_fftw_solver(resolution3d,
+                                               myDomain3d.length)
+print "FFT solver local resolution/offset: ", localres, localoffset
+##topofft = poisson.getTopology()
+
+# ---- Cartesian topology -----
+topoCart = pp.CartesianTopology(domain=myDomain3d,
+                               resolution=resolution3d, dim=3, ghosts=[nbGhosts,nbGhosts,nbGhosts])
+
+# ---- JB's topology ----------
+nbProcs = MPI.COMM_WORLD.Get_size()
+topodims = [1, 1, nbProcs]  # fits with fftw topology
+
+
+## Fields declaration
+def computeVel(x, y, z):
+#    vx = 2. / np.sqrt(3) * np.sin(2. * pi / 3.) * np.sin(x) \
+#        * np.cos(y) * np.cos(z)
+#    vy = 2. / np.sqrt(3) * np.sin(-2. * pi / 3.) * np.cos(x) \
+#        * np.sin(y) * np.cos(z)
+#    vz = 0.
+#------------------
+    vx = 0.
+    module = (x-obstOx)**2 + (y-obstOy)**2 + (z-obstOz)**2
+    if (module >= obstRadius**2):
+        vx = uinf * (1-(obstRadius**2/module))
+    vy = 0.
+    vz = 0.
+#-----------------
+#    vx = 0.
+#    vy = 0.
+#    vz = 0.
+    return vx, vy, vz
+
+
+def computeVort(x, y, z):
+#    wx = - np.cos(x) * np.sin(y) * np.sin(z)
+#    wy = - np.sin(x) * np.cos(y) * np.sin(z)
+#    wz = 2. * np.sin(x) * np.sin(y) * np.cos(z)
+#-----------------
+    wx = 0.
+    wy = 0.
+    wz = 0.
+#-----------------
+#    xc = 6./2.
+#    yc = 6./2.
+#    zc = 6./6.
+#    R = 1.5
+#    sigma = R / 3.
+#    Gamma = 0.0075
+#    dist = m.sqrt((x-xc)**2 + (y-yc)**2)
+#    s2 = (z - zc)**2 + (dist - R)**2
+#    wx = 0.
+#    wy = 0.
+#    wz = 0.
+#    if (dist != 0.):
+#        cosTheta = (x-xc) / dist
+#        sinTheta = (y-yc) / dist
+#        wTheta = Gamma / (pi * sigma**2) * m.exp(-(s2 / sigma**2))
+#        wx = - wTheta * sinTheta
+#        wy = wTheta * cosTheta
+#        wz = 0.
+    return wx, wy, wz
+
+def computeVortNorm(x, y, z):
+    norm = 0.
+#----------------
+#    xc = 6./2.
+#    yc = 6./2.
+#    zc = 6./6.
+#    R = 1.5
+#    sigma = R / 3.
+#    Gamma = 0.0075
+#    dist = m.sqrt((x-xc)**2 + (y-yc)**2)
+#    s2 = (z - zc)**2 + (dist - R)**2
+#    norm = 0.
+#    if (dist != 0.):
+#        cosTheta = (x-xc) / dist
+#        sinTheta = (y-yc) / dist
+#        wTheta = Gamma / (pi * sigma**2) * m.exp(-(s2 / sigma**2))
+#        wx = - wTheta * sinTheta
+#        wy = wTheta * cosTheta
+#        wz = 0.
+#        norm = np.sqrt(wx ** 2 + wy ** 2 + wz ** 2)
+    return norm
+
+#def computeDensity(x, y, z):
+#    if (x>=0):
+#        rho = 1000.
+#    else : 
+#        rho = 1.
+#    return rho
+
+velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
+                              name='Velocity', vector=True)
+vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
+                               name='Vorticity', vector=True)
+vortNorm = pp.AnalyticalField(domain=myDomain3d, formula=computeVortNorm,
+                             name='VortNorm', vector=False)
+#density = pp.AnalyticalField(domain=myDomain3d, formula=computeDensity,
+#                             name='Density', vector=False)
+
+
+
+## 2 - Advection solver initialization. See testScales for a working example.
+# Based on scales JB solver
+# Warning : fields input for scales should be of size (ncells), not localres.
+
+scalesres, scalesoffset, stab_coeff = \
+    scales.init_advection_solver(ncells, myDomain3d.length,
+                                 topodims, order='p_O2')
+
+print 'advection stability coeff', stab_coeff
+advecVort = pp.Advec_scales(stab_coeff, velocity, vorticity, vortNorm)
+#advecDensity = pp.Advec_scales(stab_coeff, velocity, density) 
+
+# ----> Change TOPO from Scales to Cartesian
+
+## 3 - Stretching
+stretch = pp.Stretching(velocity, vorticity)
+
+# ----> Change TOPO from Cartesian to FFT
+
+## 4 - Diffusion
+diffusion = pp.Diffusion(velocity, vorticity, viscosity=visco)
+
+## 4 - Poisson
+poisson = pp.Poisson(velocity, vorticity)
+
+# ----> Change TOPO from FFT to Cartesian
+
+## 6 - Penalization
+penal = pp.Penalization(velocity, vorticity, sphere, lambd)
+
+# ----> Change TOPO from Cartesian to Scales
+
+## Define the problem to solve
+navierStokes = pp.Problem(topoCart, [advecVort, stretch, diffusion, poisson, penal])
+#navierStokes = pp.Problem(topoCart, [advecVort, stretch, poisson])
+
+printer = pp.Printer(fields=[vorticity, velocity, vortNorm], frequency=outputModulo,
+                     outputPrefix=outputFilePrefix)
+
+navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
+
+## Problem => ParticularSover = basic.initialize()
+navierStokes.initSolver()
+
+penal.apply(0., timeStep)
+#poisson.apply(0., timeStep)
+
+## solve stretching
+navierStokes.solve()
+
+## end of time loop ##
+
+# Clean memory buffers
+ppfft.clean_fftw_solver(myDomain3d.dimension)

Examples/input.dat

-0.0
-0.0
-0.0
-6.0
-6.0
-6.0
-257
-257
-257
-1.
-10000.
-50
+-6.0
+-6.0
+-6.0
+12.0
+12.0
+12.0
+65
+65
+65
+0.02
+20000.
+10
 ./res/NS_
 2
-0.000001
+0.005
+1.
+sphere
+0.
+10.
+10E8
+1.0
+0.
+0.
+0.
+West
+0.
+0.05

HySoP/hysop/operator/diffusion_d.py

@@ -8,7 +8,6 @@ from discrete import DiscreteOperator
 from parmepy.constants import ORDER, PARMES_REAL
 import numpy as np
 import time
-import numpy.linalg as la
 
 
 class Diffusion_d(DiscreteOperator):
@@ -39,7 +38,6 @@ class Diffusion_d(DiscreteOperator):
         Apply operator.
         """
         self.compute_time = time.time()
-        print "Norm de vz pre ... ", la.norm(self.velocity[2])
 
         ind0a = self.ghosts[0]
         ind0b = self.resolution[0] - self.ghosts[0]
@@ -60,22 +58,27 @@ class Diffusion_d(DiscreteOperator):
             velocityNoG[i][...] = self.velocity[i][ind0a:ind0b,
                                                    ind1a:ind1b, ind2a:ind2b]
 
-        # Vorticity diffusion
+        # Curl + Vorticity diffusion
+#        vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
+#            fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+#                                       velocityNoG[0][...],
+#                                       velocityNoG[1][...],
+#                                       velocityNoG[2][...],
+#                                       vorticityNoG[0][...],
+#                                       vorticityNoG[1][...],
+#                                       vorticityNoG[2][...])
+
+        # Pure Vorticity diffusion
         vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
             fftw2py.solve_diffusion_3d(self.viscosity * dt,
-                                       velocityNoG[0][...],
-                                       velocityNoG[1][...],
-                                       velocityNoG[2][...],
-                                       vorticityNoG[0][...],
-                                       vorticityNoG[1][...],
-                                       vorticityNoG[2][...])
+                                                  vorticityNoG[0][...],
+                                                  vorticityNoG[1][...],
+                                                  vorticityNoG[2][...])
 
         for i in xrange(self.dim):
             self.vorticity[i][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
                 vorticityNoG[i][...]
 
-        print "Norm de vz post ... ", la.norm(self.velocity[2])
-
         self.compute_time = time.time() - self.compute_time
         self.total_time += self.compute_time
         return self.compute_time