Commit b04167d9278d6ed12fea00f78e54264423b9a45e

GNU Octave compatibility; progress indicator; visualize only nt_vis timesteps (averaging the actual solution); support both surf and contourf for plotting; configurable colormap; clarified comments
Harjoitustyo/corrosive_darcy_flow.m
(217 / 47)
  
101101
102102 s.mu = 1; % fluidin viskositeetti
103103 s.beta = 1; % fluidin kompressoituvuuskerroin
104 s.p_infty = 20; % ympäristön paine
104 s.p_infty = 20; % ympäristön paine
105105 s.v_0 = 1; % vapaa virtausnopeus (ennen huokoiseen aineeseen törmäystä)
106106
107107 s.cA0 = 1; % syövyttävän aineen konsentraatio ennen törmäystä huokoiseen aineeseen
108108
109109 s.tmax = 10.0; % simulaation loppuaika
110 s.dt = 1e-3; % aika-askeleen koko samoissa yksiköissä kuin tmax
110 s.dt = 1e-2; % aika-askeleen koko samoissa yksiköissä kuin tmax
111111 % (vakio vain yksinkertaisuuden vuoksi; RK4 ei rajoita)
112 s.nx = 51; % paikkadiskretoinnin solmujen lukumäärä
112 s.nx = 11; % paikkadiskretoinnin solmujen lukumäärä
113113
114 num_eps = 1e-8; % Numeerinen pyöristys, tätä pienemmät luvut = 0
114 nt_vis = 101; % visualisointipisteiden määrä aika-akselilla (kannattaa pitää kohtuullisena)
115115
116 num_eps = 1e-4; % Numeerinen pyöristys, tätä pienemmät luvut = 0
117
118 cm = 'jet'; % Värikartta plottia varten
119
120 % Piirtotapa. 'contourf' tai 'surf'.
121 %
122 % Tasa-arvokäyrien piirto voi olla hidasta (etenkin Octavessa, jos nt_vis on suuri),
123 % joten toinen mahdollisuus on käyttää ylhäältä päin katsottua, värikoodattua 3D-pintaa.
124 % Huomaa, että tällöin kunkin mesh-elementin väri määräytyy vasemman yläkulman (eikä keskipisteen)
125 % funktioarvon mukaisesti.
126 %
127 plot_method = 'contourf';
128% plot_method = 'surf';
129
130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 % GNU Octave -yhteensopivuus
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133
134 % MATLABin ja Octaven APIt (sekä tuettu toiminnallisuus) ovat hiukan erilaisia.
135 % Ratkaisu tähän on hiukan yhteensopivuusliimaa, jotta tämä ohjelma toimii molemmissa.
136 %
137 % Tämän osion kommentit ovat englanniksi, koska koodi on pääosin peräisin
138 % aikaisemmista tutkimusskripteistäni.
139
140 % Struct for MATLAB/Octave compatibility stuff.
141 compat = [];
142
143 % MATLAB/Octave?
144 % http://stackoverflow.com/questions/2246579/how-do-i-detect-if-im-running-matlab-or-octave
145 compat.isOctave = exist('OCTAVE_VERSION') ~= 0;
146
147 % The timing functions work differently in MATLAB/Octave.
148 %
149 % Note that "t = tic;" works the same, but Octave's toc() does not support
150 % input arguments like MATLAB's does.
151 %
152 if compat.isOctave
153 % Create a toc() function that takes an input argument and behaves like MATLAB's.
154 compat.toc = @(t) (double (tic ()) - double (t))*1e-6; % see "help tic" in Octave
155
156 % Octave's colorbar function takes the bare location as argument.
157 compat.colorbar = @(loc) colorbar(loc);
158 else
159 % Use the built-in toc() function.
160 compat.toc = @toc;
161
162 % MATLAB's colorbar function takes the keyword 'Location' and then the location as argument.
163 compat.colorbar = @(loc) colorbar('Location', loc);
164 end
165
116166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117167 % Laskentaskripti (särkyvää)
118168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119169
170 fprintf(1, 'Alustetaan ratkaisija...\n');
171
172 if strcmp(plot_method, 'contourf')
173 f_plot = @(X,Y,Z) contourf(X, Y, Z);
174 f_postplot = []; % tyhjä = ei postplot-funktiota
175 elseif strcmp(plot_method, 'surf')
176 f_plot = @(X,Y,Z) surf(X, Y, Z, 'EdgeColor', 'none');
177 f_postplot = @() view(0,90);
178 else
179 error('Tuntematon piirtotapa "%s". Tuettuja ovat "contourf" ja "surf".', plot_method);
180 end
181
120182 s.xx = transpose(linspace(0, s.L, s.nx)); % pystyvektori, helpompaa tämän koodin kannalta...
121183 s.h = s.xx(2) - s.xx(1);
122184
203203 % Aika-askelten määrä (mukaanlukien alkuhetki t = 0)
204204 s.nt = 1 + floor(s.tmax/s.dt);
205205
206 % Visualisoitavan aikavälin pituus
207 dt_vis = s.tmax / (nt_vis - 1);
208 % Kuinka monta aika-askelta otetaan per yksi visualisoitava askel
209 steps_per_vis = floor(dt_vis / s.dt);
210
206211 % Varataan muisti tuloksille. NaN on kätevä arvo alustamattomalle osalle,
207212 % koska MATLAB ei piirrä sitä kuviin.
213 %
214 % Huomaa, että vain visualisoitava määrä aika-askelia tallennetaan.
208215 %
209216 % HUOM: aika pystyakselilla (rivi), suurevektorin y arvo vaaka-akselilla (sarake).
210217 %
211 ys = nan(s.nt,s.ny);
218 ys = nan(nt_vis,s.ny);
212219
213220 % Ladotaan alkutila paikalleen:
214221 %
225225 y(s.I_cB) = s.cB0;
226226 ys(1,:) = y'; % tulokset varastoidaan aika-askeleittain vaakavektoreina
227227
228 fprintf(1, 'Ratkaistaan...\n');
228 % Laskennan edistymismittaria varten (progress indicator)
229 num_computed_ts = s.nt - 1; % montako aika-askelta oikeasti lasketaan (huom. lähtötilanne jo tunnetaan!)
230 perc = 0; % valmistunut prosenttiosuus (percentage complete)
231 perc_prev = -1; % valmistunut prosenttiosuus edellisellä aika-askeleella
232 % (käytetään muutoksen havaitsemiseen)
233 printouts_per_line = 5; % montako viestiä laitetaan per rivi
234 printouts_on_this_line = 0; % montako viestiä tällä hetkellä nykyisellä rivillä
229235
236 fprintf(1, 'Alustus valmis. Vapausasteita per muuttuja %d, laskettavia aika-askelia %d.\n', s.nx, num_computed_ts);
237
238 t_all_start = tic;
239 t_start = tic;
240 fprintf(1, 'Ratkaistaan. Laskennan tilanne:\n');
241
230242 % Integroidaan aika-askeleittain.
231243 t = 0;
244 vis_step = 2;
245 y_vis = zeros(s.ny,1);
232246 for k = 2:s.nt
233 % Reunaehdot paikalleen.
247 % Edistymismittari (progress indicator). Päivittyy prosentin välein.
248 % Huomaa k-1, koska k alkaa 2:sta ja alkutila tunnetaan.
249 %
250 perc = floor( 100 * (k-1)/num_computed_ts );
251 if perc > perc_prev
252 fprintf('%d%% (%d / %d)', perc, k-1, num_computed_ts);
253 printouts_on_this_line = printouts_on_this_line + 1;
254 if printouts_on_this_line < printouts_per_line
255 fprintf(' ');
256 else
257 fprintf('\n');
258 printouts_on_this_line = 0;
259 end
260 end
261
262 % Asetetaan reunaehdot paikalleen.
234263 %
235264 % Vasemman reunan Dirichlet-ehdot ovat helppoja.
236265 % Pakotetaan arvot ennen integrointia:
268268
269269 % Oikealla reunalla joutuu hiukan temppuilemaan, koska tarvitaan p:n arvo,
270270 % ja käytettävissä on vain p_x(L,t). Toisaalta, lisäksi käytettävissä
271 % on arvo edellisestä aika-askeleesta (tai 1. askeleella alkuarvo),
272 % ja yhtälöstä (p_t) saadaan p, kun tunnetaan p_x ja p_xx. Jälkimmäinen
273 % seuraa ensimmäisestä derivoimalla.
271 % on p:n arvo edellisestä aika-askeleesta (tai 1. askeleella alkuarvo),
272 % ja yhtälöstä (p_t) saadaan p_t, kun tunnetaan p_x ja p_xx. Jälkimmäinen
273 % seuraa ensimmäisestä derivoimalla. Aikaderivaatasta puolestaan voidaan
274 % integroida p.
274275 %
275276 % Sijoitetaan siis näin saatava p_t(L,t) paikalleen funktioon f
276277 % samalla tavalla kuin sisäosassakin, ja annetaan rk4:n hoitaa loput.
281281 y = rk4(y, t, s.dt, @(y,t) f(s,y,t));
282282
283283 % Pyöristys, jottei numeerinen epätarkkuus kaada laskua.
284 % (Muutoin voi tulla hyvin pieniä negatiivisia arvoja konsentraatioille
285 % siellä, missä ne ovat lähellä nollaa.)
284 % Konsentraatioille voi numeerisesti tulla hyvin pieniä negatiivisia
285 % arvoja siellä, missä ne ovat lähellä nollaa.
286286 %
287287 y( abs(y) < num_eps ) = 0;
288288
290290 % Kaikki suureet p, cA, cB ovat 1) numeroita, 2) epänegatiivisia.
291291 if sum(isnan(y)) || sum(y < 0)
292292 % DEBUG
293 fprintf(1, '\nRatkaisuvirhe (sanity check fail). Muuttujien arvot:\n');
293294 p = y(s.I_p)
294295 cA = y(s.I_cA)
295 cB = y(s.I_cA)
296 cB = y(s.I_cB)
296297 k
297298 t
298 error('kaboom');
299 break; % visualisoidaan kuitenkin siihen asti kuin saatiin järkevää dataa.
299300 end
300301
301 % Tallennetaan tulos piirtoa varten.
302 % Tulosten tallennus.
303 % Vain nt_vis kappaletta aika-askelia tallennetaan.
302304 %
303 % TODO: Tähän voisi virittää mekanismin, joka tallentaa vain joka n:nnen askeleen
304 % (tai yhteensä vain m kpl askelia, laskien automaattisesti n:n sopivasti).
305 %
306 ys(k,:) = y';
305 y_vis = y_vis + y; % summataan välitulos keskiarvoa varten
306 if mod(k - 1, steps_per_vis) == 0
307 % Keskiarvoistetaan tämän visualisointiaskeleen kohdalla lasketut aika-askeleet.
308 %
309 % Kannattaa huomata, että yleensä tämä ehkä menettelee, muttei ole kovin hyvä lähestymistapa.
310 % Ongelmaa ei voida välttää visualisointiresoluutiota lisäämällä, koska näytöllä on käytettävissä
311 % vain tietty määrä pikseleitä (tässä aika-akselin eli pystysuunnassa).
312 %
313 % Rehellisempää* olisi tallentaa minimi ja maksimi, ja visualisoida tämä vaihteluvälidata jotenkin.
314 % Tietysti 1+1D:ssä tämä on jo hankalaa eikä onnistu MATLABin standardiplottereilla. Yksi vaihtoehto
315 % olisi toteuttaa oma plotteri bittikarttatoiminnoilla (kuvatoiminnoilla), ja koodata esim. keskiarvo
316 % ja keskihajonta erikseen värin kirkkaus- ja saturaatiokanaviin.
317 %
318 % (* Ks. esim.
319 % Richard Fateman. Honest plotting, global extrema, and interval arithmetic. ISSAC'92, s.216--
320 % Ron Avitzur, Olaf Bachmann, Norbert Kajler. From honest to intelligent plotting. ISSAC'95, s.32--
321 % )
322 %
323 y_vis = y_vis / steps_per_vis;
307324
325 % Tallennetaan tulos piirtoa varten.
326 ys(vis_step,:) = y_vis';
327
328 % Aloitetaan uusi keskiarvoistus.
329 y_vis = zeros(s.ny,1);
330 vis_step = vis_step + 1;
331 end
332
308333 t = t + s.dt;
334 perc_prev = perc;
309335 end
336 t_elapsed = compat.toc(t_start);
337 fprintf(1, 'Ratkaisu valmis, aikaa kului %0.3gs.\n', t_elapsed);
310338
311339 fprintf(1, 'Visualisoidaan...\n');
340 t_start = tic;
312341
313342 figure(1);
314343 clf;
315344
316 [X,T] = meshgrid(s.xx, 0:s.dt:s.tmax);
345 [X,T] = meshgrid(s.xx, 0:dt_vis:s.tmax);
317346
318347 subplot(2,2, 1);
319% surf(X,T, ys(:,s.I_p), 'EdgeColor', 'none');
320% view(0,90);
321 contourf(X,T, ys(:,s.I_p));
348 f_plot(X, T, ys(:,s.I_p));
349 if ~isempty(f_postplot)
350 f_postplot();
351 end
322352 xlabel('x');
323353 ylabel('t');
324354 title('p(x,t)');
325 colorbar('Location', 'SouthOutside');
326 colormap gray;
355 compat.colorbar('SouthOutside');
327356 axis tight;
328357
358 % Huom. MATLABin/Octaven rajoitus: vain yksi colormap per figure, vaikka olisi subplottejakin.
359 % Asetamme tämän siis vain kerran.
360 colormap(cm);
361
329362 subplot(2,2, 2);
330363 phi_ys = s.phi_0 + ( (s.cB0 - ys(:,s.I_cB)) / s.cB0 ) * ( s.phi_infty - s.phi_0 );
331 contourf(X,T, phi_ys);
364 f_plot(X, T, phi_ys);
365 if ~isempty(f_postplot)
366 f_postplot();
367 end
332368 xlabel('x');
333369 ylabel('t');
334370 title('\phi(x,t)');
335 colorbar('Location', 'SouthOutside');
371 compat.colorbar('SouthOutside');
336372 axis tight;
373 caxis([0 1]); % ehkä havainnollisempi phi:lle kuin pelkkä toteutunut alue
337374
338375 subplot(2,2, 3);
339% surf(X,T, ys(:,s.I_cA), 'EdgeColor', 'none');
340% view(0,90);
341 contourf(X,T, ys(:,s.I_cA));
376 f_plot(X, T, ys(:,s.I_cA));
377 if ~isempty(f_postplot)
378 f_postplot();
379 end
342380 xlabel('x');
343381 ylabel('t');
344382 title('c_A(x,t)');
345 colorbar('Location', 'SouthOutside');
383 compat.colorbar('SouthOutside');
346384 axis tight;
347385
348386 subplot(2,2, 4);
349% surf(X,T, ys(:,s.I_cB), 'EdgeColor', 'none');
350% view(0,90);
351 contourf(X,T, ys(:,s.I_cB));
387 f_plot(X, T, ys(:,s.I_cB));
388 if ~isempty(f_postplot)
389 f_postplot();
390 end
352391 xlabel('x');
353392 ylabel('t');
354393 title('c_B(x,t)');
355 colorbar('Location', 'SouthOutside');
394 compat.colorbar('SouthOutside');
356395 axis tight;
357396
358 fprintf(1, 'Valmis.\n');
397 t_elapsed = compat.toc(t_start);
398 t_all_elapsed = compat.toc(t_all_start);
399 fprintf(1, 'Visualisointi valmis, aikaa kului %0.3gs.\n', t_elapsed);
400 fprintf(1, 'Yhteensä aikaa käytetty %0.3gs.\n', t_all_elapsed);
401
402 fprintf(1, '\nValmista tuli.\n');
359403end
360404
405%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406% Tehtävään liittyvät apufunktiot
407%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
408
361409% Yhtälön y_t = f(y, t) oikea puoli suurevektorille y.
362410%
363411% s = tehtäväparametrit (struct)
456456 y_t(s.I_p(end)) = (g(end)/s.beta) * p_xx_at_L_t;
457457end
458458
459%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
460% Yleiskäyttöiset apufunktiot
461%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
462
459463% Standardi Runge-Kutta 4. Integroi yhtälöä
460464%
461465% y_t = f(y, t)
490490% y = derivoitava kenttä (pystyvektoriformaatissa)
491491% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
492492%
493% Oletukset:
494% numel(y) >= 3
495% h > 0
496%
493497function y_x = ddx(y, h)
494498 % Keskiosa menee keskeisdifferenssillä vektoroidusti:
495499 yp = [ y(2:end) ; NaN ];
501501 y_x = (yp - ym) / (2*h);
502502
503503 % Käsitellään päät toisen kertaluvun forward- ja backward-differensseillä.
504 % (Backward: http://www.geometrictools.com/Documentation/FiniteDifferences.pdf
505 % Forward samasta kaavasta etumerkit vaihtamalla.)
504 %
505 % Backward: http://www.geometrictools.com/Documentation/FiniteDifferences.pdf
506 % Forward samasta kaavasta vaihtamalla etumerkit ja kääntämällä indeksointi
507 % takaperin, forward- ja backward-differenssien geometrian nojalla.
508 %
506509 y_x(1) = (-3*y(1) + 4*y(2) - y(3)) / (2*h);
507510 y_x(end) = ( 3*y(end) - 4*y(end-1) + y(end-2)) / (2*h);
508511end
517517% y = derivoitava kenttä (pystyvektoriformaatissa)
518518% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
519519%
520% HUOM! Tämä funktio ei käsittele kahta ensimmäistä alkiota. Tässä ratkaisijassa
521% tätä käytetään vain paineen derivaattareunaehdon kanssa loppupäässä
522% ( p_xx(L,t) = d_dx( p_x(L,t) ), missä oikean puolen p_x saadaan suoraan reunaehdosta
523% ja p_xx lasketaan tällä funktiolla yhdistäen reunaehto sekä pari sisempää pistettä
524% alueen sisältä (joiden arvot tunnetaan edelliseltä aika-askeleelta) ).
520% HUOM! Koska kyseessä on backward-differenssi, tämä funktio ei laske vektorin
521% kahden ensimmäisen alkion derivaattaa, vaan niiden derivaataksi asetetaan NaN!
525522%
523% Tässä ratkaisijassa tätä käytetään vain paineen derivaattareunaehdon kanssa loppupäässä.
524% p_xx(L,t) = d_dx( p_x(L,t), missä oikean puolen p_x saadaan suoraan reunaehdosta,
525% ja p_xx lasketaan tällä funktiolla yhdistäen reunaehto sekä pari sisempää pistettä
526% alueen sisältä. Sisäpisteiden arvot tunnetaan edelliseltä aika-askeleelta.
527%
528% Oletukset:
529% numel(y) >= 3
530% h > 0
531%
526532function y = ddx_bwd(y, h)
527533 ym1 = [ NaN ; y(1:end-1) ];
528534 ym2 = [ NaN ; NaN; y(1:end-2) ];
543543% y = derivoitava kenttä (pystyvektoriformaatissa)
544544% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
545545%
546% Oletukset:
547% numel(y) >= 4
548% h > 0
549%
546550function y_xx = d2dx2(y, h)
547551 % Sisäosa on jälleen helppo:
548552 yp = [ y(2:end) ; NaN ];
555555
556556 % Toisen kertaluvun forward- ja backward-differenssit päille.
557557 %
558 % Miten tämä tuotettiin:
558 % Miten nämä tuotettiin:
559559 % http://www.geometrictools.com/Documentation/FiniteDifferences.pdf
560560 %
561561 % Yhtälö (12) (vrt. example 1). d = 2, p = 2.
564564 % A = [ 1 1 1 1 ; 0 1 2 3 ; 0 1 4 9 ; 0 1 8 27 ];
565565 % b = [ 0 ; 0 ; 1 ; 0 ];
566566 % C = A\b
567 % 2*C
568567 % => C = 1/2 * (2, -5, 4, -1)
568 %
569 % Backward vastaavasti, tai vaihtamalla tuohon etumerkit ja kääntämällä
570 % indeksoinnin takaperin (forward- ja backward-differenssien geometrian nojalla).
569571 %
570 % Päiden käsittely tarvitaan tässä ratkaisijassa, koska muutoin
571 % NaNit propagoituisivat cB:hen (rk4:n sisällä, kun väliaskeleita lasketaan),
572 % Päiden käsittely tarvitaan tässä ratkaisijassa, koska muutoin NaNit
573 % propagoituisivat cB:hen (rk4:n sisällä, kun väliaskeleita lasketaan),
572574 % ja saataisiin phi = NaN NaN ...
573575 %
574576 y_xx(1) = 0.5*y(1) - 2.5*y(2) + 2*y(3) - 0.5*y(4);