Skip to content

Commit 83ac84f

Browse files
AndreyAPopovAndreyAPopov
authored
Merge pull request #16 from ComputationalScienceLaboratory/QG-Sylvester
QG Rhs optimization
2 parents c487b3d + 1796e6f commit 83ac84f

File tree

9 files changed

+140
-152
lines changed

9 files changed

+140
-152
lines changed

‎src/+otp/+qg/+presets/Canonical.m‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -9,22 +9,17 @@
99
addParameter(p, 'ReynoldsNumber', Re);
1010
addParameter(p, 'RossbyNumber', Ro);
1111
addParameter(p, 'Size', 'huge');
12-
addParameter(p, 'ADLambda', 0.4);
13-
addParameter(p, 'ADPasses', 4);
1412

1513
parse(p, varargin{:});
1614

1715
s = p.Results;
1816

1917
[nx, ny] = otp.qg.QuasiGeostrophicProblem.name2size(s.Size);
2018

21-
les = struct('lambda', s.ADLambda, 'passes', s.ADPasses);
22-
2319
params.nx = nx;
2420
params.ny = ny;
2521
params.Re = s.ReynoldsNumber;
2622
params.Ro = s.RossbyNumber;
27-
params.les = les;
2823

2924
%% Construct initial conditions
3025

‎src/+otp/+qg/+presets/PopovMouIliescuSandu.m‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@
22
methods
33
function obj = PopovMouIliescuSandu(varargin)
44

5-
65
defaultsize = 'huge';
76

87
Re = 450;
@@ -12,22 +11,17 @@
1211
addParameter(p, 'ReynoldsNumber', Re);
1312
addParameter(p, 'RossbyNumber', Ro);
1413
addParameter(p, 'Size', defaultsize);
15-
addParameter(p, 'ADLambda', 0.4);
16-
addParameter(p, 'ADPasses', 4);
1714

1815
parse(p, varargin{:});
1916

2017
s = p.Results;
2118

2219
[nx, ny] = otp.qg.QuasiGeostrophicProblem.name2size(s.Size);
2320

24-
les = struct('lambda', s.ADLambda, 'passes', s.ADPasses);
25-
2621
params.nx = nx;
2722
params.ny = ny;
2823
params.Re = s.ReynoldsNumber;
2924
params.Ro = s.RossbyNumber;
30-
params.les = les;
3125

3226
%% Load initial conditions
3327

‎src/+otp/+qg/QuasiGeostrophicProblem.m‎

Lines changed: 34 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,10 @@
99
end
1010
end
1111

12-
properties
12+
properties (SetAccess = private)
1313
DistanceFunction
1414
FlowVelocityMagnitude
1515
JacobianFlowVelocityMagnitudeVectorProduct
16-
RhsAD
1716
end
1817

1918
methods (Static)
@@ -110,72 +109,63 @@ function onSettingsChanged(obj)
110109
Re = obj.Parameters.Re;
111110
Ro = obj.Parameters.Ro;
112111

113-
lambda = obj.Parameters.les.lambda;
114-
passes = obj.Parameters.les.passes;
115-
%dfiltertype = obj.Parameters.les.filtertype;
116-
117112
n = [nx, ny];
118113

119-
h = 1/(nx + 1);
114+
hx = 1/(nx + 1);
115+
hy = 2/(ny + 1);
120116

121117
xdomain = [0, 1];
122118
ydomain = [0, 2];
123119

124-
domain = [0, 1; 0, 2];
125-
diffc = [1, 1];
126120
bc = 'DD';
127121

128-
% create laplacian
129-
L = otp.utils.pde.laplacian(n, domain, diffc, bc);
130-
Ddx = otp.utils.pde.Dd(n, xdomain, 1, 2, bc(1));
131-
Ddy = otp.utils.pde.Dd(n, ydomain, 2, 2, bc(2));
122+
% create the spatial derivatives
123+
Dx = otp.utils.pde.Dd(nx, xdomain, 1, 1, bc(1));
124+
DxT = Dx.';
125+
126+
Dy = otp.utils.pde.Dd(ny, ydomain, 1, 1, bc(2));
127+
DyT = Dy.';
128+
129+
% create the x and y Laplacians
130+
Lx = otp.utils.pde.laplacian(nx, xdomain, 1, bc(1));
131+
Ly = otp.utils.pde.laplacian(ny, ydomain, 1, bc(2));
132132

133-
% do a Cholesky decomposition on the negative laplacian
134-
[RdnL, ~, PdnL] = chol(-L);
135-
RdnLT = RdnL.';
136-
PdnLT = PdnL.';
133+
% Do decompositions for the eigenvalue sylvester method. See
134+
%
135+
% Kirsten, Gerhardus Petrus. Comparison of methods for solving Sylvester systems. Stellenbosch University, 2018.
136+
%
137+
% for a detailed method, the "Eigenvalue Method" which makes this particularly efficient.
138+
139+
hx2 = hx^2;
140+
hy2 = hy^2;
141+
cx = (1:nx)/(nx + 1);
142+
cy = (1:ny)/(ny + 1);
143+
L12 = -(hx2*hy2./(2*(-hx2 - hy2 + hy2*cospi(cx).' + hx2*cospi(cy))));
144+
P1 = sqrt(2/(nx + 1))*sinpi((1:nx).'*(1:nx)/(nx + 1));
145+
P2 = sqrt(2/(ny + 1))*sinpi((1:ny).'*(1:ny)/(ny + 1));
137146

138147
ys = linspace(ydomain(1), ydomain(end), ny + 2);
139148
ys = ys(2:end-1);
140149

141150
ymat = repmat(ys.', 1, nx);
142-
ymat = reshape(ymat.', prod(n), 1);
143-
F = sin(pi*(ymat - 1));
151+
F = sin(pi*(ymat.' - 1));
144152

145153
obj.Rhs = otp.Rhs(@(t, psi) ...
146-
otp.qg.f(psi, L, RdnL, RdnLT, PdnL, PdnLT, Ddx, Ddy, F, Re, Ro), ...
154+
otp.qg.f(psi, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro), ...
147155
...
148-
otp.Rhs.FieldNames.JacobianVectorProduct, @(t, psi, u) ...
149-
otp.qg.jvp(psi, u, L, RdnL, PdnL, Ddx, Ddy, Re, Ro), ...
156+
otp.Rhs.FieldNames.JacobianVectorProduct, @(t, psi, v) ...
157+
otp.qg.jvp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro), ...
150158
...
151-
otp.Rhs.FieldNames.JacobianAdjointVectorProduct, @(t, psi, u) ...
152-
otp.qg.javp(psi, u, L, RdnL, PdnL, Ddx, Ddy, Re, Ro));
159+
otp.Rhs.FieldNames.JacobianAdjointVectorProduct, @(t, psi, v) ...
160+
otp.qg.javp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro));
153161

154162

155-
% AD LES
156-
fmat = speye(prod(n)) - ((lambda*h)^2)*L;
157-
158-
[Rfmat, ~, Pfmat] = chol(fmat);
159-
RfmatT = Rfmat.';
160-
PfmatT = Pfmat.';
161-
162-
if ~isfield(obj.Parameters, 'filter') || isempty(obj.Parameters.filter)
163-
filter = @(u) Pfmat*(Rfmat\(RfmatT\(PfmatT*u)));
164-
else
165-
filter = obj.Parameters.filter;
166-
end
167-
168-
Fbar = filter(F);
169-
obj.RhsAD = otp.Rhs(@(t, psi) ...
170-
otp.qg.fAD(psi, L, RdnL, RdnLT, PdnL, PdnLT, Ddx, Ddy, Fbar, Re, Ro, filter, passes));
171-
172-
173163
%% Distance function, and flow velocity
174164
obj.DistanceFunction = @(t, y, i, j) otp.qg.distfn(t, y, i, j, nx, ny);
175165

176-
obj.FlowVelocityMagnitude = @(psi) otp.qg.flowvelmag(psi, Ddx, Ddy);
166+
obj.FlowVelocityMagnitude = @(psi) otp.qg.flowvelmag(psi, Dx, Dy);
177167

178-
obj.JacobianFlowVelocityMagnitudeVectorProduct = @(psi, u) otp.qg.jacflowvelmagvp(psi, u, Ddx, Ddy);
168+
obj.JacobianFlowVelocityMagnitudeVectorProduct = @(psi, u) otp.qg.jacflowvelmagvp(psi, u, Dx, Dy);
179169

180170
end
181171

‎src/+otp/+qg/f.m‎

Lines changed: 16 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,32 @@
1-
function dpsit = f(psi, L, RdnL, RdnLT, PdnL, PdnLT, Ddx, Ddy, F, Re, Ro)
1+
function dpsit = f(psi, Lx, Ly, P1, P2, L12, Dx, ~, ~, DyT, F, Re, Ro)
2+
3+
[nx, ny] = size(L12);
4+
5+
psi = reshape(psi, nx, ny);
26

37
% Calculate the vorticity
4-
q = -(L*psi);
8+
q = -(Lx*psi + psi*Ly);
59

610
% calculate Arakawa
7-
dpsix = Ddx*psi;
8-
dpsiy = Ddy*psi;
11+
dpsix = Dx*psi;
12+
dpsiy = psi*DyT;
913

10-
dqx = Ddx*q;
11-
dqy = Ddy*q;
14+
dqx = Dx*q;
15+
dqy = q*DyT;
1216

1317
J1 = dpsix.*dqy - dpsiy.*dqx;
14-
J2 = Ddx*(psi.*dqy) - Ddy*(psi.*dqx);
15-
J3 = Ddy*(q.*dpsix) - Ddx*(q.*dpsiy);
18+
J2 = Dx*(psi.*dqy) - (psi.*dqx)*DyT;
19+
J3 = (q.*dpsix)*DyT - Dx*(q.*dpsiy);
1620

1721
mJ = (J1 + J2 + J3)/3;
1822

1923
% almost vorticity form of the rhs
2024
dqtmq = mJ + (1/Ro)*(dpsix) + (1/Ro)*F;
2125

26+
% solve the sylvester equation
27+
nLidqtmq = P1*(L12.*(P1*dqtmq*P2))*P2;
28+
2229
% solve into stream form of the rhs
23-
dpsit = PdnL*(RdnL\(RdnLT\(PdnLT*dqtmq))) - (1/Re)*(q);
30+
dpsit = reshape(nLidqtmq - (1/Re)*(q), nx*ny, 1);
2431

2532
end

‎src/+otp/+qg/fAD.m‎

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

‎src/+otp/+qg/flowvelmag.m‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,11 @@
1-
function m = flowvelmag(psi, Ddx, Ddy)
1+
function m = flowvelmag(psi, Dx, Dy)
22

3-
m = sqrt((Ddx*psi).^2 + (Ddy*psi).^2);
3+
nx = size(Dx, 1);
4+
ny = size(Dy, 1);
5+
6+
Dxpsi = reshape(Dx*reshape(psi, nx, []), nx, ny, []);
7+
Dypsi = permute(reshape(Dy*reshape(permute(reshape(psi, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
8+
9+
m = reshape(sqrt((Dxpsi).^2 + (Dypsi).^2), nx*ny, []);
410

511
end

‎src/+otp/+qg/jacflowvelmagvp.m‎

Lines changed: 15 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,20 @@
1-
function J = jacflowvelmagvp(psi, u, Ddx, Ddy)
1+
function J = jacflowvelmagvp(psi, v, Dx, Dy)
22

3-
m = sqrt((Ddx*psi).^2 + (Ddy*psi).^2);
3+
nx = size(Dx, 1);
4+
ny = size(Dy, 1);
45

5-
dpsix = Ddx*psi;
6-
dpsiy = Ddy*psi;
6+
N = size(v, 2);
77

8-
J = (1./m).*(dpsix.*(Ddx*u) + dpsiy.*(Ddy*u));
8+
psi = repmat(psi, 1, N);
9+
10+
dpsix = reshape(Dx*reshape(psi, nx, []), nx, ny, []);
11+
dpsiy = permute(reshape(Dy*reshape(permute(reshape(psi, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
12+
13+
dvx = reshape(Dx*reshape(v, nx, []), nx, ny, []);
14+
dvy = permute(reshape(Dy*reshape(permute(reshape(v, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
15+
16+
m = sqrt((dpsix).^2 + (dpsiy).^2);
17+
18+
J = reshape((1./m).*(dpsix.*(dvx) + dpsiy.*(dvy)), nx*ny, []);
919

1020
end

‎src/+otp/+qg/javp.m‎

Lines changed: 38 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,36 +1,53 @@
1-
function Jacobianavp = javp(psi, u, L, RdnL, PdnL, Ddx, Ddy, Re, Ro)
1+
function javp = javp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, ~, Re, Ro)
22

3-
% supporting multiple mat-vecs at the same time
4-
N = size(u, 2);
5-
psi = repmat(psi, 1, N);
3+
[nx, ny] = size(L12);
64

7-
dpsix = Ddx*psi;
8-
dpsiy = Ddy*psi;
5+
psi = reshape(psi, nx, ny);
6+
v = reshape(v, nx, ny);
97

10-
q = -L*psi;
11-
dqx = Ddx*q;
12-
dqy = Ddy*q;
8+
% Calculate the vorticity
9+
q = -(Lx*psi + psi*Ly);
1310

14-
dnLu = PdnL*(RdnL\(RdnL'\(PdnL'*u)));
11+
dpsix = Dx*psi;
12+
dpsiy = psi*DyT;
13+
14+
dqx = Dx*q;
15+
dqy = q*DyT;
16+
17+
% solve the sylvester equation
18+
dnLv = P1*(L12.*(P1*v*P2))*P2;
19+
20+
DxadnLv = DxT*dnLv;
21+
DyadnLv = dnLv*Dy;
22+
23+
dpsixdnLvDy = (dpsix.*dnLv)*Dy;
24+
DxTdpsiydnLv = DxT*(dpsiy.*dnLv);
25+
26+
psiDxadnLvDy = (psi.*DxadnLv)*Dy;
27+
DxTpsiDyadnLv = DxT*(psi.*DyadnLv);
28+
29+
dpsixDyadnLv = dpsix.*DyadnLv
30+
dpsiyDxadnLv = dpsiy.*DxadnLv;
1531

16-
DdxadnLu = Ddx'*dnLu;
17-
DdyadnLu = Ddy'*dnLu;
1832

1933
% Arakawa approximation
20-
dJpsi1u = -L*(Ddy'*(dpsix.*dnLu)) + Ddx'*(dqy.*dnLu) ...
21-
+ L*(Ddx'*(dpsiy.*dnLu)) - Ddy'*(dqx.*dnLu);
34+
dJpsi1v = -(Lx*dpsixdnLvDy + dpsixdnLvDy*Ly) + DxT*(dqy.*dnLv) ...
35+
+ (Lx*DxTdpsiydnLv + DxTdpsiydnLv*Ly) - (dqx.*dnLv)*Dy;
2236

23-
dJpsi2u = -L*(Ddy'*(psi.*DdxadnLu)) + dqy.*DdxadnLu ...
24-
+ L*(Ddx'*(psi.*DdyadnLu)) - dqx.*DdyadnLu;
37+
dJpsi2v = -(Lx*psiDxadnLvDy + psiDxadnLvDy*Ly) + dqy.*DxadnLv ...
38+
+ (Lx*DxTpsiDyadnLv + DxTpsiDyadnLv*Ly) - dqx.*DyadnLv;
2539

26-
dJpsi3u = -L*(dpsix.*DdyadnLu) + Ddx'*(q.*DdyadnLu) ...
27-
+ L*(dpsiy.*DdxadnLu) - Ddy'*(q.*DdxadnLu);
40+
dJpsi3v = -(Lx*dpsixDyadnLv + dpsixDyadnLv*Ly) + DxT*(q.*DyadnLv) ...
41+
+ (Lx*dpsiyDxadnLv + dpsiyDxadnLv*Ly) - (q.*DxadnLv)*Dy;
2842

2943

30-
dJpsiu = -(dJpsi1u + dJpsi2u + dJpsi3u)/3;
44+
dJpsiv = -(dJpsi1v + dJpsi2v + dJpsi3v)/3;
3145

32-
ddqtpsiau = - dJpsiu + (1/Ro)*DdxadnLu + (1/Re)*(-L'*(L*dnLu));
46+
nLv = -(Lx*v + v*Ly);
3347

34-
Jacobianavp = ddqtpsiau;
48+
ddqtpsiav = -dJpsiv + (1/Ro)*DxadnLv - (1/Re)*(nLv);
49+
50+
javp = ddqtpsiav;
3551

3652
end
53+

0 commit comments

Comments
 (0)