diff --git a/src/+otp/+mfshallowwatersphere/+presets/Canonical.m b/src/+otp/+mfshallowwatersphere/+presets/Canonical.m new file mode 100644 index 00000000..8efcb0fc --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/+presets/Canonical.m @@ -0,0 +1,91 @@ +classdef Canonical < otp.mfshallowwatersphere.MFShallowWaterSphereProblem + methods + function obj = Canonical(varargin) + + %load('nodes100.mat', 'x', 'y', 'z'); + + load('nodes250.mat', 'x', 'y', 'z'); + + %load('nodes784.mat', 'x', 'y', 'z'); + + %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; + %u0 = 1; + % Angular velocity of the earth + Omega = otp.utils.PhysicalConstants.EarthAngularVelocity; + + coriolisForce = 2*Omega*z; + + 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 = 2.0; + %params.rbf = @otp.utils.rbf.buhmann3; + + params.rbfradius = 2; + params.rbf = @otp.utils.rbf.buhmann3; + %params.rbf = @otp.utils.rbf.quadricspline; + %params.rbf = @otp.utils.rbf.wendlandWeC2PD1; + + % 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 = 4; + [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.1; + alpha = 0; + + % 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 = obj@otp.mfshallowwatersphere.MFShallowWaterSphereProblem(tspan, ... + huv0, params); + + end + end +end diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/getrossbyhaurwitzwave.m b/src/+otp/+mfshallowwatersphere/+presets/private/getrossbyhaurwitzwave.m new file mode 100644 index 00000000..e5d967ff --- /dev/null +++ b/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 diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/nodes100.mat b/src/+otp/+mfshallowwatersphere/+presets/private/nodes100.mat new file mode 100644 index 00000000..b213aaca Binary files /dev/null and b/src/+otp/+mfshallowwatersphere/+presets/private/nodes100.mat differ diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/nodes2000.mat b/src/+otp/+mfshallowwatersphere/+presets/private/nodes2000.mat new file mode 100644 index 00000000..cfa155a1 Binary files /dev/null and b/src/+otp/+mfshallowwatersphere/+presets/private/nodes2000.mat differ diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/nodes250.mat b/src/+otp/+mfshallowwatersphere/+presets/private/nodes250.mat new file mode 100644 index 00000000..3081a1e2 Binary files /dev/null and b/src/+otp/+mfshallowwatersphere/+presets/private/nodes250.mat differ diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/nodes500.mat b/src/+otp/+mfshallowwatersphere/+presets/private/nodes500.mat new file mode 100644 index 00000000..45cc5538 Binary files /dev/null and b/src/+otp/+mfshallowwatersphere/+presets/private/nodes500.mat differ diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/nodes784.mat b/src/+otp/+mfshallowwatersphere/+presets/private/nodes784.mat new file mode 100644 index 00000000..abf975a2 Binary files /dev/null and b/src/+otp/+mfshallowwatersphere/+presets/private/nodes784.mat differ diff --git a/src/+otp/+mfshallowwatersphere/+presets/private/velocitytocartesian.m b/src/+otp/+mfshallowwatersphere/+presets/private/velocitytocartesian.m new file mode 100644 index 00000000..96c2090d --- /dev/null +++ b/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 diff --git a/src/+otp/+mfshallowwatersphere/MFShallowWaterSphereProblem.m b/src/+otp/+mfshallowwatersphere/MFShallowWaterSphereProblem.m new file mode 100644 index 00000000..ff2e430d --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/MFShallowWaterSphereProblem.m @@ -0,0 +1,210 @@ +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 + Constraints + ConstraintsDerivative + end + + methods + + function plotSphere(obj, huv, projection) + if nargin < 3 + %projection = 'eqaazim'; + projection = 'eckert4'; + 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, d2Wdx2, d2Wdy2, d2Wdz2] = 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; + + Bx2 = d2Wdx2.*(1 - x.^2) + d2Wdy2.*(-x.*y) + d2Wdz2.*(-x.*z); + By2 = d2Wdx2.*(-x.*y) + d2Wdy2.*(1 - y.^2) + d2Wdz2.*(-y.*z); + Bz2 = d2Wdx2.*(-x.*z) + d2Wdy2.*(-y.*z) + d2Wdz2.*(1 - z.^2); + + Bx2 = Bx2/a; + By2 = By2/a; + Bz2 = Bz2/a; + + + Wscale = full(sum(W, 2)); + + try + Wdecomp = decomposition(W, 'chol'); + catch + error('The selected compination of RBF, radius, and nodes did not result in a SPD collocation matrix.') + end + + + + n = size(obj.Y0, 1)/4; + + h = obj.Y0(1:n, :); + + sch0 = sum(Wdecomp\h); + + % 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, Bx2, By2, Bz2, f, g, x, y, z, a, sch0)); + + + + H0Z0E0 = otp.mfshallowwatersphere.massenstrophyenergy(obj.Y0, S, Wdecomp, Bx, By, Bz, f, g, x, y, z, a, Wscale); + + % H0Z0E0(3) = 1; + %eq = [1; 0; 0; 1]; + + obj.Constraints = @(t, huvw) ... + otp.mfshallowwatersphere.massenstrophyenergy(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z, a, Wscale)./H0Z0E0 - 1; + + + obj.ConstraintsDerivative = @(t, huvw) ... + otp.mfshallowwatersphere.jacobianmassenstrophyenergy(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z, a, Wscale)./H0Z0E0; + + + %% 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) + + 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 diff --git a/src/+otp/+mfshallowwatersphere/distfn.m b/src/+otp/+mfshallowwatersphere/distfn.m new file mode 100644 index 00000000..908a73e0 --- /dev/null +++ b/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 diff --git a/src/+otp/+mfshallowwatersphere/f.m b/src/+otp/+mfshallowwatersphere/f.m new file mode 100644 index 00000000..32453287 --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/f.m @@ -0,0 +1,93 @@ +% 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, Bx2, By2, Bz2, f, g, x, y, z, a, sch0) + +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; + +k = 0; +nu = 1e-1; + +rhsDx = u.*(Bxcu) + v.*(By*cu) + w.*(Bz*cu) + f.*(y.*w - z.*v) + g*Bxch - k*u + nu*(Bx2*u); +rhsDy = u.*(Bx*cv) + v.*(Bycv) + w.*(Bz*cv) + f.*(z.*u - x.*w) + g*Bych - k*v + nu*(By2*v); +rhsDz = u.*(Bx*cw) + v.*(By*cw) + w.*(Bzcw) + f.*(x.*v - y.*u) + g*Bzch - k*w + nu*(Bz2*w); + +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) ); + + + +% MASS CORRECTION + +% g = sum(ch) - sch0; +% +% cones = Wdecomp\ones(n, 1); +% +% J = sum(cones.^2); +% +% gamma = 1e-1; +% %gamma = 0; +% +% dhdt = dhdt - gamma*cones*(J\g); + +% mass +%mean(h); + + +%theta = atan2(z, sqrt(x.^2 + y.^2)); +%lambda = atan2(y, x); + +%m = w./cos(theta); +%z = (v + m.*sin(lambda).*sin(theta))./(cos(lambda).*cos(theta)); + + +% Enstropy +% aa = By*cw - Bz*cv; +% bb = Bx*cw - Bz*cu; +% cc = Bx*cv - By*cu + f; +% 6370000*mean((aa + bb + cc)./h); +% +% % Energy +% aa = 0.5*(h.*(u.*u + v.*v + w.*w) + g*(h.*h)); +% mean(aa); + + +%mean( (m.^2 + z.^2 + g*h).*h ) + +%mean( (m.^2 + z.^2 + g*h).*h ) + + + +% Shepard interpolation for stability +% dhdt = S*dhdt; +% dudt = S*dudt; +% dvdt = S*dvdt; +% dwdt = S*dwdt; + +dhuvwdt = [dhdt; dudt; dvdt; dwdt]; + +end + diff --git a/src/+otp/+mfshallowwatersphere/jacobianmassenstrophyenergy.m b/src/+otp/+mfshallowwatersphere/jacobianmassenstrophyenergy.m new file mode 100644 index 00000000..5297b071 --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/jacobianmassenstrophyenergy.m @@ -0,0 +1,110 @@ + +function dhuvwdt = jacobianmassenstrophyenergy(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z, a, Wscale) + +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; + + +cones = Wdecomp\ones(n, 1); + +% mass +Hh = cones.'; +Hu = zeros(1, n); +Hv = zeros(1, n); +Hw = zeros(1, n); + + + +% vorticity +% % % Qh = zeros(1, n); +% % % Qu = -(Wdecomp\(Bz.'*cones + By.'*cones)).'; +% % % Qv = (Wdecomp\(Bx.'*cones - Bz.'*cones)).'; +% % % Qw = (Wdecomp\(By.'*cones + Bx.'*cones)).'; + + +% Enstropy + +Zinner = (By*(cw - cu) + Bx*(cw + cv) - Bz*(cv + cu) + f); + + +t0 = cones.*(Zinner./h); + +Zh = -(cones.*(((Zinner.^2)./(h.^2)))).'; + +Zu = -2*(Wdecomp\(Bz.'*t0 + By.'*t0)).'; +Zv = 2*(Wdecomp\(Bx.'*t0 - Bz.'*t0)).'; +Zw = 2*(Wdecomp\(By.'*t0 + Bx.'*t0)).'; + + +% % % Zix = By*cw - Bz*cv; +% % % Ziy = Bx*cw - Bz*cu; +% % % Ziz = Bx*cv - By*cu + f; +% % % ZixP = ( Zix.*(1 - x.^2) + Ziy.*(-x.*y) + Ziz.*(-x.*z) ); +% % % ZiyP = ( Zix.*(-x.*y) + Ziy.*(1 - y.^2) + Ziz.*(-y.*z) ); +% % % ZizP = ( Zix.*(-x.*z) + Ziy.*(-y.*z) + Ziz.*(1 - z.^2) ); +% % % Zinner = (ZixP + ZiyP + ZizP); +% % % +% % % +% % % coZihi = (cones.*Zinner)./h; +% % % +% % % Zh = -(cones.*(((Zinner.^2)./(h.^2)))).'; +% % % Zu = -2*(Wdecomp\( ... +% % % By.'*(coZihi.*(1 - z.^2 - x.*z - y.*z)) ... +% % % + Bz.'*(coZihi.*(1 - y.^2 - x.*y - y.*z)) ... +% % % )).'; +% % % Zv = 2*(Wdecomp\( ... +% % % Bx.'*(coZihi.*(1 - z.^2 - x.*z - y.*z)) ... +% % % - Bz.'*(coZihi.*(1 - x.^2 - x.*y - x.*z)) ... +% % % )).'; +% % % Zw = 2*(Wdecomp\( ... +% % % Bx.'*(coZihi.*(1 - y.^2 - x.*y - y.*z)) ... +% % % + By.'*(coZihi.*(1 - x.^2 - x.*y - x.*z)) ... +% % % )).'; + + +% Zh = +% Zu = -(Wdecomp\(Bz.'*(cones./h) + By.'*(cones./h))).'; +% Zv = (Wdecomp\(Bx.'*(cones./h) - Bz.'*(cones./h))).'; +% Zw = (Wdecomp\(By.'*(cones./h) + Bx.'*(cones./h))).'; + + +%Eh = (0.5*((u.*u + v.*v + w.*w) + 2*g*h)).'/n; + +Eh = (cones.*(0.5*((u.*u + v.*v + w.*w) + 2*g*h))).'; +Eu = (cones.*(h.*u)).'; +Ev = (cones.*(h.*v)).'; +Ew = (cones.*(h.*w)).'; + + + +% Eh = (g*ch + Wdecomp\(0.5*(u.*u + v.*v + w.*w) + g*h)).'; +% Eu = (ch.*u).'; +% Ev = (ch.*v).'; +% Ew = (ch.*w).'; + + + + +%dhuvwdt = [Hh, Hu, Hv, Hw; Qh, Qu, Qv, Qw; Zh, Zu, Zv, Zw; Eh, Eu, Ev,Ew]; + +%dhuvwdt = [Hh, Hu, Hv, Hw]; + +%dhuvwdt = [Zh, Zu, Zv, Zw]; + +% dhuvwdt = [Eh, Eu, Ev, Ew]; + +%dhuvwdt = [Hh, Hu, Hv, Hw; Eh, Eu, Ev, Ew]; + +dhuvwdt = [Hh, Hu, Hv, Hw; Zh, Zu, Zv, Zw; Eh, Eu, Ev, Ew]; + +end + diff --git a/src/+otp/+mfshallowwatersphere/massenstrophyenergy.m b/src/+otp/+mfshallowwatersphere/massenstrophyenergy.m new file mode 100644 index 00000000..30991eb0 --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/massenstrophyenergy.m @@ -0,0 +1,79 @@ + +function dhuvwdt = massenstrophyenergy(huvw, S, Wdecomp, Bx, By, Bz, f, g, x, y, z, a, Wscale) + +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; + +% mass +H = sum(ch); + + + +%Qinner = (By*(cw - cu) + Bx*(cw + cv) - Bz*(cv + cu) + f); + +%lQ = Wdecomp\Qinner; + +%Q = sum(lQ); + +% Enstropy +%aa = By*cw - Bz*cv; +%bb = Bx*cw - Bz*cu; +%cc = Bx*cv - By*cu + f; +% % % % Zix = By*cw - Bz*cv; +% % % % Ziy = Bx*cw - Bz*cu; +% % % % Ziz = Bx*cv - By*cu + f; +% % % % +% % % % +% % % % ZixP = ( Zix.*(1 - x.^2) + Ziy.*(-x.*y) + Ziz.*(-x.*z) ); +% % % % ZiyP = ( Zix.*(-x.*y) + Ziy.*(1 - y.^2) + Ziz.*(-y.*z) ); +% % % % ZizP = ( Zix.*(-x.*z) + Ziy.*(-y.*z) + Ziz.*(1 - z.^2) ); +% % % % +% % % % Zinner = (ZixP + ZiyP + ZizP).^2; + +Zinner = (By*(cw - cu) + Bx*(cw + cv) - Bz*(cv + cu) + f).^2; + +lZ = Wdecomp\(Zinner./h); + +Z = sum(lZ); + +% Energy +% u2 = u.*u; +% v2 = v.*v; +% z2 = z.*z; +% EixP = ( u2.*(1 - x.^2) + v2.*(-x.*y) + z2.*(-x.*z) ); +% EiyP = ( u2.*(-x.*y) + v2.*(1 - y.^2) + z2.*(-y.*z) ); +% EizP = ( u2.*(-x.*z) + v2.*(-y.*z) + z2.*(1 - z.^2) ); +% lE = Wdecomp\(0.5*(h.*(EixP + EiyP + EizP) + g*(h.*h))); + + +%lE = ch.*(0.5*(u.*u + v.*v + w.*w) + g*h); + + +lE = Wdecomp\(0.5*(h.*(u.*u + v.*v + w.*w) + g*(h.*h))); + + + +E = sum(lE); + + +%dhuvwdt = [H; Q; Z; E]; + +dhuvwdt = [H; Z; E]; + +% dhuvwdt = E; + + +%dhuvwdt = [H; Z]; + + +end + diff --git a/src/+otp/+mfshallowwatersphere/private/aapcmap.m b/src/+otp/+mfshallowwatersphere/private/aapcmap.m new file mode 100644 index 00000000..c8c24e70 --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/private/aapcmap.m @@ -0,0 +1,262 @@ +function C = aapcmap + +C = [0,0,0; +3,1,4; +6,2,8; +9,4,11; +12,5,15; +15,6,19; +18,7,23; +21,8,26; +24,10,30; +27,11,34; +30,12,38; +33,13,41; +36,15,45; +39,16,49; +42,17,53; +45,18,57; +48,19,60; +51,21,64; +54,22,68; +57,23,72; +60,24,75; +63,25,79; +66,27,83; +69,28,87; +72,29,90; +75,30,94; +78,32,98; +81,33,102; +84,34,106; +87,35,109; +90,36,113; +93,38,117; +96,39,121; +99,40,124; +102,41,128; +105,42,132; +108,44,136; +111,45,139; +114,46,143; +117,47,147; +120,48,151; +123,50,155; +126,51,158; +129,52,162; +132,53,166; +135,55,170; +138,56,173; +141,57,177; +144,58,181; +147,59,185; +150,61,188; +153,62,192; +156,63,196; +159,64,200; +162,65,204; +165,67,207; +168,68,211; +171,69,215; +174,70,219; +177,72,222; +180,73,226; +183,74,230; +186,75,234; +189,76,237; +190,78,240; +188,79,241; +186,81,241; +183,83,241; +181,85,241; +179,86,242; +176,88,242; +174,90,242; +172,91,242; +169,93,242; +167,95,243; +165,97,243; +162,98,243; +160,100,243; +158,102,244; +155,104,244; +153,105,244; +151,107,244; +148,109,245; +146,110,245; +144,112,245; +141,114,245; +139,116,245; +137,117,246; +134,119,246; +132,121,246; +130,122,246; +127,124,247; +125,126,247; +123,128,247; +120,129,247; +118,131,248; +116,133,248; +113,135,248; +111,136,248; +109,138,248; +107,140,249; +104,141,249; +102,143,249; +100,145,249; +97,147,250; +95,148,250; +93,150,250; +90,152,250; +88,153,251; +86,155,251; +83,157,251; +81,159,251; +79,160,251; +76,162,252; +74,164,252; +72,165,252; +69,167,252; +67,169,253; +65,171,253; +62,172,253; +60,174,253; +58,176,254; +55,178,254; +53,179,254; +51,181,254; +48,183,254; +46,184,255; +44,186,255; +44,187,253; +46,188,249; +48,189,246; +51,190,242; +53,191,238; +55,192,234; +58,193,231; +60,194,227; +62,195,223; +65,196,219; +67,197,216; +69,198,212; +72,199,208; +74,200,204; +76,201,201; +79,202,197; +81,203,193; +83,204,189; +86,205,186; +88,206,182; +90,207,178; +93,208,174; +95,209,171; +97,210,167; +100,210,163; +102,211,159; +104,212,156; +107,213,152; +109,214,148; +112,215,145; +114,216,141; +116,217,137; +119,218,133; +121,219,130; +123,220,126; +126,221,122; +128,222,118; +130,223,115; +133,224,111; +135,225,107; +137,226,103; +140,227,100; +142,228,96; +144,229,92; +147,230,88; +149,231,85; +151,232,81; +154,232,77; +156,233,73; +158,234,70; +161,235,66; +163,236,62; +165,237,58; +168,238,55; +170,239,51; +172,240,47; +175,241,43; +177,242,40; +179,243,36; +182,244,32; +184,245,28; +186,246,25; +189,247,21; +191,248,17; +192,248,19; +193,248,23; +194,248,27; +195,248,30; +196,249,34; +197,249,38; +198,249,42; +199,249,45; +200,249,49; +201,249,53; +202,249,56; +203,249,60; +204,249,64; +205,250,68; +206,250,71; +207,250,75; +208,250,79; +209,250,83; +210,250,86; +211,250,90; +212,250,94; +213,250,98; +214,251,101; +215,251,105; +216,251,109; +217,251,113; +218,251,116; +219,251,120; +220,251,124; +221,251,128; +222,251,131; +223,252,135; +224,252,139; +225,252,143; +226,252,146; +227,252,150; +228,252,154; +229,252,158; +230,252,161; +231,252,165; +232,252,169; +233,253,173; +234,253,176; +235,253,180; +236,253,184; +237,253,188; +238,253,191; +239,253,195; +240,253,199; +241,253,203; +242,254,206; +243,254,210; +244,254,214; +245,254,218; +246,254,221; +247,254,225; +248,254,229; +249,254,233; +250,254,236; +251,255,240; +252,255,244; +253,255,248; +254,255,251; +255,255,255]; + +C = C/255; + +end diff --git a/src/+otp/+mfshallowwatersphere/private/rbfinterp.m b/src/+otp/+mfshallowwatersphere/private/rbfinterp.m new file mode 100644 index 00000000..af649e9e --- /dev/null +++ b/src/+otp/+mfshallowwatersphere/private/rbfinterp.m @@ -0,0 +1,55 @@ +function [W, dWdx, dWdy, dWdz, d2Wdx2, d2Wdy2, d2Wdz2] = rbfinterp(x1, y1, z1, x2, y2, z2, r, rbf) + +dx = x2 - x1.'; +dy = y2 - y1.'; +dz = z2 - z1.'; + +d = sqrt(dx.^2 + dy.^2 + dz.^2); + +n = size(d, 1); +m = size(d, 2); + +k = d/r; + +% we only support compactly supported RBFs +I = find(k <= 1); + +[Is, Js] = ind2sub(size(k), I); + +k = k(I); +d = d(I); +dx = dx(I); +dy = dy(I); +dz = dz(I); + +[Wv, dWdkv, d2Wdk2v] = rbf(k); +Wv(isnan(Wv)) = 0; +dWdkv(isnan(dWdkv)) = 0; +d2Wdk2v(isnan(d2Wdk2v)) = 0; +d2Wdk2v(isinf(d2Wdk2v)) = 0; + +dkdx = dx./(r*d); +dkdy = dy./(r*d); +dkdz = dz./(r*d); +dkdx(isnan(dkdx)) = 0; +dkdy(isnan(dkdy)) = 0; +dkdz(isnan(dkdz)) = 0; + +d2kdx2 = (dy.^2 + dz.^2)./(r*(d.^3)); +d2kdy2 = (dx.^2 + dz.^2)./(r*(d.^3)); +d2kdz2 = (dx.^2 + dy.^2)./(r*(d.^3)); +d2kdx2(isnan(d2kdx2)) = 0; +d2kdy2(isnan(d2kdy2)) = 0; +d2kdz2(isnan(d2kdz2)) = 0; + + +W = sparse(Is, Js, Wv, n, m); +dWdx = sparse(Is, Js, dWdkv.*dkdx, n, m); +dWdy = sparse(Is, Js, dWdkv.*dkdy, n, m); +dWdz = sparse(Is, Js, dWdkv.*dkdz, n, m); + +d2Wdx2 = sparse(Is, Js, d2Wdk2v.*(dkdx.^2) + dWdkv.*d2kdx2, n, m); +d2Wdy2 = sparse(Is, Js, d2Wdk2v.*(dkdy.^2) + dWdkv.*d2kdy2, n, m); +d2Wdz2 = sparse(Is, Js, d2Wdk2v.*(dkdz.^2) + dWdkv.*d2kdz2, n, m); + +end diff --git a/src/+otp/+utils/+rbf/buhmann3.m b/src/+otp/+utils/+rbf/buhmann3.m new file mode 100644 index 00000000..666cecf7 --- /dev/null +++ b/src/+otp/+utils/+rbf/buhmann3.m @@ -0,0 +1,16 @@ +function [f, df, d2f] = buhmann3(k) + +maskk0 = k == 0; + +f = (1/3) + k.^2 - (4/3)*(k.^3) + 2*(k.^2).*log(k); +f(maskk0) = (1/3) + k(maskk0).^2 - (4/3)*(k(maskk0).^3); + + +df = 4*k - 4*(k.^2) + 4*k.*log(k); +df(maskk0) = 4*k(maskk0) - 4*(k(maskk0).^2); + +d2f = 8 - 8*k + 4*log(k); +d2f(maskk0) = -Inf; + + +end diff --git a/src/+otp/+utils/+rbf/buhmann6.m b/src/+otp/+utils/+rbf/buhmann6.m new file mode 100644 index 00000000..6876322a --- /dev/null +++ b/src/+otp/+utils/+rbf/buhmann6.m @@ -0,0 +1,12 @@ +function [f, df] = buhmann6(k) + +maskk0 = k == 0; + +f = (1/15) + (19/6)*(k.^2) - (16/3)*(k.^3) + 3*(k.^4) - (16/15)*(k.^5) + (1/6)*(k.^6) + 2*(k.^2).*log(k); +f(maskk0) = (1/15) + (19/6)*(k(maskk0).^2) - (16/3)*(k(maskk0).^3) + 3*(k(maskk0).^4) - (16/15)*(k(maskk0).^5) + (1/6)*(k(maskk0).^6); + + +df = (25/3)*(k) - (16)*(k.^2) + 12*(k.^3) - (16/3)*(k.^4) + k.^5 + 4*k.*log(k); +df(maskk0) = (25/3)*(k(maskk0)) - (16)*(k(maskk0).^2) + 12*(k(maskk0).^3) - (16/3)*(k(maskk0).^4) + k(maskk0).^5; + +end diff --git a/src/+otp/+utils/+rbf/doublecos.m b/src/+otp/+utils/+rbf/doublecos.m new file mode 100644 index 00000000..5dd840d2 --- /dev/null +++ b/src/+otp/+utils/+rbf/doublecos.m @@ -0,0 +1,7 @@ +function [f, df] = doublecos(k) + +f = 4*cospi(k) + cospi(2*k) + 3; + +df = -2*pi*(2*sinpi(k) + sinpi(2*k)); + +end diff --git a/src/+otp/+utils/+rbf/euclidhat2.m b/src/+otp/+utils/+rbf/euclidhat2.m new file mode 100644 index 00000000..464e1d67 --- /dev/null +++ b/src/+otp/+utils/+rbf/euclidhat2.m @@ -0,0 +1,7 @@ +function [f, df] = euclidhat2(k) + +f = (2/pi)*(acos(k) - k.*sqrt(1 - k.^2)); + +df = -(4/pi)*sqrt(1 - k.^2); + +end diff --git a/src/+otp/+utils/+rbf/expdecay.m b/src/+otp/+utils/+rbf/expdecay.m new file mode 100644 index 00000000..c8c60d44 --- /dev/null +++ b/src/+otp/+utils/+rbf/expdecay.m @@ -0,0 +1,7 @@ +function [f, df] = expdecay(k) + +f = exp(-k) - exp(-1); + +df = -exp(-k); + +end diff --git a/src/+otp/+utils/+rbf/firstinversefn.m b/src/+otp/+utils/+rbf/firstinversefn.m new file mode 100644 index 00000000..f7643e57 --- /dev/null +++ b/src/+otp/+utils/+rbf/firstinversefn.m @@ -0,0 +1,6 @@ +function [f, df] = firstinversefn(k) + +f = 2./(1 + k) - 1; +df = -2./((1+k).^2); + +end diff --git a/src/+otp/+utils/+rbf/fourthinversefn.m b/src/+otp/+utils/+rbf/fourthinversefn.m new file mode 100644 index 00000000..f96f4e3d --- /dev/null +++ b/src/+otp/+utils/+rbf/fourthinversefn.m @@ -0,0 +1,6 @@ +function [f, df] = fourthinversefn(k) + +f = -6./(1 + k.^4) + 8./(1 + k.^3) - 1; +df = 24*(k.^2).*(-1./((1 + k.^3).^2) + k./((1 + k.^4).^2)); + +end diff --git a/src/+otp/+utils/+rbf/gc.m b/src/+otp/+utils/+rbf/gc.m new file mode 100644 index 00000000..ee8ad7c9 --- /dev/null +++ b/src/+otp/+utils/+rbf/gc.m @@ -0,0 +1,23 @@ +function [f, df] = gc(k) + +f = zeros(size(k)); +df = zeros(size(k)); + +theta = 1/2; + +mask1 = k <= theta; +mask2 = and((k > theta), (k <= 2*theta)); + +km1 = k(mask1)/theta; +km2 = k(mask2)/theta; + +f(mask1) = ones(size(km1)) - (5/3)*(km1.^2) + (5/8) * (km1.^3) + (1/2)*(km1.^4) - (1/4)*(km1.^5); +f(mask2) = 4*ones(size(km2)) - 5*km2 + (5/3)*(km2.^2) + (5/8)*(km2.^3) - (1/2)*(km2.^4) + (1/12)*(km2.^5) - (2./(3*km2)); + +df(mask1) = - (10/3)*(km1/theta) + (15/8) * (km1.^2)/theta + (2)*(km1.^3)/theta - (5/4)*(km1.^4)/theta; +df(mask2) = - 5/theta + (10/3)*(km2)/theta + (15/8)*(km2.^2)/theta - (2)*(km2.^3)/theta + (5/12)*(km2.^4)/theta + (2./(3*(km2.^2)))/theta; + +f(isinf(f)) = 0; +df(isinf(df)) = 0; + +end diff --git a/src/+otp/+utils/+rbf/hanning.m b/src/+otp/+utils/+rbf/hanning.m new file mode 100644 index 00000000..c5cf1489 --- /dev/null +++ b/src/+otp/+utils/+rbf/hanning.m @@ -0,0 +1,6 @@ +function [f, df] = hanning(k) + +f = cospi(k) + 1; +df = -sinpi(k); + +end diff --git a/src/+otp/+utils/+rbf/quadricspline.m b/src/+otp/+utils/+rbf/quadricspline.m new file mode 100644 index 00000000..29e8946d --- /dev/null +++ b/src/+otp/+utils/+rbf/quadricspline.m @@ -0,0 +1,11 @@ +function [f, df, d2f] = quadricspline(k) + +k2 = k.*k; +k3 = k.*k2; +k4 = k.*k3; + +f = (5/2)*(1 - 6*k2 + 8*k3 - 3*k4); +df = (5/2)*(-12*k + 24*k2 - 12*k3); +d2f = (5/2)*(-12 + 48*k - 36*k2); + +end diff --git a/src/+otp/+utils/+rbf/wendlandWeC2PD1.m b/src/+otp/+utils/+rbf/wendlandWeC2PD1.m new file mode 100644 index 00000000..fb5ad575 --- /dev/null +++ b/src/+otp/+utils/+rbf/wendlandWeC2PD1.m @@ -0,0 +1,7 @@ +function [f, df, d2f] = wendlandWeC2PD1(k) + +f = ((1 - k).^3).*(3*k + 1); +df = -12*((1 - k).^2).*k; +d2f = 12*(k - 1).*(3*k - 1); + +end diff --git a/src/+otp/+utils/PhysicalConstants.m b/src/+otp/+utils/PhysicalConstants.m index 2fbeec9b..e003e6cf 100644 --- a/src/+otp/+utils/PhysicalConstants.m +++ b/src/+otp/+utils/PhysicalConstants.m @@ -21,6 +21,11 @@ DaysPerYear = 365.25; % d EarthGravity = 9.80665; % m / s^2 + + SecondsInEarthDay = 24*60*60; % s + + EarthAngularVelocity = 7.292e-5; % rad / s + EarthRadius = 6.370e6; % m end end diff --git a/src/HalfExplicit_RK_NonLin2FAST.m b/src/HalfExplicit_RK_NonLin2FAST.m new file mode 100644 index 00000000..ad8ba106 --- /dev/null +++ b/src/HalfExplicit_RK_NonLin2FAST.m @@ -0,0 +1,100 @@ +classdef HalfExplicit_RK_NonLin2FAST + %HALFEXPLICIT_RK_LIN Summary of this class goes here + % Detailed explanation goes here + + properties + a + ah + b + c + ch + s + nmax + iter_tol + end + + methods + function obj = HalfExplicit_RK_NonLin2FAST(a, ah, b, c, ch, nm, tol) + obj.a = a; + obj.ah = ah; + obj.b = b; + obj.c = c; + obj.ch = ch; + obj.s = length(a); + obj.nmax = nm; + obj.iter_tol = tol; + end + + function [y_f, z_f] = time_step(obj, fx, f_z, g, dt, t, y_0, z_0) + k = zeros(length(y_0), obj.s); + + %k(:, 1) = f(t,y_0,z_0); + f_zy_0 = f_z(t,y_0); + k(:, 1) = fx(t,y_0) - f_zy_0*z_0; + %iter_jac = (dt * obj.ah(2, 2)) * jac(t, y_0, z_0); + + iter_jac = -(dt * obj.ah(2, 2)) * (f_zy_0.'*f_zy_0); + + %i = 2; + + %while true + for i = 2:obj.s + + inter_y = y_0; + inter_y_star = y_0; + for j = 1:i - 1 + inter_y = inter_y + dt * obj.a(i, j) * k(:, j); + inter_y_star = inter_y_star + dt * obj.ah(i, j)*k(:,j); + end + + %iter_func = @(x) g(t + dt * obj.ch(i), inter_y_star + dt* obj.ah(i,i) * f(t + dt * obj.c(i), inter_y ,x) ); + + fxinter = fx(t + dt * obj.c(i), inter_y); + f_zinter = f_z(t + dt * obj.c(i), inter_y); + iter_func = @(x) g(t + dt * obj.ch(i), inter_y_star + dt* obj.ah(i,i) * (fxinter - f_zinter*x ) ); + + inter_z = obj.chord_iteration(iter_func, iter_jac, z_0); + + k(:, i) = fxinter - f_zinter*inter_z; + + %k(:,i) = f(t + dt * obj.c(i), inter_y, inter_z); + + %if i >= obj.s + % break + %end + + %iter_jac = (dt * obj.ah(i + 1, i + 1)) * jac(t + dt*obj.c(i), inter_y, inter_z); + + %i = i + 1; + end + + y_f = y_0 + dt*k * obj.b'; + z_f = inter_z; + + %u_f = [;inter_z]; + end + + %Assume jac is precomputed + function y_f = chord_iteration(obj, func, jac, y_0) + i = 0; + iJ = pinv(jac); + %dJ = decomposition(-jac, 'chol'); + while i < obj.nmax + x_i = iJ * -func(y_0); + %x_i = dJ\func(y_0); + %x_i = jac\-func(y_0); + y_f = x_i + y_0; + + i = i + 1; + + if norm(x_i) < obj.iter_tol + return; + end + y_0 = y_f; + end + + + end + end +end + diff --git a/src/test.m b/src/test.m new file mode 100644 index 00000000..b87c0505 --- /dev/null +++ b/src/test.m @@ -0,0 +1,121 @@ +close all; +m = otp.mfshallowwatersphere.presets.Canonical; + +y0 = m.Y0; + +n = numel(y0); + + +% f = @(y) m.Constraints(0, y); +% f0 = m.Constraints(0, m.Y0); +% G = zeros(3, n); +% Ga = m.ConstraintsDerivative(0, m.Y0); +% h = 1e-4; +% for i = 1:n +% e = zeros(n, 1); +% e(i) = 1; +% G(:, i) = (f(y0 + h*e) - f0)/h; +% end +% abs(G - Ga)./abs(G) +% norm(G - Ga)/norm(G) +% imagesc(G - Ga) +% return + +dt = 250; + +f = m.Rhs.F; +g = m.Constraints; +dg = @(t, y) m.ConstraintsDerivative(t, y).'; + + + a_coe = [[0, 0]; + [1,0]]; + +ah_coe = [[0,0]; + [1/2,1/2]]; + +c_coe = [0; 1]; +ch_coe = [0; 1]; +b_coe = [1/2,1/2]; + +tss = HalfExplicit_RK_NonLin2FAST(a_coe, ah_coe, b_coe, c_coe, ch_coe, 100, 1e-8); + +steps = 1e6; + +fig1 = figure; + +figure; +tiledlayout(3, 1); +nexttile +fig2 = plot(nan(steps, 1), 'LineWidth', 1.2); +title('Mass') +%nexttile +%fig3 = plot(nan(steps, 1), 'LineWidth', 1.2); +%title('Vorticity') +nexttile +fig3 = plot(nan(steps, 1), 'LineWidth', 1.2); +title('Enstrophy') +nexttile +fig4 = plot(nan(steps, 1), 'LineWidth', 1.2); +title('Energy') + + +nz = numel(g(0, m.Y0)); + +cons = zeros(nz, 1); + + +Z = m.Y0; +z0 = zeros(nz, 1); +for i = 1:steps + + %[Z, z0] = tss.time_step(func, phi, Jac, dt, i*dt, Z, z0); + + [Z, z0] = tss.time_step(f, dg, g, dt, i*dt, Z, z0); + + + %sol = ode45(f, [0, dt], Z, odeset('AbsTol', 1e-1)); + %Z = sol.y(:, end); + + % [~, Z] = RK65SSP33(func, [0 dt], Z, atol, rtol); + phii = g(0, Z); + + cons(:, i) = phii; + + if rem(i, 500) == 0 + + g(0, Z) + + figure(fig1); + m.plotSphere(Z); + drawnow; + +% [hn, un, vn] = rs(Z); +% +% pv = hn; +% +% fig1.CData = [pv; pv(1, :)].'; +% fig1.FaceColor = 'interp'; +% tt1.String = "Time = "+ i*dt/3600/24 +" days"; +% drawnow; +% + fig2.YData(1:i) = cons(1, 1:i); + fig3.YData(1:i) = cons(2, 1:i); + %fig4.YData(1:i) = cons(3, 1:i); + %fig5.YData(1:i) = cons(4, 1:i); + disp("Time = "+i*dt/3600/24+" days"); + drawnow; + + end + + + +end + + + +%abs(G - Ga)./abs(G) + +%imagesc(Ga - G) + +