tfftw_utilities.cc - 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.cc (3558B)
---
1 /* Copyright (C) 2018, 2019 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 #include <cstring> // memcpy
21
22 #include "fftw_utilities.hh"
23
24 #include "pism/util/petscwrappers/Vec.hh"
25
26 namespace pism {
27
28 // Access the central part of an array "a" of size My*Mx, using offsets i_offset and
29 // j_offset specifying the corner of the part to be accessed.
30 FFTWArray::FFTWArray(fftw_complex *a, int Mx, int My, int i_offset, int j_offset)
31 : m_My(My), m_i_offset(i_offset), m_j_offset(j_offset),
32 m_array(reinterpret_cast<std::complex<double>*>(a)) {
33 (void) Mx;
34 }
35
36 // Access the whole array "a" of size My*My.
37 FFTWArray::FFTWArray(fftw_complex *a, int Mx, int My)
38 : m_My(My), m_i_offset(0), m_j_offset(0),
39 m_array(reinterpret_cast<std::complex<double>*>(a)) {
40 (void) Mx;
41 }
42
43
44 /*!
45 * Return the Discrete Fourier Transform sample frequencies.
46 */
47 std::vector<double> fftfreq(int M, double d) {
48 std::vector<double> result(M);
49
50 int N = M % 2 ? M / 2 : M / 2 - 1;
51
52 for (int i = 0; i <= N; i++) {
53 result[i] = i;
54 }
55
56 for (int i = N + 1; i < M; i++) {
57 result[i] = i - M;
58 }
59
60 // normalize
61 for (int i = 0; i < M; i++) {
62 result[i] /= (M * d);
63 }
64
65 return result;
66 }
67
68 //! \brief Fill `input` with zeros.
69 void clear_fftw_array(fftw_complex *input, int Nx, int Ny) {
70 FFTWArray fftw_in(input, Nx, Ny);
71 for (int i = 0; i < Nx; ++i) {
72 for (int j = 0; j < Ny; ++j) {
73 fftw_in(i, j) = 0.0;
74 }
75 }
76 }
77
78 //! @brief Copy `source` to `destination`.
79 void copy_fftw_array(fftw_complex *source, fftw_complex *destination, int Nx, int Ny) {
80 memcpy(destination, source, Nx * Ny * sizeof(fftw_complex));
81 }
82
83 //! Set the real part of output to input. Input has the size of My*Mx, embedded in the
84 //! bigger (output) grid of size Ny*Nx. Offsets i0 and j0 specify the location of the
85 //! subset to set.
86 /*!
87 * Sets the imaginary part to zero.
88 */
89 void set_real_part(Vec input,
90 double normalization,
91 int Mx, int My,
92 int Nx, int Ny,
93 int i0, int j0,
94 fftw_complex *output) {
95 petsc::VecArray2D in(input, Mx, My);
96 FFTWArray out(output, Nx, Ny, i0, j0);
97
98 for (int j = 0; j < My; ++j) {
99 for (int i = 0; i < Mx; ++i) {
100 out(i, j) = in(i, j) * normalization;
101 }
102 }
103 }
104
105
106 //! @brief Get the real part of input and put it in output.
107 /*!
108 * See set_real_part for details.
109 */
110 void get_real_part(fftw_complex *input,
111 double normalization,
112 int Mx, int My,
113 int Nx, int Ny,
114 int i0, int j0,
115 Vec output) {
116 petsc::VecArray2D out(output, Mx, My);
117 FFTWArray in(input, Nx, Ny, i0, j0);
118 for (int j = 0; j < My; ++j) {
119 for (int i = 0; i < Mx; ++i) {
120 out(i, j) = in(i, j).real() * normalization;
121 }
122 }
123 }
124
125 } // end of namespace pism