Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
27namespace pism {
28
29namespace petsc {
30class 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 */
40class FFTWArray {
41public:
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
53private:
55 std::complex<double> *m_array;
56};
57
58std::vector<double> fftfreq(int M, double normalization);
59
60//! Fill `input` with zeros.
61void clear_fftw_array(fftw_complex *input, int Nx, int Ny);
62
63//! Copy `source` to `destination`.
64void 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 */
72void 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 */
83void 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
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)