diff --git a/docs/make.jl b/docs/make.jl index ad513590..74f2a65f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -94,10 +94,11 @@ sitename = "GeophysicalFlows.jl", "modules/multilayerqg.md", "modules/surfaceqg.md" ], - "Forcing" => "forcing.md", + "Stochastic Forcing" => "stochastic_forcing.md", "DocStrings" => Any[ "man/types.md", - "man/functions.md"] + "man/functions.md" + ] ] ) diff --git a/docs/src/assets/README.md b/docs/src/assets/README.md new file mode 100644 index 00000000..8321fac6 --- /dev/null +++ b/docs/src/assets/README.md @@ -0,0 +1 @@ +### a placeholder directory for output generated by Docs \ No newline at end of file diff --git a/docs/src/assets/energy_budgets_Ito.png b/docs/src/assets/energy_budgets_Ito.png deleted file mode 100644 index d56f5e90..00000000 Binary files a/docs/src/assets/energy_budgets_Ito.png and /dev/null differ diff --git a/docs/src/assets/energy_budgets_SPDE_Stratonovich.png b/docs/src/assets/energy_budgets_SPDE_Stratonovich.png deleted file mode 100644 index e3d7621a..00000000 Binary files a/docs/src/assets/energy_budgets_SPDE_Stratonovich.png and /dev/null differ diff --git a/docs/src/assets/energy_budgets_Stratonovich.png b/docs/src/assets/energy_budgets_Stratonovich.png deleted file mode 100644 index ac7ebbaf..00000000 Binary files a/docs/src/assets/energy_budgets_Stratonovich.png and /dev/null differ diff --git a/docs/src/assets/energy_comparison.png b/docs/src/assets/energy_comparison.png deleted file mode 100644 index e4db1d92..00000000 Binary files a/docs/src/assets/energy_comparison.png and /dev/null differ diff --git a/docs/src/forcing.md b/docs/src/forcing.md deleted file mode 100644 index ef52614f..00000000 --- a/docs/src/forcing.md +++ /dev/null @@ -1,359 +0,0 @@ -# Forcing - -Forcing of the equations is implemented in various modules. Forcing can be either -deterministic or stochastic (random). For deterministic forcing the implementation is -straightforward; for stochastic forcing there are two main train of thoughts: -Itô calculus and Stratonovich calculus. - -Both stochastic calculi give the same results. But once we decide to use one of -the two calculi we have to remain consistent and use that calculus for everywhere. -There is a lot of confusion and mostly the confusion stems from not using the same -stochastic calculus consistently throughout the computation but rather interchanging -between the two. - -`FourierFlows.jl` uses Stratonovich calculus throughout. This choise was made -because Stratonovich calculus works the same with both stochastic and deterministic -forcing, i.e. with Stratonovich calculus we have the same chain rules for -differentiation for stochastic functions as the chain rules we learn in -normal-deterministic calculus). Therefore, the code written as is does not really -"care" of what forcing the user implements. - -If you are interested in learning more regarding the two stochastic calculi and -how they are numerically implemented then read on; otherwise you can skip this -section of the documentation. - -## Stochastic Differential Equations (SDEs) - -A differential equation in the form: - -```math - \frac{\mathrm{d} x}{\mathrm{d} t} = f(x),\quad x(t_0)=0, -``` - -can also be written in an integral form: - -```math - x(t) = \int_{t_0}^{t} f(x(s))\,\mathrm{d} s. -``` - -In a similar manner, a stochastic differential equation - -```math - \mathrm{d} x = f(x)\,\mathrm{d} t + g(x)\,\mathrm{d} W_t,\quad x(t_0)=0, -``` - -with $\mathrm{d} W_t$ a white-noise process, can be written in an integral form as: - -```math - x(t) = \int_{t_0}^{t} f(x(s))\,\mathrm{d} s + \int_{t_0}^{t} g(x(s))\,\mathrm{d} W_s . -``` - -Of course now, the last integral is a stochastic integral and there is not a single -straight-forward way of computing it --- there are a lot of different ways we can -approximate it as a Riemannian sum and each of them leads to a different answer. -The two most popular ways for computing such stochastic integrals are: - -```math -{\color{Green}\text{Itô}: \int_{t_0}^{t} g(x(s))\,\mathrm{d} W_s\approx\sum_{j} g\left(x(t_j)\right)(W_{j+1}-W_j)},\\ -{\color{Magenta}\text{Stratonovich}: \int_{t_0}^{t} g(x(s))\,\mathrm{d} W_s \approx \sum_{j} g\left(x\left(\tfrac{1}{2}(t_j+t_{j+1})\right)\right)(W_{j+1}-W_j)}. -``` - -Because the white noise process is not continuous the two definitions do not converge to the same result; the two -definitions give thoroughly different results. And to overcome that they come along with different chain rules, -i.e., chain rules that are not necessarily the same as those in plain old calculus. - -An SDE can be written also in differential form. Because we cannot formally form $\mathrm{d} W/\mathrm{d} t$, since $W$ is -nowhere differentiable, we write an SDE in differential form as: - -```math -{\color{Green}\text{Itô}: \mathrm{d} x_t = f(x_t)\mathrm{d} t + g(x_t)\mathrm{d} W_t},\\ -{\color{Magenta}\text{Stratonovich}: \mathrm{d} x_t = f(x_t) \mathrm{d} t + g(x_t) \circ \mathrm{d} W_t}. -``` - -The circle in $g(x_t)\circ\mathrm{d} W_t$ is used to differentiate between Itô or Stratonovich calculus. - -A variable change $y=G(x)$ is done as follows according to the two different calculi: - -```math -{\color{Green}\text{Itô}: \mathrm{d} y_t = \frac{\mathrm{d} G}{\mathrm{d} x}\mathrm{d} x_t + \tfrac{1}{2} g(x_t)^2 \frac{\mathrm{d}^2 G}{\mathrm{d} x^2}\mathrm{d} t =\left[ \frac{\mathrm{d} G}{\mathrm{d} x}f(x_t) + \tfrac{1}{2} g(x_t)^2 \frac{\mathrm{d}^2 G}{\mathrm{d} x^2}\right]\mathrm{d} t + \frac{\mathrm{d} G}{\mathrm{d} x}g(x_t)\mathrm{d} W_t},\\ -{\color{Magenta}\text{Stratonovich}: \mathrm{d} y_t = \frac{\mathrm{d} G}{\mathrm{d} x}\mathrm{d} x_t =\frac{\mathrm{d} G}{\mathrm{d} x} f(x_t) \mathrm{d} t + \frac{\mathrm{d} G}{\mathrm{d} x}g(x_t)\circ\mathrm{d} W_t}. -``` - -The above are the so called *stochastic chain rules*. All derivatives of $G$ are evaluated at $x_t$. - -It's easy to see that the extra drift-term in Itô's interpretation of the stochastic integral, -i.e., ${\color{Green}\tfrac{1}{2} g^2\, \mathrm{d}^2G/\mathrm{d} x^2}$ is *exactly* equal to the ensemble mean of the -Stratonovich stochastic integral. That's because the Itô stochastic integral has, by construction, -zero ensemble mean since at every instant the noise is multiplied with $g$ evaluated before the action of the -noise; $g$ and $\mathrm{d} W$ are uncorrelated and thus: - -```math -{\color{Green}\left\langle g(x_t)\mathrm{d} W_t \right\rangle =0}\quad\text{while}\quad {\color{Magenta}\left\langle g(x_t)\circ\mathrm{d} W_t \right\rangle \ne 0}. -``` - -The above is demonstrated by evaluating the simple stochastic integral: - -```math -{\color{Green}\text{Itô}: \left\langle \int_{t_0}^{t} W_s\,\mathrm{d} W_s \right\rangle \approx\sum_{j} \left\langle W_j(W_{j+1}-W_j)\right\rangle}\\ -{\color{Green}\hspace{7.3em} = \sum_{j} \left\langle W_j W_{j+1}\right\rangle - \left\langle W_jW_j\right\rangle \sim \sum_{j} t_j - t_j = 0} ,\\ -{\color{Magenta}\text{Stratonovich}: \left\langle\int_{t_0}^{t} W_s\circ\mathrm{d} W_s\right\rangle \approx \sum_{j} \left\langle \frac1{2}(W_{j} + W_{j+1}) (W_{j+1}-W_j)\right\rangle }\\ -{\color{Magenta}\hspace{7.3em} = \frac1{2}\sum_{j} \left\langle W_{j+1} W_{j+1}\right\rangle - \left\langle W_{j} W_{j}\right\rangle \sim \frac1{2}\sum_{j} t_{j+1} - t_j = \frac{t}{2}}. -``` - -SDEs rarely can be solved in closed form; most often numerical solution of SDEs is brought to the rescue. -Itô calculus has the advantage that is very easily implemented numerically. -On the other hand, the chain rule in Stratonovich calculus coincides with that in normal calculus. This stems from the fact -that in the Stratonovich interpretation the white noise process is as a series of colored noise processes with the de-correlation time tending to zero. This made Stratonovich calculus more popular in the physics community. -A nice discussion on the differences and similarities between the two calculi is given by [van Kampen](https://doi.org/10.1007/BF01007642). - -## A simple Stochastic Differential Equation (SDE): the Ornstein--Uhlenbeck process - -One of the simpler SDEs is the Ornstein--Uhlenbeck process. A variation of which is: - -```math -x(t) = \int_{t_0}^{t} -\mu x(s)\,\mathrm{d} s + \int_{t_0}^{t} \sqrt{\sigma}\,\mathrm{d} W_s . \tag{eq:OU} -``` - -Note that in differential form this is: - -```math -\mathrm{d} x_t = -\mu x_t \,\mathrm{d} t + \sqrt{\sigma}\,\mathrm{d} W_s . \tag{eq:1} -``` - -Luckily, here there is no need to distinguish between Itô and Stratonovich. This is because $g$ is independent of $x(t)$. But we stress that this is often not the case; it is only a fortuitous coincident here. - -How do we time-step this SDE numerically? Let us assume a discretization of time into time-steps -of $\tau$: $t_j=(j-1)\tau$. (What follows can be easily carried on for non-uniform time discretization.) -With that, we denote $x_j\equiv x(t_j)$. Then the Euler--Mayorama time-step scheme for \eqref{eq:1} is - -```math - x_{j+1} = x_j + (-\mu x_j)\tau + \sqrt{\sigma}(W_{j+1}-W_j). -``` - -Now let us ask the following question: How can we compute the work done by the noise? -In other words, if we are interested in the evolution of the "energy" $E\equiv \tfrac{1}{2} x^2$, how is the noise term -attributing in the growth of $E$? To answer that we first have to find the SDE that energy $E$ obeys. -But, in doing so, it is important to adopt a single interpretation for computing stochastic integrals as now a transformation of variables is needed. That is, depending on whether we choose to interpret the stochastic integrals according to Itô or to Stratonovich calculus, $E$ evolves as: - -```math -{\color{Green} \text{Itô}: \mathrm{d} E_t = \left ( -2 \mu E_t + \tfrac{1}{2} \sigma \right ) \mathrm{d} t + x_t \sqrt{\sigma} \mathrm{d} W_t } , \tag{eq:Eito} -``` - -```math -{\color{Magenta} \text{Stratonovich}: \mathrm{d} E_t = -2 \mu E_t \mathrm{d} t + x_t \circ \sqrt{\sigma} \mathrm{d} W_t} . \tag{eq:Estr} -``` - -How do we compute the work $P$ done by the noise? It is: - -```math -{\color{Green} \text{Itô}: P_t = \tfrac{1}{2} \sigma \mathrm{d} t + \sqrt{\sigma} x_t \mathrm{d} W_t \approx \tfrac{1}{2} \sigma + \sqrt{\sigma} x_j (W_{j+1} - W_j),}\\ -{\color{Magenta} \text{Stratonovich}: P_t = x_t \circ \sqrt{\sigma} \mathrm{d} W_t \approx \sqrt{\sigma} x \left ( \tfrac{1}{2} (t_j + t_{j+1}) \right ) (W_{j+1}-W_j)}. -``` - -Say we didn't know the rules for transforming Stratonovich to Itô and we were wondering what is the extra drift term -we have to include in the Itô formulations, i.e. the $\tfrac{1}{2}\sigma$ term. We can compute the Itô's drift-term using -that it is exactly equal to $\langle x_t \circ \sqrt{\sigma} \mathrm{d} W_t \rangle$; and for the latter we can use the "usual" calculus. -That is, rewrite \eqref{eq:OU} as: - -```math -\dot{x} = -\mu x + \xi, \tag{eq:OUcont} -``` - -where $\xi(t)$ is understood to be the "continuous" version of the white-noise process which is formally only -understood in terms of distributions. The forcing $\xi$ has the properties: - -```math -\left \langle \xi(t) \right \rangle = 0 \quad \text{and} \quad \left \langle \xi(t) \xi(t') \right \rangle = \sigma \delta(t-t'). -``` - -Thus we need to compute $\langle P_t \rangle = \langle x(t) \xi(t) \rangle$. But \eqref{eq:OUcont} has formally the solution: - -```math -x(t) = \mathrm{e}^{-\mu t} x(0) + \int_0^t \mathrm{e}^{-\mu(t-s)} \xi(s) \, \mathrm{d} s . -``` - -and utilizing the above we get - -```math -\langle P_t \rangle = \langle x(t) \xi(t) \rangle -= \mathrm{e}^{-\mu t} \underbrace{\langle x(0) \xi(t) \rangle}_{=0} + \int_0^t \mathrm{e}^{-\mu(t-s)} \langle \xi(t)\xi(s) \rangle \, \mathrm{d} s -= \sigma \int_0^t \mathrm{e}^{-\mu(t-s)} \delta(t-s) \, \mathrm{d} s = \frac{\sigma}{2} . -``` - -Above we used that $\int_0^t\delta(t-s)\mathrm{d} s = \tfrac{1}{2}$, which is consistent with Stratonovich symmetric interpretation -of stochastic integrals. - -### Numerical implementation - -How do we time-step \eqref{eq:Eito}? We use the Euler--Maruyama time-stepping scheme: - -```math - E_{j+1} = E_j + \left(-2\mu E_j + \frac{\sigma}{2}\right)\tau + \sqrt{\sigma}x_j(W_{j+1}-W_j). -``` -However, we cannot use Euler--Maruyama for time-stepping \eqref{eq:Estr} since the Euler--Maruyama is "Itô"-thinking. -To time-step \eqref{eq:Estr} we have to approximate $g$ in the middle of the time-step. There are many ways to do -that, one of which is the, so called, Euler--Heun method: - -```math - \widetilde{E}_{j+1} = E_j + (-2\mu E_j)\tau + \sqrt{\sigma}x_j(W_{j+1}-W_j),\\ - E_{j+1} = E_j + \left(-2\mu \frac{E_j+\widetilde{E}_{j+1}}{2}\right)\tau + \sqrt{\sigma}\frac{x_j+x_{j+1}}{2}(W_{j+1}-W_j). -``` - -![energy_comparison](assets/energy_comparison.png) - -Figure above shows a comparison of the energy evolution as done from: -- direct computation as $\tfrac{1}{2} x_t^2$, -- time-integration of \ref{eq:Eito}, and -- time-integration of \eqref{eq:Estr}. - -Figures below show the ensemble mean energy budgets (using 1000 ensemble members) as computed using Itô and -Stratonovich. For the energy budget to close we have to be consistent: if we time-step the energy equation based -on Stratonovich calculus then we must compute the work also according to Stratonovich. -(For more details see `examples/forcing/simpleSDEItoStratonovich.jl`. - -![energy_budgets_Ito](assets/energy_budgets_Ito.png) -![energy_budgets_Stratonovich](assets/energy_budgets_Stratonovich.png) - -## A simple Stochastic Partial Differential Equation (SPDE) - -We want now to transfer all the knowledge we got from the previous sections to PDEs. In particular we'll focus on the simple SPDE: - -```math -\partial_t \nabla^2\psi(\boldsymbol{x}, t) = -\mu \nabla^2\psi(\boldsymbol{x}, t) + \xi(\boldsymbol{x},t),\tag{eq:PDEcont} -``` - -which is also equivalently written as: - -```math -\mathrm{d} \nabla^2\psi_{t}(\boldsymbol{x}) = -\mu \nabla^2 \psi_{t} (\boldsymbol{x}) \mathrm{d} t + \sqrt{\sigma} \mathrm{d} W_{t} (\boldsymbol{x}) . -``` - -The form \eqref{eq:PDEcont} is the continuous version understood in the Stratonovich interpretation -(similar to \eqref{eq:OUcont}). Thus, forcing $\xi$ obeys now: - -```math -\langle \xi(\boldsymbol{x},t)\rangle = 0 \quad\text{and}\quad \langle \xi(\boldsymbol{x},t)\xi(\boldsymbol{x}',t') \rangle= Q(\boldsymbol{x}-\boldsymbol{x}')\delta(t-t'), -``` - -that is the forcing is white in time but spatially correlated; its spatial correlation is prescribed by the -function $Q$ which is, necessarily, homogeneous in all its arguments -(see discussion by [Constantinou](http://arxiv.org/abs/1503.07644); Appendix A). - -The above describes the vorticity evolution of a two-dimensional fluid $\nabla^2\psi$ which is stochastically -forced while dissipated by linear drag $\mu$. The energy of the fluid is: - -```math -E = \tfrac{1}{2}\overline{|\boldsymbol{\nabla}\psi|^2}^{x,y} = -\tfrac{1}{2}\overline{\psi\nabla^2\psi}^{x,y}, -``` - -where the overbar denotes average over $x$ and $y$. To obtain the energy equation we multiply -\eqref{eq:PDEcont} with $-\psi$ and average over the whole domain. Thus, the work done by the forcing is given by the term: - -```math -P = -\,\overline{\psi\,\xi}^{x,y}, -``` - -but the above is a stochastic integral and it is meaningless without a rule for computing the stochastic integral. - -Numerically, the work done by the forcing can be obtained Stratonovich-wise as: - -```math -\begin{aligned} -P_j = -\,\overline{\frac{\psi(\boldsymbol{x},t_j)+\psi(\boldsymbol{x},t_{j+1})}{2} \xi(\boldsymbol{x},t_{j+1}) }^{x,y}, -\end{aligned} -``` - -or Itô-wise - -```math -\begin{aligned} -P_j = -\,\overline{ \psi(\boldsymbol{x},t_j) \xi(\boldsymbol{x},t_{j+1}) }^{x,y} + \text{drift}, -\end{aligned} -``` - -But how much is the Itô drift term in this case? As in the previous section, the drift is *precisely* the ensemble -mean of the Stratonovich work, i.e.: - -```math -\textrm{Ito drift}= -\overline{\langle \underbrace{\psi(\boldsymbol{x},t)\circ \xi(\boldsymbol{x},t)}_{\textrm{Stratonovich}} \rangle}^{x,y}, -``` - -But again the above can be computed relatively easy if we use the "formal" solution of \eqref{eq:PDEcont}: - -```math -\psi(\boldsymbol{x},t) = \mathrm{e}^{-\mu t}\psi(\boldsymbol{x},0) + \int_0^t \mathrm{e}^{-\mu(t-s)}\nabla^{-2}\xi(\boldsymbol{x},s)\,\mathrm{d} s, -``` - -which implies - -```math -\text{drift} = -\overline{\mathrm{e}^{-\mu t} \underbrace{\left \langle \psi(\boldsymbol{x}, 0) \xi(\boldsymbol{x}, t) \right \rangle}_{=0}}^{x,y} - \int_0^t \mathrm{e}^{-\mu(t-s)} \overline{\nabla^{-2} \left \langle \xi(\boldsymbol{x}, s) \xi(\boldsymbol{x}, t) \right\rangle}^{x,y}\,\mathrm{d} s \\ -= -\int_0^t \mathrm{e}^{-\mu(t-s)} \overline{\underbrace{\left [ \nabla^{-2} Q (\boldsymbol{x}) \right ] \big|_{\boldsymbol{x}=0}}_{\text{independent of }x,y} \, \delta(t-s)}^{x,y} \, \mathrm{d} s = -\frac1{2} \nabla^{-2} Q(\boldsymbol{x}) \big|_{\boldsymbol{x}=0} \\ -= - \frac1{2} \left [ \nabla^{-2} \int \frac{\mathrm{d}^2 \boldsymbol{k}}{(2\pi)^2} \widehat{Q}(\boldsymbol{k}) \, \mathrm{e}^{\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{x}} \right] _{\boldsymbol{x}=0} -= \int \frac{\mathrm{d}^2 \boldsymbol{k}}{(2\pi)^2} \frac{\widehat{Q}(\boldsymbol{k})}{2k^2}. -``` - -Thus, the drift, or in this case the mean energy input rate by the stochastic forcing, is precisely determined -from the spatial correlation of the forcing. Let us denote: - -```math -\varepsilon \equiv \int \frac{\mathrm{d}^2 \boldsymbol{k}}{(2\pi)^2} \frac{\widehat{Q}(\boldsymbol{k})}{2k^2}. \tag{eq:defepsilon} -``` - -Therefore, work for a single forcing realization is computed numerically as: - -```math -\begin{aligned} -{\color{Green} \text{Itô}} &: {\color{Green} P_j = -\overline{\psi(\boldsymbol{x}, t_j) \xi(\boldsymbol{x}, t_{j+1})}^{x,y} + \varepsilon},\\ -{\color{Magenta} \text{Stratonovich}} &: {\color{Magenta} P_j = -\overline{\frac{\psi(\boldsymbol{x}, t_j)+\psi(\boldsymbol{x}, t_{j+1})}{2} \xi(\boldsymbol{x}, t_{j+1})}^{x,y}}. \tag{eq:PtStrat} -\end{aligned} -``` - -Remember, previously the work done by the stochastic forcing was: -```math -\mathrm{d} P_t = {\color{Green} \frac{\sigma}{2}\mathrm{d} t + \sqrt{\sigma} x_t \mathrm{d} W_t} = {\color{Magenta} \sqrt{\sigma} x_t \circ \mathrm{d} W_t}, -``` -and by sampling over various forcing realizations: -```math -\langle \mathrm{d} P_t \rangle = \frac{\sigma}{2} \mathrm{d} t = \langle \sqrt{\sigma} x_t \circ \mathrm{d} W_t \rangle -``` - -The code uses Stratonovich. For example, the work done by the forcing in the `TwoDTurb` module is computed based -on \eqref{eq:PtStrat} with the function - -```julia -@inline function work(s, v::ForcedVars, g) - @. v.Uh = g.invKKrsq * (v.prevsol + s.sol)/2.0 * conj(v.Fh) - 1/(g.Lx*g.Ly)*FourierFlows.parsevalsum(v.Uh, g) -end -``` - -## A less-simple SPDE - -It turns out that nothing changes if we include the nonlinear terms in the vorticity equation: - -```math -\partial_t \nabla^2 \psi(\boldsymbol{x}, t) + \mathsf{J}(\psi,\nabla^2\psi) = -\mu \nabla^2\psi(\boldsymbol{x}, t) + \xi(\boldsymbol{x},t). \tag{eq:PDEcont2} -``` - -The nonlinearity does not alter the Itô drift; thus the ensemble mean energy input by the stochastic forcing, -remains the same. We can easily verify this from the "formal" solution of \eqref{eq:PDEcont2}: - -```math -\psi(\boldsymbol{x}, t) = \mathrm{e}^{-\mu t} \psi(\boldsymbol{x}, 0) + \int_0^t \mathrm{e}^{-\mu(t-s)} \nabla^{-2} \xi(\boldsymbol{x}, s) \, \mathrm{d} s - \int_0^t \nabla^{-2} \mathsf{J} \left ( \psi(\boldsymbol{x}, s), \nabla^2 \psi(\boldsymbol{x}, s) \right ) \, \mathrm{d} s, -``` - -When multiplied with $\xi(\boldsymbol{x}, t)$ the last term vanishes since its only non-zero contribution comes from the point -$s=t$ which is of measure zero (in the integrated sense). - -Figure below shows the energy budgets for a numerical solution of \eqref{eq:PDEcont2} starting from rest -($\psi(\boldsymbol{x}, 0)=0$) in a doubly periodic square domain of size $L$ (`examples/twodturb/IsotropicRingForcing.jl`). -The forcing was prescribed to have power in a narrow ring in wavenumber space: - -```math -\widehat{Q}(\boldsymbol{k}) \propto \mathrm{e}^{-(|\boldsymbol{k}| - k_f)^2 / (2\delta_f^2)}, -``` - -with $k_f L / (2\pi) = 12$ and $\delta_f L / (2\pi) = 2$. The mean energy input rate was set to $\varepsilon = 0.1$. - -![energy_budgets_SPDE_Stratonovich](assets/energy_budgets_SPDE_Stratonovich.png) diff --git a/docs/src/modules/barotropicqgql.md b/docs/src/modules/barotropicqgql.md index ec09c99e..1f7d9a61 100644 --- a/docs/src/modules/barotropicqgql.md +++ b/docs/src/modules/barotropicqgql.md @@ -2,34 +2,49 @@ ### Basic Equations -This module solves the *quasi-linear* quasi-geostrophic barotropic vorticity equation on a beta-plane of variable fluid depth $H-h(x,y)$. -Quasi-linear refers to the dynamics that *neglect* the eddy--eddy interactions in the eddy evolution equation after an eddy--mean flow decomposition, e.g., +This module solves the *quasi-linear* quasi-geostrophic barotropic vorticity equation on a beta +plane of variable fluid depth ``H - h(x, y)``. Quasi-linear refers to the dynamics that *neglect* +the eddy--eddy interactions in the eddy evolution equation after an eddy--mean flow decomposition, e.g., -$$\phi(x, y, t) = \overline{\phi}(y, t) + \phi'(x,y,t) ,$$ +```math +\phi(x, y, t) = \overline{\phi}(y, t) + \phi'(x, y, t) , +``` -where overline above denotes a zonal mean, $\overline{\phi}(y, t) = \int \phi(x, y, t)\,\mathrm{d}x/L_x$, and prime denotes deviations from the zonal mean. This approximation is used in many process-model studies of zonation, e.g., +where overline above denotes a zonal mean, ``\overline{\phi}(y, t) = \int \phi(x, y, t) \, 𝖽x / L_x``, and prime denotes deviations from the zonal mean. This approximation is used in many process-model studies of zonation, e.g., - Farrell, B. F. and Ioannou, P. J. (2003). [Structural stability of turbulent jets.](http://doi.org/10.1175/1520-0469(2003)060<2101:SSOTJ>2.0.CO;2) *J. Atmos. Sci.*, **60**, 2101-2118. - Srinivasan, K. and Young, W. R. (2012). [Zonostrophic instability.](http://doi.org/10.1175/JAS-D-11-0200.1) *J. Atmos. Sci.*, **69 (5)**, 1633-1656. - Constantinou, N. C., Farrell, B. F., and Ioannou, P. J. (2014). [Emergence and equilibration of jets in beta-plane turbulence: applications of Stochastic Structural Stability Theory.](http://doi.org/10.1175/JAS-D-13-076.1) *J. Atmos. Sci.*, **71 (5)**, 1818-1842. - -As in the [SingleLayerQG module](singlelayerqg.md), the flow is obtained through a streamfunction $\psi$ as $(u, v) = (-\partial_y\psi, \partial_x\psi)$. All flow fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is - -$$\underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{(\partial_x v - - \partial_y u)}_{\text{relative vorticity}} + - \underbrace{\frac{f_0 h}{H}}_{\text{topographic PV}}.$$ - -The dynamical variable is the component of the vorticity of the flow normal to the plane of motion, $\zeta\equiv \partial_x v- \partial_y u = \nabla^2\psi$. Also, we denote the topographic PV with $\eta\equiv f_0 h/H$. After we apply the eddy-mean flow decomposition above, the QGPV dynamics are: - -$$\partial_t \overline{\zeta} + \mathsf{J}(\overline{\psi}, \underbrace{\overline{\zeta} + \overline{\eta}}_{\equiv \overline{q}}) + \overline{\mathsf{J}(\psi', \underbrace{\zeta' + \eta'}_{\equiv q'})} = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} -\right] \overline{\zeta} }_{\textrm{dissipation}} \ .$$ - -$$\partial_t \zeta' + \mathsf{J}(\psi', \overline{q}) + \mathsf{J}(\overline{\psi}, q') + \underbrace{\mathsf{J}(\psi', q') - \overline{\mathsf{J}(\psi', q')}}_{\textrm{EENL}} + -\beta\partial_x\psi' = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} -\right] \zeta'}_{\textrm{dissipation}} + f\ .$$ - -where $\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)$. On the right hand side, $f(x,y,t)$ is forcing (which is assumed to have zero mean, $\overline{f}=0$), $\mu$ is linear drag, and $\nu$ is hyperviscosity. Plain old viscosity corresponds to $n_{\nu}=1$. The sum of relative vorticity and topographic PV is denoted with $q\equiv\zeta+\eta$. +As in the [SingleLayerQG module](singlelayerqg.md), the flow is obtained through a +streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow fields +can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here, the QGPV is + +```math +\underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{\partial_x v + - \partial_y u}_{\text{relative vorticity}} + \underbrace{\frac{f_0 h}{H}}_{\text{topographic PV}} . +``` + +The dynamical variable is the component of the vorticity of the flow normal to the plane of +motion, ``\zeta \equiv \partial_x v - \partial_y u = \nabla^2 \psi``. Also, we denote the +topographic PV with ``\eta \equiv f_0 h/H``. After we apply the eddy-mean flow decomposition +above, the QGPV dynamics are: + +```math +\begin{aligned} +\partial_t \overline{\zeta} & + \mathsf{J}(\overline{\psi}, \underbrace{\overline{\zeta} + \overline{\eta}}_{\equiv \overline{q}}) + \overline{\mathsf{J}(\psi', \underbrace{\zeta' + \eta'}_{\equiv q'})} = \underbrace{- \left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} +\right] \overline{\zeta} }_{\textrm{dissipation}} ,\\ +\partial_t \zeta' &+ \mathsf{J}(\psi', \overline{q}) + \mathsf{J}(\overline{\psi}, q') + \underbrace{\mathsf{J}(\psi', q') - \overline{\mathsf{J}(\psi', q')}}_{\textrm{EENL}} + +\beta \partial_x \psi' = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} +\right] \zeta'}_{\textrm{dissipation}} + F . +\end{aligned} +``` + +where ``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b) - (\partial_y a)(\partial_x b)``. On +the right hand side, ``F(x, y, t)`` is forcing (which is assumed to have zero zonal mean, +``\overline{F} = 0``), ``\mu`` is linear drag, and ``\nu`` is hyperviscosity. Plain old +viscosity corresponds to ``n_{\nu} = 1``. The sum of relative vorticity and topographic PV is +denoted with ``q \equiv \zeta + \eta``. *Quasi-linear* dynamics **neglect the term eddy-eddy nonlinearity (EENL) term** above. @@ -37,19 +52,25 @@ where $\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x The equation is time-stepped forward in Fourier space: -$$\partial_t \widehat{\zeta} = - \widehat{\mathsf{J}(\psi, q)}^{\textrm{QL}} +\beta\frac{\mathrm{i}k_x}{k^2}\widehat{\zeta} -\left(\mu -+\nu k^{2n_\nu}\right) \widehat{\zeta} + \widehat{f}\ .$$ +```math +\partial_t \widehat{\zeta} = - \widehat{\mathsf{J}(\psi, q)}^{\textrm{QL}} + \beta \frac{i k_x}{|𝐤|^2} \widehat{\zeta} - \left ( \mu + \nu |𝐤|^{2n_\nu} \right ) \widehat{\zeta} + \widehat{F} . +``` -In doing so the Jacobian is computed in the conservative form: $\mathsf{J}(f,g) = -\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]$. The superscript QL in the Jacobian term above denotes that remove triad interactions that correspond to the EENL term. +In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = +\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]``. The superscript QL on the +Jacobian term above denotes that triad interactions that correspond to the EENL term are +removed. Thus: -$$\mathcal{L} = \beta\frac{\mathrm{i}k_x}{k^2} - \mu - \nu k^{2n_\nu}\ ,$$ -$$\mathcal{N}(\widehat{\zeta}) = - \mathrm{i}k_x \mathrm{FFT}(u q)^{\textrm{QL}}- - \mathrm{i}k_y \mathrm{FFT}(v q)^{\textrm{QL}}\ .$$ +```math +\begin{aligned} +L & = \beta \frac{i k_x}{|𝐤|^2} - \mu - \nu |𝐤|^{2n_\nu} ,\\ +N(\widehat{\zeta}) & = - i k_x \mathrm{FFT}(u q)^{\textrm{QL}} - i k_y \mathrm{FFT}(v q)^{\textrm{QL}} + \widehat{F}. +\end{aligned} +``` ## Examples -- `examples/barotropicqgql_betaforced.jl`: A script that simulates forced-dissipative quasi-linear quasi-geostrophic flow on a beta-plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius `kf` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation. +- `examples/barotropicqgql_betaforced.jl`: A script that simulates forced-dissipative quasi-linear quasi-geostrophic flow on a beta-plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius ``k_f`` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation. diff --git a/docs/src/modules/multilayerqg.md b/docs/src/modules/multilayerqg.md index df612b31..1f339af3 100644 --- a/docs/src/modules/multilayerqg.md +++ b/docs/src/modules/multilayerqg.md @@ -10,10 +10,10 @@ is the number of fluid layers. The QGPV in each layer is ```math -\mathrm{QGPV}_j = q_j + \underbrace{f_0 + \beta y}_{\textrm{planetary PV}} + \delta_{j,n} \underbrace{\frac{f_0 h}{H_n}}_{\textrm{topographic PV}}, \quad j = 1, \dots, n \ . +\mathrm{QGPV}_j = q_j + \underbrace{f_0 + \beta y}_{\textrm{planetary PV}} + \delta_{j, n} \underbrace{\frac{f_0 h}{H_n}}_{\textrm{topographic PV}}, \quad j = 1, \dots, n \ . ``` -where ``q_j`` incorporates the relative vorticity in each layer ``\nabla^2\psi_j`` and the +where ``q_j`` incorporates the relative vorticity in each layer ``\nabla^2 \psi_j`` and the vortex stretching terms: ```math @@ -33,9 +33,9 @@ In view of the relationships above, when we convert to Fourier space ``q``'s and related via the matrix equation ```math -\begin{pmatrix} \widehat{q}_{\boldsymbol{k}, 1}\\\vdots\\\widehat{q}_{\boldsymbol{k}, n} \end{pmatrix} = -\underbrace{\left(-|\boldsymbol{k}|^2 \mathbb{1} + \mathbb{F} \right)}_{\equiv \mathbb{S}_{\boldsymbol{k}}} -\begin{pmatrix} \widehat{\psi}_{\boldsymbol{k}, 1}\\\vdots\\\widehat{\psi}_{\boldsymbol{k}, n} \end{pmatrix} +\begin{pmatrix} \widehat{q}_{𝐤, 1}\\\vdots\\\widehat{q}_{𝐤, n} \end{pmatrix} = +\underbrace{\left(-|𝐤|^2 \mathbb{1} + \mathbb{F} \right)}_{\equiv \mathbb{S}_{𝐤}} +\begin{pmatrix} \widehat{\psi}_{𝐤, 1}\\\vdots\\\widehat{\psi}_{𝐤, n} \end{pmatrix} ``` where @@ -53,16 +53,26 @@ where Including an imposed zonal flow ``U_j(y)`` in each layer, the equations of motion are: ```math -\partial_t q_j + \mathsf{J}(\psi_j, q_j ) + (U_j - \partial_y\psi_j) \partial_x Q_j + U_j \partial_x q_j + (\partial_y Q_j)(\partial_x \psi_j) = -\delta_{j, n} \mu \nabla^2 \psi_n - \nu (-1)^{n_\nu} \nabla^{2n_\nu} q_j, +\partial_t q_j + \mathsf{J}(\psi_j, q_j ) + (U_j - \partial_y\psi_j) \partial_x Q_j + U_j \partial_x q_j + (\partial_y Q_j)(\partial_x \psi_j) = -\delta_{j, n} \mu \nabla^2 \psi_n - \nu (-1)^{n_\nu} \nabla^{2 n_\nu} q_j , ``` with ```math \partial_y Q_j \equiv \beta - \partial_y^2 U_j - (1-\delta_{j,1}) F_{j-1/2, j} (U_{j-1} - U_j) - (1 - \delta_{j,n}) F_{j+1/2, j} (U_{j+1} - U_j) + \delta_{j,n} \partial_y \eta \ , \\ -\partial_x Q_j \equiv \delta_{j,n} \partial_x \eta \ . +\partial_x Q_j \equiv \delta_{j, n} \partial_x \eta \ . ``` +### Helper functions + +```@docs +GeophysicalFlows.MultiLayerQG.set_q! +GeophysicalFlows.MultiLayerQG.set_ψ! +GeophysicalFlows.MultiLayerQG.updatevars! +``` + +### Diagnostics + The eddy kinetic energy in each layer and the eddy potential energy that corresponds to each fluid interface is computed via `energies()`: @@ -79,37 +89,50 @@ GeophysicalFlows.MultiLayerQG.fluxes ### Implementation -Matrices ``\mathbb{S}_{\boldsymbol{k}}`` as well as ``\mathbb{S}^{-1}_{\boldsymbol{k}}`` are included -in `params` as `params.S` and `params.S⁻¹` respectively. Additionally, the background PV gradients -``\partial_x Q`` and ``\partial_y Q`` are also included in the `params` as `params.Qx` and `params.Qy`. - -You can get ``\widehat{\psi}_j`` from ``\widehat{q}_j`` with `streamfunctionfrompv!(psih, qh, params, grid)`, -while to get ``\widehat{q}_j`` from ``\widehat{\psi}_j`` you need to call `pvfromstreamfunction!(qh, psih, params, grid)`. +Matrices ``\mathbb{S}_{𝐤}`` as well as ``\mathbb{S}^{-1}_{𝐤}`` are included in `params` as +`params.S` and `params.S⁻¹` respectively. Additionally, the background PV gradients +``\partial_x Q`` and ``\partial_y Q`` are also included in the `params` as `params.Qx` and +`params.Qy`. +One can get the ``\widehat{\psi}_j`` from ``\widehat{q}_j`` via +`streamfunctionfrompv!(psih, qh, params, grid)`, while the inverse, i.e. obtain ``\widehat{q}_j`` from ``\widehat{\psi}_j``, is done via `pvfromstreamfunction!(qh, psih, params, grid)`. The equations of motion are time-stepped forward in Fourier space: ```math \partial_t \widehat{q}_j = - \widehat{\mathsf{J}(\psi_j, q_j)} - \widehat{U_j \partial_x Q_j} - \widehat{U_j \partial_x q_j} -+ \widehat{(\partial_y \psi_j) \partial_x Q_j} - \widehat{(\partial_x\psi_j)(\partial_y Q_j)} + \delta_{j,n} \mu k^{2} \widehat{\psi}_n - \nu k^{2n_\nu} \widehat{q}_j \ . ++ \widehat{(\partial_y \psi_j) \partial_x Q_j} - \widehat{(\partial_x \psi_j)(\partial_y Q_j)} + \delta_{j, n} \mu |𝐤|^{2} \widehat{\psi}_n - \nu |𝐤|^{2n_\nu} \widehat{q}_j \ . ``` In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = -\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]``. +\partial_y [ (\partial_x f) g] - \partial_x[ (\partial_y f) g]``. Equations are formulated using $q_j$ as the state variables, i.e., `sol = qh`. -Thus: +The linear operator is constructed in `Equation` -```math -\begin{aligned} -\mathcal{L} & = - \nu k^{2n_\nu} \ , \\ -\mathcal{N}(\widehat{q}_j) & = - \widehat{\mathsf{J}(\psi_j, q_j)} - \widehat{U_j \partial_x Q_j} - \widehat{U_j \partial_x q_j} - + \widehat{(\partial_y \psi_j)(\partial_x Q_j)} - \widehat{(\partial_x \psi_j)(\partial_y Q_j)} + \delta_{j,n} \mu k^{2} \widehat{\psi}_n\ . -\end{aligned} +```@docs +GeophysicalFlows.MultiLayerQG.Equation +GeophysicalFlows.MultiLayerQG.hyperviscosity ``` - +The nonlinear terms is computed via + +```@docs +GeophysicalFlows.MultiLayerQG.calcN! +``` +which, in turn, calls + +```@docs +GeophysicalFlows.MultiLayerQG.calcN_advection! +``` +and + +```@docs +GeophysicalFlows.MultiLayerQG.addforcing! +``` + + ## Examples - `examples/multilayerqg_2layer.jl`: A script that simulates the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model. diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index 35324155..1fd74819 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -3,8 +3,9 @@ ### Basic Equations This module solves the barotropic or equivalent barotropic quasi-geostrophic vorticity equation -on a beta-plane of variable fluid depth ``H - h(x, y)``. The flow is obtained through a streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow -fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is +on a beta-plane of variable fluid depth ``H - h(x, y)``. The flow is obtained through a +streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow fields +can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is ```math \underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{\partial_x v @@ -16,19 +17,19 @@ fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). He where ``\ell`` is the Rossby radius of deformation. Purely barotropic dynamics corresponds to infinite Rossby radius of deformation (``\ell = \infty``), while a flow with a finite Rossby radius follows is said to obey equivalent-barotropic dynamics. We denote the sum of the relative -vorticity and the vortex stretching contributions to the QGPV with ``q = \nabla^2 \psi - \psi / \ell^2``. +vorticity and the vortex stretching contributions to the QGPV with ``q \equiv \nabla^2 \psi - \psi / \ell^2``. Also, we denote the topographic PV with ``\eta \equiv f_0 h / H``. The dynamical variable is ``q``. Thus, the equation solved by the module is: ```math -\partial_t q + \mathsf{J}(\psi, q + \eta) + -\beta \partial_x \psi = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} -\right] q}_{\textrm{dissipation}} + f \ . +\partial_t q + \mathsf{J}(\psi, q + \eta) + \beta \partial_x \psi = +\underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F \ . ``` where ``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)``. On -the right hand side, ``f(x, y, t)`` is forcing, ``\mu`` is linear drag, and ``\nu`` is hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``. +the right hand side, ``F(x, y, t)`` is forcing, ``\mu`` is linear drag, and ``\nu`` is +hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``. The kinetic energy of the fluid is computed via: @@ -53,7 +54,7 @@ GeophysicalFlows.SingleLayerQG.energy The equation is time-stepped forward in Fourier space: ```math -\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q + \eta)} + \beta \frac{\mathrm{i} k_x}{k^2 + 1/\ell^2} \widehat{q} - \left(\mu + \nu k^{2n_\nu} \right) \widehat{q} + \widehat{f} \ . +\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q + \eta)} + \beta \frac{i k_x}{|𝐤|^2 + 1/\ell^2} \widehat{q} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} \ . ``` In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = @@ -63,8 +64,8 @@ Thus: ```math \begin{aligned} -\mathcal{L} & = \beta \frac{\mathrm{i} k_x}{k^2 + 1/\ell^2} - \mu - \nu k^{2n_\nu} \ , \\ -\mathcal{N}(\widehat{q}) & = - \mathrm{i} k_x \mathrm{FFT}[u (q+\eta)] - \mathrm{i} k_y \mathrm{FFT}[v (q+\eta)] \ . +L & = \beta \frac{i k_x}{|𝐤|^2 + 1/\ell^2} - \mu - \nu |𝐤|^{2n_\nu} \ , \\ +N(\widehat{q}) & = - i k_x \mathrm{FFT}[u (q + \eta)] - i k_y \mathrm{FFT}[v (q + \eta)] + \widehat{F} \ . \end{aligned} ``` diff --git a/docs/src/modules/surfaceqg.md b/docs/src/modules/surfaceqg.md index c7f95335..aec51cb9 100644 --- a/docs/src/modules/surfaceqg.md +++ b/docs/src/modules/surfaceqg.md @@ -3,8 +3,8 @@ ### Basic Equations This module solves the non-dimensional surface quasi-geostrophic (SQG) equation for surface -buoyancy $b_s = b(x, y, z=0)$, as described in Capet et al., 2008. The buoyancy and the fluid -velocity at the surface are related through a streamfunction $\psi$ via: +buoyancy ``b_s = b(x, y, z=0)``, as described in Capet et al., 2008. The buoyancy and the fluid +velocity at the surface are related through a streamfunction ``\psi`` via: ```math (u_s, v_s, b_s) = (-\partial_y \psi, \partial_x \psi, -\partial_z \psi) . @@ -13,42 +13,42 @@ velocity at the surface are related through a streamfunction $\psi$ via: The SQG model evolves the surface buoyancy, ```math -\partial_t b_s + \mathsf{J}(\psi, b_s) = \underbrace{-\nu(-1)^{n_\nu} \nabla^{2n_\nu} b_s}_{\textrm{buoyancy diffusion}} + \underbrace{f}_{\textrm{forcing}} . +\partial_t b_s + \mathsf{J}(\psi, b_s) = \underbrace{-\nu(-1)^{n_\nu} \nabla^{2n_\nu} b_s}_{\textrm{buoyancy diffusion}} + \underbrace{F}_{\textrm{forcing}} . ``` -The evolution of buoyancy is only solved for the surface layer, but $b_s$ is a function of the vertical gradient of $\psi$. In the SQG system, the potential vorticity in the interior of the flow is identically zero. That is, relative vorticity is identical and opposite to the vertical stretching of buoyancy layers, +The evolution of buoyancy is only solved for the surface layer, but ``b_s`` is a function of the vertical gradient of ``\psi``. In the SQG system, the potential vorticity in the interior of the flow is identically zero. That is, relative vorticity is identical and opposite to the vertical stretching of buoyancy layers, ```math \underbrace{\left(\partial_x^2 + \partial_y^2 \right) \psi}_{\textrm{relative vorticity}} + \underbrace{\partial_z^2 \psi}_{\textrm{stretching term}} = 0 , ``` -with the boundary conditions $b_s = -\partial_z\psi|_{z=0}$ and $\psi \rightarrow 0$ as $z \rightarrow -\infty$. (We take here the oceanographic convention: $z \le 0$.) +with the boundary conditions ``b_s = -\partial_z\psi|_{z=0}`` and ``\psi \rightarrow 0`` as ``z \rightarrow -\infty``. (We take here the oceanographic convention: ``z \le 0``.) -These equations describe a system where the streamfunction (and hence the dynamics) at all depths is prescribed entirely by the surface buoyancy. By taking the Fourier transform in the horizontal ($x$ and $y$), the streamfunction-buoyancy relation is: +These equations describe a system where the streamfunction (and hence the dynamics) at all depths is prescribed entirely by the surface buoyancy. By taking the Fourier transform in the horizontal (``x`` and ``y``), the streamfunction-buoyancy relation is: ```math -\widehat{\psi}(k_x, k_y, z, t) = -\frac{\widehat{b_s}}{k} \,\mathrm{e}^{kz}, +\widehat{\psi}(k_x, k_y, z, t) = - \frac{\widehat{b_s}}{|𝐤|} \, e^{|𝐤|z}, ``` -where $k = \sqrt{k_x^2 + k_y^2}$ is the total horizontal wavenumber. +where ``|𝐤| = \sqrt{k_x^2 + k_y^2}`` is the total horizontal wavenumber. ### Implementation The buoyancy equation is time-stepped forward in Fourier space: ```math -\partial_t \widehat{b_s} = \underbrace{- \widehat{\mathsf{J}(\psi, b_s)}}_{\textrm{nonlinear term}} - \underbrace{\nu k^{2n_\nu} \widehat{b_s}}_{\textrm{linear term}} + \widehat{f} . +\partial_t \widehat{b_s} = - \widehat{\mathsf{J}(\psi, b_s)} - \nu |𝐤|^{2 n_\nu} \widehat{b_s} + \widehat{f} . ``` -In doing so the Jacobian is computed in the conservative form: $\mathsf{J}(f,g) = -\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]$. +In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = +\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]``. Thus: ```math \begin{aligned} -\widehat{u} &= \frac{\mathrm{i} k_y}{k} \widehat{b_s}, \qquad \widehat{v} = -\frac{\mathrm{i} k_x}{k} \widehat{b_s}, \\ -\mathcal{L} & = - \nu k^{2n_\nu},\\ -\mathcal{N}(\widehat{b_s}) & = - \mathrm{i} k_x \mathrm{FFT}(u b) - \mathrm{i} k_y \mathrm{FFT}(v b) . +\widehat{u} &= \frac{i k_y}{|𝐤|} \widehat{b_s}, \qquad \widehat{v} = -\frac{i k_x}{|𝐤|} \widehat{b_s}, \\ +L & = - \nu |𝐤|^{2n_\nu},\\ +N(\widehat{b_s}) & = - i k_x \mathrm{FFT}(u b) - i k_y \mathrm{FFT}(v b) . \end{aligned} ``` diff --git a/docs/src/stochastic_forcing.md b/docs/src/stochastic_forcing.md new file mode 100644 index 00000000..eb0e4987 --- /dev/null +++ b/docs/src/stochastic_forcing.md @@ -0,0 +1,536 @@ +# Stochastic Forcing + +Forcing terms are implemented in various modules. Forcing can be either deterministic or +stochastic (random). For deterministic forcing the implementation is straightforward; for +stochastic forcing there are two main train of thoughts: Itô calculus and Stratonovich calculus. + +Both stochastic calculi give the same results. But once we decide to use one of the two calculi +we have to remain consistent and use that calculus throughout. There can be a lot of confusion +and oftentimes confusion stems from mixing the two different stochastic calculi in a single +computation instead of using one of the two consistently all along. + +!!! note "Itô or Stratonovich in GeophysicalFlows.jl?" + All modules included in GeophysicalFlows.jl use **Stratonovich calculus**. + +The choice of Stratonovich calculus for GeophysicalFlows.jl was made since this calculus "works +the same" with both stochastic and deterministic forcing, i.e. with Stratonovich calculus we +have the same chain rules for differentiation for stochastic functions as the chain rules we +learn in normal-deterministic calculus). Therefore, with the Stratonovich calculus the code +does not really "care" whether the user implement deterministic or stochastic forcing. + +If you are interested in learning more regarding the two stochastic calculi and how they are +numerically implemented then read on; otherwise you can skip this section of the documentation. + +## Stochastic Differential Equations (SDEs) + +A differential equation in the form: + +```math + \frac{\mathrm{d} x}{\mathrm{d} t} = f(x) , \quad x(t_0) = 0, +``` + +can also be written in an integral form: + +```math + x(t) = \int_{t_0}^t f(x(s)) \, \mathrm{d} s. +``` + +In a similar manner, a stochastic differential equation + +```math + \mathrm{d} x = f(x) \, \mathrm{d} t + g(x) \, \mathrm{d} W_t , \quad x(t_0) = 0 , +``` + +with ``\mathrm{d} W_t`` a white-noise process, can be written in an integral form as: + +```math + x(t) = \int_{t_0}^{t} f(x(s)) \, \mathrm{d} s + \int_{t_0}^{t} g(x(s)) \, \mathrm{d} W_s . +``` + +Of course now, the last integral is a stochastic integral and there is not a single +straight-forward way of computing it -- there are a lot of different ways we can +approximate it as a Riemannian sum and each of them leads to a different answer. +The two most popular ways for computing such stochastic integrals are: + +```math +\begin{aligned} +{\color{Green} \text{Itô}}&: {\color{Green}\int_{t_0}^{t} g(x(s)) \, \mathrm{d} W_s \approx \sum_{j} g \left ( x(t_j) \right )(W_{j+1} - W_j)} , \\ +{\color{Magenta} \text{Stratonovich}}&: {\color{Magenta} \int_{t_0}^{t} g(x(s)) \, \mathrm{d} W_s \approx \sum_{j} g \left (x \left (\tfrac{1}{2}(t_j + t_{j+1}) \right ) \right)(W_{j+1} - W_j)} . +\end{aligned} +``` + +The difference in the two calculi above lies in the point at which we choose to evaluate ``g(x)``: +we take the start of the time-interval for ${\color{Green} \text{Itô, } t_j}$, while we use +the mid-point for ${\color{Magenta}\text{Stratonovich, } (t_j+t_{j+1})/2}$. In the case the +stochastic noise is additive, i.e., its prefactor ``g`` does not depend on the state ``x_t``, +then the two interpretations coincide. When the noise does depend on the state of the system, +i.e., ``g=g(x_t)``, then the two interpretations above give thoroughly different results. This +happens because the white noise process is not continuous and, therefore, the two interpretations +of the stochastic integrals above do not converge to the same result. + +To overcome this apparent inconsistency, the two choices above come together with different +chain rules, i.e., chain rules that are not necessarily the same as those in plain old calculus. +Let us see how different choices for computing the stochastic integrals bring about the need f +or different chain rules. Let's see how the two different interpretations for the stochastic +integral impose the necessity for different chain rules. + +An SDE can be written also in differential form. Because we cannot formally form the derivative +``\mathrm{d} W / \mathrm{d} t``, since ``W`` is nowhere differentiable, we write an SDE in +differential form as: + +```math +\begin{aligned} +{\color{Green}\text{Itô}}&: {\color{Green}\mathrm{d} x_t = f(x_t) \mathrm{d} t + g(x_t) \mathrm{d} W_t} , \\ +{\color{Magenta}\text{Stratonovich}}&: {\color{Magenta}\mathrm{d} x_t = f(x_t) \mathrm{d} t + g(x_t) \circ \mathrm{d} W_t} . +\end{aligned} +``` + +The circle in ``g(x_t) \circ \mathrm{d} W_t`` is used to differentiate between Itô or +Stratonovich calculus. + +Let's now assume we perform a variable change ``y = G(x)``. It turns out that according to +which interpretation of the stochastic integral one chooses to use, then the following chain +rule must be used: + +```math +\begin{aligned} +{\color{Green}\text{Itô}}&: {\color{Green}\mathrm{d} y_t = \frac{\mathrm{d} G}{\mathrm{d} x} \mathrm{d} x_t + \frac{1}{2} g(x_t)^2 \frac{\mathrm{d}^2 G}{\mathrm{d} x^2} \mathrm{d} t = \left[ \frac{\mathrm{d} G}{\mathrm{d} x} f(x_t) + \frac{1}{2} g(x_t)^2 \frac{\mathrm{d}^2 G}{\mathrm{d} x^2} \right] \mathrm{d} t + \frac{\mathrm{d} G}{\mathrm{d} x} g(x_t) \mathrm{d} W_t} , \\ +{\color{Magenta}\text{Stratonovich}}&: {\color{Magenta}\mathrm{d} y_t = \frac{\mathrm{d} G}{\mathrm{d} x} \mathrm{d} x_t = \frac{\mathrm{d} G}{\mathrm{d} x} f(x_t) \mathrm{d} t + \frac{\mathrm{d} G}{\mathrm{d} x} g(x_t) \circ \mathrm{d} W_t} . +\end{aligned} +``` + +The above are the so called *stochastic chain rules*. All derivatives of ``G`` are evaluated +at ``x_t``. For Stratonovich calculus, the chain rule resembles the usual chain rule one learns +in calculus; for Itô there exists an additional term, often referred to as the "drift-term": +``{\color{Green}\tfrac{1}{2} g^2 \, \mathrm{d}^2G / \mathrm{d} x^2}``. + +It's easy to see that the extra drift-term in Itô's interpretation of the stochastic integral, +is *exactly* equal to the ensemble mean over forcing realizations of the Stratonovich +stochastic integral. That's because the Itô stochastic integral has, by construction, zero +ensemble mean since at every instant the noise is multiplied with ``g`` which is evaluated +at a time instance before the action of the noise; ``g`` and ``\mathrm{d} W`` are uncorrelated +and thus: + +```math +{\color{Green} \left \langle g(x_t) \mathrm{d} W_t \right \rangle = 0} \quad \text{while} \quad {\color{Magenta} \left \langle g(x_t) \circ \mathrm{d} W_t \right \rangle \ne 0} . +``` + +The above is demonstrated by evaluating the simple stochastic integral: + +```math +\begin{aligned} +{\color{Green} \text{Itô}}&: {\color{Green} \left \langle \int_{t_0}^{t} W_s \, \mathrm{d} W_s \right \rangle \approx \sum_{j} \left \langle W_j (W_{j+1} - W_j) \right \rangle} \\ +&\hspace{7.3em} {\color{Green} = \sum_j \left \langle W_j W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \sum_{j} t_j - t_j = 0} , \\ +{\color{Magenta}\text{Stratonovich}}&: {\color{Magenta}\left \langle \int_{t_0}^{t} W_s \circ \mathrm{d} W_s \right \rangle \approx \sum_j \left \langle \frac1{2}(W_j + W_{j+1}) (W_{j+1} - W_j)\right \rangle} \\ +&\hspace{7.3em} {\color{Magenta} = \frac1{2} \sum_j \left \langle W_{j+1} W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \frac1{2} \sum_j t_{j+1} - t_j = \frac{t}{2}} . +\end{aligned} +``` + +SDEs rarely can be solved in closed form; most often numerical solution of SDEs is brought to +the rescue. Itô calculus has the advantage that is very easily implemented numerically. On +the other hand, the chain rule in Stratonovich calculus coincides with that in normal calculus. +This stems from the fact that in the Stratonovich interpretation the white noise process is as +a series of colored noise processes with the de-correlation time tending to zero. This made +Stratonovich calculus more popular in the physics community. A nice discussion on the differences +and similarities between the two calculi is given by [van Kampen](https://doi.org/10.1007/BF01007642). + +## A simple Stochastic Differential Equation (SDE): the Ornstein--Uhlenbeck process + +One of the simplest SDEs is the Ornstein--Uhlenbeck process, a variation of which is: + +```math +x(t) = - \int_{t_0}^{t} \mu x(s) \, \mathrm{d} s + \int_{t_0}^{t} \sqrt{\sigma} \, \mathrm{d} W_s . \tag{1} +``` + +Note that in differential form (1) is written as: + +```math +\mathrm{d} x_t = - \mu x_t \, \mathrm{d} t + \sqrt{\sigma} \, \mathrm{d} W_t . \tag{2} +``` + +Luckily, for (2) we don't need to distinguish between Itô and Stratonovich, since ``g`` is +independent of ``x(t)``. But note that oftentimes this is not the case; that ``g`` is independent +of ``x(t)`` is only a fortuitous coincident for this particular SDE. + +How do we time-step SDE (2) numerically? Let us assume a discretization of time into time-steps +of duration ``\tau``, i.e., ``t_j = (j-1) \tau``, ``j=1, 2, \dots``. (What follows is easily +generalized to non-uniform time discretizations.) With that in mind, we denote ``x_j \equiv x(t_j)``. +Then the Euler--Mayorama time-stepping scheme for (2) is + +```math + x_{j+1} = x_j + (-\mu x_j) \tau + \sqrt{\sigma} (W_{j+1} - W_j) . +``` + +Now let us ask the following question: How can we compute the work done by the noise? +In other words, if we are interested in the evolution of the "energy", defined as +``E \equiv \tfrac{1}{2} x^2``, then how does the noise term attribute in the growth of ``E``? +To answer that we first have to find the SDE that energy ``E`` obeys. But, in doing so, it +is important to adopt a single interpretation for computing stochastic integrals as now a +transformation of variables is needed. That is, depending on whether we choose to interpret +the stochastic integrals according to Itô or to Stratonovich calculus, ``E`` evolves according +to: + +```math +\hspace{3.35em} {\color{Green} \text{Itô}} : {\color{Green} \mathrm{d} E_t = \left ( -2 \mu E_t + \tfrac{1}{2} \sigma \right ) \mathrm{d} t + x_t \sqrt{\sigma} \mathrm{d} W_t} , \tag{3} +``` +```math +\hspace{-3.35em} {\color{Magenta} \text{Stratonovich}} : {\color{Magenta} \mathrm{d} E_t = -2 \mu E_t \mathrm{d} t + x_t \circ \sqrt{\sigma} \mathrm{d} W_t} . \tag{4} +``` + +The term ``-2 \mu E_t`` in both cases is the dissipation of energy by the ``\mu`` term; the +rest of the terms involve the noise. How do we compute the work ``P`` done by the noise? +Well, it follows that: + +```math +\begin{aligned} +{\color{Green} \text{Itô}} &: {\color{Green} P_t = \tfrac{1}{2} \sigma \mathrm{d} t + \sqrt{\sigma} x_t \mathrm{d} W_t \approx \tfrac{1}{2} \sigma \, \mathrm{d}t + \sqrt{\sigma} x_j (W_{j+1} - W_j)} , \\ +{\color{Magenta} \text{Stratonovich}} &: {\color{Magenta} P_t = x_t \circ \sqrt{\sigma} \mathrm{d} W_t \approx \sqrt{\sigma} x \left ( \tfrac{1}{2} (t_j + t_{j+1}) \right ) (W_{j+1} - W_j)} . +\end{aligned} +``` + +Now let's assume for a moment that we didn't know the rules for transforming Stratonovich to +Itô and we were wondering what is the extra drift term we have to include in the Itô formulations, +i.e., the ``\tfrac{1}{2} \sigma`` term. We can compute the Itô's drift-term using the fact that +it is exactly equal to ``\langle x_t \circ \sqrt{\sigma} \mathrm{d} W_t \rangle``; and for the +latter we can use the "usual" calculus. That is, we rewrite (1) as: + +```math +\dot{x} = -\mu x + \xi , \tag{5} +``` + +where ``\xi(t)`` is understood to be the "continuous" version of the white-noise process (which +is formally only understood in terms of distributions). The forcing ``\xi`` has the properties: + +```math +\left \langle \xi(t) \right \rangle = 0 \quad \text{and} \quad \left \langle \xi(t) \xi(t') \right \rangle = \sigma \delta(t - t') . +``` + +Thus we need to compute ``\langle P_t \rangle = \langle x(t) \xi(t) \rangle``. But (5) formally +has the solution: + +```math +x(t) = e^{-\mu t} x(0) + \int_0^t e^{-\mu (t - s)} \xi(s) \, \mathrm{d} s . +``` + +and using this solution we get + +```math +\langle P_t \rangle = \langle x(t) \xi(t) \rangle = e^{-\mu t} \underbrace{\langle x(0) \xi(t) \rangle}_{=0} + \int_0^t e^{-\mu (t - s)} \langle \xi(t) \xi(s) \rangle \, \mathrm{d} s = \sigma \int_0^t e^{- \mu (t - s)} \delta(t - s) \, \mathrm{d} s = \frac{\sigma}{2} . +``` + +Above we used that ``\int_0^t \delta(t - s) \mathrm{d} s = \tfrac{1}{2}``, which is consistent +with the Stratonovich symmetric interpretation of stochastic integrals. + +### Numerical implementation + +How do we time-step the equation for ``E``? In the case of Itô's interpretation, (3), we use +the Euler--Maruyama time-stepping scheme: + +```math + E_{j+1} = E_j + \left ( -2 \mu E_j + \frac{\sigma}{2} \right ) \tau + \sqrt{\sigma} x_j (W_{j+1} - W_j). +``` +However, we cannot use Euler--Maruyama for time-stepping the corresponding Stratonovich +version of (4), since the Euler--Maruyama scheme involves "Itô"-thinking. To time-step (4) we +have to approximate ``g`` in the middle of the time-step. There are many ways to do that, one +of which is the, so called, Euler--Heun method: + +```math +\begin{aligned} +\widetilde{E}_{j+1} &= E_j + (- 2\mu E_j) \tau + \sqrt{\sigma} x_j (W_{j+1} - W_j), \\ +E_{j+1} &= E_j + \left( -2 \mu \frac{E_j + \widetilde{E}_{j + 1}}{2} \right)\tau + \sqrt{\sigma}\frac{x_j + x_{j+1}}{2} (W_{j+1} - W_j) . +\end{aligned} +``` + +Let's apply not Euler--Maruyama and Euler--Heun schemes to time-step (3) and (4) respectively +and compare the results with those obtained from time-stepping (2) and computing ``E`` a +posteriori. + +Figure below compares the energy evolution as predicted by: +- direct computation from the ``x_t`` time-series: ``\tfrac{1}{2} x_t^2``, +- time-integration of (3) using Euler--Maruyama, and +- time-integration of (4) using Euler--Heun. + + +```@setup 1 +using Plots +Plots.default(lw=2) +``` + +```@example 1 +using Plots +using Statistics: mean +using Random: randn, seed! +seed!(1234) # for reproducing the same plots + + μ = 0.2 # drag + σ = 0.2 # noise strength + dt = 0.01 # timestep + nsteps = 2001 # total timesteps +n_realizations = 1000 # how many forcing realizations + +t = 0:dt:(nsteps-1)*dt # time + +ΔW = randn(nsteps, n_realizations) * sqrt(dt) # noise + +# Numerical calculation +x = zeros(size(ΔW)) +E_ito = zeros(size(ΔW)) +E_str = zeros(size(ΔW)) +E_numerical = zeros(size(ΔW)) + +for j = 2:nsteps # time step the equation + + # time-step dxₜ = -μ xₜ + √σ dWₜ + @. x[j, :] = x[j-1, :] + -μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :] + + # time-step dEₜ = (-2μ Eₜ + ½σ) dt + √σ xₜ dWₜ + @. E_ito[j, :] = E_ito[j-1, :] + (-2μ * E_ito[j-1, :] + + σ/2) * dt + sqrt(σ) * x[j-1, :] * ΔW[j-1, :] + + # time-step dEₜ = -2μ Eₜ dt + √σ xₜ ∘ dWₜ + xbar = @. x[j-1, :] - μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :] + Ebar = @. E_str[j-1, :] + (-2μ * E_str[j-1, :]) * dt + sqrt(σ) * x[j-1, :] * ΔW[j-1, :] + @. E_str[j, :] = E_str[j-1, :] + (-2μ * (E_str[j-1, :] + + Ebar) / 2) * dt + sqrt(σ) * (x[j-1, :] + xbar) / 2 * ΔW[j-1, :] +end + +# direct computation of E from xₜ +@. E_numerical = 0.5 * x^2 + +# compare the three E(t) solutions +plot(μ * t, [E_numerical[:, 1] E_ito[:, 1] E_str[:, 1]], + linewidth = [3 2 1], + label = ["½ xₜ²" "Eₜ (Ito)" "Eₜ (Stratonovich) "], + linestyle = [:solid :dash :dashdot], + xlabel = "μ t", + ylabel = "E", + legend = :topleft, + title = "comparison of E(t) for single realization") + +savefig("assets/energy_comparison.svg"); nothing # hide +``` + +![energy_comparison](assets/energy_comparison.svg) + +Now we can further compute the "energy" budgets, i.e., the work done by the noise versus the +energy loss by the ``\mu`` term, using Itô and Stratonovich formalisms. Figures below show +the ensemble mean energy budgets (using 1000 ensemble members) as computed using Itô and +Stratonovich calculus. For the energy budget to close we have to be consistent: if we time-step +the energy equation based on Stratonovich calculus then we must compute the work also according +to Stratonovich and vice versa. + +```@example 1 +# theoretical results for ⟨E⟩ and d⟨E⟩/dt + E_theory = @. σ/4μ * (1 - exp(-2μ * t)) +dEdt_theory = @. σ/2 * exp(-2μ * t) + +# compute d⟨E⟩/dt numerically +dEdt_ito = mean(E_ito[2:nsteps, :] .- E_ito[1:nsteps-1, :], dims=2) / dt + +# compute the work and dissipation +work_ito = mean(sqrt(σ) * ΔW[1:nsteps-1, :] / dt .* x[1:nsteps-1, :] .+ σ/2, dims=2) +diss_ito = 2*μ * (mean(E_ito[1:nsteps-1, :], dims=2)) + +# Ensemble mean energy budgets from the Itô integration + +plot_E = plot(μ * t, [E_theory mean(E_ito, dims=2)], + linewidth = [3 2], + label = ["theoretical ⟨E⟩" "⟨E⟩ from $n_realizations ensemble member(s)"], + xlabel = "μ t", + ylabel = "E", + legend = :bottomright, + title = "Ito: 𝖽Eₜ = (-2μ Eₜ + ½σ) 𝖽t + xₜ √σ 𝖽Wₜ") + +plot_Ebudget = plot(μ * t[1:nsteps-1], [dEdt_ito work_ito.-diss_ito dEdt_theory[1:nsteps-1]], + linestyle = [:dash :dashdot :solid], + linewidth = [2 1 3], + label = ["numerical 𝖽⟨E⟩/𝖽t" "⟨work - dissipation⟩" "theoretical 𝖽⟨E⟩/𝖽t"], + legend = :topright, + xlabel = "μ t") + +plot(plot_E, plot_Ebudget, layout=(2, 1)) + +savefig("assets/energy_budgets_Ito.svg"); nothing # hide +``` + +![energy_budgets_Ito](assets/energy_budgets_Ito.svg) + + +```@example 1 +# compute d⟨E⟩/dt numerically +dEdt_str = mean(E_str[2:nsteps, :] .- E_str[1:nsteps-1, :], dims=2) / dt + +# compute the work and dissipation +work_str = mean(sqrt(σ) * ΔW[1:nsteps-1, :] / dt .* (x[1:nsteps-1, :] .+ x[2:nsteps, :])/2, dims=2) +diss_str = 2*μ * (mean(E_str[1:nsteps-1, :], dims=2)) + +plot_E = plot(μ * t, [E_theory mean(E_str, dims=2)], + linewidth = [3 2], + label = ["theoretical ⟨E⟩" "⟨E⟩ from $n_realizations ensemble member(s)"], + xlabel = "μ t", + ylabel = "E", + legend = :bottomright, + title = "Stratonovich: 𝖽Eₜ = -2μ Eₜ 𝖽t + xₜ ∘ √σ 𝖽Wₜ") + +plot_Ebudget = plot(μ * t[1:nsteps-1], [dEdt_str[1:nsteps-1] work_str[1:nsteps-1].-diss_str[1:nsteps-1] dEdt_theory[1:nsteps-1]], + linestyle = [:dash :dashdot :solid], + linewidth = [2 1 3], + label = ["numerical 𝖽⟨E⟩/𝖽t" "⟨work - dissipation⟩" "theoretical 𝖽⟨E⟩/𝖽t"], + legend = :bottomleft, + xlabel = "μ t") + +plot(plot_E, plot_Ebudget, layout=(2, 1)) + +savefig("assets/energy_budgets_Stratonovich.svg"); nothing # hide +``` + +![energy_budgets_Stratonovich](assets/energy_budgets_Stratonovich.svg) + + +## A simple Stochastic Partial Differential Equation (SPDE) + +We would like now to transfer all the knowledge we got from the previous sections to PDEs. +In particular we'll start by focussing on the simple SPDE: + +```math +\partial_t \nabla^2 \psi(\bm{x}, t) = - \mu \nabla^2 \psi(\bm{x}, t) + \xi(\bm{x}, t) , \tag{6} +``` + +with periodic boundary conditions in both ``x`` and ``y``. SPDE (6) is also equivalently written as: + +```math +\mathrm{d} \nabla^2 \psi_{t}(\bm{x}) = - \mu \nabla^2 \psi_{t} (\bm{x}) \mathrm{d} t + \sqrt{\sigma} \mathrm{d} W_{t} (\bm{x}) . +``` + +The form (6) is the continuous version, similar to (2). In this SPDE, since the forcing is +purely additive, i.e., it does not depend on the state of the system, both Itô and Stratonovich +interpretations coincide. + +The forcing ``\xi`` obeys: + +```math +\langle \xi(\bm{x},t) \rangle = 0 \quad \text{and} \quad \langle \xi(\bm{x}, t) \xi(\bm{x}', t') \rangle= Q(\bm{x} - \bm{x}') \delta(t-t'), +``` + +that is, the forcing is white in time but spatially correlated; its spatial correlation is +prescribed by the function ``Q`` which is, necessarily, homogeneous in all its arguments +(see discussion by [Constantinou (2015)](http://arxiv.org/abs/1503.07644); Appendix A). + +Equation (6) above describes the vorticity evolution of a two-dimensional fluid ``\nabla^2 \psi`` +that is stochastically forced while dissipated by linear drag ``\mu``. The energy of the +fluid is: + +```math +E = \tfrac{1}{2} \overline{|\bm{\nabla}\psi|^2}^{x,y} = -\tfrac{1}{2} \overline{\psi \nabla^2\psi}^{x,y}, +``` + +where the overbar denotes average over ``x`` and ``y`` and an integration-by-parts was carried +through in the last equality. To obtain the energy equation we multiply (6) with ``-\psi`` and +average over the whole domain. Thus, the work done by the forcing is given by: + +```math +P = - \, \overline{\psi \, \xi}^{x,y}, +``` + +but the above is a stochastic integral and it is meaningless without a rule for computing the stochastic integral. + +Numerically, the work done by the forcing at the ``j``-th timestep can be obtained +Stratonovich-wise via: + +```math +\begin{aligned} +P_j = - \, \overline{\frac{\psi(\bm{x}, t_j) + \psi(\bm{x}, t_{j+1})}{2} \xi(\bm{x}, t_{j+1}) }^{x,y} , +\end{aligned} +``` + +or Itô-wise as + +```math +\begin{aligned} +P_j = -\, \overline{ \psi(\bm{x}, t_j) \xi(\bm{x}, t_{j+1}) }^{x,y} + \text{drift} . +\end{aligned} +``` + +But how much is the Itô drift term in this case? As in the previous section, the drift is +*precisely* the ensemble mean of the Stratonovich work, i.e.: + +```math +\textrm{Ito drift}= - \overline{\langle \underbrace{\psi(\bm{x}, t) \circ \xi(\bm{x}, t)}_{\textrm{Stratonovich}} \rangle}^{x,y} . +``` + +But again, the above can be computed using the "formal" solution of (6): + +```math +\psi(\bm{x}, t) = e^{-\mu t} \psi(\bm{x}, 0) + \int_0^t e^{- \mu (t - s)} \nabla^{-2} \xi(\bm{x}, s) \, \mathrm{d} s , +``` + +which implies + +```math +\begin{aligned} +\text{drift} & = -\overline{e^{-\mu t} \underbrace{\left \langle \psi(\bm{x}, 0) \xi(\bm{x}, t) \right \rangle}_{=0}}^{x,y} - \int_0^t e^{- \mu (t - s)} \overline{\nabla^{-2} \left \langle \xi(\bm{x}, s) \xi(\bm{x}, t) \right\rangle}^{x,y} \, \mathrm{d} s \\ +& = -\int_0^t e^{-\mu(t - s)} \overline{\underbrace{\left [ \nabla^{-2} Q (\bm{x}) \right ] \big|_{\bm{x}=0}}_{\text{independent of }x,y} \, \delta(t - s)}^{x,y} \, \mathrm{d} s \\ +& = -\frac1{2} \nabla^{-2} Q(\bm{x}) \big|_{\bm{x}=0} \\ +& = - \frac1{2} \left [ \nabla^{-2} \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \widehat{Q}(\bm{k}) \, e^{i \bm{k} \bm{\cdot}\bm{x}} \right ]_{\bm{x}=0} \\ +& = \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \frac{\widehat{Q}(\bm{k})}{2 |\bm{k}|^2} . +\end{aligned} +``` + +Thus, the drift, or in this case the mean energy input rate by the stochastic forcing, is +precisely determined from the spatial correlation of the forcing, ``Q``. Let us denote the +drift as: + +```math +\varepsilon \equiv \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \frac{\widehat{Q}(\bm{k})}{2 |\bm{k}|^2} . \tag{7} +``` + +Using the above, the work for a single forcing realization at the ``j``-th timestep is numerically +computed as: + +```math +{\color{Green} \text{Itô}} : {\color{Green} P_j = -\overline{\psi(\bm{x}, t_j) \xi(\bm{x}, t_{j+1})}^{x,y} + \varepsilon} , \tag{8} +``` +```math +{\color{Magenta} \text{Stratonovich}} : {\color{Magenta} P_j = -\overline{\frac{\psi(\bm{x}, t_j) + \psi(\bm{x}, t_{j+1})}{2} \xi(\bm{x}, t_{j+1})}^{x,y}} . \tag{9} +``` + +Remember, previously the work done by the stochastic forcing was: +```math +\mathrm{d} P_t = {\color{Green} \frac{\sigma}{2}\mathrm{d} t + \sqrt{\sigma} x_t \mathrm{d} W_t} = {\color{Magenta} \sqrt{\sigma} x_t \circ \mathrm{d} W_t} , +``` +and by sampling over various forcing realizations: +```math +\langle \mathrm{d} P_t \rangle = \frac{\sigma}{2} \mathrm{d} t = \langle \sqrt{\sigma} x_t \circ \mathrm{d} W_t \rangle . +``` + +All modules in GeophysicalFlows.jl use Stratonovich calculus. For example, the domain-averaged +energy injected per unit time by the forcing in the `TwoDNavierStokes` module is computed +using (9) via + +```@docs +GeophysicalFlows.TwoDNavierStokes.energy_work +``` + +## A bit more elaborate SPDE + +It turns out everything carries through if in our SPDE above for the 2D vorticity equation we +also include the nonlinear advection terms: + +```math +\partial_t \nabla^2 \psi(\bm{x}, t) + \mathsf{J}(\psi, \nabla^2 \psi) = - \mu \nabla^2 \psi(\bm{x}, t) + \xi(\bm{x}, t) . \tag{10} +``` + +The nonlinearity does not alter the Itô drift; thus the ensemble mean energy input by the +stochastic forcing, remains the same. We can easily verify this from the "formal" solution +of (10): + +```math +\psi(\bm{x}, t) = e^{- \mu t} \psi(\bm{x}, 0) + \int_0^t e^{- \mu (t - s)} \nabla^{-2} \xi(\bm{x}, s) \, \mathrm{d} s - \int_0^t \nabla^{-2} \mathsf{J} \left ( \psi(\bm{x}, s), \nabla^2 \psi(\bm{x}, s) \right ) \mathrm{d} s , +``` + +When multiplied with ``\xi(\bm{x}, t)`` the last term vanishes since its only non-zero +contribution comes from the point ``s = t``, which is of measure zero (in the integrated sense). + +A demonstration of how the energy budgets can be computed for a case with stochastic forcing +can be found in an [example of the TwoDNavierStokes](../generated/twodnavierstokes_stochasticforcing_budgets/) +module. diff --git a/src/barotropicqgql.jl b/src/barotropicqgql.jl index 1966fb21..2b89d2c5 100644 --- a/src/barotropicqgql.jl +++ b/src/barotropicqgql.jl @@ -369,8 +369,7 @@ Returns the domain-averaged rate of work of energy by the forcing, `params.Fh`. end @inline function work(sol, vars::StochasticForcedVars, grid) - @. vars.uh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich - # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + @. vars.uh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 8463bce6..5b0f419b 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -209,21 +209,42 @@ numberoflayers(::SingleLayerParams) = 1 # Equations # --------- -function hyperdissipation(dev, params, grid) +""" + hyperviscosity(dev, params, grid) +Returns the linear operator `L` that corresponds to (hyper)-viscosity of order ``n_ν`` with +coefficient ``ν`` for ``n`` fluid layers. +```math +L_j = - ν |𝐤|^{2 n_ν}, j = 1, ...,n . +``` +""" +function hyperviscosity(dev, params, grid) T = eltype(grid) L = ArrayType(dev){T}(undef, (grid.nkr, grid.nl, numberoflayers(params))) @. L = - params.ν * grid.Krsq^params.nν @views @. L[1, 1, :] = 0 + return L end +""" + LinearEquation(dev, params, grid) +Returns the `equation` for a multi-layer quasi-geostrophic problem with `params` and `grid`. +The linear opeartor ``L`` includes only (hyper)-viscosity and is computed via +`hyperviscosity(dev, params, grid)`. The nonlinear term is computed via function `calcNlinear!()`. +""" function LinearEquation(dev, params, grid) - L = hyperdissipation(dev, params, grid) + L = hyperviscosity(dev, params, grid) return FourierFlows.Equation(L, calcNlinear!, grid) end - + +""" + Equation(dev, params, grid) +Returns the `equation` for a multi-layer quasi-geostrophic problem with `params` and `grid`. +The linear opeartor ``L`` includes only (hyper)-viscosity and is computed via +`hyperviscosity(dev, params, grid)`. The nonlinear term is computed via function `calcN!()`. +""" function Equation(dev, params, grid) - L = hyperdissipation(dev, params, grid) + L = hyperviscosity(dev, params, grid) return FourierFlows.Equation(L, calcN!, grid) end @@ -319,8 +340,8 @@ end """ calcS!(S, Fp, Fm, nlayers, grid) -Constructs the array S, which consists of nlayer x nlayer static arrays S_kl that relate -the q's and ψ's at every wavenumber: q̂_{k, l} = S_kl * ψ̂_{k, l}. +Constructs the array ``𝕊``, which consists of `nlayer` x `nlayer` static arrays ``𝕊_𝐤`` that +relate the ``q̂_j``'s and ``ψ̂_j``'s at every wavenumber: ``q̂_𝐤 = 𝕊_𝐤 ψ̂_𝐤``. """ function calcS!(S, Fp, Fm, nlayers, grid) F = Matrix(Tridiagonal(Fm, -([Fp; 0] + [0; Fm]), Fp)) @@ -335,8 +356,8 @@ end """ calcS⁻¹!(S, Fp, Fm, nlayers, grid) -Constructs the array S⁻¹, which consists of nlayer x nlayer static arrays (S_kl)⁻¹ that -relate the q's and ψ's at every wavenumber: ψ̂_{k, l} = (S_kl)⁻¹ * q̂_{k, l}. +Constructs the array ``𝕊⁻¹``, which consists of `nlayer` x `nlayer` static arrays ``(𝕊_𝐤)⁻¹`` +that relate the ``q̂_j``'s and ``ψ̂_j``'s at every wavenumber: ``ψ̂_𝐤 = (𝕊_𝐤)⁻¹ q̂_𝐤``. """ function calcS⁻¹!(S⁻¹, Fp, Fm, nlayers, grid) T = eltype(grid) @@ -355,6 +376,15 @@ end # Solvers # ------- +""" + calcN!(N, sol, t, clock, vars, params, grid) +Compute the nonlinear term, that is the advection term and the forcing, +```math +N(q̂_j) = - \\widehat{𝖩(ψ_j, q_j)} - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} + δ_{j, n} μ |𝐤|^2 ψ̂_n + F̂_j , +``` +by calling `calcN_advection!` and `addforcing!`. +""" function calcN!(N, sol, t, clock, vars, params, grid) nlayers = numberoflayers(params) calcN_advection!(N, sol, vars, params, grid) @@ -374,7 +404,11 @@ end """ calcN_advection!(N, sol, vars, params, grid) -Calculates the advection term. +Compute the advection term and stores it in `N`: +```math +N(q̂_j) = - \\widehat{𝖩(ψ_j, q_j)} - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} . +``` """ function calcN_advection!(N, sol, vars, params, grid) @. vars.qh = sol @@ -418,7 +452,12 @@ end """ calcN_linearadvection!(N, sol, vars, params, grid) -Calculates the advection term of the linearized equations. +Compute the advection term of the linearized equations and stores it in `N`: +```math +N(q̂_j) = - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} . +``` + """ function calcN_linearadvection!(N, sol, vars, params, grid) @. vars.qh = sol @@ -456,11 +495,17 @@ function calcN_linearadvection!(N, sol, vars, params, grid) return nothing end +""" + addforcing!(N, sol, t, clock, vars, params, grid) +When the problem includes forcing, calculate the forcing term ``F̂_j`` and add it to the +nonlinear term ``N``. +""" addforcing!(N, sol, t, clock, vars::Vars, params, grid) = nothing function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid) params.calcFq!(vars.Fqh, sol, t, clock, vars, params, grid) @. N += vars.Fqh + return nothing end @@ -524,8 +569,8 @@ set_q!(prob, q) = set_q!(prob.sol, prob.params, prob.vars, prob.grid, q) set_ψ!(params, vars, grid, sol, ψ) set_ψ!(prob) -Set the solution `prob.sol` to correspond to the transform of streamfunction `ψ` and -updates variables. +Set the solution `prob.sol` to the transform `qh` that corresponds to streamfunction `ψ` +and updates variables. """ function set_ψ!(sol, params, vars, grid, ψ) A = typeof(vars.ψ) @@ -561,11 +606,12 @@ is the number of layers in the fluid. (When ``n=1``, only the kinetic energy is The kinetic energy at the ``j``-th fluid layer is ```math -\\textrm{KE}_j = \\frac{H_j}{H} \\int \\frac1{2} |\\boldsymbol{\\nabla} \\psi_j|^2 \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} \\ , \\quad j = 1, \\dots, n \\ , +𝖪𝖤_j = \\frac{H_j}{H} \\int \\frac1{2} |{\\bf ∇} ψ_j|^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{H_j}{H} \\sum_{𝐤} |𝐤|² |ψ̂_j|², \\ j = 1, ..., n \\ , ``` -while the potential energy that corresponds to the interface ``j+1/2`` (i.e., interface between the ``j``-th and ``(j+1)``-th fluid layer) is +while the potential energy that corresponds to the interface ``j+1/2`` (i.e., the interface +between the ``j``-th and ``(j+1)``-th fluid layer) is ```math -\\textrm{PE}_{j+1/2} = \\int \\frac1{2} \\frac{f_0^2}{g'_{j+1/2}} (\\psi_j - \\psi_{j+1})^2 \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} \\ , \\quad j = 1, \\dots, n-1 \\ . +𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} \\sum_{𝐤} |ψ_j - ψ_{j+1}|², \\ j = 1, ..., n-1 \\ . ``` """ function energies(vars, params, grid, sol) @@ -612,15 +658,14 @@ verticalfluxes``_{3/2},...,``verticalfluxes``_{n-1/2}``, where ``n`` is the tota The lateral eddy fluxes whithin the ``j``-th fluid layer are ```math -\\textrm{lateralfluxes}_j = \\frac{H_j}{H} \\int U_j \\, v_j \\, \\partial_y u_j -\\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} \\ , \\quad j = 1, \\dots, n \\ , +\\textrm{lateralfluxes}_j = \\frac{H_j}{H} \\int U_j v_j ∂_y u_j +\\frac{𝖽x 𝖽y}{L_x L_y} , \\ j = 1, ..., n \\ , ``` while the vertical eddy fluxes at the ``j+1/2``-th fluid interface (i.e., interface between the ``j``-th and ``(j+1)``-th fluid layer) are ```math -\\textrm{verticalfluxes}_{j+1/2} = \\int \\frac{f_0^2}{g'_{j+1/2} H} (U_j - U_{j+1}) \\, -v_{j+1} \\, \\psi_{j} \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} \\ , \\quad -j = 1 , \\dots , n-1. +\\textrm{verticalfluxes}_{j+1/2} = \\int \\frac{f₀²}{g'_{j+1/2} H} (U_j - U_{j+1}) \\, +v_{j+1} ψ_{j} \\frac{𝖽x 𝖽y}{L_x L_y} , \\ j = 1, ..., n-1. ``` """ function fluxes(vars, params, grid, sol) @@ -640,7 +685,7 @@ function fluxes(vars, params, grid, sol) lateralfluxes *= grid.dx * grid.dy / (grid.Lx * grid.Ly * sum(params.H)) for j = 1:nlayers-1 - CUDA.@allowscalar verticalfluxes[j] = sum(@views @. params.f₀^2 / params.g′[j] * (params.U[: ,:, j] - params.U[:, :, j+1]) * vars.v[:, :, j+1] * vars.ψ[:, :, j] ; dims=(1, 2))[1] + CUDA.@allowscalar verticalfluxes[j] = sum(@views @. params.f₀^2 / params.g′[j] * (params.U[: ,:, j] - params.U[:, :, j+1]) * vars.v[:, :, j+1] * vars.ψ[:, :, j]; dims=(1, 2))[1] CUDA.@allowscalar verticalfluxes[j] *= grid.dx * grid.dy / (grid.Lx * grid.Ly * sum(params.H)) end diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index bea14bdc..b17eae70 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -316,7 +316,7 @@ end Returns the domain-averaged kinetic energy of the fluid, ```math -\\textrm{KE} = \\int \\frac1{2} |\\boldsymbol{\\nabla} \\psi|^2 \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} = \\sum_{\\boldsymbol{k}} \\frac1{2} |\\boldsymbol{k}|^2 |\\hat{\\psi}|^2 \\ . +\\textrm{KE} = \\int \\frac1{2} |\\bm{\\nabla} \\psi|^2 \\frac{\\mathrm{d}^2 \\bm{x}}{L_x L_y} = \\sum_{\\bm{k}} \\frac1{2} |\\bm{k}|^2 |\\hat{\\psi}|^2 \\ . ``` """ function kinetic_energy(sol, vars, params, grid) @@ -334,7 +334,7 @@ kinetic_energy(prob) = kinetic_energy(prob.sol, prob.vars, prob.params, prob.gri Returns the domain-averaged potential energy of the fluid, ```math -\\textrm{PE} = \\int \\frac1{2} \\frac{\\psi^2}{\\ell^2} \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} = \\sum_{\\boldsymbol{k}} \\frac1{2} \\frac{|\\hat{\\psi}|^2}{\\ell^2} \\ . +\\textrm{PE} = \\int \\frac1{2} \\frac{\\psi^2}{\\ell^2} \\frac{\\mathrm{d}^2 \\bm{x}}{L_x L_y} = \\sum_{\\bm{k}} \\frac1{2} \\frac{|\\hat{\\psi}|^2}{\\ell^2} \\ . ``` """ function potential_energy(sol, vars, params::EquivalentBarotropicQGParams, grid) @@ -366,7 +366,7 @@ flow. Returns the domain-averaged enstrophy ```math -\\int \\frac1{2} (q + \\eta)^2 \\frac{\\mathrm{d}^2 \\boldsymbol{x}}{L_x L_y} = \\sum_{\\boldsymbol{k}} \\frac1{2} |\\hat{q} + \\hat{\\eta}|^2 \\ . +\\int \\frac1{2} (q + \\eta)^2 \\frac{\\mathrm{d}^2 \\bm{x}}{L_x L_y} = \\sum_{\\bm{k}} \\frac1{2} |\\hat{q} + \\hat{\\eta}|^2 \\ . ``` """ function enstrophy(sol, vars, params, grid) @@ -428,8 +428,7 @@ end @inline function energy_work(sol, vars::StochasticForcedVars, params::BarotropicQGParams, grid) energy_workh = vars.uh # use vars.uh as scratch variable - @. energy_workh = grid.invKrsq * (vars.prevsol + sol)/2 * conj(vars.Fh) # Stratonovich - # @. energy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + @. energy_workh = grid.invKrsq * (vars.prevsol + sol)/2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end @@ -453,8 +452,7 @@ end @inline function enstrophy_work(sol, vars::StochasticForcedVars, params::BarotropicQGParams, grid) enstrophy_workh = vars.uh # use vars.uh as scratch variable - @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich - # @. enstrophy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid) end diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index 0f44dc95..d9ea6def 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -398,8 +398,7 @@ end @inline function energy_work(sol, vars::StochasticForcedVars, grid) energy_workh = vars.uh # use vars.uh as scratch variable - @. energy_workh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich - # @. energy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + @. energy_workh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid) end @@ -424,8 +423,7 @@ end @inline function enstrophy_work(sol, vars::StochasticForcedVars, grid) enstrophy_workh = vars.uh # use vars.uh as scratch variable - @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich - # @. enstrophy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid) end diff --git a/test/test_barotropicqgql.jl b/test/test_barotropicqgql.jl index dbeba4b4..ae2d484c 100644 --- a/test/test_barotropicqgql.jl +++ b/test/test_barotropicqgql.jl @@ -96,10 +96,7 @@ function test_bqgql_stochasticforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dEdt_computed = W[2:E.i] + ε - D[1:E.i-1] - R[1:E.i-1] # Ito - dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1] # Stratonovich + dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1] return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-3) end diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 03e1e0d7..47825381 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -97,10 +97,7 @@ function test_1layerqg_stochasticforcing_energybudget(dev::Device=CPU(); n=256, dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dEdt_computed = W[2:E.i] + ε - D[1:E.i-1] - R[1:E.i-1] # Ito - dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1] # Stratonovich + dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1] return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-3) end @@ -206,10 +203,7 @@ function test_1layerqg_stochasticforcing_enstrophybudget(dev::Device=CPU(); n=25 dZdt_numerical = (Z[2:Z.i] - Z[1:Z.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dZdt_computed = W[2:Z.i] + εᶻ - D[1:Z.i-1] - R[1:Z.i-1] # Ito - dZdt_computed = W[2:Z.i] - D[1:Z.i-1] - R[1:Z.i-1] # Stratonovich + dZdt_computed = W[2:Z.i] - D[1:Z.i-1] - R[1:Z.i-1] return isapprox(dZdt_numerical, dZdt_computed, rtol=1e-3) end diff --git a/test/test_surfaceqg.jl b/test/test_surfaceqg.jl index dcef845c..a9a74542 100644 --- a/test/test_surfaceqg.jl +++ b/test/test_surfaceqg.jl @@ -208,10 +208,7 @@ function test_sqg_stochasticforcing_buoyancy_variance_budget(dev::Device=CPU(); dBdt_numerical = (B[2:B.i] - B[1:B.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dBdt_computed = W[2:B.i] + εᵇ - D[1:B.i-1] # Ito - dBdt_computed = W[2:B.i] - D[1:B.i-1] # Stratonovich + dBdt_computed = W[2:B.i] - D[1:B.i-1] return isapprox(dBdt_numerical, dBdt_computed, rtol = 1e-4) end diff --git a/test/test_twodnavierstokes.jl b/test/test_twodnavierstokes.jl index f04df809..1e07e4c3 100644 --- a/test/test_twodnavierstokes.jl +++ b/test/test_twodnavierstokes.jl @@ -69,10 +69,7 @@ function test_twodnavierstokes_stochasticforcing_energybudget(dev::Device=CPU(); dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dEdt_computed = W[2:E.i] + ε + D[1:E.i-1] + R[1:E.i-1] # Ito - dEdt_computed = W[2:E.i] + D[1:E.i-1] + R[1:E.i-1] # Stratonovich + dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1] return isapprox(dEdt_numerical, dEdt_computed, rtol = 1e-3) end @@ -124,10 +121,7 @@ function test_twodnavierstokes_stochasticforcing_enstrophybudget(dev::Device=CPU dZdt_numerical = (Z[2:Z.i] - Z[1:Z.i-1]) / prob.clock.dt - # If the Ito interpretation was used for the work - # then we need to add the drift term - # dZdt_computed = W[2:Z.i] + εᶻ + D[1:Z.i-1] + R[1:Z.i-1] # Ito - dZdt_computed = W[2:Z.i] + D[1:Z.i-1] + R[1:Z.i-1] # Stratonovich + dZdt_computed = W[2:Z.i] - D[1:Z.i-1] - R[1:Z.i-1] return isapprox(dZdt_numerical, dZdt_computed, rtol = 1e-3) end