adding Jacobian Adjoint Vector Products in the new form

Author: AndreyAPopov

src/+otp/+qg/QuasiGeostrophicProblem.m

@@ -131,10 +131,12 @@ function onSettingsChanged(obj)
             Ddx = otp.utils.pde.Dd(n, xdomain, 1, 2, bc(1));
             Ddy = otp.utils.pde.Dd(n, ydomain, 2, 2, bc(2));
 
+            % create the spatial derivatives
             Dx = otp.utils.pde.Dd(nx, xdomain, 1, 1, bc(1));
+            DxT = Dx.';
 
-            % The transpose here is important
-            DyT = otp.utils.pde.Dd(ny, ydomain, 1, 1, bc(2)).';
+            Dy  = otp.utils.pde.Dd(ny, ydomain, 1, 1, bc(2));
+            DyT = Dy.';
 
             % create the x and y Laplacians
             Lx = otp.utils.pde.laplacian(nx, xdomain, 1, bc(1));
@@ -145,12 +147,11 @@ function onSettingsChanged(obj)
             RdnLT = RdnL.';
             PdnLT = PdnL.';
 
-            % do decompositions for the eigenvalue sylvester method. see
+            % Do decompositions for the eigenvalue sylvester method. See
             %
             %       Kirsten, Gerhardus Petrus. Comparison of methods for solving Sylvester systems. Stellenbosch University, 2018.
             %
-            % for a detailed method, the "Eigenvalue Method" which makes this
-            % particularly efficient
+            % for a detailed method, the "Eigenvalue Method" which makes this particularly efficient.
 
             hx2 = hx^2;
             hy2 = hy^2;
@@ -167,13 +168,13 @@ function onSettingsChanged(obj)
             F = sin(pi*(ymat.' - 1));
 
             obj.Rhs = otp.Rhs(@(t, psi) ...
-                otp.qg.f(psi, Lx, Ly, P1, P2, L12, Dx, DyT, F, Re, Ro), ...
+                otp.qg.f(psi, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro), ...
                 ...
                 otp.Rhs.FieldNames.JacobianVectorProduct, @(t, psi, v) ...
-                otp.qg.jvp(psi, v, Lx, Ly, P1, P2, L12, Dx, DyT, F, Re, Ro), ...
+                otp.qg.jvp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro), ...
                 ...
                 otp.Rhs.FieldNames.JacobianAdjointVectorProduct, @(t, psi, v) ...
-                otp.qg.javp(psi, v, L, RdnL, PdnL, Ddx, Ddy, Re, Ro));
+                otp.qg.javp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, F, Re, Ro));
 
 
             % AD LES

src/+otp/+qg/f.m

@@ -1,4 +1,4 @@
-function dpsit = f(psi, Lx, Ly, P1, P2, L12, Dx, DyT, F, Re, Ro)
+function dpsit = f(psi, Lx, Ly, P1, P2, L12, Dx, ~, ~, DyT, F, Re, Ro)
 
 [nx, ny] = size(L12);
 

src/+otp/+qg/javp.m

@@ -1,36 +1,53 @@
-function Jacobianavp = javp(psi, u, L, RdnL, PdnL, Ddx, Ddy, Re, Ro)
+function javp = javp(psi, v, Lx, Ly, P1, P2, L12, Dx, DxT, Dy, DyT, ~, Re, Ro)
 
-% supporting multiple mat-vecs at the same time
-N = size(u, 2);
-psi = repmat(psi, 1, N);
+[nx, ny] = size(L12);
 
-dpsix = Ddx*psi;
-dpsiy = Ddy*psi;
+psi = reshape(psi, nx, ny);
+v   = reshape(v,   nx, ny);
 
-q   = -L*psi;
-dqx = Ddx*q;
-dqy = Ddy*q;
+% Calculate the vorticity
+q = -(Lx*psi + psi*Ly);
 
-dnLu = PdnL*(RdnL\(RdnL'\(PdnL'*u)));
+dpsix = Dx*psi;
+dpsiy = psi*DyT;
+
+dqx = Dx*q;
+dqy = q*DyT;
+
+% solve the sylvester equation
+dnLv = P1*(L12.*(P1*v*P2))*P2;
+
+DxadnLv = DxT*dnLv;
+DyadnLv = dnLv*Dy;
+
+dpsixdnLvDy = (dpsix.*dnLv)*Dy;
+DxTdpsiydnLv = DxT*(dpsiy.*dnLv);
+
+psiDxadnLvDy = (psi.*DxadnLv)*Dy;
+DxTpsiDyadnLv = DxT*(psi.*DyadnLv);
+
+dpsixDyadnLv = dpsix.*DyadnLv
+dpsiyDxadnLv = dpsiy.*DxadnLv;
 
-DdxadnLu = Ddx'*dnLu;
-DdyadnLu = Ddy'*dnLu;
 
 % Arakawa approximation
-dJpsi1u = -L*(Ddy'*(dpsix.*dnLu))   + Ddx'*(dqy.*dnLu) ...
-    + L*(Ddx'*(dpsiy.*dnLu))       - Ddy'*(dqx.*dnLu);
+dJpsi1v = -(Lx*dpsixdnLvDy + dpsixdnLvDy*Ly)   + DxT*(dqy.*dnLv) ...
+    + (Lx*DxTdpsiydnLv + DxTdpsiydnLv*Ly)      - (dqx.*dnLv)*Dy;
 
-dJpsi2u = -L*(Ddy'*(psi.*DdxadnLu)) + dqy.*DdxadnLu ...
-    + L*(Ddx'*(psi.*DdyadnLu))     - dqx.*DdyadnLu;
+dJpsi2v = -(Lx*psiDxadnLvDy + psiDxadnLvDy*Ly) + dqy.*DxadnLv ...
+    + (Lx*DxTpsiDyadnLv + DxTpsiDyadnLv*Ly)    - dqx.*DyadnLv;
 
-dJpsi3u = -L*(dpsix.*DdyadnLu)   + Ddx'*(q.*DdyadnLu) ...
-    + L*(dpsiy.*DdxadnLu)       - Ddy'*(q.*DdxadnLu);
+dJpsi3v = -(Lx*dpsixDyadnLv + dpsixDyadnLv*Ly) + DxT*(q.*DyadnLv) ...
+    + (Lx*dpsiyDxadnLv + dpsiyDxadnLv*Ly)      - (q.*DxadnLv)*Dy;
 
 
-dJpsiu = -(dJpsi1u + dJpsi2u + dJpsi3u)/3;
+dJpsiv = -(dJpsi1v + dJpsi2v + dJpsi3v)/3;
 
-ddqtpsiau = - dJpsiu + (1/Ro)*DdxadnLu  + (1/Re)*(-L'*(L*dnLu));
+nLv = -(Lx*v + v*Ly);
 
-Jacobianavp = ddqtpsiau; 
+ddqtpsiav = -dJpsiv + (1/Ro)*DxadnLv  - (1/Re)*(nLv);
+
+javp = ddqtpsiav; 
 
 end
+

src/+otp/+qg/jvp.m

@@ -1,4 +1,4 @@
-function jvp = jvp(psi, v, Lx, Ly, P1, P2, L12, Dx, DyT, ~, Re, Ro)
+function jvp = jvp(psi, v, Lx, Ly, P1, P2, L12, Dx, ~, ~, DyT, ~, Re, Ro)
 
 [nx, ny] = size(L12);
 
@@ -14,31 +14,31 @@
 dqx = Dx*q;
 dqy = q*DyT;
 
-nLu = -(Lx*v + v*Ly);
-Dxu = Dx*v;
-Dyu = v*DyT;
+nLv = -(Lx*v + v*Ly);
+Dxv = Dx*v;
+Dyv = v*DyT;
 
-DxnLu = Dx*nLu;
-DynLu = nLu*DyT;
+DxnLv = Dx*nLv;
+DynLv = nLv*DyT;
 
 % Arakawa approximation
-dJpsi1u = dpsix.*DynLu     + dqy.*Dxu ...
-    - dpsiy.*DxnLu         - dqx.*Dyu;
+dJpsi1v = dpsix.*DynLv     + dqy.*Dxv ...
+    - dpsiy.*DxnLv         - dqx.*Dyv;
 
-dJpsi2u = Dx*(psi.*DynLu)  + Dx*(dqy.*v) ...
-    - (psi.*DxnLu)*DyT     - (dqx.*v)*DyT;
+dJpsi2v = Dx*(psi.*DynLv)  + Dx*(dqy.*v) ...
+    - (psi.*DxnLv)*DyT     - (dqx.*v)*DyT;
 
-dJpsi3u = (dpsix.*nLu)*DyT + (q.*Dxu)*DyT ...
-    - Dx*(dpsiy.*nLu)      - Dx*(q.*Dyu);
+dJpsi3v = (dpsix.*nLv)*DyT + (q.*Dxv)*DyT ...
+    - Dx*(dpsiy.*nLv)      - Dx*(q.*Dyv);
 
-dJpsiu = -(dJpsi1u + dJpsi2u + dJpsi3u)/3;
+dJpsiv = -(dJpsi1v + dJpsi2v + dJpsi3v)/3;
 
-ddqtpsivp = -dJpsiu + (1/Ro)*Dxu;
+ddqtpsivp = -dJpsiv + (1/Ro)*Dxv;
 
 % solve the sylvester equation
 nLidqtmq = P1*(L12.*(P1*ddqtpsivp*P2))*P2;
 
 % solve into stream form of the Jacobian vp
-jvp = reshape(nLidqtmq - (1/Re)*(nLu), nx*ny, 1);
+jvp = reshape(nLidqtmq - (1/Re)*(nLv), nx*ny, 1);
 
 end