PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
fftw_utilities.hh
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2020 PISM Authors
2  *
3  * This file is part of PISM.
4  *
5  * PISM is free software; you can redistribute it and/or modify it under the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13  * details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with PISM; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // Utilities for serial models using FFTW and extended computational grids.
21 
22 #include <vector>
23 #include <complex>
24 
25 #include <fftw3.h>
26 
27 namespace pism {
28 
29 namespace petsc {
30 class Vec;
31 } // end of namespace petsc
32 
33 /*!
34  * Template class for accessing the central part of an extended grid, i.e. PISM's grid
35  * surrounded by "padding" necessary to reduce artifacts coming from interpreting model
36  * inputs as periodic in space.
37  *
38  * Allows accessing a 1D array using 2D indexing.
39  */
40 class FFTWArray {
41 public:
42  // Access the central part of an array "a" of size My*Mx, using offsets i_offset and
43  // j_offset specifying the corner of the part to be accessed.
44  FFTWArray(fftw_complex *a, int Mx, int My, int i_offset, int j_offset);
45 
46  // Access the whole array "a" of size My*My.
47  FFTWArray(fftw_complex *a, int Mx, int My);
48 
49  inline std::complex<double> &operator()(int i, int j) {
50  return m_array[(m_j_offset + j) + m_My * (m_i_offset + i)];
51  }
52 
53 private:
54  const int m_My, m_i_offset, m_j_offset;
55  std::complex<double> *m_array;
56 };
57 
58 std::vector<double> fftfreq(int M, double normalization);
59 
60 //! Fill `input` with zeros.
61 void clear_fftw_array(fftw_complex *input, int Nx, int Ny);
62 
63 //! Copy `source` to `destination`.
64 void copy_fftw_array(fftw_complex *source, fftw_complex *destination, int Nx, int Ny);
65 
66 //! Set the real part of output to input. Input has the size of My*Mx, embedded in the
67 //! bigger (output) grid of size Ny*Nx. Offsets i0 and j0 specify the location of the
68 //! subset to set.
69 /*!
70  * Sets the imaginary part to zero.
71  */
72 void set_real_part(petsc::Vec &input,
73  double normalization,
74  int Mx, int My,
75  int Nx, int Ny,
76  int i0, int j0,
77  fftw_complex *output);
78 
79 //! \brief Get the real part of input and put it in output.
80 /*!
81  * See set_real_part for details.
82  */
83 void get_real_part(fftw_complex *input,
84  double normalization,
85  int Mx, int My,
86  int Nx, int Ny,
87  int i0, int j0,
88  petsc::Vec &output);
89 
90 } // end of namespace pism
const int m_j_offset
const int m_i_offset
FFTWArray(fftw_complex *a, int Mx, int My, int i_offset, int j_offset)
std::complex< double > * m_array
std::complex< double > & operator()(int i, int j)
void clear_fftw_array(fftw_complex *input, int Nx, int Ny)
Fill input with zeros.
void copy_fftw_array(fftw_complex *source, fftw_complex *destination, int Nx, int Ny)
Copy source to destination.
void get_real_part(fftw_complex *input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, petsc::Vec &output)
Get the real part of input and put it in output.
void set_real_part(petsc::Vec &input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, fftw_complex *output)
std::vector< double > fftfreq(int M, double normalization)