From 0d9fdb9ae65451ab6c985ae47634087efd1f5461 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Fri, 15 Mar 2013 10:29:21 +0000
Subject: [PATCH] 2D poisson solver included in Parmes
---
HySoP/hysop/f2py/fftw2py.f90 | 3 +
HySoP/hysop/operator/diffusion_fft.py | 47 ++++----
HySoP/hysop/operator/poisson_fft.py | 37 +++---
.../operator/tests/test_diff_poisson_3D.py | 107 ++++++++++++++++++
HySoP/hysop/operator/tests/test_poisson_2D.py | 76 +++++++++++++
5 files changed, 238 insertions(+), 32 deletions(-)
create mode 100755 HySoP/hysop/operator/tests/test_diff_poisson_3D.py
create mode 100755 HySoP/hysop/operator/tests/test_poisson_2D.py
diff --git a/HySoP/hysop/f2py/fftw2py.f90 b/HySoP/hysop/f2py/fftw2py.f90
index 8d63c42df..91ed3efa1 100755
--- a/HySoP/hysop/f2py/fftw2py.f90
+++ b/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
diff --git a/HySoP/hysop/operator/diffusion_fft.py b/HySoP/hysop/operator/diffusion_fft.py
index b13a06f65..03cbbec63 100644
--- a/HySoP/hysop/operator/diffusion_fft.py
+++ b/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]
diff --git a/HySoP/hysop/operator/poisson_fft.py b/HySoP/hysop/operator/poisson_fft.py
index d3d85b8ac..fdbf2b62f 100644
--- a/HySoP/hysop/operator/poisson_fft.py
+++ b/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 ==============
diff --git a/HySoP/hysop/operator/tests/test_diff_poisson_3D.py b/HySoP/hysop/operator/tests/test_diff_poisson_3D.py
new file mode 100755
index 000000000..4e9698a3c
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_diff_poisson_3D.py
@@ -0,0 +1,107 @@
+# -*- 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()
diff --git a/HySoP/hysop/operator/tests/test_poisson_2D.py b/HySoP/hysop/operator/tests/test_poisson_2D.py
new file mode 100755
index 000000000..ab0c362c9
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_poisson_2D.py
@@ -0,0 +1,76 @@
+# -*- 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()
--
GitLab