diff --git a/examples/00_misc/06_fourier.py b/examples/00_misc/06_fourier.py index 5dc1c608..360f4708 100644 --- a/examples/00_misc/06_fourier.py +++ b/examples/00_misc/06_fourier.py @@ -21,16 +21,12 @@ model = gs.Gaussian(dim=2, var=1, len_scale=200) # Next, we hand the cov. model to the spatial random field class -# and set the generator to `Fourier`. We will let the class figure out the -# modes internally, by handing over `period` and `mode_rel_cutoff` which is the cutoff -# value of the spectral density, relative to the maximum spectral density at -# the origin. Simply put, we will use `mode_rel_cutoff`% of the spectral -# density for the calculations. The argument `period` is set to the domain -# size. +# and set the generator to `Fourier`. The `mode_no` argument sets the number of +# Fourier modes per dimension. The argument `period` is set to the domain size. srf = gs.SRF( model, generator="Fourier", - mode_rel_cutoff=0.99, + mode_no=[32, 32], period=L, seed=1681903, ) @@ -40,24 +36,3 @@ # GSTools has a few simple visualization methods built in. srf.plot() - -# Alternatively, we could calculate the modes ourselves and hand them over to -# GSTools. Therefore, we set the cutoff values to absolut values in Fourier -# space. But always check, if you cover enough of the spectral density to not -# run into numerical problems. -modes_cutoff = [1.0, 1.0] - -# Next, we have to compute the numerical step size in Fourier space. This choice -# influences the periodicity, which we want to set to the domain size by -modes_delta = 2 * np.pi / L - -# Now, we calculate the modes with -modes = [np.arange(0, modes_cutoff[d], modes_delta[d]) for d in range(2)] - -# And we can create a new instance of the SRF class with our own modes. -srf_modes = gs.SRF( - model, - generator="Fourier", - modes=modes, - seed=494754, -) diff --git a/examples/00_misc/07_fourier_trans.py b/examples/00_misc/07_fourier_trans.py index f1d40be3..8331f766 100644 --- a/examples/00_misc/07_fourier_trans.py +++ b/examples/00_misc/07_fourier_trans.py @@ -24,7 +24,7 @@ srf = gs.SRF( model, generator="Fourier", - mode_rel_cutoff=0.999, + mode_no=[32, 32], period=L, seed=1681903, ) diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index c6c74aee..712fbece 100755 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -599,49 +599,35 @@ class Fourier(Generator): def __init__( self, model, - mode_rel_cutoff=None, - period=None, - modes=None, + mode_no, + period, seed=None, verbose=False, **kwargs, ): if kwargs: warnings.warn("gstools.Fourier: **kwargs are ignored") - if ( - mode_rel_cutoff is not None - and period is not None - and modes is not None - ): - warnings.warn( - "gstools.Fourier: mode_rel_cutoff & period are ignored, as " - "modes is provided." - ) - if mode_rel_cutoff is None and period is None and modes is None: - raise ValueError("Fourier: No mode information provided.") - dim = model.dim - if ( - modes is None - and period is not None - and mode_rel_cutoff is not None - ): - if len(period) != dim: - raise ValueError("Fourier: Dimension mismatch.") - self._mode_rel_cutoff = mode_rel_cutoff - self._period = np.array(period) - self._delta_k = 2.0 * np.pi / self._period - modes_cutoff = self.calc_modes_cutoff(model, self._mode_rel_cutoff) - self._modes_cutoff = self._fill_to_dim(dim, modes_cutoff) - modes = [ - np.arange(-self._modes_cutoff[d], self._modes_cutoff[d], self._delta_k[d]) - for d in range(dim) - ] - elif modes is not None: - self._delta_k = np.array( - [modes[d][1] - modes[d][0] for d in range(dim)] + if len(mode_no) != dim: + raise ValueError( + "Fourier: Dimension mismatch in argument mode_no." + ) + if len(period) != dim: + raise ValueError("Fourier: Dimension mismatch in argument period.") + if (np.asarray([m % 2 for m in mode_no]) != 0).any(): + raise ValueError("Fourier: Odd mode_no not supported.") + + self._period = np.array(period) + self._delta_k = 2.0 * np.pi / self._period + anis = np.insert(model.anis.copy(), 0, 1.0) + modes = [ + np.arange( + -mode_no[d] / 2.0 * self._delta_k[d] / anis[d], + mode_no[d] / 2.0 * self._delta_k[d] / anis[d], + self._delta_k[d], ) - self._period = 2.0 * np.pi / self._delta_k + for d in range(dim) + ] # initialize attributes self._modes = generate_grid(modes) @@ -813,41 +799,6 @@ def _fill_to_dim( r = np.pad(r, (0, dim - len(r)), "edge") return r - def calc_modes_cutoff( - self, model, mode_rel_cutoff - ): # pylint: disable=R6301 - """Find the cutoff value so that `mode_rel_cutoff`% of the spectrum is kept. - - This helper function uses a least squares algorithm to determine the - cutoff value so that `mode_rel_cutoff`% of the spectrum is kept. - - Parameters - ---------- - model : :any:`CovModel` - Covariance model - mode_rel_cutoff : :class:`float` - - Returns - ------- - :class:`float` - the cutoff value - """ - norm = model.spectral_density(0) - # the first len_scale is a good enough first guess - try: - len_scale = model.len_scale[0] - except IndexError: - len_scale = model.len_scale - k_cutoff0 = np.sqrt(1.0 / len_scale) - mode_rel_cutoff_r = 1.0 - mode_rel_cutoff - res = least_squares( - lambda k, mode_rel_cutoff_r: model.spectral_density(k) - - mode_rel_cutoff_r * norm, - k_cutoff0, - args=(mode_rel_cutoff_r,), - ) - return res.x[0] - @property def seed(self): """:class:`int`: Seed of the master RNG. @@ -873,6 +824,16 @@ def model(self): def model(self, model): self.update(model) + @property + def modes(self): + """:class:`numpy.ndarray`: Modes on which the spectrum is calculated.""" + return self._modes + + @property + def period(self): + """:class:`numpy.ndarray`: Period length of the spatial random field.""" + return self._period + @property def verbose(self): """:class:`bool`: Verbosity of the generator.""" diff --git a/tests/test_fouriergen.py b/tests/test_fouriergen.py old mode 100644 new mode 100755 index 0a1c3cb9..23ea7267 --- a/tests/test_fouriergen.py +++ b/tests/test_fouriergen.py @@ -2,13 +2,11 @@ This is the unittest of the Fourier class. """ -import copy import unittest import numpy as np import gstools as gs -from gstools.field.generator import Fourier class TestFourier(unittest.TestCase): @@ -17,56 +15,70 @@ def setUp(self): self.cov_model_1d = gs.Gaussian(dim=1, var=0.5, len_scale=10.0) self.cov_model_2d = gs.Gaussian(dim=2, var=2.0, len_scale=30.0) self.cov_model_3d = gs.Gaussian(dim=3, var=2.1, len_scale=21.0) - L = [80, 30, 91] - self.x = np.linspace(0, L[0], 11) - self.y = np.linspace(0, L[1], 31) - self.z = np.linspace(0, L[2], 13) + self.L = [80, 30, 91] + self.x = np.linspace(0, self.L[0], 11) + self.y = np.linspace(0, self.L[1], 31) + self.z = np.linspace(0, self.L[2], 13) - cutoff_rel = 0.999 - cutoff_abs = 1 - dk = [2 * np.pi / l for l in L] - - self.modes_1d = [np.arange(0, cutoff_abs, dk[0])] - self.modes_2d = self.modes_1d + [np.arange(0, cutoff_abs, dk[1])] - self.modes_3d = self.modes_2d + [np.arange(0, cutoff_abs, dk[2])] + self.mode_no = [12, 6, 14] self.srf_1d = gs.SRF( self.cov_model_1d, generator="Fourier", - modes=self.modes_1d, + mode_no=[self.mode_no[0]], + period=[self.L[0]], seed=self.seed, ) self.srf_2d = gs.SRF( self.cov_model_2d, generator="Fourier", - modes=self.modes_2d, + mode_no=self.mode_no[:2], + period=self.L[:2], seed=self.seed, ) self.srf_3d = gs.SRF( self.cov_model_3d, generator="Fourier", - mode_rel_cutoff=cutoff_rel, - period=L, + mode_no=self.mode_no, + period=self.L, seed=self.seed, ) def test_1d(self): field = self.srf_1d((self.x,), mesh_type="structured") - self.assertAlmostEqual(field[0], 0.40939877176695477) + self.assertAlmostEqual(field[0], 0.6236929351309081) def test_2d(self): field = self.srf_2d((self.x, self.y), mesh_type="structured") - self.assertAlmostEqual(field[0, 0], 1.6338790313270515) + self.assertAlmostEqual(field[0, 0], -0.1431996611581266) def test_3d(self): field = self.srf_3d((self.x, self.y, self.z), mesh_type="structured") - self.assertAlmostEqual(field[0, 0, 0], 0.3866689424599251) + self.assertAlmostEqual(field[0, 0, 0], -1.0433325279452803) + + def test_periodicity_1d(self): + field = self.srf_1d((self.x,), mesh_type="structured") + self.assertAlmostEqual(field[0], field[-1]) - def test_periodicity(self): + def test_periodicity_2d(self): field = self.srf_2d((self.x, self.y), mesh_type="structured") self.assertAlmostEqual( field[0, len(self.y) // 2], field[-1, len(self.y) // 2] ) + self.assertAlmostEqual( + field[len(self.x) // 2, 0], field[len(self.x) // 2, -1] + ) + + def test_periodicity_3d(self): + field = self.srf_3d((self.x, self.y, self.z), mesh_type="structured") + self.assertAlmostEqual( + field[0, len(self.y) // 2, 0], field[-1, len(self.y) // 2, 0] + ) + self.assertAlmostEqual(field[0, 0, 0], field[0, -1, 0]) + self.assertAlmostEqual( + field[len(self.x) // 2, len(self.y) // 2, 0], + field[len(self.x) // 2, len(self.y) // 2, -1], + ) def test_assertions(self): # unstructured grids not supported @@ -74,5 +86,12 @@ def test_assertions(self): self.assertRaises( ValueError, self.srf_2d, (self.x, self.y), mesh_type="unstructured" ) - with self.assertRaises(ValueError): - gs.SRF(self.cov_model_2d, generator="Fourier") + self.assertRaises( + ValueError, + gs.SRF, + self.cov_model_2d, + generator="Fourier", + mode_no=[13, 50], + period=self.L[:2], + seed=self.seed, + )