From 05b4239d1935c8b4599fc4f2e2e925d238573adb Mon Sep 17 00:00:00 2001 From: LSchueler Date: Tue, 9 Aug 2022 21:28:37 +0200 Subject: [PATCH] [WIP] Add support for nd vectors on md field This is a first suggestion of how to implement spatio-temporal vector fields. It is hardly tested and more a foundation for a discussion. --- src/gstools/field/base.py | 2 +- src/gstools/field/generator.py | 19 +++++++++++++++---- src/gstools/field/summator.pyx | 13 +++++++------ 3 files changed, 23 insertions(+), 11 deletions(-) diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py index 690042af..1a24ff8f 100755 --- a/src/gstools/field/base.py +++ b/src/gstools/field/base.py @@ -578,7 +578,7 @@ def pos(self, pos): self._pos, self._field_shape = format_struct_pos_dim(pos, self.dim) # prepend dimension if we have a vector field if self.value_type == "vector": - self._field_shape = (self.dim,) + self._field_shape + self._field_shape = (self.generator.vec_dim,) + self._field_shape if self.latlon: raise ValueError("Field: Vector fields not allowed for latlon") diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index 31cc7eea..3c7087a2 100644 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -339,6 +339,8 @@ class IncomprRandMeth(RandMeth): the mean velocity in x-direction mode_no : :class:`int`, optional number of Fourier modes. Default: ``1000`` + vec_dim : :class:`int`, optional + vector dimension, in case it mismatches the model dimension seed : :class:`int` or :any:`None`, optional the seed of the random number generator. If "None", a random seed is used. Default: :any:`None` @@ -391,18 +393,27 @@ def __init__( model, mean_velocity=1.0, mode_no=1000, + vec_dim=None, seed=None, verbose=False, sampling="auto", **kwargs, ): - if model.dim < 2 or model.dim > 3: + if vec_dim is None and (model.dim < 2 or model.dim > 3): raise ValueError( - "Only 2D and 3D incompressible fields can be generated." + "Only 2D and 3D incompressible vectors can be generated." + ) + if vec_dim is not None and (vec_dim < 2 or vec_dim > 3): + raise ValueError( + "Only 2D and 3D incompressible vectors can be generated." ) super().__init__(model, mode_no, seed, verbose, sampling, **kwargs) self.mean_u = mean_velocity + if vec_dim is None: + self.vec_dim = model.dim + else: + self.vec_dim = vec_dim self._value_type = "vector" def __call__(self, pos): @@ -425,7 +436,7 @@ def __call__(self, pos): """ pos = np.asarray(pos, dtype=np.double) summed_modes = summate_incompr( - self._cov_sample, self._z_1, self._z_2, pos + self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos ) nugget = self.get_nugget(summed_modes.shape) @@ -458,7 +469,7 @@ def _create_unit_vector(self, broadcast_shape, axis=0): the unit vector """ shape = np.ones(len(broadcast_shape), dtype=int) - shape[0] = self.model.dim + shape[0] = self.vec_dim e1 = np.zeros(shape) e1[axis] = 1.0 diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx index ecd7ea58..01414987 100644 --- a/src/gstools/field/summator.pyx +++ b/src/gstools/field/summator.pyx @@ -50,6 +50,7 @@ cdef (double) abs_square(const double[:] vec) nogil: def summate_incompr( + const int vec_dim, const double[:, :] cov_samples, const double[:] z_1, const double[:] z_2, @@ -58,24 +59,24 @@ def summate_incompr( cdef int i, j, d cdef double phase cdef double k_2 - cdef int dim = pos.shape[0] + cdef int field_dim = pos.shape[0] - cdef double[:] e1 = np.zeros(dim, dtype=float) + cdef double[:] e1 = np.zeros(vec_dim, dtype=float) e1[0] = 1. - cdef double[:] proj = np.empty(dim) + cdef double[:] proj = np.empty(vec_dim) cdef int X_len = pos.shape[1] cdef int N = cov_samples.shape[1] - cdef double[:, :] summed_modes = np.zeros((dim, X_len), dtype=float) + cdef double[:, :] summed_modes = np.zeros((vec_dim, X_len), dtype=float) for i in range(X_len): for j in range(N): k_2 = abs_square(cov_samples[:, j]) phase = 0. - for d in range(dim): + for d in range(field_dim): phase += cov_samples[d, j] * pos[d, i] - for d in range(dim): + for d in range(vec_dim): proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2 summed_modes[d, i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))