Fix Fourier meth. eq. and let user generate modes

src/gstools/field/generator.py

@@ -32,8 +32,8 @@
 else:
     from gstools.field.summator import (
         summate,
-        summate_incompr,
         summate_fourier,
+        summate_incompr,
     )
 
 __all__ = ["Generator", "RandMeth", "IncomprRandMeth", "Fourier"]
@@ -551,8 +551,6 @@ class Fourier(Generator):
         Number of Fourier modes per dimension.
     mode_truncation : :class:`list`
         Cut-off values of the Fourier modes.
-    period_len : :class:`float` or :class:`list`, optional
-        Period length of the field in each dim as a factor of the domain size.
     seed : :class:`int`, optional
         The seed of the random number generator.
         If "None", a random seed is used. Default: :any:`None`
@@ -594,49 +592,33 @@ class Fourier(Generator):
     def __init__(
         self,
         model,
-        modes_no,
-        modes_truncation,
-        period_len=None,
+        modes,
         seed=None,
         verbose=False,
         **kwargs,
     ):
         if kwargs:
             warnings.warn("gstools.Fourier: **kwargs are ignored")
         # initialize attributes
-        self._modes_truncation = self._fill_to_dim(
-            model.dim, modes_truncation, np.double
+        self._modes = generate_grid(modes)
+        self._modes_no = [len(m) for m in modes]
+        self._delta_k = np.array(
+            [modes[d][1] - modes[d][0] for d in range(model.dim)]
         )
-        self._modes_no = self._fill_to_dim(model.dim, modes_no, int)
-        self._modes = []
-        [
-            self._modes.append(
-                np.linspace(
-                    -self._modes_truncation[d] / 2,
-                    self._modes_truncation[d] / 2,
-                    self._modes_no[d],
-                    endpoint=False,
-                ).T
-            )
-            for d in range(model.dim)
-        ]
 
-        self._period_len = self._fill_to_dim(
-            model.dim, period_len, np.double, 1.0
-        )
         self._verbose = bool(verbose)
         # initialize private attributes
         self._model = None
         self._seed = None
         self._rng = None
         self._z_1 = None
         self._z_2 = None
-        self._spectral_density_sqrt = None
+        self._spectrum_factor = None
         self._value_type = "scalar"
         # set model and seed
         self.update(model, seed)
 
-    def __call__(self, pos, add_nugget=True, phase_factor=2.*np.pi, spec_factor=1., var_factor=1.,):
+    def __call__(self, pos, add_nugget=True):
         """Calculate the modes for the Fourier method.
 
         This method  calls the `summate_*` Cython methods, which are the
@@ -655,34 +637,17 @@ def __call__(self, pos, add_nugget=True, phase_factor=2.*np.pi, spec_factor=1.,
             the random modes
         """
         pos = np.asarray(pos, dtype=np.double)
-        domain_size = pos.max(axis=1) - pos.min(axis=1)
-        self._modes = [
-            self._modes[d] / domain_size[d] * self._period_len[d]
-            for d in range(self._model.dim)
-        ]
 
-        self._modes = generate_grid(self._modes)
-
-        # pre calc. the spectral density for all wave numbers
-        # they are handed over to Cython
-        k_norm = np.linalg.norm(self._modes, axis=0)
-        self._spectral_density_sqrt = np.sqrt(
-            self._model.spectral_density(k_norm)
-        )
         summed_modes = summate_fourier(
-            self._spectral_density_sqrt,
+            self._spectrum_factor,
             self._modes,
             self._z_1,
             self._z_2,
             pos,
-            phase_factor,
-            spec_factor,
+            config.NUM_THREADS,
         )
         nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0
-        return (
-            np.sqrt(var_factor * 2.0 * self.model.var / np.prod(domain_size)) * summed_modes
-            + nugget
-        )
+        return summed_modes + nugget
 
     def get_nugget(self, shape):
         """
@@ -782,6 +747,12 @@ def reset_seed(self, seed=np.nan):
         # normal distributed samples for randmeth
         self._z_1 = self._rng.random.normal(size=np.prod(self._modes_no))
         self._z_2 = self._rng.random.normal(size=np.prod(self._modes_no))
+        # pre calc. the spectrum for all wave numbers they are handed over to
+        # Cython, which doesn't have access to the CovModel
+        k_norm = np.linalg.norm(self._modes, axis=0)
+        self._spectrum_factor = np.sqrt(
+            2.0 * self._model.spectrum(k_norm) * np.prod(self._delta_k)
+        )
 
     def _fill_to_dim(self, dim, values, dtype, default_value=None):
         """Fill an array with last element up to len(dim)."""
@@ -824,26 +795,6 @@ def model(self):
     def model(self, model):
         self.update(model)
 
-    @property
-    def modes_truncation(self):
-        """:class:`list`: Cut-off values of the Fourier modes."""
-        return self._modes_truncation
-
-    @modes_truncation.setter
-    def modes_truncation(self, modes_truncation):
-        self._modes_truncation = modes_truncation
-
-    @property
-    def period_len(self):
-        """:class:`list`: Period length of the field in each dim."""
-        return self._period_len
-
-    @period_len.setter
-    def period_len(self, period_len):
-        self._period_len = self._fill_to_dim(
-            self._model.dim, period_len, np.double, 1.0
-        )
-
     @property
     def verbose(self):
         """:class:`bool`: Verbosity of the generator."""

src/gstools/field/srf.py

@@ -15,10 +15,10 @@
 
 from gstools.field.base import Field
 from gstools.field.generator import (
+    Fourier,
     Generator,
     IncomprRandMeth,
     RandMeth,
-    Fourier,
 )
 from gstools.field.upscaling import var_coarse_graining, var_no_scaling
 
@@ -116,9 +116,6 @@ def __call__(
         mesh_type="unstructured",
         post_process=True,
         store=True,
-        phase_factor=2.*np.pi,
-        spec_factor=1.,
-        var_factor=1.,
     ):
         """Generate the spatial random field.
 
@@ -163,18 +160,7 @@ def __call__(
         # get isometrized positions and the resulting field-shape
         iso_pos, shape = self.pre_pos(pos, mesh_type)
         # generate the field
-        try:
-            field = np.reshape(
-                self.generator(
-                    iso_pos,
-                    phase_factor=phase_factor,
-                    spec_factor=spec_factor,
-                    var_factor=var_factor
-                ),
-                shape,
-            )
-        except TypeError:
-            field = np.reshape(self.generator(iso_pos, ), shape)
+        field = np.reshape(self.generator(iso_pos), shape)
         # upscaled variance
         if not np.isscalar(point_volumes) or not np.isclose(point_volumes, 0):
             scaled_var = self.upscaling_func(self.model, point_volumes)

src/gstools/field/summator.pyx

@@ -10,7 +10,7 @@ IF OPENMP:
     cimport openmp
 
 cimport numpy as np
-from libc.math cimport pi, cos, sin, sqrt
+from libc.math cimport cos, pi, sin, sqrt
 
 
 def set_num_threads(num_threads):
@@ -100,13 +100,12 @@ def summate_incompr(
 
 
 def summate_fourier(
-    const double[:] spectral_density_sqrt,
+    const double[:] spectrum_factor,
     const double[:, :] modes,
     const double[:] z_1,
     const double[:] z_2,
     const double[:, :] pos,
-    const double phase_factor,
-    const double spec_factor,
+    num_threads=None,
     ):
     cdef int i, j, d
     cdef double phase
@@ -117,14 +116,15 @@ def summate_fourier(
 
     cdef double[:] summed_modes = np.zeros(X_len, dtype=float)
 
-    for i in prange(X_len, nogil=True):
+    cdef int num_threads_c = set_num_threads(num_threads)
+
+    for i in prange(X_len, nogil=True, num_threads=num_threads_c):
         for j in range(N):
             phase = 0.
             for d in range(dim):
                 phase += modes[d, j] * pos[d, i]
-            # OpenMP doesn't like *= after +=... seems to be a compiler specific thing
-            # phase = phase * 2. * pi
-            phase = phase * phase_factor
-            summed_modes[i] += spec_factor * spectral_density_sqrt[j] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))
+            summed_modes[i] += (
+                spectrum_factor[j] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))
+            )
 
     return np.asarray(summed_modes)

tests/test_fouriergen.py

@@ -14,57 +14,56 @@
 class TestFourier(unittest.TestCase):
     def setUp(self):
         self.seed = 19900408
-        self.cov_model_1d = gs.Gaussian(dim=1, var=0.5, len_scale=10.)
-        self.cov_model_2d = gs.Gaussian(dim=2, var=2.0, len_scale=30.)
-        self.cov_model_3d = gs.Gaussian(dim=3, var=2.1, len_scale=21.)
-        self.x = np.linspace(0, 80, 11)
-        self.y = np.linspace(0, 30, 31)
-        self.z = np.linspace(0, 91, 13)
+        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.modes_no_1d = 20
-        self.trunc_1d = 8
-        self.modes_no_2d = [16, 7]
-        self.trunc_2d = [16, 7]
-        self.modes_no_3d = [16, 7, 11]
-        self.trunc_3d = [16, 7, 12]
+        dk = [2 * np.pi / l for l in L]
+
+        self.modes_1d = [np.arange(0, 2, dk[0])]
+        self.modes_2d = self.modes_1d + [np.arange(0, 2, dk[1])]
+        self.modes_3d = self.modes_2d + [np.arange(0, 2, dk[2])]
 
         self.srf_1d = gs.SRF(
             self.cov_model_1d,
             generator="Fourier",
-            modes_no=self.modes_no_1d,
-            modes_truncation=self.trunc_1d,
+            modes=self.modes_1d,
             seed=self.seed,
         )
         self.srf_2d = gs.SRF(
             self.cov_model_2d,
             generator="Fourier",
-            modes_no=self.modes_no_2d,
-            modes_truncation=self.trunc_2d,
+            modes=self.modes_2d,
             seed=self.seed,
         )
         self.srf_3d = gs.SRF(
             self.cov_model_3d,
             generator="Fourier",
-            modes_no=self.modes_no_3d,
-            modes_truncation=self.trunc_3d,
+            modes=self.modes_3d,
             seed=self.seed,
         )
 
     def test_1d(self):
         field = self.srf_1d((self.x,), mesh_type="structured")
-        self.assertAlmostEqual(field[0], 0.9009981010688789)
+        self.assertAlmostEqual(field[0], 0.40939882638496783)
 
     def test_2d(self):
         field = self.srf_2d((self.x, self.y), mesh_type="structured")
-        self.assertAlmostEqual(field[0, 0], 1.1085370190533947)
+        self.assertAlmostEqual(field[0, 0], 0.8176311251780369)
 
     def test_3d(self):
         field = self.srf_3d((self.x, self.y, self.z), mesh_type="structured")
-        self.assertAlmostEqual(field[0, 0, 0], 1.7648407965681794)
+        self.assertAlmostEqual(field[0, 0, 0], -1.2636015063084773)
 
     def test_periodicity(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[0, len(self.y) // 2], field[-1, len(self.y) // 2]
+        )
 
     def test_assertions(self):
         # unstructured grids not supported