update example for flow past sphere

Examples/NSDebug.py

@@ -12,24 +12,26 @@ from parmepy.f2py import fftw2py
 import numpy as np
 from parmepy.fields.continuous import Field
 from parmepy.mpi.topology import Cartesian
+from parmepy.fields.variable_parameter import VariableParameter
 from parmepy.domain.obstacle.sphere import Sphere
 from parmepy.operator.advection import Advection
 from parmepy.operator.stretching import Stretching
 from parmepy.operator.poisson import Poisson
 from parmepy.operator.velocity_correction import VelocityCorrection
+#from parmepy.operator.vorticity_correction import VorticityCorrection
 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.monitors import Printer, Energy_enstrophy, DragAndLift,\
-    Reprojection_criterion
+    Reprojection_criterion, ProfilesOutput
 from parmepy.problem.simulation import Simulation
 from parmepy.constants import VTK, HDF5
-from parmepy.domain.obstacle.planes import PlaneBoundaries
+from parmepy.domain.obstacle.planes import PlaneBoundaries, SubPlane
+from parmepy.domain.obstacle.controlBox import ControlBox
 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
-from parmepy.domain.obstacle.controlBox import ControlBox
 import parmepy.tools.numpywrappers as npw
 import parmepy.tools.io_utils as io
 ## ----------- A 3d problem -----------
@@ -41,17 +43,35 @@ cos = np.cos
 sin = np.sin
 
 ## Domain
-boxlength = npw.realarray([2.0 * pi, 2 * pi, 2 * pi])
-boxorigin = npw.realarray([0., 0., 0.])
+#boxlength = npw.realarray([8.0 * pi, 2.0 * pi, 2.0 * pi])
+#boxorigin = npw.realarray([-2.0, -pi, -pi])
+boxlength = npw.realarray([10.24, 5.12, 5.12])
+boxorigin = npw.realarray([-2.0, -2.56, -2.56])
 box = pp.Box(dim, length=boxlength, origin=boxorigin)
 
 ## Global resolutionsimu
-nbElem = [129]*3
+#nbElem = [1025, 513, 513]
+nbElem = [129, 65, 65]
 
-# Upstream flow velocity
+
+# Initial upstream flow velocity
 uinf = 1.0
 
-# Function to compute potential flow past a sphere
+def computeFlowrate(simu):
+    ##Time-dependant flow rate
+    t = simu.tk
+    Tstart = 3.0
+    flowrate = np.zeros(3)
+    flowrate[0] = uinf * box.length[1] * box.length[2]
+    if t >= Tstart and t <= Tstart + 1.0:
+        flowrate[1] = sin(pi * (t - Tstart)) * \
+                      box.length[1] * box.length[2]
+    ##constant flow rate
+#    flowrate = np.zeros(3)
+#    flowrate[0] = uinf * box.length[1] * box.length[2]
+    return flowrate
+
+# Function to compute initial velocity
 def computeVel(res, x, y, z, t):
     res[0][...] = uinf
     res[1][...] = 0.
@@ -72,16 +92,23 @@ velo = Field(domain=box, formula=computeVel,
 vorti = Field(domain=box, formula=computeVort,
               name='Vorticity', isVector=True)
 
-# Topo for operators that required ghosts points (stretching)
+## 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
-topostr = Cartesian(box, box.dimension, nbElem, ghosts=ghosts)
+topo = Cartesian(box, box.dimension, nbElem,
+                 ghosts=ghosts)
 
 ## === Obstacles (sphere + up and down plates) ===
-sphere_pos = (boxlength - boxorigin) * 0.5
+#sphere_pos = (boxlength - boxorigin) * 0.5
+sphere_pos = [0., 0., 0.]
 sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
 
-bc = PlaneBoundaries(box, 2, thickness=0.1)
-
 ## === Computational Operators ===
 op = {}
 op['advection'] = Advection(velo, vorti,
@@ -90,7 +117,7 @@ op['advection'] = Advection(velo, vorti,
 
 op['stretching'] = Stretching(velo, vorti,
                               resolutions={velo: nbElem, vorti: nbElem},
-                              topo=topostr)
+                              topo=topo)
 
 op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
 
@@ -108,15 +135,20 @@ opfft = op['poisson']
 topofft = opfft.discreteFields[opfft.vorticity].topology
 topocurl = op['curl'].discreteFields[op['curl'].invar].topology
 topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
+# space step
+ref_step = topofft.mesh.space_step
+step_dir = ref_step[0]
 
-op['penalization'] = Penalization(velo, [sphere, bc], coeff=[1e8],
+## Penalization velocity
+op['penalization'] = Penalization(velo, [sphere], coeff=[1e8],
                                   topo=topofft, resolutions={velo: nbElem})
 
-ref_rate = uinf * box.length[1] * box.length[2]
+## Time-dependant required-flowrate (Variable Parameter)
+req_flowrate = VariableParameter(formula=computeFlowrate)
 op['correction'] = VelocityCorrection(velo, vorti,
                                       resolutions={velo: nbElem,
                                                    vorti: nbElem},
-                                      req_flowrate=ref_rate, topo=topofft)
+                                      req_flowrate=req_flowrate, topo=topofft)
 
 op['penalization'].discretize()
 op['correction'].discretize()
@@ -126,37 +158,73 @@ distr = {}
 
 distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
 distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
+
+# 1 - Curl to stretching (velocity only)
+#distr['curl2str'] = Redistribute([velo], op['curl'], op['stretching'])
+# 2 - Advection to stretching
+# velocity only
+#distr['adv2str_v'] = Redistribute([velo], op['advection'], op['stretching'])
+# vorticity only
+#distr['adv2str_w'] = Redistribute([vorti], op['advection'], op['stretching'])
+# Both
 distr['adv2str'] = Redistribute([velo, vorti],
                                 op['advection'], op['stretching'])
+# 3 - Stretching to Diffusion
 distr['str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
 distr['curl2adv'] = Redistribute([velo, vorti], op['curl'], op['advection'])
 distr['curl2str'] = Redistribute([velo, vorti], op['curl'], op['stretching'])
 
+#distr['adv2corr'] = Redistribute([velo, vorti],
+#                                op['advection'], op['correc_vorti'])
+
 for ope in distr.values():
     ope.discretize()
 
 # === Simulation setup ===
-simu = Simulation(tinit=0.0, tend=10., nbIter=2000, iterMax=1000000)
+simu = Simulation(tinit=0.0, tend=75., timeStep=0.0125, iterMax=1000000)
 
 ## === Monitoring operators ===
 monitors = {}
 # Save velo and vorti values on topofft
 monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
-                                 formattype=HDF5, prefix='fft_io')
+                                 formattype=HDF5)
 
 monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
-                                  formattype=HDF5, prefix='curl_io')
-monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topoadvec,
-                                   formattype=HDF5, prefix='advec_io',
+                                  formattype=HDF5, prefix='curl_io_vorti')
+
+monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
+                                   frequency=1, formattype=HDF5, prefix='advec_io',
                                    xmfalways=True)
-monitors['printerStr'] = Printer(variables=[velo, vorti], topo=topostr,
-                                 formattype=HDF5, prefix='str_io',
-                                 xmfalways=True)
+
+thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir], 
+                        lengths=[10.24, 5.12, 4.0 * step_dir])
+#thickSliceXY = ControlBox(box, origin=[-2.0, -pi, -2.0 * step_dir], 
+#                        lengths=[8.0*pi, 2.0*pi, 4.0 * step_dir])
+monitors['printerSliceXY'] = Printer(variables=[velo, vorti], topo=topofft,
+                                   frequency=100, formattype=HDF5, prefix='sliceXY', 
+                                   subset=thickSliceXY, xmfalways=True)
+
+thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56], 
+                        lengths=[10.24, 4.0 * step_dir, 5.12])
+monitors['printerSliceXZ'] = Printer(variables=[velo, vorti], topo=topofft,
+                                   frequency=100, formattype=HDF5, prefix='sliceXZ', 
+                                   subset=thickSliceXZ, xmfalways=True)
+
+subBox = ControlBox(box, origin=[-0.7, -1.25, -1.25], lengths=[7.0, 2.5, 2.5])
+monitors['printerSubBox'] = Printer(variables=[velo, vorti], topo=topofft,
+                                   frequency=100, formattype=HDF5, prefix='subBox', 
+                                   subset=subBox, xmfalways=True)
 # Compute and save enstrophy/Energy, on topofft
 io_default = {}
-monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topoadvec,
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
                                       io_params={},
                                       viscosity=VISCOSITY, isNormalized=False)
+
+# Compute velocity and vorticity magnitude on topofft and print it in data file
+io_default = {}
+monitors['profiles'] = ProfilesOutput(velo, vorti, obstacles=[sphere],
+                                      topo=topofft, io_params={"frequency":1000})
+
 # Setup for reprojection
 reprojection_constant = 0.04
 reprojection_rate = 1
@@ -172,19 +240,17 @@ vdfft = velo.discreteFields[topofft].data
 
 # Compute and save drag and lift on topofft
 # 3D Control box
-ref_step = topofft.mesh.space_step
 cbpos = boxorigin + 10 * ref_step
-cblength = boxlength - 20 * ref_step
+cblength = boxlength - 25 * ref_step
 cb = ControlBox(box, cbpos, cblength)
 coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
 
-print coeffForce
+#print coeffForce
 # Monitor for forces
-## monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-##                                  topo, cb,
-##                                  obstacles=[sphere], io_params={})
 monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-                                 topostr, cb, io_params={})
+                                  topo, cb, obstacles=[sphere], io_params={})
+#monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+#                                 topo, cb, io_params={})
 
 # === Setup for all declared operators ===
 for ope in op.values():
@@ -205,21 +271,18 @@ allops = dict(op.items() + distr.items() + monitors.items())
 # - compute vorticity with curl
 def initFields_mode1():
     # The velocity is initialized on the fftw topology
-    # penalized and distributed to curl topo
+    # penalized and distributed on curl topo
     velo.initialize(topo=topofft)
     op['penalization'].apply(simu)
     distr['fft2curl'].apply(simu)
     # Vorticity is initialized after penalization of velocity
     op['curl'].apply(simu)
-    # Distribute vorticity and velocity on topostr
-    distr['curl2str'].apply(simu)
-    # Warning : distribute operators do not update ghost points. This is done
-    # by computational operators. So a synchronize call may be required
-    # to monitor proper values.
-    #op['stretching'].discreteOperator._synchronize(op['stretching'].discreteOperator.velocity.data + op['stretching'].discreteOperator.vorticity.data)
-    op['stretching'].updateGhosts()
+    op['advection'].apply(simu)
+    # Distribute vorticity on topofft
+    distr['adv2str'].apply(simu)
+#    distr['curl2adv'].apply(simu)
     # From this point both velocity and vorticity are initialized
-    # on topostr and topocurl. velocity is also init. on topofft.
+    # on topocurl and topoadvec. velocity is also init. on topofft.
 
 ## Call initialization
 initFields_mode1()
@@ -227,73 +290,70 @@ initFields_mode1()
 ##### Check initialize process results #####
 norm_vel_fft = velo.norm(topofft)
 norm_vel_curl = velo.norm(topocurl)
-norm_vel_str = velo.norm(topostr)
-
+norm_vort_fft = vorti.norm(topofft)
 norm_vort_curl = vorti.norm(topocurl)
-norm_vort_str = vorti.norm(topostr)
-
-print norm_vort_curl, norm_vort_str
-print norm_vel_curl, norm_vel_fft
-assert np.allclose(norm_vort_curl, norm_vort_str)
-assert np.allclose(norm_vel_curl, norm_vel_fft)
-assert np.allclose(norm_vel_curl, norm_vel_str)
-
-wref = np.asarray([0., 4354.98713193, 612.8061093], dtype=np.float64)
-vref = np.asarray([1418.04337028, 0., 0.], dtype=np.float64)
-assert np.allclose(norm_vel_curl, vref)
-assert np.allclose(norm_vort_str, wref)
+
+#print norm_vort_curl, norm_vort_fft
+#print norm_vel_curl, norm_vel_fft
+#assert np.allclose(norm_vort_curl, norm_vort_fft)
+#assert np.allclose(norm_vel_curl, norm_vel_fft)
+
 # Print initial state
-for mon in monitors.values():
-    mon.apply(simu)
+#for mon in monitors.values():
+#    mon.apply(simu)
 
 ## Note Franck : tests init ok, same values on both topologies, for v and w.
 # === After init ===
 #
-# v,w init on topostr
+# v init on topocurl (which may be equal to topofft)
 #
 # === Definition of the 'one-step' sequence of operators ===
 #
-#  A - w = stretching(v,w), topostr
-#  B
-#   1 - w = diffusion(w), topofft
-#   2 - w = reprojection(w), topofft
-#   3 - v = poisson(w), topofft
-#   4 - v = correction(v, w), topofft
-#   5 - v = penalization(v), topofft
-#  C - w = curl(v), topocurl
-#  D - w = advection(v,w), topoadvec
+#  1 - w = curl(v), topocurl
+#  2 - w = advection(v,w), topo-advec
+#  3 - w = stretching(v,w), topostr
+#  4 - w = diffusion(w), topofft
+#  5 - v = poisson(w), topofft
+#  6 - v = correction(v, w), topofft
+#  7 - v = penalization(v), topofft
 #
-
-#  Required bridges and in/out var:
-#  A, IN : v, w on topostr
-#     OUT : w on topostr
-#  A - B bridge topostr-topofft for w
-#  B, IN : w
-#     OUT : w, v on topofft
-#  B - C : bridge topofft - topocurl for v
-#  C, IN : v on topocurl
-#     OUT : w on topocurl
-#  C - D : bridge topocurl - topoadvec for v and w
-#  D, IN : v,w on topoadvec
-#     OUT : w on topoadvec
-#  D - A : bridge topoadvec - topostr for v and w
-
-
-# Monitors :
-#  - energy(v,w), works for any topo
-#  - reprojection(w), needs ghost points
-#  - printer(v,w), works for any topo
-#  - forces(v,w), needs ghost points
-
-fullseq = ['stretching',
-           'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
+#  Required bridges :
+# 1 --> 2 : (v,w) on topo-advec
+# 2 --> 3 : (v,w) on topostr
+# 3 --> 4 : w on topofft
+# 7 --> 1 : v on topocurl
+
+#fullseq = ['advection',  # in: v,w out: w, on topoadvec
+#            'adv2str', 'stretching',  # in: v,w out: w on topostr
+#            # in: w, out: v, w on topofft
+#            'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
+#            'penalization',
+#            'fft2curl', 'curl', 'forces',
+#            'curl2adv', 'energy', 'profiles', 'printerAdvec',
+#            ]
+
+fullseq = ['stretching', 
+           'str2diff', 'diffusion', 'poisson', 'correction',
            'penalization',
-           'fft2curl', 'curl',
-           'curl2adv',
+           'fft2curl', 'curl', 
+           'curl2fft', 'poisson', 'correction', 
            'advection',
-           'energy', 'printerAdvec',
+#           'printerSliceXY', 'printerSliceXZ', 'printerSubBox',
+           'energy', 
            'adv2str', 'forces']
 
+#fullseq = ['advection', 
+#           'diffusion',
+#           'poisson', 
+#           'correction',
+#           'penalization',
+#           'fft2curl',
+#           'curl', 
+#           'curl2fft', 
+#           'poisson',
+#           'correction',
+#           'printerSlice']
+
 
 def run(sequence):
     for name in sequence:
@@ -302,12 +362,12 @@ def run(sequence):
         allops[name].apply(simu)
 
 seq = fullseq
-seq = []
+
 simu.initialize()
 
 while not simu.isOver:
-    #if topocurl.rank == 0:
-    #    simu.printState()
+    if topocurl.rank == 0:
+        simu.printState()
     run(seq)
     simu.advance()
 
@@ -343,21 +403,3 @@ for ope in distr.values():
     ope.finalize()
 for monit in monitors.values():
     monit.finalize()
-
-from abc import ABCMeta, abstractmethod
-
-
-class A(object):
-
-    def toto(self):
-        print "toto A "
-
-    def titi(self):
-        self.toto()
-
-
-class B(A):
-
-    def toto(self):
-        print "toto B"
-