NSdebug update

Examples/NSDebug.py

@@ -24,7 +24,7 @@ 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, ProfilesOutput
+    Reprojection_criterion
 from parmepy.problem.simulation import Simulation
 from parmepy.constants import VTK, HDF5
 from parmepy.domain.obstacle.planes import PlaneBoundaries, SubPlane
@@ -49,9 +49,11 @@ 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)
 
+nbElem = [129, 65, 65]
+
+
 ## Global resolutionsimu
 #nbElem = [1025, 513, 513]
-nbElem = [129, 65, 65]
 
 
 # Initial upstream flow velocity
@@ -101,8 +103,8 @@ vorti = Field(domain=box, formula=computeVort,
 # - 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)
+topostr = Cartesian(box, box.dimension, nbElem,
+                    ghosts=ghosts)
 
 ## === Obstacles (sphere + up and down plates) ===
 #sphere_pos = (boxlength - boxorigin) * 0.5
@@ -117,7 +119,7 @@ op['advection'] = Advection(velo, vorti,
 
 op['stretching'] = Stretching(velo, vorti,
                               resolutions={velo: nbElem, vorti: nbElem},
-                              topo=topo)
+                              topo=topostr)
 
 op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
 
@@ -158,18 +160,8 @@ 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'])
@@ -193,11 +185,11 @@ monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
                                   formattype=HDF5, prefix='curl_io_vorti')
 
 monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
-                                   frequency=1, formattype=HDF5, prefix='advec_io',
+                                   formattype=HDF5, prefix='advec_io',
                                    xmfalways=True)
 
 thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir], 
-                        lengths=[10.24, 5.12, 4.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,
@@ -221,9 +213,9 @@ monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
                                       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})
+#io_default = {}
+#monitors['profiles'] = ProfilesOutput(velo, vorti, obstacles=[sphere],
+#                                 topo=topofft, io_params={"frequency":1000})
 
 # Setup for reprojection
 reprojection_constant = 0.04
@@ -248,7 +240,7 @@ coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
 #print coeffForce
 # Monitor for forces
 monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-                                  topo, cb, obstacles=[sphere], io_params={})
+                                  topostr, cb, obstacles=[sphere], io_params={})
 #monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
 #                                 topo, cb, io_params={})
 
@@ -277,10 +269,10 @@ def initFields_mode1():
     distr['fft2curl'].apply(simu)
     # Vorticity is initialized after penalization of velocity
     op['curl'].apply(simu)
+    distr['curl2adv'].apply(simu)
     op['advection'].apply(simu)
-    # Distribute vorticity on topofft
+    # Distribute vorticity on topostr
     distr['adv2str'].apply(simu)
-#    distr['curl2adv'].apply(simu)
     # From this point both velocity and vorticity are initialized
     # on topocurl and topoadvec. velocity is also init. on topofft.
 
@@ -293,14 +285,12 @@ norm_vel_curl = velo.norm(topocurl)
 norm_vort_fft = vorti.norm(topofft)
 norm_vort_curl = vorti.norm(topocurl)
 
-#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)
+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 ===
@@ -332,27 +322,14 @@ norm_vort_curl = vorti.norm(topocurl)
 #            'curl2adv', 'energy', 'profiles', 'printerAdvec',
 #            ]
 
-fullseq = ['stretching', 
+fullseq = ['stretching',
            'str2diff', 'diffusion', 'poisson', 'correction',
            'penalization',
-           'fft2curl', 'curl', 
-           'curl2fft', 'poisson', 'correction', 
+           'fft2curl', 'curl',
+           'curl2fft', 'poisson', 'correction',
            'advection',
-#           'printerSliceXY', 'printerSliceXZ', 'printerSubBox',
-           'energy', 
-           'adv2str', 'forces']
-
-#fullseq = ['advection', 
-#           'diffusion',
-#           'poisson', 
-#           'correction',
-#           'penalization',
-#           'fft2curl',
-#           'curl', 
-#           'curl2fft', 
-#           'poisson',
-#           'correction',
-#           'printerSlice']
+           # 'printerSliceXY', 'printerSliceXZ', 'printerSubBox',
+           'energy', 'adv2str', 'forces']
 
 
 def run(sequence):
@@ -374,23 +351,6 @@ while not simu.isOver:
 
 ## print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
 
-
-## print 'aaaa', velo.norm(topofft)
-## op['correction'].apply(simu)
-## print 'iiii', velo.norm(topofft)
-
-## sref = op['correction'].discreteOperator.surfRef
-
-## flowrate = velo.integrateOnSurface(sref, topofft)
-## print flowrate
-## assert (flowrate - ref_rate) < 1e-10
-
-## flowratey = velo.integrateOnSurface(sref, topofft, component=1)
-## print flowratey
-## flowratez = velo.integrateOnSurface(sref, topofft, component=2)
-## print flowratez
-
-
 # === Finalize for all declared operators ===
 
 # Note FP : bug in fft finalize. To be checked ...

Examples/testVisu.py

@@ -43,7 +43,7 @@ def func3D(res, x, y, z, t):
     return res
 
 
-scal3D = Field(domain=dom, formula=func3D, name='Toto')
+scal3D = Field(domain=dom, formula=func3D, name='Scal')
 scalRef = Field(domain=dom, formula=func3D, name='Ref')
 resolution3D = np.ones(3) * nb
 ghosts = np.zeros(3) + 2
@@ -127,3 +127,8 @@ for pref in filelist[1:]:
 ##     #print scal3D.norm(topo)
 ##     #assert np.allclose(sc3D.data[0][topo.mesh.iCompute], scRef.data[0][topo.mesh.iCompute])
 
+
+velo = Field(domain=dom, name='Velocity', isVector=True)
+vorti = Field(domain=dom, name='vty', isVector=True)
+reader = Reader(variables=[velo, vorti], topo=topo, prefix='curl_io_00000')
+