From 430cf288629cc60b07606c13f211ae70504a1445 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Thu, 22 Jan 2015 18:54:53 +0100
Subject: [PATCH] poisson/pressure (not really tested ...)
---
HySoP/hysop/default_methods.py | 3 +-
HySoP/hysop/f2py/fftw2py.f90 | 46 ++++--
.../hysop/numerics/differential_operations.py | 59 ++++++++
HySoP/hysop/operator/differential.py | 24 ++-
HySoP/hysop/operator/discrete/differential.py | 26 +++-
HySoP/hysop/operator/discrete/poisson_fft.py | 139 +++++++++++-------
HySoP/hysop/operator/poisson.py | 49 +++---
HySoP/src/fftw/fft3d.f90 | 73 +++++++--
8 files changed, 323 insertions(+), 96 deletions(-)
diff --git a/HySoP/hysop/default_methods.py b/HySoP/hysop/default_methods.py
index 9e17fbf0e..a18e8c6df 100644
--- a/HySoP/hysop/default_methods.py
+++ b/HySoP/hysop/default_methods.py
@@ -28,7 +28,8 @@ BAROCLINIC = {SpaceDiscretisation: FD_C_4}
DIFFUSION = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
-POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True,
+ Formulation: 'velocity'}
STRETCHING = {TimeIntegrator: RK3, Formulation: "Conservative",
SpaceDiscretisation: FD_C_4}
diff --git a/HySoP/hysop/f2py/fftw2py.f90 b/HySoP/hysop/f2py/fftw2py.f90
index cfa244b94..721bcda06 100755
--- a/HySoP/hysop/f2py/fftw2py.f90
+++ b/HySoP/hysop/f2py/fftw2py.f90
@@ -64,6 +64,38 @@ contains
end if
end subroutine init_fftw_solver
+
+ !> Initialisation of fftw context : create plans and memory buffers
+ !! @param[in] resolution global resolution of the discrete domain
+ !! @param[in] lengths width of each side of the domain
+ !! @param[in] comm MPI communicator
+ !! @param[out] datashape local dimension of the input/output field
+ !! @param[out] offset absolute index of the first component of the local field
+ subroutine init_fftw_solver_scalar(resolution,lengths,comm,datashape,offset,dim,fftw_type_real)
+
+ integer, intent(in) :: dim
+ integer, dimension(dim),intent(in) :: resolution
+ real(pk),dimension(dim), intent(in) :: lengths
+ integer(ik), dimension(dim), intent(out) :: datashape
+ integer(ik), dimension(dim), intent(out) :: offset
+ integer, intent(in) :: comm
+ logical, optional :: fftw_type_real
+ !f2py optional :: dim=len(resolution)
+ !f2py intent(hide) dim
+ !f2py logical optional, intent(in) :: fftw_type_real = 1
+
+ integer :: ierr
+
+ ! Duplicate comm into client_data::main_comm (used later in fft2d and fft3d)
+ call mpi_comm_dup(comm, main_comm, ierr)
+
+ print*, "Init fftw/poisson solver for a 3d problem"
+ call init_r2c_scalar_3d(resolution,lengths)
+
+ call getParamatersTopologyFFTW3d(datashape,offset)
+
+ end subroutine init_fftw_solver_scalar
+
!> Free memory allocated for fftw-related objects (plans and buffers)
subroutine clean_fftw_solver(dim)
@@ -225,18 +257,14 @@ contains
!> Compute the pressure from the velocity field, solving a Poisson equation.
!! \f{eqnarray*} \Delta p ' &=& rhs \f}
!! with rhs depending on the first derivatives of the velocity field
- !! @param[in] rhs
- !! @param[out] pressure
- !! in the following, rhs is the right hand side 3D-scalar field, and pressure is
- !! a 3D scalar field.
- subroutine pressure_3d(rhs, pressure, ghosts)
- real(pk),dimension(:,:,:), intent(in) :: rhs
+ !! @param[in, out] pressure
+ !! in the following, pressure is used as inout parameter. It must contains the rhs of poisson equation.
+ subroutine pressure_3d(pressure, ghosts)
integer, dimension(3), intent(in) :: ghosts
real(pk),dimension(:,:,:),intent(inout):: pressure
+ !f2py intent(in,out) :: pressure
- !f2py intent(in,out) :: omega_x,omega_y,omega_z
-
- call r2c_scalar_3d(rhs, ghosts)
+ call r2c_scalar_3d(pressure, ghosts)
call filter_pressure_3d()
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index ab4767efd..5dab84270 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -481,3 +481,62 @@ class GradVxW(DifferentialOperation):
npw.add(result[comp], self._work[0], result[comp])
diagnostics[1] = max(diagnostics[1], np.max(self._work[1]))
return result, diagnostics
+
+
+class DivAdvection(DifferentialOperation):
+ """
+ Computes \f$ - \nabla .(V . \nabla V) \f$, with
+ \f$V\f$ a vector field.
+ """
+ def __init__(self, topo, method=FD_C_4):
+ DifferentialOperation.__init__(self, topo)
+ if method is FD_C_4:
+ self.fcall = self.FDCentral
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
+ assert (topo.ghosts() >= 2).all(),\
+ 'you need a ghost layer for FD4 scheme.'
+
+ elif method is FD_C_2:
+ self.fcall = self.FDCentral
+ self.fd_scheme = FD_C_2(topo.mesh.space_step)
+ assert (topo.ghosts() >= 1).all(),\
+ 'you need a ghost layer for FD2 scheme.'
+
+ else:
+ raise ValueError("FD scheme Not yet implemented")
+ self._indices = topo.mesh.iCompute
+ self.fd_scheme.computeIndices(self._indices)
+ # dimension of the fields
+ self._dim = topo.domain.dimension
+
+ def __call__(self, var1, result):
+ return self.fcall(var1, result)
+
+ def FDCentral(self, var1, result):
+
+
+ work = npw.zeros_like(rhs)
+
+ j1 = np.arange(self.velocity.nbComponents)
+ rhs[...] = 0.0
+
+ for i in xrange(self.velocity.nbComponents):
+ j2 = np.where(j1 != i)[0]
+ self._fd_scheme.compute(vdata[j2[0]], j2[1], work)
+ self._fd_scheme.compute_and_mult(vdata[j2[1]], j2[0], work)
+ rhs[...] += work[...]
+ work[...] = 0.0
+ rhs[...] *= -1.0
+
+ for i in xrange(self.velocity.nbComponents):
+ j2 = np.where(j1 != i)[0]
+ self._fd_scheme.compute(vdata[j2[0]], j2[0], work)
+ self._fd_scheme.compute_and_mult(vdata[j2[1]], j2[1], work)
+ rhs[...] += work[...]
+ work[...] = 0.0
+
+ rhs[...] *= 2.0
+
+
+
+ return result
diff --git a/HySoP/hysop/operator/differential.py b/HySoP/hysop/operator/differential.py
index b10cb26db..060a6bd0f 100644
--- a/HySoP/hysop/operator/differential.py
+++ b/HySoP/hysop/operator/differential.py
@@ -5,7 +5,8 @@ Differential operators
"""
from hysop.constants import debug
from hysop.operator.computational import Computational
-from hysop.operator.discrete.differential import CurlFFT, CurlFD, GradFD
+from hysop.operator.discrete.differential import CurlFFT, CurlFD,\
+ GradFD, DivAdvectionFD
from hysop.methods_keys import SpaceDiscretisation
from hysop.operator.continuous import opsetup
from hysop.numerics.finite_differences import FD_C_4,\
@@ -130,3 +131,24 @@ class Grad(Differential):
outvar=self.discreteFields[self.outvar],
method=self.method, rwork=rwork)
self._is_uptodate = True
+
+
+class DivAdvection(Differential):
+ """
+ Computes \f$ outVar = -\nabla .(invar . \nabla invar) \f$
+ """
+
+ @debug
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ if self.method[SpaceDiscretisation].mro()[1] is not FiniteDifference:
+ raise ValueError("DivAdvection operator only\
+ implemented with finite differences. Please change\
+ method[SpaceDiscretisation] value.")
+ # Create a discrete operator, according to the chosen method
+ # (only finite differences at the time).
+ self.discrete_op = DivAdvectionFD(
+ invar=self.discreteFields[self.invar],
+ outvar=self.discreteFields[self.outvar],
+ method=self.method, rwork=rwork)
+ self._is_uptodate = True
diff --git a/HySoP/hysop/operator/discrete/differential.py b/HySoP/hysop/operator/discrete/differential.py
index b00d6a074..4b7436052 100644
--- a/HySoP/hysop/operator/discrete/differential.py
+++ b/HySoP/hysop/operator/discrete/differential.py
@@ -6,7 +6,8 @@ Discretisation of the differential operators (curl, grad ...)
"""
from hysop.constants import debug
from hysop.operator.discrete.discrete import DiscreteOperator
-from hysop.numerics.differential_operations import Curl, GradV
+from hysop.numerics.differential_operations import Curl, GradV,\
+ DivAdvection
import hysop.tools.numpywrappers as npw
from abc import ABCMeta, abstractmethod
from hysop.numerics.update_ghosts import UpdateGhosts
@@ -149,3 +150,26 @@ class GradFD(Differential):
def apply(self, simulation=None):
self._synchronize(self.invar.data)
self.outvar.data = self._function(self.invar.data, self.outvar.data)
+
+
+class DivAdvectionFD(Differential):
+ """
+ Compute the grad of a discrete field, using finite differences.
+ """
+
+ def __init__(self, **kwds):
+
+ super(DivAdvectionFD, self).__init__(**kwds)
+ # prepare ghost points synchro for velocity
+ self._synchronize = UpdateGhosts(self.invar.topology,
+ self.invar.nb_components)
+ dim = self.domain.dimension
+ assert dim * self.outvar.nb_components == self.invar.nb_components
+ self._function = DivAdvection(self.invar.topology,
+ self.method[SpaceDiscretisation])
+
+ @debug
+ @profile
+ def apply(self, simulation=None):
+ self._synchronize(self.invar.data)
+ self.outvar.data = self._function(self.invar.data, self.outvar.data)
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index aa031458e..1b1be1875 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -22,42 +22,42 @@ class PoissonFFT(DiscreteOperator):
"""
@debug
- def __init__(self, velocity, vorticity, projection=None,
+ def __init__(self, output_field, input_field, projection=None,
filterSize=None, correction=None, **kwds):
"""
Constructor.
- @param[out] velocity : discretization of the solution field
- @param[in] vorticity : discretization of the RHS (mind the minus rhs!)
+ @param[out] output_field : discretization of the solution field
+ @param[in] input_field : discretization of the RHS (mind the minus rhs!)
@param projection : if None, no projection. Else:
- either the value of the frequency of reprojection, never updated.
- or Reprojection discrete operator. In that case, a criterion
- depending on the vorticity will be computed at each time step, if
+ depending on the input_field will be computed at each time step, if
criterion > threshold, then frequency projection is active.
@param filterSize :
- @param correction : operator used to shift velocity according
+ @param correction : operator used to shift output_field according
to a given input (fixed) flowrate.
See hysop.operator.velocity_correction.
Default = None.
"""
# Base class initialisation
assert 'variables' not in kwds, 'variables parameter is useless.'
- super(PoissonFFT, self).__init__(variables=[velocity, vorticity],
+ super(PoissonFFT, self).__init__(variables=[output_field, input_field],
**kwds)
## Solution field
- self.velocity = velocity
+ self.output_field = output_field
## RHS field
- self.vorticity = vorticity
- ## Solenoidal projection of vorticity ?
+ self.input_field = input_field
+ ## Solenoidal projection of input_field ?
self.projection = projection
## Filter size array = domainLength/(CoarseRes-1)
self.filterSize = filterSize
- # If 2D problem, vorticity must be a scalar
- self.dim = self.velocity.domain.dimension
+ # If 2D problem, input_field must be a scalar
+ self.dim = self.output_field.domain.dimension
if self.dim == 2:
- assert self.vorticity.nb_components == 1
+ assert self.input_field.nb_components == 1
self.correction = correction
- self.input = [self.vorticity]
- self.output = [self.velocity]
+ self.input = [self.input_field]
+ self.output = [self.output_field]
## The function called during apply
self.solve = None
@@ -67,16 +67,22 @@ class PoissonFFT(DiscreteOperator):
self._select_solve()
def _select_solve(self):
+
+ """
+ TODO : add pressure solver selection
+ f(output.nbComponents) + pb type 'pressure-poisson'
+ """
+
## Multiresolution ?
- multires = self.velocity.topology.mesh != self.vorticity.topology.mesh
+ multires = self.output_field.topology.mesh != self.input_field.topology.mesh
# connexion to the required apply function
if self.dim == 2:
self._solve = self._solve2D
elif self.dim == 3:
- # If there is a projection, vorticity is also an output
+ # If there is a projection, input_field is also an output
if self.projection is not None:
- self.output.append(self.vorticity)
+ self.output.append(self.input_field)
if multires:
self._solve = self._solve3D_proj_multires
else:
@@ -95,7 +101,7 @@ class PoissonFFT(DiscreteOperator):
else:
raise AttributeError('Not implemented for 1D problems.')
- # Operator to shift velocity according to an input required flowrate
+ # Operator to shift output_field according to an input required flowrate
if self.correction is not None:
self.solve = self._solve_and_correct
else:
@@ -114,43 +120,43 @@ class PoissonFFT(DiscreteOperator):
"""
Solve 2D poisson problem
"""
- self.velocity.data[0], self.velocity.data[1] =\
- fftw2py.solve_poisson_2d(self.vorticity.data[0],
- self.velocity.data[0],
- self.velocity.data[1])
+ self.output_field.data[0], self.output_field.data[1] =\
+ fftw2py.solve_poisson_2d(self.input_field.data[0],
+ self.output_field.data[0],
+ self.output_field.data[1])
def _project(self):
"""
- apply projection onto vorticity
+ apply projection onto input_field
"""
- ghosts_w = self.vorticity.topology.ghosts()
- self.vorticity.data[0], self.vorticity.data[1], \
- self.vorticity.data[2] = \
- fftw2py.projection_om_3d(self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2], ghosts_w)
+ ghosts_w = self.input_field.topology.ghosts()
+ self.input_field.data[0], self.input_field.data[1], \
+ self.input_field.data[2] = \
+ fftw2py.projection_om_3d(self.input_field.data[0],
+ self.input_field.data[1],
+ self.input_field.data[2], ghosts_w)
def _solve3D_multires(self, simu=None):
"""
3D, multiresolution
"""
- # Projects vorticity values from fine to coarse grid
+ # Projects input_field values from fine to coarse grid
# in frequencies space by nullifying the smallest modes
- vortFilter = npw.copy(self.vorticity.data)
+ vortFilter = npw.copy(self.input_field.data)
vortFilter[0], vortFilter[1], vortFilter[2] = \
fftw2py.multires_om_3d(self.filterSize[0], self.filterSize[1],
- self.filterSize[2], self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2])
-
- # Solves Poisson equation using filter vorticity
- ghosts_v = self.velocity.topology.ghosts()
- ghosts_w = self.vorticity.topology.ghosts()
- self.velocity.data[0], self.velocity.data[1], self.velocity.data[2] = \
+ self.filterSize[2], self.input_field.data[0],
+ self.input_field.data[1],
+ self.input_field.data[2])
+
+ # Solves Poisson equation using filter input_field
+ ghosts_v = self.output_field.topology.ghosts()
+ ghosts_w = self.input_field.topology.ghosts()
+ self.output_field.data[0], self.output_field.data[1], self.output_field.data[2] = \
fftw2py.solve_poisson_3d(vortFilter[0], vortFilter[1],
- vortFilter[2], self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2], ghosts_w, ghosts_v)
+ vortFilter[2], self.output_field.data[0],
+ self.output_field.data[1],
+ self.output_field.data[2], ghosts_w, ghosts_v)
def _solve3D_proj_multires(self, simu):
"""
@@ -172,21 +178,48 @@ class PoissonFFT(DiscreteOperator):
"""
Basic solve
"""
- # Solves Poisson equation using usual vorticity
- ghosts_v = self.velocity.topology.ghosts()
- ghosts_w = self.vorticity.topology.ghosts()
- 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], ghosts_w, ghosts_v)
+ # Solves Poisson equation using usual input_field
+ ghosts_v = self.output_field.topology.ghosts()
+ ghosts_w = self.input_field.topology.ghosts()
+ self.output_field.data[0], self.output_field.data[1], self.output_field.data[2] =\
+ fftw2py.solve_poisson_3d(self.input_field.data[0],
+ self.input_field.data[1],
+ self.input_field.data[2],
+ self.output_field.data[0],
+ self.output_field.data[1],
+ self.output_field.data[2], ghosts_w, ghosts_v)
def _solve_and_correct(self, simu):
self._solve(simu.currentIteration)
self.correction.apply(simu)
+
+ def _solve_3d_scalar_fd(self, simu=None):
+ """solve poisson-pressure like problem
+ input = 3D vector field
+ output = 3D scalar field
+ """
+ # Compute rhs = f(input) inplace
+ # --> output == rhs
+ # Call fftw filter
+ # !!! pressure3d use the same arg for input and output
+ # ---> input_field will be overwritten
+ ghosts = self.output_field.topology.ghosts()
+ self.output_field.data[0] = fftw2py.pressure_3d(
+ self.input_field.data[0], ghosts)
+
+ def _solve_3d_scalar(self, simu=None):
+ """solve poisson-pressure like problem
+ input = 3D vector field
+ output = 3D scalar field
+ """
+ # # Call fftw filter
+ # self._output_field.data[0] = fftw2py.solve_poisson_3d_pressure(
+ # self._input_field.data[0],
+ # self._input_field.data[1],
+ # self._input_field.data[2])
+ pass
+
@debug
@profile
def apply(self, simulation=None):
@@ -197,4 +230,4 @@ class PoissonFFT(DiscreteOperator):
Clean memory (fftw plans and so on)
"""
pass
- #fftw2py.clean_fftw_solver(self.velocity.dimension)
+ #fftw2py.clean_fftw_solver(self.output_field.dimension)
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index 6043dedc6..9a237ae88 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -27,39 +27,46 @@ class Poisson(Computational):
"""
@debug
- def __init__(self, velocity, vorticity, flowrate=None,
+ def __init__(self, output_field, input_field, flowrate=None,
projection=None, **kwds):
"""
Constructor for the Poisson problem.
- @param[out] velocity : solution field
- @param[in] vorticity : rhs field
- @param[in] flowrate : a flow rate value (through input surf,
+ @param[out] output_field : solution field
+ @param[in] input_field : rhs field
+ @param[in] flowrate : a flow rate value (through input_field surf,
normal to xdir) used to compute a correction of the solution field.
- Default = 0 (no correction). See hysop.operator.velocity_correction.
+ Default = 0 (no correction). See hysop.operator.output_field_correction.
@param projection : if None, no projection. Else:
- either the value of the frequency of reprojection, never update.
- or a tuple = (frequency, threshold).
In that case, a criterion
- depending on the vorticity will be computed at each time step, if
+ depending on the input_field will be computed at each time step, if
criterion > threshold, then frequency projection is active.
+
+ Note about method:
+ - SpaceDiscretisation == fftw
+ - Formulation = 'velocity' or 'pressure'
+ velocity : laplacian(phi) = -w and v = nabla X psi, in = vorticity, out = velo
+ pressure : laplacian(p) = -nabla.(u.nabla u, in = velo, out = pressure
"""
# Warning : for fftw all variables must have
# the same resolution.
assert 'variables' not in kwds, 'variables parameter is useless.'
- super(Poisson, self).__init__(variables=[velocity, vorticity], **kwds)
+ super(Poisson, self).__init__(variables=[output_field, input_field],
+ **kwds)
## solution of the problem
- self.velocity = velocity
+ self.output_field = output_field
## -(right-hand side)
- self.vorticity = vorticity
+ self.input_field = input_field
if self.method is None:
self.method = default.POISSON
if self.method[SpaceDiscretisation] is not 'fftw':
raise AttributeError("Method not yet implemented.")
- self.input = [self.vorticity]
- self.output = [self.velocity]
+ self.input = [self.input_field]
+ self.output = [self.output_field]
if flowrate is not None:
self.withCorrection = True
self._flowrate = flowrate
@@ -70,25 +77,25 @@ class Poisson(Computational):
self._config = kwds
if projection is not None:
- self.output.append(self.vorticity)
+ self.output_field.append(self.input_field)
def discretize(self):
# Poisson solver based on fftw
if self.method[SpaceDiscretisation] is 'fftw':
super(Poisson, self)._fftw_discretize()
if self.withCorrection:
- toporef = self.discreteFields[self.velocity].topology
+ toporef = self.discreteFields[self.output_field].topology
if 'discretization' in self._config:
self._config['discretization'] = toporef
self.correction = VelocityCorrection(
- self.velocity, self.vorticity,
+ self.output_field, self.input_field,
req_flowrate=self._flowrate, **self._config)
self.correction.discretize()
if isinstance(self.projection, tuple):
freq = self.projection[0]
threshold = self.projection[1]
- self.projection = Reprojection(self.vorticity,
+ self.projection = Reprojection(self.input_field,
threshold, freq,
**self._config)
self.projection.discretize()
@@ -108,16 +115,16 @@ class Poisson(Computational):
# Activate projection, if required
if isinstance(self.projection, Reprojection):
# Projection frequency is updated at each
- # time step, and depends on the vorticity
+ # time step, and depends on the input_field
self.projection.setup(rwork=rwork)
projection_discr = self.projection.discrete_op
else:
projection_discr = self.projection
- self.discrete_op = PoissonFFT(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- correction=cd,
- rwork=rwork, iwork=iwork,
- projection=projection_discr)
+ self.discrete_op = PoissonFFT(self.discreteFields[self.output_field],
+ self.discreteFields[self.input_field],
+ correction=cd,
+ rwork=rwork, iwork=iwork,
+ projection=projection_discr)
self._is_uptodate = True
diff --git a/HySoP/src/fftw/fft3d.f90 b/HySoP/src/fftw/fft3d.f90
index 89071c81d..6bb0c32a4 100755
--- a/HySoP/src/fftw/fft3d.f90
+++ b/HySoP/src/fftw/fft3d.f90
@@ -22,7 +22,7 @@ module fft3d
private
- public :: init_r2c_3d,init_c2c_3d,r2c_3d,r2c_scalar_3d,c2c_3d,c2r_3d,c2r_scalar_3d,cleanFFTW_3d,&
+ public :: init_r2c_3d,init_c2c_3d,init_r2c_scalar_3d, r2c_3d,r2c_scalar_3d,c2c_3d,c2r_3d,c2r_scalar_3d,cleanFFTW_3d,&
getParamatersTopologyFFTW3d,filter_poisson_3d,filter_curl_diffusion_3d, &
init_r2c_3d_many, r2c_3d_many, c2r_3d_many, filter_diffusion_3d_many,&
filter_poisson_3d_many, filter_diffusion_3d, filter_curl_3d, filter_projection_om_3d,&
@@ -270,6 +270,60 @@ contains
end subroutine init_r2c_3d
+ !> Initialisation of the fftw context for real to complex transforms (forward and backward)
+ !! @param[in] resolution global domain resolution
+ !! @param[in] lengths width of each side of the domain
+ subroutine init_r2c_scalar_3d(resolution,lengths)
+
+ integer, dimension(3), intent(in) :: resolution
+ real(mk),dimension(3), intent(in) :: lengths
+ !! Size of the local memory required for fftw (cbuffer)
+ integer(C_INTPTR_T) :: alloc_local,halfLength
+
+ if(is3DUpToDate) return
+
+ ! init fftw mpi context
+ call fftw_mpi_init()
+
+ if(rank==0) open(unit=21,file=filename,form="formatted")
+ allocate(fft_resolution(3))
+ fft_resolution(:) = resolution(:)-1
+ halfLength = fft_resolution(c_X)/2+1
+
+ ! compute "optimal" size (according to fftw) for local data (warning : dimension reversal)
+ alloc_local = fftw_mpi_local_size_3d_transposed(fft_resolution(c_Z),fft_resolution(c_Y),halfLength,&
+ main_comm,local_resolution(c_Z),local_offset(c_Z),local_resolution(c_Y),local_offset(c_Y));
+
+ ! init c_X part. This is required to compute kx with the same function in 2d and 3d cases.
+ local_offset(c_X) = 0
+ local_resolution(c_X) = halfLength
+
+ ! allocate local buffer (used to save datain/dataout ==> in-place transform!!)
+ cbuffer1 = fftw_alloc_complex(alloc_local)
+
+ ! link rdatain and dataout to cbuffer, setting the right dimensions for each
+ call c_f_pointer(cbuffer1, rdatain1, [2*halfLength,fft_resolution(c_Y),local_resolution(c_Z)])
+ call c_f_pointer(cbuffer1, dataout1, [halfLength, fft_resolution(c_Z), local_resolution(c_Y)])
+
+ rdatain1 = 0.0
+
+ ! create MPI plans for in-place forward/backward DFT (note dimension reversal)
+ plan_forward1 = fftw_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain1, dataout1, &
+ main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
+ plan_backward1 = fftw_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
+ main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
+
+ call computeKx(lengths(c_X))
+ call computeKy(lengths(c_Y))
+ call computeKz(lengths(c_Z))
+
+ normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
+ !! call fft3d_diagnostics(alloc_local)
+ manycase = .false.
+ is3DUpToDate = .true.
+
+ end subroutine init_r2c_scalar_3d
+
!> forward transform - The result is saved in local buffers
!! @param[in] omega_x 3d scalar field, x-component of the input vector field
!! @param[in] omega_y 3d scalar field, y-component of the input vector field
@@ -342,9 +396,9 @@ contains
!> forward transform - The result is saved in a local buffer
!! @param[in] omega 3d scalar field, x-component of the input vector field
!! @param[in] ghosts, number of points in the ghost layer of input field.
- subroutine r2c_scalar_3d(omega, ghosts)
+ subroutine r2c_scalar_3d(scalar, ghosts)
- real(mk),dimension(:,:,:),intent(in) :: omega
+ real(mk),dimension(:,:,:),intent(in) :: scalar
integer, dimension(3), intent(in) :: ghosts
real(mk) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
@@ -359,7 +413,7 @@ contains
jg = j + ghosts(c_Y)
do i = 1, fft_resolution(c_X)
ig = i + ghosts(c_X)
- rdatain1(i,j,k) = omega(ig,jg,kg)
+ rdatain1(i,j,k) = scalar(ig,jg,kg)
end do
end do
end do
@@ -372,10 +426,10 @@ contains
end subroutine r2c_scalar_3d
!> Backward fftw transform
- !! @param[in,out] omega 3d scalar field
- !! @param[in] ghosts, number of points in the ghost layer of in/out omega scalar field.
- subroutine c2r_scalar_3d(omega, ghosts)
- real(mk),dimension(:,:,:),intent(inout) :: omega
+ !! @param[in,out] scalar 3d scalar field
+ !! @param[in] ghosts, number of points in the ghost layer of in/out scalar field.
+ subroutine c2r_scalar_3d(scalar, ghosts)
+ real(mk),dimension(:,:,:),intent(inout) :: scalar
integer, dimension(3), intent(in) :: ghosts
real(mk) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
@@ -389,7 +443,7 @@ contains
jg = j + ghosts(c_Y)
do i = 1, fft_resolution(c_X)
ig = i + ghosts(c_X)
- omega(ig,jg,kg) = rdatain1(i,j,k)*normFFT
+ scalar(ig,jg,kg) = rdatain1(i,j,k)*normFFT
end do
end do
end do
@@ -845,7 +899,6 @@ contains
end do
end subroutine filter_pressure_3d
-
!> Solve Poisson problem in the Fourier space :
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ v &=& \nabla\times\psi \f}
subroutine filter_poisson_3d_many()
--
GitLab