testCurl pour Chloe

Examples/testCurl.py

+import numpy as np
+from parmepy.domain.obstacle.controlBox import ControlBox
+import parmepy as pp
+import parmepy.mpi as mpi
+from parmepy.domain.obstacle.planes import SubSpace, SubPlane
+from parmepy.operator.monitors.printer import Printer
+from parmepy.problem.simulation import Simulation
+from parmepy.constants import HDF5
+import numpy as np
+## pi constant
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+nb = 129
+
+Lx = Ly = Lz = 2 * pi
+dom = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[-pi, -pi, -pi])
+resol3D = [nb, nb, nb]
+
+from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
+    Remesh, Support, Splitting, dtAdvecCrit, SpaceDiscretisation, GhostUpdate
+from parmepy.methods import FD_C_4
+
+## 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
+
+velo = pp.Field(domain=dom, formula=computeVel, name='Velocity', isVector=True)
+velo2 = pp.Field(domain=dom, formula=computeVel, name='Veloty2', isVector=True)
+vortifft = pp.Field(domain=dom, name='vortifft', isVector=True)
+vortiFD = pp.Field(domain=dom, isVector=True)
+vortiref = pp.Field(domain=dom, name='vortiref', formula=computeVort, isVector=True)
+vortiref2 = pp.Field(domain=dom, name='vortiref2', formula=computeVort, isVector=True)
+
+from parmepy.operator.differential import Curl
+from parmepy.f2py import fftw2py
+topo3 = mpi.topology.Cartesian(dom, 3, resol3D, ghosts=[2, 2, 2])
+
+curlfft = Curl(velo, vortifft, resolutions={velo: resol3D, vortifft: resol3D},
+               method={SpaceDiscretisation: fftw2py, GhostUpdate: False})
+
+curlFD = Curl(velo2, vortiFD, resolutions={velo2: resol3D, vortiFD: resol3D},
+              method={SpaceDiscretisation: FD_C_4, GhostUpdate: True},
+              topo=topo3)
+
+curlfft.discretize()
+curlFD.discretize()
+
+topofft = curlfft.discreteFields[curlfft.outvar].topology
+vortiref.discretize(topofft)
+vortiref2.discretize(topo3)
+velo.initialize(topofft)
+velo2.initialize(topo3)
+vortiref.initialize(topofft)
+vortiref2.initialize(topo3)
+curlfft.setUp()
+curlFD.setUp()
+
+wdfft = vortifft.discreteFields[topofft].data
+wFD = vortiFD.discreteFields[topo3].data
+wREF1 = vortiref.discreteFields[topofft].data
+wREF2 = vortiref2.discreteFields[topo3].data
+
+vd1 = velo.discreteFields[topofft].data
+vd2 = velo2.discreteFields[topo3].data
+
+from parmepy.constants import VTK
+printer = Printer(variables=[vortiref, vortifft],
+                  topo=topofft,
+                  frequency=1,
+                  prefix='./testCurl',
+                  formattype=VTK)
+printer.setUp()
+printer2 = Printer(variables=[vortiFD, vortiref2],
+                   topo=topo3,
+                   frequency=1,
+                   prefix='./testCurl2',
+                   formattype=VTK)
+
+printer2.setUp()
+simulation = Simulation()
+
+curlfft.apply(simulation)
+curlFD.apply(simulation)
+
+printer.apply(simulation)
+printer2.apply(simulation)
+
+ind = topo3.mesh.iCompute
+indfft = topofft.mesh.iCompute
+
+
+for i in xrange(3):
+    print vortifft.norm()
+    print vortiFD.norm()
+    print vortiref.norm()
+    print vortiref2.norm()
+    ## print np.allclose(wdfft[i][...], wFD[i][ind])
+    ## print np.allclose(wdfft[i][...], wREF1[i])
+    ## print np.allclose(wFD[i][ind], wREF2[i][ind])
+    ## print np.allclose(wREF1[i][indfft], wREF2[i][ind])
+    ## print np.sum(wdfft[i][...] - wFD[i][ind])
+    ## print 'max :', np.abs((wdfft[i][indfft] - wFD[i][ind])).max()
+    ## print 'max fft/ref1 :', np.abs((wdfft[i] - wREF1[i])).max()
+    ## print 'max fft/ref11 :', np.abs((wdfft[i][indfft] - wREF1[i][indfft])).max()
+    ## print 'max fft/ref11 :', np.abs((wdfft[i][...] - wREF1[i][...])).max()
+    ## print 'max fft/ref2 :', np.abs((wdfft[i][indfft] - wREF2[i][ind])).max()
+    ## print 'max FD/ref:', np.abs((wFD[i][ind] - wREF2[i][ind])).max()
+    ## print 'max ref/ref:', np.abs((wREF1[i][indfft] - wREF2[i][ind])).max()
+