Skip to content

Commit 1d7a942

Browse files
author
Michael J. Schmidt
committed
initial commit
0 parents  commit 1d7a942

25 files changed

+2948
-0
lines changed

‎1D_2domain/DTerrorAnalysis_1D_2omega.m‎

Lines changed: 233 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,233 @@
1+
%==========================================================================
2+
% This script conducts a convergence analysis in terms of time step length,
3+
% dt, for the discontinuous D(x) MTPT method employing a semi-analytical solution.
4+
% This script, will generate Figure 10 in:
5+
% "A Mass-transfer Particle-tracking Method for Simulating Transport
6+
% with Discontinuous Diffusion Coefficients,"
7+
% Advances in Water Resources, 2020.
8+
%==========================================================================
9+
10+
% close all
11+
clear variables
12+
13+
% position of point-mass source
14+
x0 = -6.5;
15+
% length of 1D domain
16+
L = 5e1;
17+
% number of particles for mass-transfer (MT) simulation
18+
N = 5e3;
19+
% total simulation time
20+
maxT = 6;
21+
22+
% number of iterations to normalize using the Sinkhorn-Knopp Algorithm
23+
nNorm = 1e3;
24+
25+
% number of members in ensemble (i.e., the number of D2 values to consider)
26+
numDEns = 3;
27+
28+
% subdomain boundary location
29+
gamma = -1;
30+
31+
% number of refinements in dt
32+
numDT = 6;
33+
NvecStart = 1;
34+
35+
% half refinements in dt
36+
Nvec = linspace(NvecStart, numDT + NvecStart - 1, numDT);
37+
Nvec = 2.^Nvec;
38+
dtVec = 1 ./ Nvec;
39+
40+
% print 1 / dt and dt vectors
41+
inv_dt = Nvec
42+
dtVec
43+
44+
% diffusion coefficients (D2vec is length numDEns)
45+
D1 = 5;
46+
D2vec = [2.5 0.5 0.05];
47+
48+
% locations of stationary MT particles
49+
X = linspace(x0 - L / 2, x0 + L / 2, N)';
50+
51+
% calculate analytical solution, using Eqs. (12)-(16)
52+
analyticalSoln = zeros(N, 3);
53+
for i = 1 : 3
54+
analyticalSoln(:, i) = MT_CnJ_arb_dCut(X, x0, D1, D2vec(i), gamma, maxT);
55+
end
56+
57+
% calculated height to make plots pretty
58+
height = 1.1 * max(max(analyticalSoln));
59+
60+
% dirac IC
61+
[~, sourceX] = min(abs(X - x0));
62+
X(sourceX) = x0;
63+
mass = zeros(N, 1);
64+
mass(sourceX) = 1 / (L / N);
65+
mass0 = mass;
66+
67+
%% MT simulation
68+
69+
MTsoldtVec = zeros(N, 3);
70+
71+
% dimension of array is: # of error metrics x dt ensemble x D2 ensemble
72+
errVec = zeros(4, numDT, 3);
73+
74+
% dt ensemble loop
75+
for dtEns = 1 : numDT
76+
77+
% time step length and calculated number of time steps to take
78+
dt = dtVec(dtEns);
79+
nSteps = ceil(maxT / dt);
80+
81+
% D2 ensemble loop
82+
for Dens = 1 : numDEns
83+
84+
D2 = D2vec(Dens);
85+
86+
% calculate the max interaction distance for rangesearch()
87+
dist = 3 * sqrt(4 * max([D1 D2]) * dt);
88+
89+
% conduct the rangesearch to find nearby particles
90+
[idx, r] = rangesearch(X, X, dist, 'BucketSize', ceil(1e-2 * N));
91+
92+
% determine how many particles are nearby and preallocate the vectors
93+
% to build the sparse weight matrix
94+
Nclose = sum(cellfun('length', idx));
95+
row = zeros(Nclose, 1);
96+
col = zeros(Nclose, 1);
97+
val = zeros(Nclose, 1);
98+
99+
% calculate the entries of the weight matrix
100+
start = 1;
101+
for i = 1 : N
102+
finish = start - 1 + length(idx{i});
103+
104+
% this builds the weight matrix with the predictor-corrector solution
105+
% according to line 3 in Algorithm 2
106+
row(start : finish) = i;
107+
col(start : finish) = idx{i};
108+
val(start : finish) = PrCo_1D_2omega(X(idx{i}), X(i), D1, D2, gamma, dt);
109+
start = finish + 1;
110+
end
111+
112+
% when N is super large the sparse methods and clearing large arrays saves memory
113+
clear idx r
114+
115+
% create the sparse weight matrix
116+
Wmat = sparse(row, col, val);
117+
clear row col val
118+
119+
% normalize via SK, with nNorm iterations, ending with a column
120+
% normalization
121+
Wmat = sinkhornKnoppCol(Wmat, 'MaxIter', nNorm);
122+
123+
% build the transfer matrix according to Algorithm 3
124+
Tmat = Wmat;
125+
126+
% initialize the mass vector
127+
mass = mass0;
128+
129+
% conduct the mass transfers
130+
for i = 1 : nSteps
131+
mass = Tmat * mass;
132+
end
133+
134+
% save the final mass vector for plotting
135+
MTsoldtVec(:, Dens) = mass;
136+
137+
% calculate the error in RMSE, then inf, 2, and 1 norms
138+
errVec(1, dtEns, Dens) = sqrt(mean((MTsoldtVec(:, Dens) - analyticalSoln(:, Dens)).^2));
139+
errVec(2, dtEns, Dens) = norm(MTsoldtVec(:, Dens) - analyticalSoln(:, Dens), inf);
140+
errVec(3, dtEns, Dens) = norm(MTsoldtVec(:, Dens) - analyticalSoln(:, Dens), 2);
141+
errVec(4, dtEns, Dens) = norm(MTsoldtVec(:, Dens) - analyticalSoln(:, Dens), 1);
142+
143+
end
144+
145+
end
146+
147+
%% print out final mass and percent mass change
148+
149+
sumM0 = sum(mass0);
150+
151+
fprintf('MT Initial Mass = %f\n', sumM0)
152+
153+
for i = 1 : numDEns
154+
ensMsum = sum(MTsoldtVec(:, i));
155+
fprintf('MT Final Mass (Ens = %i) = %f\n', i, ensMsum)
156+
fprintf(' Percent Mass Change = %.2f\n', 100 * abs(sumM0 - ensMsum) / sumM0)
157+
end
158+
159+
160+
%% Final time solution plots
161+
162+
% get color order for plotting
163+
color = lines(7);
164+
165+
% spatial mass plot for final time step
166+
fig = 2;
167+
figure(fig)
168+
clf
169+
fill([X(1) gamma gamma X(1) X(1)],height.*[0 0 1 1 0],[0 0.25 1],'FaceAlpha',0.05); grid on; hold on;
170+
fill([X(end) gamma gamma X(end) X(end)],height.*[0 0 1 1 0],[1 0.25 0],'FaceAlpha',0.05);
171+
172+
p1 = plot(X, analyticalSoln(:, 1), 'k', 'linewidth', 17);
173+
p2 = plot(X(:), MTsoldtVec(:, 1), '-.+', 'color', color(1, :), 'linewidth', 2,...
174+
'MarkerSize', 10);
175+
176+
plot(X, analyticalSoln(:, 2), 'k', 'linewidth', 17)
177+
plot(X(:), MTsoldtVec(:, 2), '-.+', 'color', color(2, :), 'linewidth', 2,...
178+
'MarkerSize', 10);
179+
180+
plot(X, analyticalSoln(:, 3), 'k', 'linewidth', 17)
181+
plot(X(:), MTsoldtVec(:, 3), '-.+', 'color', color(3, :), 'linewidth', 2,...
182+
'MarkerSize', 10);
183+
184+
p4 = line([x0 x0], [0 height], 'color', 'b', 'linewidth', 5);
185+
p5 = line([gamma gamma], [0 height], 'color', 'r', 'linewidth', 5);
186+
xlim([X(1) X(end)]); ylim([0 height]);
187+
[~, hobj, ~, ~] = legend([p1 p2 p4 p5], {'\textbf{Analytic Solution}',...
188+
'\textbf{Mass-Transfer}', '\textbf{Source Location} {\boldmath$(x_0)$}',...
189+
'\textbf{Subdomain Boundary} {\boldmath($\gamma$)}'}, 'Interpreter',...
190+
'latex', 'FontSize', 32, 'Location', 'northwest');
191+
set(hobj, 'linewidth', 4);
192+
rectangle('Position', [X(end) - 6.8, 0.48 * height, 5.9, 0.0116], 'FaceColor', 'w',...
193+
'EdgeColor', 'k', 'LineWidth', 1)
194+
for i = 1 : 3
195+
text(X(end) - 1, 0.6 * height - 0.004 * (i - 1), ['{\boldmath $D_{2}=\ $}\bf' num2str(D2vec(i), '%1.2f')],...
196+
'FontSize', 32, 'color', color(i, :), 'HorizontalAlignment',...
197+
'right', 'VerticalAlignment', 'middle', 'Interpreter', 'latex');
198+
end
199+
rectangle('Position', [X(1) + 0.9, 0.58 * height, 5.9, 0.0036], 'FaceColor', 'w',...
200+
'EdgeColor', 'k', 'LineWidth', 1)
201+
text(X(1) + 1, 0.6 * height, ['{\boldmath $D_{1}=\ $}\bf' num2str(D1, '%1.2f')], 'FontSize', 32, 'color', 'k',...
202+
'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', 'Interpreter', 'latex');
203+
xlabel('\textbf{\emph{x}-position [L]}', 'Interpreter', 'latex', 'FontSize', 32)
204+
ylabel('\textbf{Concentration [mol L{\boldmath$^{-1}$}]}', 'Interpreter', 'latex', 'FontSize', 32)
205+
figure(fig)
206+
set(gcf, 'Position', [0, 100, 1800, 900])
207+
208+
%% Convergence plots
209+
210+
fig = 27;
211+
figure(fig)
212+
clf
213+
for i = 1 : Dens
214+
subplot(1, 3, i)
215+
loglog(dtVec, errVec(2, :, i), '-o', 'LineWidth', 5.5)
216+
hold on
217+
loglog(dtVec, errVec(3, :, i), '-o', 'LineWidth', 5.5)
218+
p2 = polyfit(log10(dtVec), log10(errVec(2, :, i)), 1);
219+
plot(dtVec, 10.^(p2(1) * log10(dtVec) + p2(2)), '--d', 'LineWidth', 3.5)
220+
p3 = polyfit(log10(dtVec), log10(errVec(3, :, i)), 1);
221+
plot(dtVec, 10.^(p3(1) * log10(dtVec) + p3(2)), '--d', 'LineWidth', 3.5)
222+
xlabel('\boldmath{$\Delta t$}','Interpreter','latex','FontSize', 22)
223+
ylabel('\textbf{Error}','Interpreter','latex','FontSize', 22)
224+
legend({'\boldmath{$\ell^\infty$} \textbf{Norm}', '\boldmath{$\ell^2$} \textbf{Norm}', ['\boldmath{$\mathcal{O}(\Delta t^p), p = $} \bf',...
225+
num2str(p2(1))], ['\boldmath{$\mathcal{O}(\Delta t^p), p = $} \bf',...
226+
num2str(p3(1))]}, 'Interpreter', 'latex','FontSize', 22, 'Location', 'northwest')
227+
title(['\boldmath{$D_2 = $} \bf', num2str(D2vec(i))], 'Interpreter', 'latex','FontSize', 22)
228+
ax = gca;
229+
ax.YLim(2) = ax.YLim(2) * 10;
230+
end
231+
ylim([1e-3 1e0])
232+
figure(fig)
233+
set(gcf, 'Position', [0, 100, 1400, 450])

‎1D_2domain/MT_CnJ_arb_dCut.m‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
% analytical solution for composite diffusion, adapted from Carslaw and
2+
% Jaeger, "Conduction of Heat in Solids", p. 363
3+
4+
% Inputs/Assumptions:
5+
% 1. X, the 1D position vector, must be in ascending order and contain zero
6+
% 2. x0 is the x-location of the Dirac pulse source
7+
% 3. D1 is diffusion coefficient for x < gamma
8+
% 4. D2 is diffusion coefficient for x > gamma
9+
% 5. maxT is the final simulation time of interest
10+
% 6. the initial condition is a pulse of mass 1, centered at x = x0
11+
% 7. gamma is the location of the interface between the subdomains with
12+
% distinct diffusion coefficients
13+
14+
function solution = MT_CnJ_arb_dCut(X, x0, D1, D2, gamma, maxT)
15+
16+
% Eq. (12)
17+
if x0 >= gamma
18+
19+
% solution for x < 0
20+
c1 = @(x) (D2 * D1 * D2^(-1 / 2)) /...
21+
((D2 * sqrt(D1) + D1 * sqrt(D2)) * sqrt(pi * D1 * maxT)) *...
22+
exp(-(x - gamma - (x0 - gamma) * sqrt(D1 / D2)).^2 / (4 * D1 * maxT));
23+
24+
% solution for x > 0
25+
c2 = @(x) 1 / (2 * sqrt(pi * D2 * maxT)) *...
26+
exp(-(x - x0).^2 / (4 * D2 * maxT)) + ...
27+
(D2 * sqrt(D1) - D1 * sqrt(D2)) /...
28+
(2 * (D2 * sqrt(D1) + D1 * sqrt(D2)) * sqrt(pi * D2 * maxT)) *...
29+
exp(-(x + x0 - 2 * gamma).^2 / (4 * D2 * maxT));
30+
31+
negIndex = find(X <= gamma);
32+
posIndex = find(X > gamma);
33+
34+
solution = zeros(size(X));
35+
36+
solution(negIndex) = c1(X(negIndex));
37+
solution(posIndex) = c2(X(posIndex));
38+
39+
% Eq. (16)
40+
else
41+
42+
temp = D1;
43+
D1 = D2;
44+
D2 = temp;
45+
46+
% solution for x > 0
47+
c1 = @(x) (D2 * D1 * D2^(-1 / 2)) /...
48+
((D2 * sqrt(D1) + D1 * sqrt(D2)) * sqrt(pi * D1 * maxT)) *...
49+
exp(-(x - gamma - (x0 - gamma) * sqrt(D1 / D2)).^2 / (4 * D1 * maxT));
50+
51+
% solution for x < 0
52+
c2 = @(x) 1 / (2 * sqrt(pi * D2 * maxT)) *...
53+
exp(-(x - x0).^2 / (4 * D2 * maxT)) + ...
54+
(D2 * sqrt(D1) - D1 * sqrt(D2)) /...
55+
(2 * (D2 * sqrt(D1) + D1 * sqrt(D2)) * sqrt(pi * D2 * maxT)) *...
56+
exp(-(x + x0 - 2 * gamma).^2 / (4 * D2 * maxT));
57+
58+
negIndex = find(X <= gamma);
59+
posIndex = find(X > gamma);
60+
61+
solution = zeros(size(X));
62+
63+
solution(negIndex) = c2(X(negIndex));
64+
solution(posIndex) = c1(X(posIndex));
65+
66+
end
67+

0 commit comments

Comments
 (0)