Flow velocity improvements in the new system

Author: AndreyAPopov

src/+otp/+qg/QuasiGeostrophicProblem.m

@@ -119,10 +119,6 @@ function onSettingsChanged(obj)
 
             bc = 'DD';
 
-            % create Laplacian
-            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.';
@@ -167,9 +163,9 @@ function onSettingsChanged(obj)
             %% Distance function, and flow velocity
             obj.DistanceFunction = @(t, y, i, j) otp.qg.distfn(t, y, i, j, nx, ny);
 
-            obj.FlowVelocityMagnitude = @(psi) otp.qg.flowvelmag(psi, Ddx, Ddy);
+            obj.FlowVelocityMagnitude = @(psi) otp.qg.flowvelmag(psi, Dx, Dy);
 
-            obj.JacobianFlowVelocityMagnitudeVectorProduct = @(psi, u) otp.qg.jacflowvelmagvp(psi, u, Ddx, Ddy);
+            obj.JacobianFlowVelocityMagnitudeVectorProduct = @(psi, u) otp.qg.jacflowvelmagvp(psi, u, Dx, Dy);
 
         end
 

src/+otp/+qg/flowvelmag.m

@@ -1,5 +1,11 @@
-function m = flowvelmag(psi, Ddx, Ddy)
+function m = flowvelmag(psi, Dx, Dy)
 
-m = sqrt((Ddx*psi).^2 + (Ddy*psi).^2);
+nx = size(Dx, 1);
+ny = size(Dy, 1);
+
+Dxpsi = reshape(Dx*reshape(psi, nx, []), nx, ny, []);
+Dypsi = permute(reshape(Dy*reshape(permute(reshape(psi, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
+
+m = reshape(sqrt((Dxpsi).^2 + (Dypsi).^2), nx*ny, []);
 
 end

src/+otp/+qg/jacflowvelmagvp.m

@@ -1,10 +1,20 @@
-function J = jacflowvelmagvp(psi, u, Ddx, Ddy)
+function J = jacflowvelmagvp(psi, v, Dx, Dy)
 
-m = sqrt((Ddx*psi).^2 + (Ddy*psi).^2);
+nx = size(Dx, 1);
+ny = size(Dy, 1);
 
-dpsix = Ddx*psi;
-dpsiy = Ddy*psi;
+N = size(v, 2);
 
-J = (1./m).*(dpsix.*(Ddx*u) + dpsiy.*(Ddy*u));
+psi = repmat(psi, 1, N);
+
+dpsix = reshape(Dx*reshape(psi, nx, []), nx, ny, []);
+dpsiy = permute(reshape(Dy*reshape(permute(reshape(psi, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
+
+dvx = reshape(Dx*reshape(v, nx, []), nx, ny, []);
+dvy = permute(reshape(Dy*reshape(permute(reshape(v, nx, ny, []), [2, 1, 3]), ny, []), ny, nx, []), [2, 1, 3]);
+
+m = sqrt((dpsix).^2 + (dpsiy).^2);
+
+J = reshape((1./m).*(dpsix.*(dvx) + dpsiy.*(dvy)), nx*ny, []);
 
 end