Skip to content

Commit 0725884

Browse files
AndreyAPopovAndreyAPopov
authored
Merge pull request #121 from ComputationalScienceLaboratory/QG-docs
Quasi-geostrophic problem documentation
2 parents 9e151ca + a2a13d5 commit 0725884

7 files changed

Lines changed: 197 additions & 37 deletions

File tree

‎docs/references.bib‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -179,6 +179,28 @@ @article{Pet86
179179
doi = {10.1137/0723054}
180180
}
181181

182+
@article{PMSI21,
183+
title={A multifidelity ensemble Kalman filter with reduced order control variates},
184+
author={Popov, Andrey A and Mou, Changhong and Sandu, Adrian and Iliescu, Traian},
185+
journal={SIAM Journal on Scientific Computing},
186+
volume={43},
187+
number={2},
188+
pages={A1134--A1162},
189+
year={2021},
190+
publisher={SIAM}
191+
}
192+
193+
@article{SSWZI11,
194+
title={Approximate deconvolution large eddy simulation of a barotropic ocean circulation model},
195+
author={San, Omer and Staples, Anne E and Wang, Zhu and Iliescu, Traian},
196+
journal={Ocean Modelling},
197+
volume={40},
198+
number={2},
199+
pages={120--132},
200+
year={2011},
201+
publisher={Elsevier}
202+
}
203+
182204
@article{Sch75,
183205
title = {A new approach to explain the ``high irradiance responses'' of photomorphogenesis on the basis of phytochrome},
184206
volume = {2},

‎tests/+problemtests/validateqg.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
fprintf(' Testing Quasi-Geostrophic Equations\n');
44

5-
model = otp.quasigeostrophic.presets.PopovMouIliescuSandu('Nx', 16, 'Ny', 32);
5+
model = otp.quasigeostrophic.presets.PopovMouSanduIliescu('Nx', 16, 'Ny', 32);
66
model.TimeSpan = [0, 0.109];
77

88
[~] = otp.utils.Solver.Nonstiff(model.RHSADLES.F, model.TimeSpan, model.Y0);

‎toolbox/+otp/+quasigeostrophic/+presets/Canonical.m‎

Lines changed: 29 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,32 @@
11
classdef Canonical < otp.quasigeostrophic.QuasiGeostrophicProblem
2+
% A preset of the Quasi-geostrophic equations starting at the unstable
3+
% state $ψ_0 = 0$.
4+
%
5+
% The value of the Reynolds number is $Re= 450$, the Rossby number is
6+
% $Ro=0.0036$, and the spatial discretization is $255$ interior units in
7+
% the $x$ direction and $511$ interior units in the $y$ direction.
8+
%
9+
% The timespan starts at $t=0$ and ends at $t=100$, which is
10+
% roughly equivalent to $25.13$ years.
11+
%
12+
213
methods
314
function obj = Canonical(varargin)
4-
tspan = [0, 100];
15+
% Create the Canonical Quasi-geostrophic problem object.
16+
%
17+
% Parameters
18+
% ----------
19+
% varargin
20+
% A variable number of name-value pairs. The accepted names are
21+
%
22+
% - ``ReynoldsNumber`` – Value of $Re$.
23+
% - ``RossbyNumber`` – Value of $Ro$.
24+
% - ``Nx`` – Spatial discretization in $x$.
25+
% - ``Ny`` – Spatial discretization in $y$.
26+
% - ``ADLambda`` – Scaling factor for approximate deconvolution RHS
27+
% - ``ADPasses`` – Number of AD passes
28+
%
29+
530
params = otp.quasigeostrophic.QuasiGeostrophicParameters( ...
631
'Nx', 255, ...
732
'Ny', 511, ...
@@ -10,6 +35,9 @@
1035
'ADLambda', 1, ...
1136
'ADPasses', 4, ...
1237
varargin{:});
38+
39+
%% Do the rest
40+
tspan = [0, 100];
1341
psi0 = zeros(params.Nx * params.Ny, 1);
1442
obj = obj@otp.quasigeostrophic.QuasiGeostrophicProblem(tspan, psi0, params);
1543
end

‎toolbox/+otp/+quasigeostrophic/+presets/PopovMouIliescuSandu.m‎

Lines changed: 0 additions & 27 deletions
This file was deleted.

‎toolbox/+otp/+quasigeostrophic/+presets/PopovMouSanduIliescu.m‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
classdef PopovMouSanduIliescu < otp.quasigeostrophic.QuasiGeostrophicProblem
2+
% A preset for the Quasi-geostrophic equations created for
3+
% :cite:p:`PMSI21`.
4+
%
5+
% The initial condition is given by an internal file and was created by
6+
% integrating the canonical preset until time $t=100$.
7+
%
8+
% The value of the Reynolds number is $Re= 450$, the Rossby number is
9+
% $Ro=0.0036$, and the spatial discretization is $63$ interior units in
10+
% the $x$ direction and $127$ interior units in the $y$ direction.
11+
%
12+
% The timespan starts at $t=100$ and ends at $t=100.0109$, which is
13+
% roughly equivalent to one day in model time.
14+
15+
methods
16+
function obj = PopovMouSanduIliescu(varargin)
17+
% Create the PopovMouSanduIliescu Quasi-geostrophic problem object.
18+
%
19+
% Parameters
20+
% ----------
21+
% varargin
22+
% A variable number of name-value pairs. The accepted names are
23+
%
24+
% - ``ReynoldsNumber`` – Value of $Re$.
25+
% - ``RossbyNumber`` – Value of $Ro$.
26+
% - ``Nx`` – Spatial discretization in $x$.
27+
% - ``Ny`` – Spatial discretization in $y$.
28+
% - ``ADLambda`` – Scaling factor for approximate deconvolution RHS
29+
% - ``ADPasses`` – Number of AD passes
30+
31+
params = otp.quasigeostrophic.QuasiGeostrophicParameters( ...
32+
'Nx', 255, ...
33+
'Ny', 511, ...
34+
'ReynoldsNumber', 450, ...
35+
'RossbyNumber', 0.0036, ...
36+
'ADLambda', 1, ...
37+
'ADPasses', 4, ...
38+
varargin{:});
39+
40+
% OCTAVE BUG: Octave gives a warning on loading the data, even
41+
% though MATLAB supports this type of private folder loading
42+
spy0s = load('PMISQGICsp.mat');
43+
psi0 = reshape(double(spy0s.y0), 255, 511);
44+
45+
psi0 = otp.quasigeostrophic.QuasiGeostrophicProblem.resize(psi0, [params.Nx, params.Ny]);
46+
psi0 = reshape(psi0, [], 1);
47+
48+
%% Do the rest
49+
oneday = 24*80/176251.2;
50+
51+
tspan = [100, 100 + oneday];
52+
53+
obj = obj@otp.quasigeostrophic.QuasiGeostrophicProblem(tspan, ...
54+
psi0, params);
55+
end
56+
end
57+
end

‎toolbox/+otp/+quasigeostrophic/QuasiGeostrophicParameters.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
classdef QuasiGeostrophicParameters < otp.Parameters
2-
% Parameters for the quasi-geostrophic problem.
3-
2+
% Parameters for the Quasi-geostrophic problem.
3+
44
properties
55
%Nx is the number of grid points in the x direction
66
Nx %MATLAB ONLY: (1,1) {mustBeInteger, mustBePositive} = 63

‎toolbox/+otp/+quasigeostrophic/QuasiGeostrophicProblem.m‎

Lines changed: 86 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,77 @@
11
classdef QuasiGeostrophicProblem < otp.Problem
2-
2+
% A chaotic PDE modeling the flow of a fluid on the earth.
3+
%
4+
% The governing partial differential equation that is discretized is,
5+
%
6+
% $$
7+
% Δψ_t = -J(ψ,ω) - {Ro}^{-1} \partial_x ψ -{Re}^{-1} Δω - {Ro}^{-1} F,
8+
% $$
9+
% where the Jacobian term is a quadratic function,
10+
% $$
11+
% J(ψ,ω) \equiv \partial_x ψ \partial_x ω - \partial_x ψ \partial_x ω,
12+
% $$
13+
% the relationship between the vorticity $ω$ and the stream function $ψ$ is
14+
% $$
15+
% ω = -Δψ,
16+
% $$
17+
% the term $Δ$ is the two dimensional Laplacian over the
18+
% discretization, the terms $\partial_x$ and $\partial_y$ are the first derivatives
19+
% in the $x$ and $y$ directions respectively, $Ro$ is the Rossby number, $Re$ is the Reynolds
20+
% number, and $F$ is a forcing term.
21+
% The spatial domain is fixed to $x ∈ [0, 1]$ and $y ∈ [0, 2]$, and
22+
% the boundary conditions of the PDE are assumed to be zero dirichlet
23+
% everywhere.
24+
%
25+
% A second order finite difference approximation is performed on the
26+
% grid to create the first derivative operators and the Laplacian
27+
% operator.
28+
%
29+
% The Laplacian is defined using the 5-point stencil
30+
% $$
31+
% Δ = \begin{bmatrix} & 1 & \\ 1 & -4 & 1\\ & 1 & \end{bmatrix},
32+
% $$
33+
% which is scaled with respect to the square of the step size in each
34+
% respective direction.
35+
% The first derivatives,
36+
% $$
37+
% \partial_x &= \begin{bmatrix} & 0 & \\ 1/2 & 0 & 1/2\\ & 0 & \end{bmatrix},\\
38+
% \partial_y &= \begin{bmatrix} & 1/2 & \\ 0 & 0 & 0\\ & 1/2 & \end{bmatrix},
39+
% $$
40+
% are the standard second order central finite difference operators in
41+
% the $x$ and $y$ directions.
42+
%
43+
% The Jacobian is discretized using the Arakawa approximation,
44+
%
45+
% $$
46+
% J(ψ,ω) = \frac{1}{3}[ψ_x ω_y - ψ_y ω_x + (ψ ω_y)_x - (ψ ω_x)_y + (ψ_x ω)_y - (ψ_y ω)_x],
47+
% $$
48+
%
49+
% in order for the system to not become unstable.
50+
%
51+
% The Poisson equation is solved by the eigenvalue sylvester method for
52+
% computational efficiency.
53+
%
54+
% An Approximate deconvolution large eddy simulation closure model from
55+
% :cite:p:`SSWZI11` is also implemented to have the same level of
56+
% accuracy with a coarser grid size.
57+
%
58+
% Notes
59+
% -----
60+
% +---------------------+-----------------------------------------------------------+
61+
% | Type | ODE |
62+
% +---------------------+-----------------------------------------------------------+
63+
% | Number of Variables | $Nx \times Ny$ |
64+
% +---------------------+-----------------------------------------------------------+
65+
% | Stiff | not typically, depending on $Re$, $Ro$, $Nx$, and $Ny$ |
66+
% +---------------------+-----------------------------------------------------------+
67+
%
68+
% Example
69+
% -------
70+
% >>> problem = otp.quasigeostrophic.presets.PopovMouSanduIliescu;
71+
% >>> sol = problem.solve();
72+
% >>> problem.movie(sol);
73+
%
74+
375
methods
476
function obj = QuasiGeostrophicProblem(timeSpan, y0, parameters)
577

@@ -28,10 +100,18 @@
28100

29101
methods (Static)
30102

31-
function u = resize(u, newsize)
32-
% resize uses interpolation to resize states
33-
34-
s = size(u);
103+
function psi = resize(psi, newsize)
104+
% Resizes the state onto a new grid by performing interpolation
105+
%
106+
% Parameters
107+
% ----------
108+
% psi : numeric(nx, ny)
109+
% the old state on the $x \times y$ grid.
110+
% newsize : numeric(1, 2)
111+
% the new size as a two-tuple $[nx, ny]$ indicating the new state of the system.
112+
%
113+
114+
s = size(psi);
35115

36116
X = linspace(0, 1, s(1) + 2);
37117
Y = linspace(0, 2, s(2) + 2).';
@@ -43,7 +123,7 @@
43123
Xnew = Xnew(2:end-1);
44124
Ynew = Ynew(2:end-1);
45125

46-
u = interp2(Y, X, u, Ynew, Xnew);
126+
psi = interp2(Y, X, psi, Ynew, Xnew);
47127

48128
end
49129

0 commit comments

Comments
 (0)