Source code for fluidsim.base.forcing.anisotropic
# -*- coding: utf-8 -*-
""" Anisotropic (:mod:`fluidsim.base.forcing.anisotropic`)
==========================================================
.. autoclass:: TimeCorrelatedRandomPseudoSpectralAnisotropic
:members:
:private-members:
"""
from math import degrees
from math import pi
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from fluiddyn.calcul.easypyfft import fftw_grid_size
from fluidsim.base.forcing.specific import TimeCorrelatedRandomPseudoSpectral
from fluidsim.util import ensure_radians
[docs]
class TimeCorrelatedRandomPseudoSpectralAnisotropic(
TimeCorrelatedRandomPseudoSpectral
):
"""Random normalized anisotropic forcing.
.. inheritance-diagram:: TimeCorrelatedRandomPseudoSpectralAnisotropic
"""
tag = "tcrandom_anisotropic"
[docs]
@classmethod
def _complete_params_with_default(cls, params):
"""This static method is used to complete the *params* container."""
super(
TimeCorrelatedRandomPseudoSpectral, cls
)._complete_params_with_default(params)
params.forcing._set_child(
"tcrandom_anisotropic",
{"angle": "45°", "delta_angle": None, "kz_negative_enable": False},
)
params.forcing.tcrandom_anisotropic._set_doc(
"""
See :mod:`fluidsim.base.forcing.anisotropic`.
angle: float or str
Angle between the wavevector and the horizontal characterising the forcing
region.
delta_angle: float or None
Control the shape of the forcing region in k-space. If None, ``(khmin,
khmax, kvmin, kvmax)`` are computed from the angle,
``params.forcing.nkmin_forcing`` and ``params.forcing.nkmax_forcing``.
kz_negative_enable: bool
If True, modes with negative kz are also forced.
"""
)
def __init__(self, sim):
super().__init__(sim)
if self.params.forcing.normalized.type == "particular_k":
raise NotImplementedError
def _create_params_coarse(self, fft_size):
params_coarse = super()._create_params_coarse(fft_size)
self.angle = angle = ensure_radians(self.params.forcing[self.tag].angle)
tmp = self.params.forcing.tcrandom_anisotropic
try:
delta_angle = tmp.delta_angle
except AttributeError:
# loading old simul with delta_angle
delta_angle = None
else:
delta_angle = ensure_radians(delta_angle)
if delta_angle is None:
self.khmax_forcing = np.sin(angle) * self.kmax_forcing
self.kvmax_forcing = np.cos(angle) * self.kmax_forcing
else:
self.khmin_forcing = (
np.sin(angle - 0.5 * delta_angle) * self.kmin_forcing
)
self.kvmin_forcing = (
np.cos(angle + 0.5 * delta_angle) * self.kmin_forcing
)
self.khmax_forcing = (
np.sin(angle + 0.5 * delta_angle) * self.kmax_forcing
)
self.kvmax_forcing = (
np.cos(angle - 0.5 * delta_angle) * self.kmax_forcing
)
if hasattr(params_coarse.oper, "nz"):
# 3d
kymax_forcing = self.khmax_forcing
else:
# 2d
kymax_forcing = self.kvmax_forcing
# The "+ 1" aims to give some gap between the kxmax and
# the boundary of the oper_coarse.
try:
params_coarse.oper.nx = 2 * fftw_grid_size(
int(self.khmax_forcing / self.oper.deltakx) + 1
)
except AttributeError:
pass
try:
params_coarse.oper.ny
except AttributeError:
pass
else:
params_coarse.oper.ny = 2 * fftw_grid_size(
int(kymax_forcing / self.oper.deltaky) + 1
)
try:
params_coarse.oper.nz
except AttributeError:
pass
else:
params_coarse.oper.nz = 2 * fftw_grid_size(
int(self.kvmax_forcing / self.oper.deltakz) + 1
)
return params_coarse
[docs]
def _compute_cond_no_forcing(self):
"""Computes condition no forcing of the anisotropic case."""
angle = self.angle
tmp = self.params.forcing.tcrandom_anisotropic
try:
delta_angle = tmp.delta_angle
except AttributeError:
# loading old simul with delta_angle
delta_angle = None
else:
delta_angle = ensure_radians(delta_angle)
kf_min = self.kmin_forcing
kf_max = self.kmax_forcing
try:
self.params.oper.nz
except AttributeError:
ndim = 2
else:
ndim = 3
if delta_angle is None:
self.khmin_forcing = np.sin(angle) * self.kmin_forcing
self.kvmin_forcing = np.cos(angle) * self.kmin_forcing
if ndim == 2:
Kh = self.oper_coarse.KX
Kv = self.oper_coarse.KY
else:
Kh = np.sqrt(self.oper_coarse.Kx**2 + self.oper_coarse.Ky**2)
Kv = self.oper_coarse.Kz
COND_NO_F_KH = np.logical_or(
Kh > self.khmax_forcing, Kh < self.khmin_forcing
)
COND_NO_F_KV = np.logical_or(
Kv > self.kvmax_forcing, Kv < self.kvmin_forcing
)
if self.params.forcing.tcrandom_anisotropic.kz_negative_enable:
COND_NO_F_KV = np.logical_and(
COND_NO_F_KV,
np.logical_or(
Kv < -self.kvmax_forcing, Kv > -self.kvmin_forcing
),
)
COND_NO_F = np.logical_or(COND_NO_F_KH, COND_NO_F_KV)
COND_NO_F[self.oper_coarse.shapeK_loc[0] // 2] = True
COND_NO_F[:, self.oper_coarse.shapeK_loc[1] - 1] = True
else:
if ndim == 2:
K = np.sqrt(self.oper_coarse.KX**2 + self.oper_coarse.KY**2)
Kv = self.oper_coarse.KY
else:
K = np.sqrt(
self.oper_coarse.Kx**2
+ self.oper_coarse.Ky**2
+ self.oper_coarse.Kz**2
)
Kv = self.oper_coarse.Kz
K_nozero = K.copy()
K_nozero[K_nozero == 0] = 1e-14
theta = np.arccos(Kv / K_nozero)
del K_nozero
COND_NO_F_K = np.logical_or(K > kf_max, K < kf_min)
COND_NO_F_THETA = np.logical_or(
theta > angle + 0.5 * delta_angle,
theta < angle - 0.5 * delta_angle,
)
if self.params.forcing.tcrandom_anisotropic.kz_negative_enable:
COND_NO_F_THETA = np.logical_and(
COND_NO_F_THETA,
np.logical_or(
theta < pi - angle - 0.5 * delta_angle,
theta > pi - angle + 0.5 * delta_angle,
),
)
COND_NO_F = np.logical_or(COND_NO_F_K, COND_NO_F_THETA)
COND_NO_F[self.oper_coarse.shapeK_loc[0] // 2] = True
COND_NO_F[:, self.oper_coarse.shapeK_loc[1] - 1] = True
return COND_NO_F
[docs]
def plot_forcing_region(self):
"""Plots the forcing region"""
pforcing = self.params.forcing
khmin_forcing = self.khmin_forcing
khmax_forcing = self.khmax_forcing
kvmin_forcing = self.kvmin_forcing
kvmax_forcing = self.kvmax_forcing
kf_min = self.kmin_forcing
kf_max = self.kmax_forcing
tmp = self.params.forcing.tcrandom_anisotropic
try:
delta_angle = tmp.delta_angle
except AttributeError:
# loading old simul with delta_angle
delta_angle = None
try:
self.params.oper.nz
except AttributeError:
ndim = 2
else:
ndim = 3
# Define forcing region
coord_x = khmin_forcing
coord_y = kvmin_forcing
width = khmax_forcing - khmin_forcing
height = kvmax_forcing - kvmin_forcing
if ndim == 2:
Kh = self.oper_coarse.KX
Kv = self.oper_coarse.KY
deltakh = self.oper.deltakx
deltakv = self.oper.deltaky
else:
Kh = np.sqrt(self.oper_coarse.Kx**2 + self.oper_coarse.Ky**2)
Kv = self.oper_coarse.Kz
deltakh = self.oper.deltakx
deltakv = self.oper.deltakz
fig, ax = plt.subplots()
ax.set_aspect("equal")
title = (
pforcing.type
+ "; "
+ rf"$nk_{{min}} = {pforcing.nkmin_forcing} \delta k_v$; "
+ rf"$nk_{{max}} = {pforcing.nkmax_forcing} \delta k_v$; "
+ "\n"
+ r"$\theta_f = {:.0f}^\circ$; ".format(degrees(self.angle))
+ rf"Forced modes = {self.nb_forced_modes}"
)
ax.set_title(title)
ax.set_xlabel(r"$k_h$")
ax.set_ylabel(r"$k_v$")
# Parameters figure
# Set limits to 125% of the kf_max
factor = 1.2
ax.set_xlim([0.0, factor * kf_max])
ax.set_ylim([0.0, factor * kf_max])
xticks = np.arange(0.0, factor * kf_max, deltakv)
yticks = np.arange(0.0, factor * kf_max, deltakv)
ax.set_xticks(xticks)
ax.set_yticks(yticks)
if delta_angle is None:
# Plot forcing region
ax.add_patch(
patches.Rectangle(
xy=(coord_x, coord_y), width=width, height=height, fill=False
)
)
# Plot lines forcing region
ax.plot(
[khmin_forcing, khmin_forcing],
[0, kvmin_forcing],
"k--",
linewidth=0.8,
)
ax.plot(
[khmax_forcing, khmax_forcing],
[0, kvmin_forcing],
"k--",
linewidth=0.8,
)
ax.plot(
[0, khmin_forcing],
[kvmin_forcing, kvmin_forcing],
"k--",
linewidth=0.8,
)
ax.plot(
[0, khmin_forcing],
[kvmax_forcing, kvmax_forcing],
"k--",
linewidth=0.8,
)
# Location labels 0.8% the length of the axis
factor = 0.008
loc_label_y = abs(Kv).max() * factor
loc_label_x = abs(Kh).max() * factor
ax.text(loc_label_x + khmin_forcing, loc_label_y, r"$k_{h,min}$")
ax.text(loc_label_x + khmax_forcing, loc_label_y, r"$k_{h,max}$")
ax.text(loc_label_x, kvmin_forcing + loc_label_y, r"$k_{v,min}$")
ax.text(loc_label_x, kvmax_forcing + loc_label_y, r"$k_{v,max}$")
else:
# Plot forcing region
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=(kf_min + kf_max),
height=(kf_min + kf_max),
angle=0,
theta1=90.0 - degrees(self.angle),
theta2=90.0,
linestyle="dotted",
)
)
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=2.1 * kf_max,
height=2.1 * kf_max,
angle=0,
theta1=90.0
- degrees(self.angle)
- 0.5 * degrees(delta_angle),
theta2=90.0
- degrees(self.angle)
+ 0.5 * degrees(delta_angle),
linestyle="--",
)
)
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=2 * kf_min,
height=2 * kf_min,
angle=0,
theta1=90.0
- degrees(self.angle)
- 0.5 * degrees(delta_angle),
theta2=90.0
- degrees(self.angle)
+ 0.5 * degrees(delta_angle),
linestyle="-",
)
)
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=2 * kf_max,
height=2 * kf_max,
angle=0,
theta1=90.0
- degrees(self.angle)
- 0.5 * degrees(delta_angle),
theta2=90.0
- degrees(self.angle)
+ 0.5 * degrees(delta_angle),
linestyle="-",
)
)
# Plot arc kmin and kmax
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=2 * kf_min,
height=2 * kf_min,
angle=0,
theta1=0.0,
theta2=90.0,
linestyle="-.",
)
)
ax.add_patch(
patches.Arc(
xy=(0, 0),
width=2 * kf_max,
height=2 * kf_max,
angle=0,
theta1=0.0,
theta2=90.0,
linestyle="-.",
)
)
# Plot lines angle & lines forcing region
xmin = khmin_forcing
xmax = self.kmax_forcing * np.sin(self.angle - 0.5 * delta_angle)
ymin = self.kmin_forcing * np.cos(self.angle - 0.5 * delta_angle)
ymax = kvmax_forcing
ax.plot([xmin, xmax], [ymin, ymax], color="k", linewidth=1)
xmin = self.kmin_forcing * np.sin(self.angle + 0.5 * delta_angle)
xmax = khmax_forcing
ymin = kvmin_forcing
ymax = self.kmax_forcing * np.cos(self.angle + 0.5 * delta_angle)
ax.plot([xmin, xmax], [ymin, ymax], color="k", linewidth=1)
# Location labels kmin and kmax
factor = 0.015
loc_label_y = abs(Kv).max() * factor
loc_label_x = abs(Kh).max() * factor
ax.text(loc_label_y + self.kmin_forcing, loc_label_y, r"$k_{f,min}$")
ax.text(loc_label_x + self.kmax_forcing, loc_label_y, r"$k_{f,max}$")
# Location label angle \theta
factor = 1.1
loc_label_y = (
(kf_min + kf_max) * 0.5 * np.cos(self.angle * 0.5) * factor
)
loc_label_x = (
(kf_min + kf_max) * 0.5 * np.sin(self.angle * 0.5) * factor
)
ax.text(loc_label_x, loc_label_y, r"$\theta_f$")
# Location label delta_angle \delta \theta
factor = 1.1
loc_label_y = kf_max * np.cos(self.angle) * factor
loc_label_x = kf_max * np.sin(self.angle) * factor
ax.text(loc_label_x, loc_label_y, r"$\delta \theta_f$")
# Plot forced modes in red
indices_forcing = np.argwhere(self.COND_NO_F == False)
for i, index in enumerate(indices_forcing):
if ndim == 2:
ax.plot(
Kh[0, index[1]],
Kv[index[0], 0],
"ro",
label="Forced mode" if i == 0 else "",
)
else:
ax.plot(
Kh[0, index[1], index[2]],
Kv[index[0], 0, 0],
"ro",
label="Forced mode" if i == 0 else "",
)
ax.grid(linestyle="--", alpha=0.4)
ax.legend()
class TimeCorrelatedRandomPseudoSpectralAnisotropic3D(
TimeCorrelatedRandomPseudoSpectralAnisotropic
):
"""Random normalized anisotropic forcing.
.. inheritance-diagram:: TimeCorrelatedRandomPseudoSpectralAnisotropic3D
"""
tag = "tcrandom_anisotropic"
@classmethod
def _complete_params_with_default(cls, params):
super()._complete_params_with_default(params)
params.forcing.tcrandom_anisotropic.delta_angle = "10°"