tfftw_utilities.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tfftw_utilities.hh (2867B)
---
1 /* Copyright (C) 2018 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 #include "pism/util/petscwrappers/Vec.hh"
28
29 namespace pism {
30
31 /*!
32 * Template class for accessing the central part of an extended grid, i.e. PISM's grid
33 * surrounded by "padding" necessary to reduce artifacts coming from interpreting model
34 * inputs as periodic in space.
35 *
36 * Allows accessing a 1D array using 2D indexing.
37 */
38 class FFTWArray {
39 public:
40 // Access the central part of an array "a" of size My*Mx, using offsets i_offset and
41 // j_offset specifying the corner of the part to be accessed.
42 FFTWArray(fftw_complex *a, int Mx, int My, int i_offset, int j_offset);
43
44 // Access the whole array "a" of size My*My.
45 FFTWArray(fftw_complex *a, int Mx, int My);
46
47 inline std::complex<double> &operator()(int i, int j) {
48 return m_array[(m_j_offset + j) + m_My * (m_i_offset + i)];
49 }
50
51 private:
52 const int m_My, m_i_offset, m_j_offset;
53 std::complex<double> *m_array;
54 };
55
56 std::vector<double> fftfreq(int M, double normalization);
57
58 //! Fill `input` with zeros.
59 void clear_fftw_array(fftw_complex *input, int Nx, int Ny);
60
61 //! Copy `source` to `destination`.
62 void copy_fftw_array(fftw_complex *source, fftw_complex *destination, int Nx, int Ny);
63
64 //! Set the real part of output to input. Input has the size of My*Mx, embedded in the
65 //! bigger (output) grid of size Ny*Nx. Offsets i0 and j0 specify the location of the
66 //! subset to set.
67 /*!
68 * Sets the imaginary part to zero.
69 */
70 void set_real_part(Vec input,
71 double normalization,
72 int Mx, int My,
73 int Nx, int Ny,
74 int i0, int j0,
75 fftw_complex *output);
76
77 //! \brief Get the real part of input and put it in output.
78 /*!
79 * See set_real_part for details.
80 */
81 void get_real_part(fftw_complex *input,
82 double normalization,
83 int Mx, int My,
84 int Nx, int Ny,
85 int i0, int j0,
86 Vec output);
87
88 } // end of namespace pism