Skip to content

Quasi-geostrophic problem documentation #121

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 3 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions toolbox/+otp/+quasigeostrophic/+presets/Canonical.m
Original file line number Diff line number Diff line change
@@ -1,6 +1,28 @@
classdef Canonical < otp.quasigeostrophic.QuasiGeostrophicProblem
% A preset of the Quasi-geostrophic equations starting at the unstable
% state $ψ_0 = 0$.
%
% The value of the Reynolds number is $Re= 450$, the Rossby number is
% $Ro=0.0036$, and the spatial discretization is $255$ interior units in
% the $x$ direction and $511$ interior units in the $y$ direction.
%
% The timespan starts at $t=0$ and ends at $t=100$, which is
% roughly equivalent to $25.13$ years.
%

methods
function obj = Canonical(varargin)
% Create the Canonical Quasi-geostrophic problem object.
%
% Parameters
% ----------
% varargin
% A variable number of name-value pairs. The accepted names are
%
% - ``ReynoldsNumber`` – Value of $Re$.
% - ``RossbyNumber`` – Value of $Ro$.
% - ``Size`` – Two-tuple of the spatial discretization, $[nx, ny]$.
%

Re = 450;
Ro = 0.0036;
Expand Down
Original file line number Diff line number Diff line change
@@ -1,8 +1,32 @@
classdef PopovMouIliescuSandu < otp.quasigeostrophic.QuasiGeostrophicProblem
classdef PopovMouSanduIliescu < otp.quasigeostrophic.QuasiGeostrophicProblem
% A preset for the Quasi-geostrophic equations created for
% CITE ME HERE.
%
% The initial condition is given by an internal file and was created by
% integrating the canonical preset until time $t=100$.
%
% The value of the Reynolds number is $Re= 450$, the Rossby number is
% $Ro=0.0036$, and the spatial discretization is $63$ interior units in
% the $x$ direction and $127$ interior units in the $y$ direction.
%
% The timespan starts at $t=100$ and ends at $t=100.0109$, which is
% roughly equivalent to one day in model time.

methods
function obj = PopovMouIliescuSandu(varargin)
function obj = PopovMouSanduIliescu(varargin)
% Create the PopovMouSanduIliescu Quasi-geostrophic problem object.
%
% Parameters
% ----------
% varargin
% A variable number of name-value pairs. The accepted names are
%
% - ``ReynoldsNumber`` – Value of $Re$.
% - ``RossbyNumber`` – Value of $Ro$.
% - ``Size`` – Two-tuple of the spatial discretization, $[nx, ny]$.
%

defaultsize = [255, 511];
defaultsize = [63, 127];

Re = 450;
Ro = 0.0036;
Expand Down Expand Up @@ -40,9 +64,9 @@

%% Do the rest

sixHours = 6*80/176251.2;
oneday = 24*80/176251.2;

tspan = [100, 100 + sixHours];
tspan = [100, 100 + oneday];

obj = [email protected](tspan, ...
psi0, params);
Expand Down
89 changes: 83 additions & 6 deletions toolbox/+otp/+quasigeostrophic/QuasiGeostrophicProblem.m
Original file line number Diff line number Diff line change
@@ -1,5 +1,74 @@
classdef QuasiGeostrophicProblem < otp.Problem

% A chaotic PDE modeling the flow of a fluid on the earth.
%
% The governing partial differential equation that is discretized is,
%
% $$
% Δψ_t = -J(ψ,ω) - {Ro}^{-1} \partial_x ψ -{Re}^{-1} Δω - {Ro}^{-1} F,
% $$
% where the Jacobian term is a quadratic function,
% $$
% J(ψ,ω) \equiv \partial_x ψ \partial_x ω - \partial_x ψ \partial_x ω,
% $$
% the relationship between the vorticity $ω$ and the stream function $ψ$ is
% $$
% ω = -Δψ,
% $$
% the term $Δ$ is the two dimensional Laplacian over the
% discretization, the terms $\partial_x$ and $\partial_y$ are the first derivatives
% in the $x$ and $y$ directions respectively, $Ro$ is the Rossby number, $Re$ is the Reynolds
% number, and $F$ is a forcing term.
% The spatial domain is fixed to $x ∈ [0, 1]$ and $y ∈ [0, 2]$, and
% the boundary conditions of the PDE are assumed to be zero dirichlet
% everywhere.
%
% A second order finite difference approximation is performed on the
% grid to create the first derivative operators and the Laplacian
% operator.
%
% The Laplacian is defined using the 5-point stencil
% $$
% Δ = \begin{bmatrix} & 1 & \\ 1 & -4 & 1\\ & 1 & \end{bmatrix},
% $$
% which is scaled with respect to the square of the step size in each
% respective direction.
% The first derivatives,
% $$
% \partial_x &= \begin{bmatrix} & 0 & \\ 1/2 & 0 & 1/2\\ & 0 & \end{bmatrix},\\
% \partial_y &= \begin{bmatrix} & 1/2 & \\ 0 & 0 & 0\\ & 1/2 & \end{bmatrix},
% $$
% are the standard second order central finite difference operators in
% the $x$ and $y$ directions.
%
% The Jacobian is discretized using the Arakawa approximation
% CITEME
% $$
% J(ψ,ω) = \frac{1}{3}[ψ_x ω_y - ψ_y ω_x + (ψ ω_y)_x - (ψ ω_x)_y + (ψ_x ω)_y - (ψ_y ω)_x],
% $$
% in order for the system to not become unstable.
%
% The Poisson equation is solved by the eigenvalue sylvester method for
% computational efficiency.
%
% A ADLES MORE HERE.
%
% Notes
% -----
% +---------------------+-----------------------------------------------------------+
% | Type | ODE |
% +---------------------+-----------------------------------------------------------+
% | Number of Variables | $nx \times ny$ |
% +---------------------+-----------------------------------------------------------+
% | Stiff | not typically, depending on $Re$, $Ro$, $Nx$, and $Ny$ |
% +---------------------+-----------------------------------------------------------+
%
% Example
% -------
% >>> problem = otp.quasigeostrophic.presets.PopovMouSanduIliescu;
% >>> sol = problem.solve();
% >>> problem.movie(sol);
%

methods
function obj = QuasiGeostrophicProblem(timeSpan, y0, parameters)

Expand All @@ -19,10 +88,18 @@

methods (Static)

function u = resize(u, newsize)
% resize uses interpolation to resize states

s = size(u);
function psi = resize(psi, newsize)
% Resizes the state onto a new grid by performing interpolation
%
% Parameters
% ----------
% psi : numeric(nx, ny)
% the old state on the $x \times y$ grid.
% newsize : numeric(1, 2)
% the new size as a two-tuple $[nx, ny]$ indicating the new state of the system.
%

s = size(psi);

X = linspace(0, 1, s(1) + 2);
Y = linspace(0, 2, s(2) + 2).';
Expand All @@ -34,7 +111,7 @@
Xnew = Xnew(2:end-1);
Ynew = Ynew(2:end-1);

u = interp2(Y, X, u, Ynew, Xnew);
psi = interp2(Y, X, psi, Ynew, Xnew);

end

Expand Down