Better 2D Curvilinear Operators (MATLAB)

examples/matlab_octave/curvConvergenceTest.m

@@ -0,0 +1,218 @@
+%% Convergence Test for Cuvilinear Operators
+% 
+% In 2D space, MOLE operators can only calculate the xi and eta derivatives
+% On a curvilinear domain, the x and xi, and y and eta directions are
+% not always the same. We can assume a transformation of x = x(xi, eta)
+% and y = y(xi, eta) with Jacobian J = x_xi y_eta - x_eta y_xi.
+% 
+% If we have a function u(x,y), via chain rule we get:
+%   u_xi  = u_x x_xi  + u_y y_xi
+%   u_eta = u_x x_eta + u_y y_eta
+% 
+% Assuming a non-zero Jacobian, we can invert the mapping to obtain:
+%   u_x = (u_xi y_eta - u_eta y_xi) / J
+%   u_y = (u_eta x_xi - u_xi x_eta) / J
+% 
+% Everything on the right is easily calculable using MOLE's operators
+% 
+
+%% Comparison between current implementation and new implementation
+% 
+% Current Implementation:
+%   - Nodal Jacobian derivatives
+%   - MATLAB interpolation (2nd order)
+%   - Hard to understand implementation (DI2 and GI2 functions: It is 
+%       obvious what they are doing, but not how they are doing it.)
+%   - Cannot handle periodic domains
+% 
+% New Implementation
+%   - MOLE gradient used to calculate Jacobian derivatives
+%   - MOLE interpolators
+%   - Easy to understand implementation
+%   - Can handle periodic domains
+% 
+
+%% Problem
+% 
+% ∆u = 0
+% 
+% 1 < R < 2
+% 0 < θ < pi
+% 
+% x = R cos(θ)
+% y = R sin(θ)
+% 
+% Dirichlet Boundary Conditions
+% 
+% Exact Solution: (many options)
+%   u(x,y) = sin(x) sinh(y)
+%   u(x,y) = x^2 - y^2
+%   u(x,y) = ln(x^2 + y^2)
+%   u(x,y) = x^3 - 3 x y^2
+% 
+close all; clear; clc;
+
+
+%% Fixed mesh, varying k Comparison
+% Parameters
+k = [2 4 6 8];
+m = 100;
+n =  25;
+numCenters = (m+2) * (n+2);
+dx = pi / m;
+dy = 1  / n;
+
+% ue = @(X,Y) sin(X) .* sinh(Y);
+% ue = @(X,Y) X.^2 - Y.^2;
+% ue = @(X,Y) log(X.^2 + Y.^2);
+ue = @(X,Y) X.^3 - 3 * X .* Y.^2;
+
+errorCur = zeros(size(k));
+errorNew = zeros(size(k));
+
+dc = [1;1;1;1]; nc = [0;0;0;0];
+
+for i = 1:numel(k)
+
+    fprintf("Starting k = %d\n", k(i))
+    % Nodes to make current operators
+    rsN = 1:dy:2; tsN = pi:-dx:0;
+    [TSN,RSN] = meshgrid(tsN,rsN);
+    xn = RSN .* cos(TSN);
+    yn = RSN .* sin(TSN);
+
+    % Centers to make new operators
+    rsC = [1 1+dy/2:dy:2-dy/2 2];
+    tsC = [pi pi-dx/2:-dx:dx/2 0];
+    [TSC,RSC] = meshgrid(tsC,rsC);
+    X = RSC .* cos(TSC);
+    Y = RSC .* sin(TSC);
+
+    % Build operators
+    curG = grad2DCurv(k(i),xn,yn);
+    curD =  div2DCurv(k(i),xn,yn);
+    curL = curD * curG;
+
+    newG = grad2DCurv(k(i),X,Y,dc,nc);
+    newD =  div2DCurv(k(i),X,Y,dc,nc);
+    newL = newD * newG;
+
+    % Boundary Conditions
+    u = ue(X,Y);
+
+    l = u(:,1);  r = u(:,end);
+    b = u(1,:)'; t = u(end,:)';
+    v = {l(2:end-1);r(2:end-1);b;t};
+    B = zeros(numCenters,1);
+    [curL,B] = addScalarBC2D(curL,B,k(i),m,dx,n,dy,dc,nc,v);
+    [newL,B] = addScalarBC2D(newL,B,k(i),m,dx,n,dy,dc,nc,v);
+
+    curU = curL \ B;
+    newU = newL \ B;
+
+    curU = reshape(curU,m+2,n+2)';
+    newU = reshape(newU,m+2,n+2)';
+
+    % Maximum Absolute Error
+    errorCur(i) = max(max(abs(u-curU)));
+    errorNew(i) = max(max(abs(u-newU)));
+
+end
+
+% Plot results
+figure
+hold on
+plot(k, errorCur, 'LineWidth', 2)
+plot(k, errorNew, 'LineWidth', 2)
+hold off
+yscale log
+xlim([k(1), k(end)])
+xlabel('k')
+ylabel('Maximum Absolute Error')
+legend('Current Implementation', 'New Implementation','Location','southwest')
+grid on
+title("Maximum Absolute Error vs k")
+subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, m = 100, n = 25")
+
+
+%% Fixed k, varying mesh Comparison
+k = 4;
+minCells = 2*k+1;
+m = [4 8 16 32 64] * minCells;
+n = [1 2  4  8 16] * minCells;
+numCenters = (m+2) .* (n+2);
+dx = pi ./ m;
+dy = 1  ./ n;
+
+errorCur = zeros(size(numCenters));
+errorNew = zeros(size(numCenters));
+
+dc = [1;1;1;1]; nc = [0;0;0;0];
+
+for i = 1:numel(numCenters)
+
+    fprintf("Starting m = %d, n = %d\n", m(i),n(i))
+    % Nodes for current operators
+    rsN = 1:dy(i):2; tsN = pi:-dx(i):0;
+    [TSN,RSN] = meshgrid(tsN,rsN);
+    xn = RSN .* cos(TSN);
+    yn = RSN .* sin(TSN);
+
+    % Centers for new operators
+    rsC = [1 1+dy(i)/2:dy(i):2-dy(i)/2 2];
+    tsC = [pi pi-dx(i)/2:-dx(i):dx(i)/2 0];
+    [TSC,RSC] = meshgrid(tsC,rsC);
+    X = RSC .* cos(TSC);
+    Y = RSC .* sin(TSC);
+
+    % Build operators
+    curG = grad2DCurv(k,xn,yn);
+    curD =  div2DCurv(k,xn,yn);
+    curL = curD * curG;
+
+    newG = grad2DCurv(k,X,Y,dc,nc);
+    newD =  div2DCurv(k,X,Y,dc,nc);
+    newL = newD * newG;
+
+    % Boundary Conditions
+    u = ue(X,Y);
+
+    l = u(:,1);  r = u(:,end);
+    b = u(1,:)'; t = u(end,:)';
+    v = {l(2:end-1);r(2:end-1);b;t};
+    B = zeros(numCenters(i),1);
+    [curL,B] = addScalarBC2D(curL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v);
+    [newL,B] = addScalarBC2D(newL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v);
+
+    curU = curL \ B;
+    newU = newL \ B;
+
+    curU = reshape(curU,m(i)+2,n(i)+2)';
+    newU = reshape(newU,m(i)+2,n(i)+2)';
+
+    % Maximum Absolute Error
+    errorCur(i) = max(max(abs(u-curU)));
+    errorNew(i) = max(max(abs(u-newU)));
+
+end
+
+% Plot results
+h = numCenters.^-0.5;
+refCur = errorCur(1) * (h / h(1)).^k;
+refNew = errorNew(1) * (h / h(1)).^k;
+figure
+hold on
+plot(h, errorCur, 'LineWidth', 2)
+plot(h, errorNew, 'LineWidth', 2)
+plot(h, refCur, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--')
+plot(h, refNew, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--')
+hold off
+xscale log
+yscale log
+xlim([min(h) max(h)])
+xlabel('1 / sqrt(Number of Centers)')
+ylabel('Maximum Absolute Error')
+legend('Current Implementation', 'New Implementation',"k =  " + k + " Reference Line",'Location','southeast')
+grid on
+title("Maximum Absolute Error vs 1 / sqrt(Number of Centers)")
+subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, k = " + k)

examples/matlab_octave/elliptic2DCurvPeriodic.m

@@ -0,0 +1,88 @@
+clear; clc; close all;
+addpath("../../src/matlab_octave/")
+% 
+% ∆u = 0
+% 
+% 1 < r < 2
+% 0 < θ < 2 pi
+% 
+% x = r cos(θ)
+% y = r sin(θ)
+% 
+% Dirichlet Boundary Conditions
+% 
+% Exact solution
+% u(x, y) = x^2 - y^2
+% 
+
+% Parameters
+k = 4;
+m = 100;
+n = 25;
+dx = 2 * pi / m;
+dy = 1 / n;
+dc = [0;0;1;1];
+nc = [0;0;0;0];
+
+% Exact solution
+ue = @(X, Y) X.^2 - Y.^2;
+
+% Create Mesh
+rs = [1 (1 + dy / 2) : dy : (2 - dy / 2) 2];
+ts = 0 : dx : (2*pi - dx);
+[T,R] = meshgrid(ts, rs);
+X = R .* cos(T);
+Y = R .* sin(T);
+
+% Build Operators
+G = grad2DCurv(k, X, Y, dc, nc);
+D = div2DCurv(k, X, Y, dc, nc);
+L = D * G;
+
+% Boundary Conditions
+u = ue(X, Y);
+
+l = 0; % Left and right boundaries don't exist
+r = 0;
+b = u(1, :)';
+t = u(end, :)';
+v = {l; r; b; t};
+B = zeros(numel(X), 1);
+[L0, B0] = addScalarBC2D(L, B, k, m, dx, n, dy, dc, nc, v);
+
+ua = L0 \ B0;
+ua = reshape(ua, m, n + 2)';
+
+% Plot results
+% Join left and right edges for nicer plotting
+X = [X X(:,1)];
+Y = [Y Y(:,1)];
+u = [u u(:,1)];
+ua = [ua ua(:,1)];
+
+figure
+surf(X, Y, ua, "EdgeColor", "none")
+title("Approximate Solution")
+xlabel("X")
+ylabel("Y")
+axis equal
+view([0 90])
+cb = colorbar;
+cb.Label.String = "u(X,Y)";
+colormap("jet")
+
+figure
+surf(X, Y, u, "EdgeColor", "none")
+title("Exact Solution")
+xlabel("X")
+ylabel("Y")
+axis equal
+view([0 90])
+cb = colorbar;
+cb.Label.String = "u(X,Y)";
+colormap("jet")
+
+max_err = max(max(abs(ua - u)));
+disp("Maximum Absolute Error: " + max_err)
+rel_err = 100 * max_err / (max(max(u)) - min(min(u)));
+disp("Maximum Relative Error: " + rel_err)