from pprint import pprint
import numpy as np
from numpy.fft import fft, ifft, fftshift

class FourierFourierSpatialDiscretization:

    def __init__(self, config):
        self.length_x = config['length_x']
        self.length_y = config['length_y']
        self.num_points_x = config['num_points_x']
        self.num_points_y = config['num_points_y']

        self.__build_grid__()
        self.__build_wavenumbers__()

    def __build_grid__(self):
        x_1d, dx = np.linspace(0, self.length_x, self.num_points_x, False, True)
        y_1d, dy = np.linspace(0, self.length_y, self.num_points_y, False, True)
        self.x, self.y = np.meshgrid(x_1d, y_1d)

        self.min_grid_spacing = min(dx,dy)

    def __build_wavenumbers__(self):
        k_1d = fftshift(np.array(range(0, self.num_points_x))*
                                               ((2*np.pi)/self.length_x))
        l_1d = fftshift(np.array(range(0, self.num_points_y))*
                                               ((2*np.pi)/self.length_y))
        self.k, self.l = np.meshgrid(k_1d, l_1d)

    def differentiate_x(self, field):
        return np.real(ifft((1.j*self.k)*fft(field, axis=1), axis=1))

    def differentiate_y(self, field):
        return np.real(ifft((1.j*self.l)*fft(field, axis=0), axis=0))


@@ -412,16 +412,48 @@
 bers__()
+%0A        self.__build_filter__()
 %0A%0A    de
@@ -1047,56 +1047,991 @@
 y))%0A
-        self.k, self.l = np.meshgrid(k_1d, l_1d)
+%0A        self.kmax = max(k_1d)%0A        self.lmax = max(l_1d)%0A        self.k, self.l = np.meshgrid(k_1d, l_1d)%0A%0A%0A    def __build_filter__(self):%0A        cutoff = 0.65 # 0.65 typically,  CUTOFF OF 0 CORRESPONDS TO HYPERVISCOSITY%0A        epsf = 1e-15 # FILTER STRENGTH AT HIGH WAVENUMBERS%0A%0A        kcrit = self.kmax*cutoff%0A        lcrit = self.lmax*cutoff%0A        filter_order = 4%0A%0A        self.filter = np.ones(%5Bself.num_points_y, self.num_points_x%5D)%0A        # Filter in x.%0A        mask = np.abs(self.k) %3C kcrit%0A        self.filter *= (mask + (1-mask)*np.exp(np.log(epsf)*(np.power((self.k-kcrit)/(self.kmax-kcrit), filter_order))))%0A%0A        # Filter in y.%0A        mask = np.abs(self.l) %3C lcrit%0A        self.filter *= (mask + (1-mask)*np.exp(np.log(epsf)*(np.power((self.l-lcrit)/(self.lmax-lcrit), filter_order))))%0A%0A        # Remove the Nyquist frequency entirely%0A        self.filter%5B:, np.floor(self.num_points_x/2+1)%5D = 0%0A        self.filter%5Bnp.floor(self.num_points_y/2+1), :%5D = 0%0A
 %0A%0A  

Write implementation for spectral filter (Fourier-Fourier).