Skip to content
Snippets Groups Projects
Commit 17379d1a authored by Jean-Baptiste Keck's avatar Jean-Baptiste Keck
Browse files

np.fft -> pyfftw, removed debug dumps

parent 4f9088cf
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
hysop/backend/device/opencl/operator/poisson_rotational.py
+ 4
− 26
View file @ 17379d1a
...@@ -130,7 +130,6 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic): ...@@ -130,7 +130,6 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
self.kernel_buffers[0][i].bind_memory_object(K1d[start:end]) self.kernel_buffers[0][i].bind_memory_object(K1d[start:end])
self.kernel_buffers[1][i].bind_memory_object(K2d[start:end]) self.kernel_buffers[1][i].bind_memory_object(K2d[start:end])
start = end start = end
self.K1d = K1d self.K1d = K1d
self.K2d = K2d self.K2d = K2d
self.Ksizes = Ksizes self.Ksizes = Ksizes
...@@ -225,26 +224,14 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic): ...@@ -225,26 +224,14 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
# /!\ clFFT use the destination buffer as a scratch # /!\ clFFT use the destination buffer as a scratch
# so we reverse the order of forward transforms. # so we reverse the order of forward transforms.
it = simulation.current_iteration #it = simulation.current_iteration
t = simulation.t() #t = simulation.t()
import numpy as np #debug_dumper(it, t, 'W', self.dW)
np.random.seed(42)
for i in xrange(3):
Wi = self.W_buffers[i]
wi = self.dW[i].handle
ri = np.full(shape=wi.shape, dtype=wi.dtype, fill_value=np.nan)
ri[self.dW.compute_slices] = np.random.rand(*Wi.shape).astype(Wi.dtype) + 1.0*i
wi[...] = ri
debug_dumper(it, t, 'W', self.dW)
for fp in self.forward_W_plans: for fp in self.forward_W_plans:
evt, = fp.enqueue(queue=self.queue) evt, = fp.enqueue(queue=self.queue)
debug_dumper(it, t, 'W0_hat', (self.fft_buffers[0].get().handle,))
debug_dumper(it, t, 'W1_hat', (self.fft_buffers[1].get().handle,))
debug_dumper(it, t, 'W2_hat', (self.fft_buffers[2].get().handle,))
# project and recover vorticity if required # project and recover vorticity if required
if self._do_project(simulation): if self._do_project(simulation):
assert False assert False
...@@ -256,20 +243,11 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic): ...@@ -256,20 +243,11 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
# recover velocity # recover velocity
evt = self.compute_velocity_kernel(queue=self.queue) evt = self.compute_velocity_kernel(queue=self.queue)
debug_dumper(it, t, 'U0_hat', (self.fft_buffers[0].get().handle,))
debug_dumper(it, t, 'U1_hat', (self.fft_buffers[1].get().handle,))
debug_dumper(it, t, 'U2_hat', (self.fft_buffers[2].get().handle,))
for bp in self.backward_U_plans: for bp in self.backward_U_plans:
evt, = bp.enqueue() evt, = bp.enqueue()
debug_dumper(it, t, 'Upreg', self.dU)
if (self._exchange_U_ghosts is not None): if (self._exchange_U_ghosts is not None):
evt = self._exchange_U_ghosts(queue=self.queue) evt = self._exchange_U_ghosts(queue=self.queue)
debug_dumper(it, t, 'Upostg', self.dU) #debug_dumper(it, t, 'U', self.dU)
import sys
sys.exit(1)
hysop/backend/host/python/operator/poisson_rotational.py
+ 17
− 44
View file @ 17379d1a
import pyfftw import pyfftw, multiprocessing
from hysop.tools.types import check_instance, first_not_None from hysop.tools.types import check_instance, first_not_None
from hysop.tools.decorators import debug from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw from hysop.tools.numpywrappers import npw
...@@ -38,10 +38,12 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator): ...@@ -38,10 +38,12 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator):
iK2[Z] = 0.0 iK2[Z] = 0.0
FFT = pyfftw.interfaces.numpy_fft FFT = pyfftw.interfaces.numpy_fft
nthreads = multiprocessing.cpu_count()
self.K = K1 self.K = K1
self.iK2 = iK2 self.iK2 = iK2
self.FFT = FFT self.FFT = FFT
self.nthreads = nthreads
@op_apply @op_apply
def apply(self, simulation=None, **kwds): def apply(self, simulation=None, **kwds):
...@@ -60,64 +62,46 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator): ...@@ -60,64 +62,46 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator):
def _solve2d(self, **kwds): def _solve2d(self, **kwds):
W_buffers, U_buffers = self.W_buffers, self.U_buffers W_buffers, U_buffers = self.W_buffers, self.U_buffers
iK2, (Ky,Kx) = self.iK2, self.K iK2, (Ky,Kx) = self.iK2, self.K
fft = self.FFT FFT, nthreads = self.FFT, self.nthreads
W0, = W_buffers W0, = W_buffers
U0, U1 = U_buffers U0, U1 = U_buffers
W0_hat = FFT.rfftn(W0) W0_hat = FFT.rfftn(W0, threads=nthreads)
W0_hat[...] *= iK2 W0_hat[...] *= iK2
Ui_hat = npw.empty_like(W0_hat) Ui_hat = npw.empty_like(W0_hat)
Ui_hat[...] = W0_hat Ui_hat[...] = W0_hat
Ui_hat[...] *= -Ky Ui_hat[...] *= -Ky
U0[...] = FFT.irfftn(Ui_hat, U0.shape) U0[...] = FFT.irfftn(Ui_hat, s=U0.shape, threads=nthreads)
Ui_hat[...] = W0_hat Ui_hat[...] = W0_hat
Ui_hat[...] *= +Kx Ui_hat[...] *= +Kx
U1[...] = FFT.irfftn(Ui_hat, U1.shape) U1[...] = FFT.irfftn(Ui_hat, s=U1.shape, threads=nthreads)
self.dU.exchange_ghosts() self.dU.exchange_ghosts()
def _solve3d(self, simulation, project=False, debug_dumper=None, **kwds): def _solve3d(self, simulation, project=False, debug_dumper=None, **kwds):
it = simulation.current_iteration
t = simulation.t()
W_buffers, U_buffers = self.W_buffers, self.U_buffers W_buffers, U_buffers = self.W_buffers, self.U_buffers
iK2, (Kz,Ky,Kx) = self.iK2, self.K iK2, (Kz,Ky,Kx) = self.iK2, self.K
FFT = self.FFT FFT, nthreads = self.FFT, self.nthreads
W0,W1,W2 = W_buffers W0,W1,W2 = W_buffers
U0,U1,U2 = U_buffers U0,U1,U2 = U_buffers
import numpy as np W0_hat = FFT.rfftn(W0, threads=nthreads).astype(self.ctype)
np.random.seed(42) W1_hat = FFT.rfftn(W1, threads=nthreads).astype(self.ctype)
for i in xrange(3): W2_hat = FFT.rfftn(W2, threads=nthreads).astype(self.ctype)
Wi = self.W_buffers[i]
wi = self.dW[i].handle
ri = np.full(shape=wi.shape, dtype=wi.dtype, fill_value=np.nan)
ri[self.dW.compute_slices] = np.random.rand(*Wi.shape).astype(Wi.dtype) + 1.0*i
wi[...] = ri
debug_dumper(it, t, 'W', self.dW)
W0_hat = FFT.rfftn(W0).astype(self.ctype)
W1_hat = FFT.rfftn(W1).astype(self.ctype)
W2_hat = FFT.rfftn(W2).astype(self.ctype)
debug_dumper(it, t, 'W0_hat', W0_hat )
debug_dumper(it, t, 'W1_hat', W1_hat )
debug_dumper(it, t, 'W2_hat', W2_hat )
if project: if project:
assert False assert False
divW = iK2*(Kx*W0_hat + Ky*W1_hat + Kz*W2_hat) divW = iK2*(Kx*W0_hat + Ky*W1_hat + Kz*W2_hat)
W0_hat[...] -= Kx*divW W0_hat[...] -= Kx*divW
W1_hat[...] -= Ky*divW W1_hat[...] -= Ky*divW
W2_hat[...] -= Kz*divW W2_hat[...] -= Kz*divW
W0[...] = FFT.irfftn(W0_hat, W0.shape) W0[...] = FFT.irfftn(W0_hat, s=W0.shape, threads=nthreads)
W1[...] = FFT.irfftn(W1_hat, W1.shape) W1[...] = FFT.irfftn(W1_hat, s=W1.shape, threads=nthreads)
W2[...] = FFT.irfftn(W2_hat, W2.shape) W2[...] = FFT.irfftn(W2_hat, s=W2.shape, threads=nthreads)
self.dW.exchange_ghosts() self.dW.exchange_ghosts()
W0_hat[...] *= iK2 W0_hat[...] *= iK2
...@@ -127,23 +111,12 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator): ...@@ -127,23 +111,12 @@ class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator):
Ui_hat = npw.empty_like(W0_hat) Ui_hat = npw.empty_like(W0_hat)
Ui_hat[...] = (-Ky*W2_hat + Kz*W1_hat) Ui_hat[...] = (-Ky*W2_hat + Kz*W1_hat)
debug_dumper(it, t, 'U0_hat', Ui_hat ) U0[...] = FFT.irfftn(Ui_hat, s=U0.shape, threads=nthreads)
U0[...] = FFT.irfftn(Ui_hat, U0.shape)
Ui_hat[...] = (-Kz*W0_hat + Kx*W2_hat) Ui_hat[...] = (-Kz*W0_hat + Kx*W2_hat)
debug_dumper(it, t, 'U1_hat', Ui_hat ) U1[...] = FFT.irfftn(Ui_hat, s=U1.shape, threads=nthreads)
U1[...] = FFT.irfftn(Ui_hat, U1.shape)
Ui_hat[...] = (-Kx*W1_hat + Ky*W0_hat) Ui_hat[...] = (-Kx*W1_hat + Ky*W0_hat)
debug_dumper(it, t, 'U2_hat', Ui_hat ) U2[...] = FFT.irfftn(Ui_hat, s=U2.shape, threads=nthreads)
U2[...] = FFT.irfftn(Ui_hat, U2.shape)
debug_dumper(it, t, 'Upreg', self.dU)
self.dU.exchange_ghosts() self.dU.exchange_ghosts()
debug_dumper(it, t, 'Upostg', self.dU)
import sys
sys.exit(1)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment