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.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2019, 2020, 2021, 2023 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 "pism/util/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 normalization) {
48  std::vector<double> result(M);
49 
50  int N = (M % 2) != 0 ? 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 * normalization);
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  */
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  petsc::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
FFTWArray(fftw_complex *a, int Mx, int My, int i_offset, int j_offset)
Wrapper around VecGetArray2d and VecRestoreArray2d.
Definition: Vec.hh:55
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)