2D poisson solver included in Parmes

HySoP/hysop/f2py/fftw2py.f90

@@ -71,13 +71,16 @@ contains
   subroutine solve_poisson_2d(omega,velocity_x,velocity_y)
     real(pk),dimension(:,:),intent(in):: omega
     real(pk),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
+    real(pk) :: start
     !f2py intent(in,out) :: velocity_x,velocity_y
+    start = MPI_WTime()
 
     call r2c_2d(omega)
     
     call filter_poisson_2d()
 
     call c2r_2d(velocity_x,velocity_y)
+    print *, "fortran resolution time : ", MPI_WTime() - start
   
   end subroutine solve_poisson_2d
 

HySoP/hysop/operator/diffusion_fft.py

@@ -33,6 +33,7 @@ class DiffusionFFT(DiscreteOperator):
             self.diff.discreteFieldId[self.diff.vorticity]]
         ## Viscosity.
         self.viscosity = viscosity
+        self.domain = self.diff.domain
         ## Boolean : pure vort diffusion or curl(u) + vort diffusion.
         self.with_curl = with_curl
         self.compute_time = 0.
@@ -48,27 +49,35 @@ class DiffusionFFT(DiscreteOperator):
         """
         self.compute_time = time.time()
 
-        if self.with_curl:
-
-            self.vorticity.data[0], \
-            self.vorticity.data[1], \
-            self.vorticity.data[2] = \
-                fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
-                                                self.velocity.data[0],
-                                                self.velocity.data[1],
-                                                self.velocity.data[2],
-                                                self.vorticity.data[0],
-                                                self.vorticity.data[1],
-                                                self.vorticity.data[2])
+        if (self.domain.dimension == 2):
+
+            raise ValueError("Not yet implemented")
+
+        elif (self.domain.dimension == 3):
+
+            if self.with_curl:
+                self.vorticity.data[0], \
+                self.vorticity.data[1], \
+                self.vorticity.data[2] = \
+                    fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+                                                    self.velocity.data[0],
+                                                    self.velocity.data[1],
+                                                    self.velocity.data[2],
+                                                    self.vorticity.data[0],
+                                                    self.vorticity.data[1],
+                                                    self.vorticity.data[2])
+
+            else:
+                self.vorticity.data[0], \
+                self.vorticity.data[1], \
+                self.vorticity.data[2] = \
+                    fftw2py.solve_diffusion_3d(self.viscosity * dt,
+                                               self.vorticity.data[0],
+                                               self.vorticity.data[1],
+                                               self.vorticity.data[2])
 
         else:
-            self.vorticity.data[0], \
-            self.vorticity.data[1], \
-            self.vorticity.data[2] = \
-                fftw2py.solve_diffusion_3d(self.viscosity * dt,
-                                           self.vorticity.data[0],
-                                           self.vorticity.data[1],
-                                           self.vorticity.data[2])
+            raise ValueError("invalid problem dimension")
 
 
 #        ind0a = self.topology.mesh.local_start[0]

HySoP/hysop/operator/poisson_fft.py

@@ -15,7 +15,7 @@ class PoissonFFT(DiscreteOperator):
     """
 
     @debug
-    def __init__(self, poisson, method='', domain=None):
+    def __init__(self, poisson, method=''):
         """
         Constructor.
         @param velocity : descretized velocity to update from vorticity.
@@ -31,7 +31,7 @@ class PoissonFFT(DiscreteOperator):
         self.vorticity = self.poisson.vorticity.discreteField[
             self.poisson.discreteFieldId[self.poisson.vorticity]]
         self.method = method
-        self.domain = domain
+        self.domain = self.poisson.domain
         self.compute_time = 0.
 
     @debug
@@ -51,21 +51,32 @@ class PoissonFFT(DiscreteOperator):
         self.compute_time = time.time()
 
         #=== RECOVER VELOCITY FIELD FROM VORTICITY FIELD ===
-        self.velocity.data[0],\
-        self.velocity.data[1],\
-        self.velocity.data[2] = \
-            fftw2py.solve_poisson_3d(self.vorticity.data[0],
-                                     self.vorticity.data[1],
-                                     self.vorticity.data[2],
-                                     self.velocity.data[0],
-                                     self.velocity.data[1],
-                                     self.velocity.data[2])
+        if (self.domain.dimension == 2):
+            self.velocity.data[0],\
+            self.velocity.data[1] =\
+                fftw2py.solve_poisson_2d(self.vorticity.data,
+                                        self.velocity.data[0],
+                                        self.velocity.data[1])
+
+        elif (self.domain.dimension == 3):
+            self.velocity.data[0],\
+            self.velocity.data[1],\
+            self.velocity.data[2] = \
+                fftw2py.solve_poisson_3d(self.vorticity.data[0],
+                                        self.vorticity.data[1],
+                                        self.vorticity.data[2],
+                                        self.velocity.data[0],
+                                        self.velocity.data[1],
+                                        self.velocity.data[2])
+
+        else:
+            raise ValueError("invalid problem dimension")
 
 
         #============== VELOCITY CORRECTION ================
 
-        surf = (self.domain.length[1]-self.domain.origin[1]) * \
-               (self.domain.length[2]-self.domain.origin[2])
+#        surf = (self.domain.length[1]-self.domain.origin[1]) * \
+#               (self.domain.length[2]-self.domain.origin[2])
 
         #============= CORRECTION VEL X : u_x ==============
 

HySoP/hysop/operator/tests/test_diff_poisson_3D.py

+# -*- coding: utf-8 -*-
+import time
+
+import parmepy as pp
+import numpy as np
+from parmepy.constants import PARMES_REAL, ORDER
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from math import sqrt, pi, exp
+
+
+def computeVel(x, y, z):
+    vx = 0.
+    vy = 0.
+    vz = 0.
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    xc = 1. / 2.
+    yc = 1. / 2.
+    zc = 1. / 4.
+    R = 0.2
+    sigma = R / 2.
+    Gamma = 0.0075
+    dist = 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) * \
+                 exp(-(s2 / sigma ** 2))
+        wx = - wTheta * sinTheta
+        wy = wTheta * cosTheta
+        wz = 0.
+    return wx, wy, wz
+
+def test_Diff_Poisson():
+    # Parameters
+    nb = 65
+    dim = 3
+    boxLength = [1., 1., 1.]
+    boxMin = [0., 0., 0.]
+    nbElem = [nb, nb, nb]
+
+    timeStep = 0.01
+    finalTime = 150 * timeStep
+    outputFilePrefix = './Energies'
+    outputModulo = 10
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## FFT Diffusion operators and FFT Poisson solver
+    diffusion = Diffusion(velo, vorti,
+                          resolutions={velo: nbElem,
+                                       vorti: nbElem},
+                          method='',
+                          viscosity=0.002, 
+                          with_curl=False
+                         )
+
+    poisson = Poisson(velo, vorti,
+                      resolutions={velo: nbElem,
+                                   vorti: nbElem},
+                      method=''
+                     )
+
+    ## Problem
+    pb = NSProblem([diffusion, poisson],
+                    monitors=[Energy_enstrophy(
+                                fields=[velo, vorti],
+                                frequency=outputModulo,
+                                outputPrefix=outputFilePrefix)])
+
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+    test_Diff_Poisson()

HySoP/hysop/operator/tests/test_poisson_2D.py

+# -*- coding: utf-8 -*-
+import time
+
+import parmepy as pp
+import numpy as np
+from parmepy.constants import PARMES_REAL, ORDER
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from math import sqrt, pi, exp
+
+
+def computeVel(x, y):
+    vx = 0.
+    vy = 0.
+    return vx, vy
+
+def computeVort(x, y):
+    w = np.cos(y) * np.cos(x)
+    return w
+
+def test_Diff_Poisson():
+    # Parameters
+    nb = 65
+    dim = 2
+    boxLength = [1., 1.]
+    boxMin = [0., 0.]
+    nbElem = [nb, nb]
+
+    timeStep = 0.01
+    finalTime = 150 * timeStep
+    outputFilePrefix = './Energies'
+    outputModulo = 10
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=False)
+
+    ## FFT Poisson solver
+    poisson = Poisson(velo, vorti,
+                      resolutions={velo: nbElem,
+                                   vorti: nbElem},
+                      method=''
+                     )
+
+    ## Problem
+    pb = NSProblem([poisson])
+
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+    test_Diff_Poisson()