From c3d17efa46f7df2733c5c605f0c687be38adaaee Mon Sep 17 00:00:00 2001 From: Joshua Benabou Date: Wed, 17 Aug 2022 00:20:01 -0700 Subject: [PATCH 1/2] added functionality for zero-velocity incompressible 3-component vector field (Kraichnan 1970 method) --- src/gstools/field/generator.py | 316 +++++++++++++++++++++++++++++++-- src/gstools/field/srf.py | 6 +- src/gstools/field/summator.pyx | 98 +++++++++- 3 files changed, 400 insertions(+), 20 deletions(-) diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index 6c8c905b..58e63eea 100644 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -9,6 +9,8 @@ .. autosummary:: RandMeth IncomprRandMeth + IncomprRandZeroVelMeth + GenericRandVectorFieldMeth """ # pylint: disable=C0103, W0222, C0412 import warnings @@ -24,9 +26,9 @@ # pylint: disable=E0401 from gstools_core import summate, summate_incompr else: - from gstools.field.summator import summate, summate_incompr + from gstools.field.summator import summate, summate_incompr, summate_incompr_zero_vel, summate_generic_vector_field -__all__ = ["RandMeth", "IncomprRandMeth"] +__all__ = ["RandMeth", "IncomprRandMeth", "IncomprRandZeroVelMeth", "GenericRandVectorFieldMeth", ] SAMPLING = ["auto", "inversion", "mcmc"] @@ -97,10 +99,10 @@ def __init__( ): if kwargs: warnings.warn("gstools.RandMeth: **kwargs are ignored") - # initialize attributes + # initialize atributes self._mode_no = int(mode_no) self._verbose = bool(verbose) - # initialize private attributes + # initialize private atributes self._model = None self._seed = None self._rng = None @@ -232,12 +234,14 @@ def reset_seed(self, seed=np.nan): if seed is None or not np.isnan(seed): self._seed = seed self._rng = RNG(self._seed) + # normal distributed samples for randmeth self._z_1 = self._rng.random.normal(size=self._mode_no) self._z_2 = self._rng.random.normal(size=self._mode_no) + # sample uniform on a sphere sphere_coord = self._rng.sample_sphere(self.model.dim, self._mode_no) - # sample radii according to radial spectral density of the model + # sample radii acording to radial spectral density of the model if self.sampling == "inversion" or ( self.sampling == "auto" and self.model.has_ppf ): @@ -339,6 +343,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 +397,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,14 +440,14 @@ 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) e1 = self._create_unit_vector(summed_modes.shape) return ( - self.mean_u * e1 + #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields + self.mean_u * np.sqrt(self.model.var / self._mode_no) * summed_modes @@ -458,8 +473,289 @@ 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 return e1 + +class IncomprRandZeroVelMeth(RandMeth): + r"""RandMeth for incompressible random vector fields with zero velocity, using Eq. 20 of Kraichnan (1970) + + Parameters + ---------- + model : :any:`CovModel` + covariance model + mean_velocity : :class:`float`, optional + 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` + verbose : :class:`bool`, optional + State if there should be output during the generation. + Default: :any:`False` + sampling : :class:`str`, optional + Sampling strategy. Either + + * "auto": select best strategy depending on given model + * "inversion": use inversion method + * "mcmc": use mcmc sampling + + **kwargs + Placeholder for keyword-args + + Notes + ----- + The Randomization method is used to generate isotropic + spatial incompressible random vector fields characterized + by a given covariance model. The equation is [Kraichnan1970]_: + + .. math:: + u\left(x\right)= + \bar{u}\sqrt{\frac{\sigma^{2}}{N}}\cdot + \sum_{j=1}^{N}\left( + W_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+ + W_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right) + \right) + + where: + + * :math:`\bar u` : mean velocity in :math:`e_1` direction + * :math:`N` : fourier mode number + * :math:`Z_{k,j}` : random samples from a normal distribution, each of size (N,vec_dim) + * :math:`k_j` : samples from the spectral density distribution of + the covariance model + * :math:`W_{k,j}` : cross-product of the random normal vector Z_{k,j} with the wave-vector k_j + ensuring the incompressibility + + References + ---------- + .. [Kraichnan1970] Kraichnan, R. H., + "Diffusion by a random velocity field.", + The physics of fluids, 13(1), 22-31., (1970) + """ + + def __init__( + self, + model, + mean_velocity=1.0, + mode_no=1000, + vec_dim=None, + seed=None, + verbose=False, + sampling="auto", + **kwargs, + ): + if vec_dim is None and (model.dim < 2 or model.dim > 3): + raise ValueError( + "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" + + # for using the Kraichnan method for zero-velocity fluid, z_1 and z_2 must instead contain N=self._mode_no independent realizations of normal vectors of size vec_dim, + mean = np.zeros(self.vec_dim) + cov = np.identity(self.vec_dim) + + self._z_1 = self._rng.random.multivariate_normal(mean, cov, size=self._mode_no) # shape (_mode_no, self.vec_dim) + self._z_2 = self._rng.random.multivariate_normal(mean, cov, size=self._mode_no) + print("shape of z_1: ",np.shape(self._z_1)) + + def __call__(self, pos): + """Calculate the random modes for the randomization method. + + This method calls the `summate_incompr_*` Cython methods, + which are the heart of the randomization method. + In this class the method contains a projector to + ensure the incompressibility of the vector field. + + Parameters + ---------- + pos : (d, n), :class:`numpy.ndarray` + the position tuple with d dimensions and n points. + + Returns + ------- + :class:`numpy.ndarray` + the random modes + """ + pos = np.asarray(pos, dtype=np.double) + print("\nStarting summate_incompr_zero_vel") + + summed_modes = summate_incompr_zero_vel( + self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos + ) + print("\nFinished mode summation!") + + nugget = self.get_nugget(summed_modes.shape) + + return ( + #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields + + self.mean_u + * np.sqrt(self.model.var / self._mode_no) + * summed_modes + + nugget + ) + + + +class GenericRandVectorFieldMeth(RandMeth): + r"""RandMeth for incompressible random vector fields. + + Parameters + ---------- + model : :any:`CovModel` + covariance model + mean_velocity : :class:`float`, optional + 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` + verbose : :class:`bool`, optional + State if there should be output during the generation. + Default: :any:`False` + sampling : :class:`str`, optional + Sampling strategy. Either + + * "auto": select best strategy depending on given model + * "inversion": use inversion method + * "mcmc": use mcmc sampling + + **kwargs + Placeholder for keyword-args + + Notes + ----- + The Randomization method is used to generate isotropic + spatial incompressible random vector fields characterized + by a given covariance model. The equation is [Kraichnan1970]_: + + .. math:: + u_i\left(x\right)= \bar{u_i} \delta_{i1} + + \bar{u_i}\sqrt{\frac{\sigma^{2}}{N}}\cdot + \sum_{j=1}^{N}p_i(k_{j})\left( + Z_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+ + Z_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right) + \right) + + where: + + * :math:`\bar u` : mean velocity in :math:`e_1` direction + * :math:`N` : fourier mode number + * :math:`Z_{k,j}` : random samples from a normal distribution + * :math:`k_j` : samples from the spectral density distribution of + the covariance model + * :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector + ensuring the incompressibility + + References + ---------- + .. [Kraichnan1970] Kraichnan, R. H., + "Diffusion by a random velocity field.", + The physics of fluids, 13(1), 22-31., (1970) + """ + + def __init__( + self, + model, + mean_velocity=1.0, + mode_no=1000, + vec_dim=None, + seed=None, + verbose=False, + sampling="auto", + **kwargs, + ): + if vec_dim is None and (model.dim < 2 or model.dim > 3): + raise ValueError( + "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): + """Calculate the random modes for the randomization method. + + This method calls the `summate_incompr_*` Cython methods, + which are the heart of the randomization method. + In this class the method contains a projector to + ensure the incompressibility of the vector field. + + Parameters + ---------- + pos : (d, n), :class:`numpy.ndarray` + the position tuple with d dimensions and n points. + + Returns + ------- + :class:`numpy.ndarray` + the random modes + """ + pos = np.asarray(pos, dtype=np.double) + summed_modes = summate_generic_vector_field( + self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos + ) + nugget = self.get_nugget(summed_modes.shape) + + e1 = self._create_unit_vector(summed_modes.shape) + + return ( + #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields + + self.mean_u + * np.sqrt(self.model.var / self._mode_no) + * summed_modes + + nugget + ) + + def _create_unit_vector(self, broadcast_shape, axis=0): + """Create a unit vector. + + Can be multiplied with a vector of shape broadcast_shape + + Parameters + ---------- + broadcast_shape : :class:`tuple` + the shape of the array with which + the unit vector is to be multiplied + axis : :class:`int`, optional + the direction of the unit vector. Default: ``0`` + + Returns + ------- + :class:`numpy.ndarray` + the unit vector + """ + shape = np.ones(len(broadcast_shape), dtype=int) + shape[0] = self.vec_dim + + e1 = np.zeros(shape) + e1[axis] = 1.0 + return e1 \ No newline at end of file diff --git a/src/gstools/field/srf.py b/src/gstools/field/srf.py index e8896138..88bb040c 100644 --- a/src/gstools/field/srf.py +++ b/src/gstools/field/srf.py @@ -13,7 +13,7 @@ import numpy as np from gstools.field.base import Field -from gstools.field.generator import IncomprRandMeth, RandMeth +from gstools.field.generator import IncomprRandMeth, IncomprRandZeroVelMeth, RandMeth, GenericRandVectorFieldMeth from gstools.field.upscaling import var_coarse_graining, var_no_scaling __all__ = ["SRF"] @@ -21,8 +21,10 @@ GENERATOR = { "RandMeth": RandMeth, "IncomprRandMeth": IncomprRandMeth, + "IncomprRandZeroVelMeth": IncomprRandZeroVelMeth, "VectorField": IncomprRandMeth, "VelocityField": IncomprRandMeth, + "GenericRandVectorFieldMeth": GenericRandVectorFieldMeth, } """dict: Standard generators for spatial random fields.""" @@ -119,7 +121,7 @@ def __call__( the position tuple, containing main direction and transversal directions seed : :class:`int`, optional - seed for RNG for resetting. Default: keep seed from generator + seed for RNG for reseting. Default: keep seed from generator point_volumes : :class:`float` or :class:`numpy.ndarray` If your evaluation points for the field are coming from a mesh, they are probably representing a certain element volume. diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx index ecd7ea58..2bc5047c 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,25 +59,106 @@ 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]) + k_2 = abs_square(cov_samples[:vec_dim, 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): - proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2 + for d in range(vec_dim): + proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2 #!!! this gives realizations with zero curl in x-direction summed_modes[d, i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase)) return np.asarray(summed_modes) + + +# function makes incompressible random vector field with zero velocity, unlike summate_incompr which makes the x-axis a preferential direction +# z_1 has shape (vec_dim,) such that z_1[i] for i =0,...,vec_dim-1 is a 1D array of length N = cov_samples.shape[1] +def summate_incompr_zero_vel( + const int vec_dim, + const double[:, :] cov_samples, + const double[:,:] z_1, + const double[:,:] z_2, + const double[:, :] pos + ): + cdef int i, j, d + cdef double phase + cdef int field_dim = pos.shape[0] + + cdef int X_len = pos.shape[1] + cdef int N = cov_samples.shape[1] + + cdef double[:, :] summed_modes = np.zeros((vec_dim, X_len), dtype=float) + cdef double[:] w_1 = np.empty(vec_dim, dtype=float) # !!! Joshua hopefully initialized w, v correctly, this is just for memory allocation right? + cdef double[:] w_2 = np.empty(vec_dim, dtype=float) + + for i in range(X_len): + for j in range(N): + # compute cross-product of wavevectors with random vectors z_1,z_2 of size vec_dim drawn from (vec_dim)-dimensional Gaussian with unit covariance matrix + + # !!! using the following two lines is cleaner, but code takes forever (type instability?) WHY????? + + #w_1 = np.cross(z_1[j,:], cov_samples[:vec_dim, j]) # 1D array of size vec_dim + #w_2 = np.cross(z_2[j,:], cov_samples[:vec_dim, j]) + + # Pending trying to get above two lines to work, I am hard coding the cross-product for vec_dim =3 + # !!! WARNING: this only works for vec_dim = 3 + w_1[0] = z_1[j,1] * cov_samples[2, j] - z_1[j,2] * cov_samples[1, j] + w_1[1] = z_1[j,2] * cov_samples[0, j] - z_1[j,0] * cov_samples[2, j] + w_1[2] = z_1[j,0] * cov_samples[1, j] - z_1[j,1] * cov_samples[0, j] + + w_2[0] = z_2[j,1] * cov_samples[2, j] - z_2[j,2] * cov_samples[1, j] + w_2[1] = z_2[j,2] * cov_samples[0, j] - z_2[j,0] * cov_samples[2, j] + w_2[2] = z_2[j,0] * cov_samples[1, j] - z_2[j,1] * cov_samples[0, j] + + phase = 0. + for d in range(field_dim): + phase += cov_samples[d, j] * pos[d, i] + for d in range(vec_dim): + summed_modes[d, i] += w_1[d] * cos(phase) + w_2[d] * sin(phase) #1 * cos(phase) + 1 * sin(phase) + + return np.asarray(summed_modes) + +def summate_generic_vector_field( + const int vec_dim, + const double[:, :] cov_samples, + const double[:] z_1, + const double[:] z_2, + const double[:, :] pos + ): + cdef int i, j, d + cdef double phase + cdef double k_2 + cdef int field_dim = pos.shape[0] + + cdef double[:] e1 = np.zeros(vec_dim, dtype=float) + e1[0] = 1. + 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((vec_dim, X_len), dtype=float) + + for i in range(X_len): + for j in range(N): + k_2 = abs_square(cov_samples[:vec_dim, j]) + phase = 0. + for d in range(field_dim): + phase += cov_samples[d, j] * pos[d, i] + for d in range(vec_dim): + #proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2 #!!! don't want incompressibility projector here + summed_modes[d, i] += (z_1[j] * cos(phase) + z_2[j] * sin(phase)) + + return np.asarray(summed_modes) From 48d0decb5caa71ec36142efa39721c56f21242f4 Mon Sep 17 00:00:00 2001 From: Joshua Benabou Date: Tue, 23 Aug 2022 20:19:54 -0700 Subject: [PATCH 2/2] added periodic boundary conditions for spatial dimensions, feature implemented only for generator IncomprRandZeroVelMeth --- src/gstools/field/generator.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index 58e63eea..1b157310 100644 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -175,6 +175,7 @@ def update(self, model=None, seed=np.nan): If :any:`None`, a random seed is used. If :any:`numpy.nan`, the actual seed will be kept. Default: :any:`numpy.nan` """ + # check if a new model is given if isinstance(model, CovModel): if self.model != model: @@ -234,7 +235,7 @@ def reset_seed(self, seed=np.nan): if seed is None or not np.isnan(seed): self._seed = seed self._rng = RNG(self._seed) - + # normal distributed samples for randmeth self._z_1 = self._rng.random.normal(size=self._mode_no) self._z_2 = self._rng.random.normal(size=self._mode_no) @@ -249,12 +250,14 @@ def reset_seed(self, seed=np.nan): rad = self._rng.sample_dist( size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0 ) + else: rad = self._rng.sample_ln_pdf( ln_pdf=self.model.ln_spectral_rad_pdf, size=self._mode_no, sample_around=1.0 / self.model.len_rescaled, ) + # get fully spatial samples by multiplying sphere samples and radii self._cov_sample = rad * sphere_coord @@ -548,6 +551,8 @@ def __init__( seed=None, verbose=False, sampling="auto", + periodic_bc=False, + box_len=None, **kwargs, ): if vec_dim is None and (model.dim < 2 or model.dim > 3): @@ -567,6 +572,9 @@ def __init__( self.vec_dim = vec_dim self._value_type = "vector" + self.periodic_bc=periodic_bc + self.box_len=box_len + # for using the Kraichnan method for zero-velocity fluid, z_1 and z_2 must instead contain N=self._mode_no independent realizations of normal vectors of size vec_dim, mean = np.zeros(self.vec_dim) cov = np.identity(self.vec_dim) @@ -594,8 +602,21 @@ def __call__(self, pos): the random modes """ pos = np.asarray(pos, dtype=np.double) + print("\nStarting summate_incompr_zero_vel") + if self.periodic_bc: + if not self.box_len==None: + print("\nimposing periodic boundary conditions on spatial coordinates!") + print("\nrounding first 'vec_dim' vector components in cov_sample to multiples of 2pi/box_length") + fac = 2*np.pi/np.array(self.box_len) + self._cov_sample[:self.vec_dim,:] = fac[:,None]*np.round(self._cov_sample[:self.vec_dim,:]/fac[:,None]) + else: + raise ValueError( + "For periodic boundary conditions on spatial coordinates, specify parameter box_len, an array of lengths of the box in each spatial dimension. The length of box_len must be equal to vec_dim." + ) + + summed_modes = summate_incompr_zero_vel( self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos )