tfftw_utilities.cc - pism - [fork] customized build of PISM, the parallel ice s… | |
git clone git://src.adamsgaard.dk/pism | |
Log | |
Files | |
Refs | |
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… | |
6 * terms of the GNU General Public License as published by the Free Soft… | |
7 * Foundation; either version 3 of the License, or (at your option) any … | |
8 * version. | |
9 * | |
10 * PISM is distributed in the hope that it will be useful, but WITHOUT A… | |
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FIT… | |
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… | |
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 offset… | |
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 … | |
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, in… | |
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, e… | |
84 //! bigger (output) grid of size Ny*Nx. Offsets i0 and j0 specify the lo… | |
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 |