Skip to content

Commit 5bac07f

Browse files
jbrzenskjbrzensk
authored
Merge branch 'main' into task/remove_elliptic1DaddscalarBC.m
2 parents 804c7f8 + aa229f9 commit 5bac07f

21 files changed

Lines changed: 1708 additions & 581 deletions

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
[![MATLAB File Exchange](https://www.mathworks.com/matlabcentral/images/matlab-file-exchange.svg)](https://www.mathworks.com/matlabcentral/fileexchange/124870-mole)
55
[![License](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
66
[![Build Status](https://github.com/csrc-sdsu/mole/actions/workflows/ci.yml/badge.svg)](https://github.com/csrc-sdsu/mole/actions/workflows/ci.yml)
7+
[![codecov](https://codecov.io/gh/csrc-sdsu/mole/branch/main/graph/badge.svg)](https://codecov.io/gh/csrc-sdsu/mole)
78
[![Documentation](https://readthedocs.org/projects/mole-docs/badge/?version=main)](https://mole-docs.readthedocs.io/en/main/)
89

910
## Description

‎examples/matlab_octave/curvConvergenceTest.m‎

Lines changed: 218 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
%% Convergence Test for Cuvilinear Operators
2+
%
3+
% In 2D space, MOLE operators can only calculate the xi and eta derivatives
4+
% On a curvilinear domain, the x and xi, and y and eta directions are
5+
% not always the same. We can assume a transformation of x = x(xi, eta)
6+
% and y = y(xi, eta) with Jacobian J = x_xi y_eta - x_eta y_xi.
7+
%
8+
% If we have a function u(x,y), via chain rule we get:
9+
% u_xi = u_x x_xi + u_y y_xi
10+
% u_eta = u_x x_eta + u_y y_eta
11+
%
12+
% Assuming a non-zero Jacobian, we can invert the mapping to obtain:
13+
% u_x = (u_xi y_eta - u_eta y_xi) / J
14+
% u_y = (u_eta x_xi - u_xi x_eta) / J
15+
%
16+
% Everything on the right is easily calculable using MOLE's operators
17+
%
18+
19+
%% Comparison between current implementation and new implementation
20+
%
21+
% Current Implementation:
22+
% - Nodal Jacobian derivatives
23+
% - MATLAB interpolation (2nd order)
24+
% - Hard to understand implementation (DI2 and GI2 functions: It is
25+
% obvious what they are doing, but not how they are doing it.)
26+
% - Cannot handle periodic domains
27+
%
28+
% New Implementation
29+
% - MOLE gradient used to calculate Jacobian derivatives
30+
% - MOLE interpolators
31+
% - Easy to understand implementation
32+
% - Can handle periodic domains
33+
%
34+
35+
%% Problem
36+
%
37+
% ∆u = 0
38+
%
39+
% 1 < R < 2
40+
% 0 < θ < pi
41+
%
42+
% x = R cos(θ)
43+
% y = R sin(θ)
44+
%
45+
% Dirichlet Boundary Conditions
46+
%
47+
% Exact Solution: (many options)
48+
% u(x,y) = sin(x) sinh(y)
49+
% u(x,y) = x^2 - y^2
50+
% u(x,y) = ln(x^2 + y^2)
51+
% u(x,y) = x^3 - 3 x y^2
52+
%
53+
close all; clear; clc;
54+
55+
56+
%% Fixed mesh, varying k Comparison
57+
% Parameters
58+
k = [2 4 6 8];
59+
m = 100;
60+
n = 25;
61+
numCenters = (m+2) * (n+2);
62+
dx = pi / m;
63+
dy = 1 / n;
64+
65+
% ue = @(X,Y) sin(X) .* sinh(Y);
66+
% ue = @(X,Y) X.^2 - Y.^2;
67+
% ue = @(X,Y) log(X.^2 + Y.^2);
68+
ue = @(X,Y) X.^3 - 3 * X .* Y.^2;
69+
70+
errorCur = zeros(size(k));
71+
errorNew = zeros(size(k));
72+
73+
dc = [1;1;1;1]; nc = [0;0;0;0];
74+
75+
for i = 1:numel(k)
76+
77+
fprintf("Starting k = %d\n", k(i))
78+
% Nodes to make current operators
79+
rsN = 1:dy:2; tsN = pi:-dx:0;
80+
[TSN,RSN] = meshgrid(tsN,rsN);
81+
xn = RSN .* cos(TSN);
82+
yn = RSN .* sin(TSN);
83+
84+
% Centers to make new operators
85+
rsC = [1 1+dy/2:dy:2-dy/2 2];
86+
tsC = [pi pi-dx/2:-dx:dx/2 0];
87+
[TSC,RSC] = meshgrid(tsC,rsC);
88+
X = RSC .* cos(TSC);
89+
Y = RSC .* sin(TSC);
90+
91+
% Build operators
92+
curG = grad2DCurv(k(i),xn,yn);
93+
curD = div2DCurv(k(i),xn,yn);
94+
curL = curD * curG;
95+
96+
newG = grad2DCurv(k(i),X,Y,dc,nc);
97+
newD = div2DCurv(k(i),X,Y,dc,nc);
98+
newL = newD * newG;
99+
100+
% Boundary Conditions
101+
u = ue(X,Y);
102+
103+
l = u(:,1); r = u(:,end);
104+
b = u(1,:)'; t = u(end,:)';
105+
v = {l(2:end-1);r(2:end-1);b;t};
106+
B = zeros(numCenters,1);
107+
[curL,B] = addScalarBC2D(curL,B,k(i),m,dx,n,dy,dc,nc,v);
108+
[newL,B] = addScalarBC2D(newL,B,k(i),m,dx,n,dy,dc,nc,v);
109+
110+
curU = curL \ B;
111+
newU = newL \ B;
112+
113+
curU = reshape(curU,m+2,n+2)';
114+
newU = reshape(newU,m+2,n+2)';
115+
116+
% Maximum Absolute Error
117+
errorCur(i) = max(max(abs(u-curU)));
118+
errorNew(i) = max(max(abs(u-newU)));
119+
120+
end
121+
122+
% Plot results
123+
figure
124+
hold on
125+
plot(k, errorCur, 'LineWidth', 2)
126+
plot(k, errorNew, 'LineWidth', 2)
127+
hold off
128+
yscale log
129+
xlim([k(1), k(end)])
130+
xlabel('k')
131+
ylabel('Maximum Absolute Error')
132+
legend('Current Implementation', 'New Implementation','Location','southwest')
133+
grid on
134+
title("Maximum Absolute Error vs k")
135+
subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, m = 100, n = 25")
136+
137+
138+
%% Fixed k, varying mesh Comparison
139+
k = 4;
140+
minCells = 2*k+1;
141+
m = [4 8 16 32 64] * minCells;
142+
n = [1 2 4 8 16] * minCells;
143+
numCenters = (m+2) .* (n+2);
144+
dx = pi ./ m;
145+
dy = 1 ./ n;
146+
147+
errorCur = zeros(size(numCenters));
148+
errorNew = zeros(size(numCenters));
149+
150+
dc = [1;1;1;1]; nc = [0;0;0;0];
151+
152+
for i = 1:numel(numCenters)
153+
154+
fprintf("Starting m = %d, n = %d\n", m(i),n(i))
155+
% Nodes for current operators
156+
rsN = 1:dy(i):2; tsN = pi:-dx(i):0;
157+
[TSN,RSN] = meshgrid(tsN,rsN);
158+
xn = RSN .* cos(TSN);
159+
yn = RSN .* sin(TSN);
160+
161+
% Centers for new operators
162+
rsC = [1 1+dy(i)/2:dy(i):2-dy(i)/2 2];
163+
tsC = [pi pi-dx(i)/2:-dx(i):dx(i)/2 0];
164+
[TSC,RSC] = meshgrid(tsC,rsC);
165+
X = RSC .* cos(TSC);
166+
Y = RSC .* sin(TSC);
167+
168+
% Build operators
169+
curG = grad2DCurv(k,xn,yn);
170+
curD = div2DCurv(k,xn,yn);
171+
curL = curD * curG;
172+
173+
newG = grad2DCurv(k,X,Y,dc,nc);
174+
newD = div2DCurv(k,X,Y,dc,nc);
175+
newL = newD * newG;
176+
177+
% Boundary Conditions
178+
u = ue(X,Y);
179+
180+
l = u(:,1); r = u(:,end);
181+
b = u(1,:)'; t = u(end,:)';
182+
v = {l(2:end-1);r(2:end-1);b;t};
183+
B = zeros(numCenters(i),1);
184+
[curL,B] = addScalarBC2D(curL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v);
185+
[newL,B] = addScalarBC2D(newL,B,k,m(i),dx(i),n(i),dy(i),dc,nc,v);
186+
187+
curU = curL \ B;
188+
newU = newL \ B;
189+
190+
curU = reshape(curU,m(i)+2,n(i)+2)';
191+
newU = reshape(newU,m(i)+2,n(i)+2)';
192+
193+
% Maximum Absolute Error
194+
errorCur(i) = max(max(abs(u-curU)));
195+
errorNew(i) = max(max(abs(u-newU)));
196+
197+
end
198+
199+
% Plot results
200+
h = numCenters.^-0.5;
201+
refCur = errorCur(1) * (h / h(1)).^k;
202+
refNew = errorNew(1) * (h / h(1)).^k;
203+
figure
204+
hold on
205+
plot(h, errorCur, 'LineWidth', 2)
206+
plot(h, errorNew, 'LineWidth', 2)
207+
plot(h, refCur, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--')
208+
plot(h, refNew, 'LineWidth', 1.5, 'Color', 'k', 'LineStyle', '--')
209+
hold off
210+
xscale log
211+
yscale log
212+
xlim([min(h) max(h)])
213+
xlabel('1 / sqrt(Number of Centers)')
214+
ylabel('Maximum Absolute Error')
215+
legend('Current Implementation', 'New Implementation',"k = " + k + " Reference Line",'Location','southeast')
216+
grid on
217+
title("Maximum Absolute Error vs 1 / sqrt(Number of Centers)")
218+
subtitle("∆u = 0, 1 < R < 2, 0 < theta < π, k = " + k)

‎examples/matlab_octave/elliptic2DCurvPeriodic.m‎

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
clear; clc; close all;
2+
addpath("../../src/matlab_octave/")
3+
%
4+
% ∆u = 0
5+
%
6+
% 1 < r < 2
7+
% 0 < θ < 2 pi
8+
%
9+
% x = r cos(θ)
10+
% y = r sin(θ)
11+
%
12+
% Dirichlet Boundary Conditions
13+
%
14+
% Exact solution
15+
% u(x, y) = x^2 - y^2
16+
%
17+
18+
% Parameters
19+
k = 4;
20+
m = 100;
21+
n = 25;
22+
dx = 2 * pi / m;
23+
dy = 1 / n;
24+
dc = [0;0;1;1];
25+
nc = [0;0;0;0];
26+
27+
% Exact solution
28+
ue = @(X, Y) X.^2 - Y.^2;
29+
30+
% Create Mesh
31+
rs = [1 (1 + dy / 2) : dy : (2 - dy / 2) 2];
32+
ts = 0 : dx : (2*pi - dx);
33+
[T,R] = meshgrid(ts, rs);
34+
X = R .* cos(T);
35+
Y = R .* sin(T);
36+
37+
% Build Operators
38+
G = grad2DCurv(k, X, Y, dc, nc);
39+
D = div2DCurv(k, X, Y, dc, nc);
40+
L = D * G;
41+
42+
% Boundary Conditions
43+
u = ue(X, Y);
44+
45+
l = 0; % Left and right boundaries don't exist
46+
r = 0;
47+
b = u(1, :)';
48+
t = u(end, :)';
49+
v = {l; r; b; t};
50+
B = zeros(numel(X), 1);
51+
[L0, B0] = addScalarBC2D(L, B, k, m, dx, n, dy, dc, nc, v);
52+
53+
ua = L0 \ B0;
54+
ua = reshape(ua, m, n + 2)';
55+
56+
% Plot results
57+
% Join left and right edges for nicer plotting
58+
X = [X X(:,1)];
59+
Y = [Y Y(:,1)];
60+
u = [u u(:,1)];
61+
ua = [ua ua(:,1)];
62+
63+
figure
64+
surf(X, Y, ua, "EdgeColor", "none")
65+
title("Approximate Solution")
66+
xlabel("X")
67+
ylabel("Y")
68+
axis equal
69+
view([0 90])
70+
cb = colorbar;
71+
cb.Label.String = "u(X,Y)";
72+
colormap("jet")
73+
74+
figure
75+
surf(X, Y, u, "EdgeColor", "none")
76+
title("Exact Solution")
77+
xlabel("X")
78+
ylabel("Y")
79+
axis equal
80+
view([0 90])
81+
cb = colorbar;
82+
cb.Label.String = "u(X,Y)";
83+
colormap("jet")
84+
85+
max_err = max(max(abs(ua - u)));
86+
disp("Maximum Absolute Error: " + max_err)
87+
rel_err = 100 * max_err / (max(max(u)) - min(min(u)));
88+
disp("Maximum Relative Error: " + rel_err)

0 commit comments

Comments
 (0)