From 20c75ebb5f46c60749a9b2163b1b265a3d3324fa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 24 Jan 2014 13:02:14 +0000
Subject: [PATCH] testCurl pour Chloe
---
Examples/testCurl.py | 122 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 122 insertions(+)
create mode 100644 Examples/testCurl.py
diff --git a/Examples/testCurl.py b/Examples/testCurl.py
new file mode 100644
index 000000000..b7e300c14
--- /dev/null
+++ b/Examples/testCurl.py
@@ -0,0 +1,122 @@
+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()
+
--
GitLab