Meshfree Shallow Water on a Sphere

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

@@ -0,0 +1,47 @@
+classdef Canonical < otp.mfshallowwatersphere.MFShallowWaterSphereProblem
+    methods
+        function obj = Canonical(varargin)
+
+            load('mesh4000.mat', 'lambda', 'phi');
+
+            % mean water height
+            H = 5.768e4;
+            % earth gravity
+            g = 9.8;
+            % radius of the earth
+            a = 6.370e6;
+            % initial velocity
+            u0 = 20;
+            % Angular speed of the earth
+            Omega = 7.292e-5;
+
+
+            params = struct;
+            params.gravity = otp.utils.PhysicalConstants.EarthGravity;
+            params.radius = a;
+            params.angularSpeed = Omega;
+
+            params.lambda = lambda;
+            params.phi = phi;
+
+            phiT = phi.';
+            lambdaT = lambda.';
+
+            h = (1/g)*(H + 2*Omega*a*u0*( sin(phiT).^3 ).*cos(phiT).*sin(lambdaT));
+            u = -3*u0*sin(phiT).*( cos(phiT).^2 ).*sin(lambdaT) + u0*( sin(phiT).^3 ).*sin(lambdaT);
+            v = u0*( sin(phiT).^2 ).*cos(lambdaT);
+
+            huv0 = [h; u; v];
+
+            %% Do the rest
+
+            oneday = 24*60*60;
+
+            tspan = [0, oneday];
+
+            obj = obj@otp.mfshallowwatersphere.MFShallowWaterSphereProblem(tspan, ...
+                huv0, params);
+
+        end
+    end
+end

src/+otp/+mfshallowwatersphere/MFShallowWaterSphereProblem.m

@@ -0,0 +1,156 @@
+classdef MFShallowWaterSphereProblem < otp.Problem
+
+    methods
+        function obj = MFShallowWaterSphereProblem(timeSpan, y0, parameters)
+
+            obj@otp.Problem('Meshfree Shallow Water on a Sphere', [], ...
+                timeSpan, y0, parameters);
+
+        end
+    end
+
+    properties (SetAccess = private)
+        DistanceFunction
+    end
+
+    properties (Access = private)
+        PlottingInterp
+        PlottingLatitude
+        PlottingLongitude
+    end
+
+    methods
+
+        function plotSphere(obj, huv, projection)
+            if nargin < 3
+                projection = 'eqaazim';
+            end
+
+            load('coastlines', 'coastlat', 'coastlon');
+
+            n = size(huv, 1)/3;
+            h = huv(1:n);
+            u = huv((n+1):(2*n));
+            v = huv((2*n+1):end);
+
+            Winterp = obj.PlottingInterp;
+
+            Nplot = sqrt(size(Winterp, 1));
+
+            lat = obj.PlottingLatitude;
+            lon = obj.PlottingLongitude;
+
+            cmap = interp1([0; 0.5; 1], [1, 0, 0; 1, 1, 1; 0, 0.3, 0.8], linspace(0, 1, 500));
+
+            colormap(cmap);
+
+            subplot(1, 3, 1); cla;
+            axesm(projection);
+            contourfm(lat,lon,reshape(Winterp*h, Nplot, Nplot),'LineStyle','none')
+            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),'LineStyle','none')
+            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),'LineStyle','none')
+            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)            
+            lambda = obj.Parameters.lambda;
+            phi = obj.Parameters.phi;
+            g = obj.Parameters.gravity;
+            a = obj.Parameters.radius;
+            Omega = obj.Parameters.angularSpeed;
+
+            f = 2*Omega*sin(phi.');
+
+            cosphi = cos(phi.');
+            tanphi = tan(phi.');
+
+            % create the interpolation matrix and derivatives
+            rbf = @otp.utils.rbf.quadratic;
+
+            interpolationradius = pi/3;
+            [W, dWdl, dWdp] = rbfinterp(lambda, phi, lambda, phi, interpolationradius, rbf);
+
+            % create the interploation matrix for plotting
+            Nplot = 50;
+            lambdaplot = linspace(-pi, pi, Nplot);
+            phiplot = linspace(-pi/2, pi/2, Nplot);
+            [lambdainterpgrid, phiinterpgrid] = meshgrid(lambdaplot, phiplot);
+            radiusplot = pi/6;
+            Wplot = rbfinterp(lambda, phi, lambdainterpgrid(:).', phiinterpgrid(:).', radiusplot, rbf);
+
+
+            plotlongitude = 360*(lambdainterpgrid/pi + 1)/2;
+            plotlatitude = 180*(phiinterpgrid/(pi/2))/2;
+
+            obj.PlottingInterp = Wplot;
+            obj.PlottingLatitude = plotlatitude;
+            obj.PlottingLongitude = plotlongitude;
+
+            % set the right hand side
+            obj.Rhs = otp.Rhs(@(t, huv) ...
+                otp.mfshallowwatersphere.f(huv, W, dWdl, dWdp, cosphi, tanphi, g, a, f));
+
+
+            %% Distance function
+            obj.DistanceFunction = @(t, huv, i, j) otp.mfshallowwatersphere.distfn(t, huv, i, j, lambda, phi);
+
+        end
+
+        function validateNewState(obj, newTimeSpan, newY0, newParameters)
+
+            validateNewState@otp.Problem(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 = internalSolve@otp.Problem(obj, 'Method', @ode45, varargin{:});
+        end
+
+    end
+end

src/+otp/+mfshallowwatersphere/distfn.m

@@ -0,0 +1,20 @@
+function d = distfn(~, ~, i, j, lambda, phi)
+
+n = numel(lambda);
+
+i = mod(i - 1, n) + 1;
+j = mod(j - 1, n) + 1;
+
+li = lambda(i);
+pi = phi(i);
+
+lj = lambda(j);
+pj = phi(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,32 @@
+function dhuvdt = f(huv, W, dWdl, dWdp, cosphi, tanphi, g, a, f)
+
+% unpack huv into its constituent parts
+n = size(huv, 1)/3;
+
+h = huv(1:n);
+u = huv((n+1):(2*n));
+v = huv((2*n+1):end);
+
+% calculate the derivatives with respect to phi and lambda
+dhdl = dWdl*h;
+dhdp = dWdp*h;
+dudl  = dWdl*u;
+dudp  = dWdp*u;
+dvdl = dWdl*v;
+dvdp = dWdp*v;
+
+% calculate 
+dcospvdp = dWdp*(cosphi.*v);
+
+dhdt = -(u./(a*cosphi)).*dhdl - (v/a).*dhdp - (h./(a*cosphi)).*(dudl + dcospvdp);
+dudt = -(u./(a*cosphi)).*dudl - (v/a).*dudp - (g./(a*cosphi)).*dhdl + (f + (u/a).*tanphi).*v;
+dvdt = -(u./(a*cosphi)).*dvdl - (v/a).*dvdp - (g./a).*dhdp - (f + (u/a).*tanphi).*u;
+
+% Interpolate the derivatives to smooth them out
+dhdt = W*dhdt;
+dudt = W*dudt;
+dvdt = W*dvdt;
+
+dhuvdt = [dhdt; dudt; dvdt];
+
+end