Meshfree Shallow Water on a Sphere

src/+otp/+mfshallowwatersphere/+presets/Canonical.m

@@ -0,0 +1,78 @@
+classdef Canonical < otp.mfshallowwatersphere.MFShallowWaterSphereProblem
+    methods
+        function obj = Canonical(varargin)
+
+            load('nodes500.mat', 'x', 'y', 'z');
+
+            % mean water height
+            H = 5.768e4;
+            % earth gravity
+            g = otp.utils.PhysicalConstants.EarthGravity;
+            % radius of the earth
+            a = otp.utils.PhysicalConstants.EarthRadius;
+            % initial velocity for the perturbation
+            u0 = 20;
+            % Angular velocity of the earth
+            Omega = otp.utils.PhysicalConstants.EarthAngularVelocity;
+
+            coriolisForce = 2*Omega*z;
+
+            params = struct;
+            params.gravity = otp.utils.PhysicalConstants.EarthGravity;
+            params.radius = a;
+            params.angularSpeed = Omega;
+            params.coriolisForce = coriolisForce;
+
+            params.x = x;
+            params.y = y;
+            params.z = z;
+
+            params.rbfradius = 1.5;
+            params.rbf = @otp.utils.rbf.buhmann3;
+
+            % convert from Cartesian to spherical coordinates
+            theta = atan2(z, sqrt(x.^2 + y.^2));
+            lambda = atan2(y, x);
+
+            %% Define the initial conditions
+            % first get a stableish Rossby-Haurwitz wave with a wave number
+            % of R = 7
+            R = 7;
+            [h, zonalwind, meridionalwind] = getrossbyhaurwitzwave(x, y, z, Omega, a, g, R);
+
+            % then, define perturbations to this wave in terms of the T-Z
+            % initial condition.
+            hpert              = (1/g)*(H + 2*Omega*a*u0*( sin(theta).^3 ).*cos(theta).*sin(lambda));
+            zonalwindpert      = (-3*u0*sin(theta).*( cos(theta).^2 ).*sin(lambda) + u0*( sin(theta).^3 ).*sin(lambda));
+            meridionalwindpert = (u0*( sin(theta).^2 ).*cos(lambda));
+
+            % remove the excess height and velocity
+            %hpert = hpert - mean(hpert);
+            %zonalwindpert = zonalwindpert - mean(zonalwindpert);
+            %meridionalwindpert = meridionalwindpert - mean(meridionalwindpert);
+
+            % mixing coefficient
+            alpha = 0.5;
+
+            % perturb the R-H wave
+            h              = alpha*h              + (1 - alpha)*hpert;
+            zonalwind      = alpha*zonalwind      + (1 - alpha)*zonalwindpert;
+            meridionalwind = alpha*meridionalwind + (1 - alpha)*meridionalwindpert;
+
+            % finally, convert to Cartesian coordinates
+            [u, v, w] = velocitytocartesian(x, y, z, zonalwind, meridionalwind);
+
+            huv0 = [h; u; v; w];
+
+            %% Do the rest
+
+            oneday = 24*60*60;
+
+            tspan = [0, oneday];
+
+            obj = [email protected](tspan, ...
+                huv0, params);
+
+        end
+    end
+end

src/+otp/+mfshallowwatersphere/+presets/private/getrossbyhaurwitzwave.m

@@ -0,0 +1,33 @@
+% See
+%   Williamson, David L., John B. Drake, James J. Hack, Rüdiger Jakob, and Paul N. Swarztrauber. 
+%   "A standard test set for numerical approximations to the shallow water equations in spherical geometry." 
+%   Journal of Computational Physics 102, no. 1 (1992): 211-224.
+function [h, zonalwind, meridionalwind] = getrossbyhaurwitzwave(x, y, z, Omega, a, g, wavenumber)
+
+if nargin < 7
+    wavenumber = 4;
+end
+
+omega = 7.848e-6;
+K = omega;
+h0 = 8.0e3;
+R = wavenumber;
+
+% convert from Cartesian to spherical coordinates
+theta = atan2(z, sqrt(x.^2 + y.^2));
+lambda = atan2(y, x);
+
+zonalwind = a*omega*cos(theta) + a*K*( cos(theta).^(R-1) ).*( R*(sin(theta).^2) - (cos(theta).^2)).*cos(R*lambda);
+meridionalwind = -a*K*R*( cos(theta).^(R-1) ).*sin(theta).*sin(R*lambda);
+
+At = (R + 1)*( cos(theta).^2 ) + (2*(R^2) - R - 2) - (2*(R^2))./( cos(theta).^2 );
+A = (omega/2)*(2*Omega + omega)*( cos(theta).^2 ) + ((K^2)/4)*( cos(theta).^(2*R) ).*At;
+
+Bt = R^2 + 2*R + 2 - ((R + 1).^2).*( cos(theta).^2 );
+B = (2*(Omega + omega)*K)/((R+1)*(R+2)).*( cos(theta).^R ).*Bt;
+
+C = ((K^2)/4)*( cos(theta).^(2*R) ).*( (R+1)*( cos(theta).^2 ) - (R+2));
+
+h = (1/g)*(h0 + (a^2)*( A + B.*cos(R*lambda) + C.*cos(2*R*lambda) ));
+
+end

src/+otp/+mfshallowwatersphere/+presets/private/velocitytocartesian.m

@@ -0,0 +1,12 @@
+function [u, v, w] = velocitytocartesian(x, y, z, zonalwind, meridionalwind)
+
+% convert from Cartesian to spherical coordinates
+theta = atan2(z, sqrt(x.^2 + y.^2));
+lambda = atan2(y, x);
+
+% convert from spherical to Cartesian velocity
+u = (-zonalwind.*sin(lambda).*cos(theta) - meridionalwind.*cos(lambda).*sin(theta));
+v = (zonalwind.*cos(lambda).*cos(theta) - meridionalwind.*sin(lambda).*sin(theta));
+w = (meridionalwind.*cos(theta));
+
+end

src/+otp/+mfshallowwatersphere/MFShallowWaterSphereProblem.m

@@ -0,0 +1,173 @@
+classdef MFShallowWaterSphereProblem < otp.Problem
+
+    methods
+        function obj = MFShallowWaterSphereProblem(timeSpan, y0, parameters)
+
+            [email protected]('Meshfree Shallow Water on a Sphere', [], ...
+                timeSpan, y0, parameters);
+
+        end
+    end
+
+    properties (SetAccess = private)
+        DistanceFunction
+    end
+
+    methods
+
+        function plotSphere(obj, huv, projection)
+            if nargin < 3
+                projection = 'eqaazim';
+            end
+
+
+            x = obj.Parameters.x;
+            y = obj.Parameters.y;
+            z = obj.Parameters.z;
+            rbf = obj.Parameters.rbf;
+
+            Nplot = 50;
+            lambdaplot = linspace(-pi, pi, Nplot);
+            thetaplot = linspace(-pi/2, pi/2, Nplot);
+            [lambdainterpgrid, thetainterpgrid] = meshgrid(lambdaplot, thetaplot);
+
+            % get the cartesian coordinates for the uniform mesh
+            [x2, y2, z2] = sph2cart(lambdainterpgrid(:), thetainterpgrid(:), ones(numel(lambdainterpgrid), 1));
+            radiusplot = 0.5;
+
+            %
+            Winterp = rbfinterp(x, y, z, x2, y2, z2, radiusplot, rbf);
+            Winterp = Winterp./sum(Winterp, 2);
+
+            lon = 360*(lambdainterpgrid/pi + 1)/2;
+            lat = 180*(thetainterpgrid/(pi/2))/2;
+
+
+            load('coastlines', 'coastlat', 'coastlon');
+
+            n = size(huv, 1)/4;
+            h = huv(1:n);
+            %u = huv((n+1):(2*n));
+            %v = huv((2*n+1):end);
+
+            Nplot = sqrt(size(Winterp, 1));
+
+            cmap = interp1([0; 0.5; 1], [1, 0, 0; 1, 1, 1; 0, 0.3, 0.8], linspace(0, 1, 500));
+            levels = 20;
+
+            %colormap(cmap);
+
+            %subplot(1, 3, 1); 
+            cla;
+            axesm(projection);
+            contourfm(lat, lon, reshape(Winterp*h, Nplot, Nplot), levels, 'LineStyle','none');
+            colorbar;
+            ax = gca;
+            setm(ax,'FLineWidth', 3, 'Grid','on')
+            l = plotm(coastlat, coastlon, '-k');
+            l.Color = [l.Color, 0.5];
+
+            %             subplot(1, 3, 2); cla;
+            %             axesm(projection);
+            %             contourfm(lat, lon, reshape(Winterp*u, Nplot, Nplot), levels, 'LineStyle','none');
+            %             colorbar;
+            %             ax = gca;
+            %             setm(ax,'FLineWidth', 3, 'Grid','on')
+            %             l = plotm(coastlat, coastlon, '-k');
+            %             l.Color = [l.Color, 0.5];
+            %
+            %             subplot(1, 3, 3); cla;
+            %             axesm(projection);
+            %             contourfm(lat, lon, reshape(Winterp*v, Nplot, Nplot), levels, 'LineStyle','none');
+            %             colorbar;
+            %             ax = gca;
+            %             setm(ax,'FLineWidth', 3, 'Grid','on')
+            %             l = plotm(coastlat, coastlon, '-k');
+            %             l.Color = [l.Color, 0.5];
+
+            drawnow;
+
+
+        end
+
+
+    end
+
+
+    methods (Access = protected)
+
+        function onSettingsChanged(obj)            
+            x = obj.Parameters.x;
+            y = obj.Parameters.y;
+            z = obj.Parameters.z;
+            g = obj.Parameters.gravity;
+            a = obj.Parameters.radius;
+            rbf = obj.Parameters.rbf;
+            rbfradius = obj.Parameters.rbfradius;
+            f = obj.Parameters.coriolisForce;
+
+            % create the interpolation matrix and derivatives
+            [W, dWdx, dWdy, dWdz] = rbfinterp(x, y, z, x, y, z, rbfradius, rbf);
+
+            Bx = dWdx.*(1 - x.^2) + dWdy.*(-x.*y)    + dWdz.*(-x.*z);
+            By = dWdx.*(-x.*y)    + dWdy.*(1 - y.^2) + dWdz.*(-y.*z);
+            Bz = dWdx.*(-x.*z)    + dWdy.*(-y.*z)    + dWdz.*(1 - z.^2);
+
+            Bx = Bx/a;
+            By = By/a;
+            Bz = Bz/a;
+
+            try 
+                Wdecomp = decomposition(W, 'chol');
+            catch
+                error('The selected compination of RBF, radius, and nodes did not result in a SPD collocation matrix.')
+            end
+
+            % build Shepard interpolation matrix
+            S = W./sum(W, 2);
+
+            % set the right hand side
+            obj.Rhs = otp.Rhs(@(t, huvw) ...
+                otp.mfshallowwatersphere.f(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z));
+
+
+            %% Distance function
+
+            theta = atan2(z, sqrt(x.^2 + y.^2));
+            lambda = atan2(y, x);
+
+            obj.DistanceFunction = @(t, huv, i, j) otp.mfshallowwatersphere.distfn(t, huv, i, j, theta, lambda);
+
+        end
+
+        function validateNewState(obj, newTimeSpan, newY0, newParameters)
+
+            [email protected](obj, ...
+                newTimeSpan, newY0, newParameters)
+
+            %otp.utils.StructParser(newParameters) ...
+            %    .checkField('nx', 'finite', 'scalar', 'integer', 'positive') ...
+            %    .checkField('ny', 'finite', 'scalar', 'integer', 'positive') ...
+            %    .checkField('Re', 'finite', 'scalar', 'real', 'positive') ...
+            %    .checkField('Ro', 'finite', 'scalar', 'real');
+
+        end
+
+        function label = internalIndex2label(obj, index)
+
+
+            label = [];
+
+            %[i, j] = ind2sub([obj.Parameters.nx, obj.Parameters.ny], index);
+
+            %label = sprintf('(%d, %d)', i, j);
+
+        end
+
+        function sol = internalSolve(obj, varargin)
+            % This really requires an SSP method
+            sol = [email protected](obj, 'Method', @ode45, varargin{:});
+        end
+
+    end
+end

src/+otp/+mfshallowwatersphere/distfn.m

@@ -0,0 +1,20 @@
+function d = distfn(~, ~, i, j, theta, lambda)
+
+n = numel(lambda);
+
+i = mod(i - 1, n) + 1;
+j = mod(j - 1, n) + 1;
+
+li = lambda(i);
+pi = theta(i);
+
+lj = lambda(j);
+pj = theta(j);
+
+dell = li - lj;
+
+nom = sqrt(  ( cos(pj).*sin(dell) ).^2 + ( cos(pi).*sin(pj) - sin(pi).*cos(pj).*cos(dell) ).^2 );
+den = sin(pi).*sin(pj) + cos(pi).*cos(pj).*cos(dell);
+d = atan2( nom, den ).';
+
+end

src/+otp/+mfshallowwatersphere/f.m

@@ -0,0 +1,47 @@
+% See
+%   Flyer, Natasha, and Grady B. Wright. 
+%   "A radial basis function method for the shallow water equations on a sphere." 
+%   Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 465, no. 2106 (2009): 1949-1976.
+
+function dhuvwdt = f(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z)
+
+n = size(huvw, 1)/4;
+
+h = huvw(1:n, :);
+u = huvw((n+1):(2*n), :);
+v = huvw((2*n+1):(3*n), :);
+w = huvw((3*n+1):end, :);
+
+ch = Wdecomp\h;
+cu = Wdecomp\u;
+cv = Wdecomp\v;
+cw = Wdecomp\w;
+
+Bxch = Bx*ch;
+Bych = By*ch;
+Bzch = Bz*ch;
+
+Bxcu = Bx*cu;
+Bycv = By*cv;
+Bzcw = Bz*cw;
+
+rhsDx = u.*(Bxcu) + v.*(By*cu) + w.*(Bz*cu) + f.*(y.*w - z.*v) + g*Bxch;
+rhsDy = u.*(Bx*cv) + v.*(Bycv) + w.*(Bz*cv) + f.*(z.*u - x.*w) + g*Bych;
+rhsDz = u.*(Bx*cw) + v.*(By*cw) + w.*(Bzcw) + f.*(x.*v - y.*u) + g*Bzch;
+
+dudt = -( rhsDx.*(1 - x.^2) + rhsDy.*(-x.*y) + rhsDz.*(-x.*z) );
+dvdt = -( rhsDx.*(-x.*y) + rhsDy.*(1 - y.^2) + rhsDz.*(-y.*z) );
+dwdt = -( rhsDx.*(-x.*z) + rhsDy.*(-y.*z) + rhsDz.*(1 - z.^2) );
+
+dhdt = -( u.*Bxch + v.*Bych + w.*Bzch + h.*(Bxcu + Bycv + Bzcw) );
+
+% Shepard interpolation for stability
+dhdt = S*dhdt;
+dudt = S*dudt;
+dvdt = S*dvdt;
+dwdt = S*dwdt;
+
+dhuvwdt = [dhdt; dudt; dvdt; dwdt];
+
+end
+