Commit 983dbba5dda485175ab12c9096e64f0ac3e14560

  • Tree SHA1: ee84593
  • Parent SHA1: b04167d (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)
  • raw diff | raw patch
IBVPRK4 (bwith BC callback (gives better accuracy with BCs); online visualization; swapped x and t axes in spacetime sheets; no num_eps, just discard anything below 0; clarified comments somomewhat more
Harjoitustyo/corrosive_darcy_flow.m
(196 / 82)
  
9393
9494 s.phi_0 = 0.5; % huokoinen aine = 50% tyhjää alussa
9595 s.phi_infty = 0.9; % 1 - 0.9 = 10% syöpymätöntä ainetta
96 s.kappa = 1; % materiaalin permeabiliteetti
96 s.kappa = 10; % materiaalin permeabiliteetti
9797
9898 s.k1 = 0.2; % reaktionopeus, komponentti cB
9999 s.k2 = 0.1; % reaktionopeus, komponentti cA
100100 s.cB0 = 1; % huokoisen aineen syöpyvän osan konsentraatio alussa
101101
102102 s.mu = 1; % fluidin viskositeetti
103 s.beta = 1; % fluidin kompressoituvuuskerroin
103 s.beta = 5; % fluidin kompressoituvuuskerroin
104104 s.p_infty = 20; % ympäristön paine
105 s.v_0 = 1; % vapaa virtausnopeus (ennen huokoiseen aineeseen törmäystä)
105 s.v_0 = 2; % 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
109 s.tmax = 10.0; % simulaation loppuaika
109 s.tmax = 20.0; % simulaation loppuaika
110110 s.dt = 1e-2; % aika-askeleen koko samoissa yksiköissä kuin tmax
111111 % (vakio vain yksinkertaisuuden vuoksi; RK4 ei rajoita)
112 s.nx = 11; % paikkadiskretoinnin solmujen lukumäärä
112 s.nx = 21; % paikkadiskretoinnin solmujen lukumäärä
113113
114 nt_vis = 101; % visualisointipisteiden määrä aika-akselilla (kannattaa pitää kohtuullisena)
114 % Laskentatulosta voi visualisoida laskennan edistymisen aikana.
115 %
116 % Tämä hidastaa laskentaa huomattavasti, mutta voi olla havainnollista katseltavaa.
117 %
118 online_vis = 0; % laskennanaikainen visualisointi päälle/pois (1 = päällä, 0 = pois)
119 online_vis_wait = 0.01; % kauanko odotetaan laskennanaikaisessa visualisaatiossa piirron jälkeen,
120 % ennen kuin aletaan laskea seuraavaa aika-askelta. Arvo sekunteina.
115121
116 num_eps = 1e-4; % Numeerinen pyöristys, tätä pienemmät luvut = 0
122 % Laskennan lopuksi (joka tapauksessa) piirrettävän aika-avaruuslakanan piirtoasetukset.
123 %
124 nt_vis = 101; % Visualisointipisteiden määrä aika-akselilla (kannattaa pitää kohtuullisena)
125 cm = 'jet'; % Värikartta funktioiden arvoille. MATLABin/Octaven mukana tulevan värikartan nimi,
126 % tai vektori RGB-arvoja MATLABin colorbar-formaatissa.
117127
118 cm = 'jet'; % Värikartta plottia varten
119
120 % Piirtotapa. 'contourf' tai 'surf'.
128 % Piirtotapa, 'contourf' tai 'surf'.
121129 %
122130 % Tasa-arvokäyrien piirto voi olla hidasta (etenkin Octavessa, jos nt_vis on suuri),
123131 % joten toinen mahdollisuus on käyttää ylhäältä päin katsottua, värikoodattua 3D-pintaa.
124132 % Huomaa, että tällöin kunkin mesh-elementin väri määräytyy vasemman yläkulman (eikä keskipisteen)
125133 % funktioarvon mukaisesti.
126134 %
127 plot_method = 'contourf';
128% plot_method = 'surf';
135 % Tämä voi olla hyödyllistä myös, jos ratkaisu ei konvergoi (ja näin ollen osa alueesta jää alustamatta).
136 % Ainakin Octaven contourf hidastuu huomattavasti tai jopa jumiutuu, jos datassa on NaNeja.
137 %
138% plot_method = 'contourf';
139 plot_method = 'surf';
129140
130141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131142 % GNU Octave -yhteensopivuus
187187 f_plot = @(X,Y,Z) surf(X, Y, Z, 'EdgeColor', 'none');
188188 f_postplot = @() view(0,90);
189189 else
190 error('Tuntematon piirtotapa "%s". Tuettuja ovat "contourf" ja "surf".', plot_method);
190 error('Tuntematon piirtotapa "%s". Mahdolliset valinnat ovat "contourf" ja "surf".', plot_method);
191191 end
192192
193193 s.xx = transpose(linspace(0, s.L, s.nx)); % pystyvektori, helpompaa tämän koodin kannalta...
211211 % Huomaa, että pisteet ovat järjestyksessä kasvavan x:n mukaan
212212 % (tätä käytetään reunaehtojen käsittelyssä).
213213
214 % Aika-askelten määrä (mukaanlukien alkuhetki t = 0)
214 % Aika-askelten määrä, mukaanlukien alkuhetki t = 0
215215 s.nt = 1 + floor(s.tmax/s.dt);
216216
217217 % Visualisoitavan aikavälin pituus
218218 dt_vis = s.tmax / (nt_vis - 1);
219
219220 % Kuinka monta aika-askelta otetaan per yksi visualisoitava askel
220221 steps_per_vis = floor(dt_vis / s.dt);
221222
225225 %
226226 % Huomaa, että vain visualisoitava määrä aika-askelia tallennetaan.
227227 %
228 % HUOM: aika pystyakselilla (rivi), suurevektorin y arvo vaaka-akselilla (sarake).
228 % HUOM: rivi = sama aika-askel kaikille suureille,
229 % sarake = suurevektorin y arvo tietyssä pisteessä (ja tietylle suureelle) koko historian yli.
229230 %
230 ys = nan(nt_vis,s.ny);
231 ys = nan(nt_vis, s.ny);
231232
232233 % Ladotaan alkutila paikalleen:
233234 %
236236 y(s.I_p) = s.p_infty - (s.phi_0 * s.v_0 * s.mu / s.kappa) * s.xx;
237237 y(s.I_cA) = 0;
238238 y(s.I_cB) = s.cB0;
239 ys(1,:) = y'; % tulokset varastoidaan aika-askeleittain vaakavektoreina
240239
241 % Laskennan edistymismittaria varten (progress indicator)
240 % Toteutetaan reunaehdot myös alkuhetkellä:
241 y = bc_callback(s, y, 0);
242
243 % Kuten yllä mainittiin, tulokset varastoidaan aika-askeleittain vaakavektoreina.
244 % Tallennetaan lähtötilanne.
245 ys(1,:) = y';
246
247 % Alustetaan laskennan edistymismittari (progress indicator)
242248 num_computed_ts = s.nt - 1; % montako aika-askelta oikeasti lasketaan (huom. lähtötilanne jo tunnetaan!)
243249 perc = 0; % valmistunut prosenttiosuus (percentage complete)
244250 perc_prev = -1; % valmistunut prosenttiosuus edellisellä aika-askeleella
278278 end
279279 end
280280
281 % Asetetaan reunaehdot paikalleen.
282 %
283 % Vasemman reunan Dirichlet-ehdot ovat helppoja.
284 % Pakotetaan arvot ennen integrointia:
285 y(s.I_p(1)) = s.p_infty;
286 y(s.I_cA(1)) = s.cA0;
287
288 % Oikealla reunalla joutuu hiukan temppuilemaan, koska tarvitaan p:n arvo,
289 % ja käytettävissä on vain p_x(L,t). Toisaalta, lisäksi käytettävissä
290 % on p:n arvo edellisestä aika-askeleesta (tai 1. askeleella alkuarvo),
291 % ja yhtälöstä (p_t) saadaan p_t, kun tunnetaan p_x ja p_xx. Jälkimmäinen
292 % seuraa ensimmäisestä derivoimalla. Aikaderivaatasta puolestaan voidaan
293 % integroida p.
294 %
295 % Sijoitetaan siis näin saatava p_t(L,t) paikalleen funktioon f
296 % samalla tavalla kuin sisäosassakin, ja annetaan rk4:n hoitaa loput.
297 % Ks. funktion f toteutus.
298
299281 % Aikaintegroidaan yksi aika-askel.
300 y = rk4(y, t, s.dt, @(y,t) f(s,y,t));
282 %
283 % Tehtäväparametrit annetaan funktioille f() ja bc_callback() luomalla nimettömät apufunktiot.
284 % Integraattori pitää bc_callback()-funktion avulla huolta siitä, että reunaehdot toteutuvat.
285 %
286 y = rk4_ibvp(y, t, s.dt, @(y,t) f(s,y,t), @(y,t) bc_callback(s,y,t));
301287
302 % Pyöristys, jottei numeerinen epätarkkuus kaada laskua.
303 % Konsentraatioille voi numeerisesti tulla hyvin pieniä negatiivisia
304 % arvoja siellä, missä ne ovat lähellä nollaa.
305 %
306 y( abs(y) < num_eps ) = 0;
288 % HACK: konsentraatioille voi numeerisesti tulla pieniä negatiivisia
289 % arvoja siellä, missä ne ovat lähellä nollaa. Korjataan nämä pois.
290 %
291 % Tämä aiheuttaa laskuun epätarkkuutta (aikaevoluutioyhtälö ei toteudu tarkasti),
292 % mutta mahdollistaa (approksimatiivisen) laskemisen pienemmällä määrällä vapausasteita.
293 %
294 % (Huomaa, että koska yhtälö (p_t) on lähinnä parabolista tyyppiä, tiedetään von Neumann
295 % -stabiilisuusanalyysistä, että dt on oltava verrannollinen dx^2:een. Aika-askeleita
296 % tarvitaan siis hyvin äkkiä todella paljon, kun nx kasvaa. Suhteellisen pienellä nx:llä
297 % laskeminen on siis tälle ohjelmalle hyvinkin hyödyllinen ominaisuus.)
298 %
299 y( y < 0 ) = 0;
307300
308 % Tehdään pari sanity checkiä.
301 % Välituloksen sanity check.
302 %
309303 % Kaikki suureet p, cA, cB ovat 1) numeroita, 2) epänegatiivisia.
310 if sum(isnan(y)) || sum(y < 0)
311 % DEBUG
304 % Jälkimmäinen on pakotettu toteutumaan yllä, tarkistetaan vain ensimmäinen.
305 if sum(isnan(y))
312306 fprintf(1, '\nRatkaisuvirhe (sanity check fail). Muuttujien arvot:\n');
313307 p = y(s.I_p)
314308 cA = y(s.I_cA)
315309 cB = y(s.I_cB)
316310 k
317311 t
318 break; % visualisoidaan kuitenkin siihen asti kuin saatiin järkevää dataa.
312
313 % Visualisoidaan kuitenkin siihen asti kuin saatiin järkevää dataa.
314 % Mutta vain, jos plotteri on sellainen, joka ei jumiudu NaNeista...
315 if strcmp(plot_method, 'surf')
316 fprintf(1, 'Ratkaisu ei konvergoinut. Visualisoidaan ratkennut osa.\n');
317 break;
318 else
319 error('Ratkaisu ei konvergoinut.');
320 end
319321 end
320322
321323 % Tulosten tallennus.
329329 %
330330 % Kannattaa huomata, että yleensä tämä ehkä menettelee, muttei ole kovin hyvä lähestymistapa.
331331 % Ongelmaa ei voida välttää visualisointiresoluutiota lisäämällä, koska näytöllä on käytettävissä
332 % vain tietty määrä pikseleitä (tässä aika-akselin eli pystysuunnassa).
332 % vain tietty määrä pikseleitä (tässä aika-akselin eli vaakasuunnassa).
333333 %
334334 % Rehellisempää* olisi tallentaa minimi ja maksimi, ja visualisoida tämä vaihteluvälidata jotenkin.
335335 % Tietysti 1+1D:ssä tämä on jo hankalaa eikä onnistu MATLABin standardiplottereilla. Yksi vaihtoehto
343343 %
344344 y_vis = y_vis / steps_per_vis;
345345
346 if online_vis
347 figure(2);
348 clf;
349
350 subplot(2,2, 1);
351 plot(s.xx, y_vis(s.I_p));
352 xlabel('x');
353 % t+dt, koska visualisoitava arvo kuvaa aika-askeleen lopun tilannetta.
354 my_label = sprintf('p(x,t), t = %0.3g', t+s.dt);
355 ylabel(my_label);
356 axis tight;
357
358 subplot(2,2, 2);
359 phi_ys = s.phi_0 + ( (s.cB0 - y_vis(s.I_cB)) / s.cB0 ) * ( s.phi_infty - s.phi_0 );
360 plot(s.xx, phi_ys);
361 xlabel('x');
362 my_label = sprintf('\\phi(x,t), t = %0.3g', t+s.dt);
363 ylabel(my_label);
364 axis([s.xx(1) s.xx(end) 0 1]); % ehkä havainnollisempi phi:lle kuin pelkkä toteutunut alue
365
366 subplot(2,2, 3);
367 plot(s.xx, y_vis(s.I_cA));
368 xlabel('x');
369 my_label = sprintf('c_A(x,t), t = %0.3g', t+s.dt);
370 ylabel(my_label);
371 axis tight;
372
373 subplot(2,2, 4);
374 plot(s.xx, y_vis(s.I_cB));
375 xlabel('x');
376 my_label = sprintf('c_B(x,t), t = %0.3g', t+s.dt);
377 ylabel(my_label);
378 axis tight;
379
380 drawnow;
381 pause(online_vis_wait);
382 end % online_vis
383
346384 % Tallennetaan tulos piirtoa varten.
347385 ys(vis_step,:) = y_vis';
348386
404404 [X,T] = meshgrid(s.xx, 0:dt_vis:s.tmax);
405405
406406 subplot(2,2, 1);
407 f_plot(X, T, ys(:,s.I_p));
407 f_plot(T, X, ys(:,s.I_p));
408408 if ~isempty(f_postplot)
409409 f_postplot();
410410 end
411 xlabel('x');
412 ylabel('t');
411 xlabel('t');
412 ylabel('x');
413413 title('p(x,t)');
414414 compat.colorbar('SouthOutside');
415415 axis tight;
420420
421421 subplot(2,2, 2);
422422 phi_ys = s.phi_0 + ( (s.cB0 - ys(:,s.I_cB)) / s.cB0 ) * ( s.phi_infty - s.phi_0 );
423 f_plot(X, T, phi_ys);
423 f_plot(T, X, phi_ys);
424424 if ~isempty(f_postplot)
425425 f_postplot();
426426 end
427 xlabel('x');
428 ylabel('t');
427 xlabel('t');
428 ylabel('x');
429429 title('\phi(x,t)');
430 caxis([0 1]); % ehkä havainnollisempi phi:lle kuin pelkkä toteutunut alue
430431 compat.colorbar('SouthOutside');
431432 axis tight;
432 caxis([0 1]); % ehkä havainnollisempi phi:lle kuin pelkkä toteutunut alue
433433
434434 subplot(2,2, 3);
435 f_plot(X, T, ys(:,s.I_cA));
435 f_plot(T, X, ys(:,s.I_cA));
436436 if ~isempty(f_postplot)
437437 f_postplot();
438438 end
439 xlabel('x');
440 ylabel('t');
439 xlabel('t');
440 ylabel('x');
441441 title('c_A(x,t)');
442442 compat.colorbar('SouthOutside');
443443 axis tight;
444444
445445 subplot(2,2, 4);
446 f_plot(X, T, ys(:,s.I_cB));
446 f_plot(T, X, ys(:,s.I_cB));
447447 if ~isempty(f_postplot)
448448 f_postplot();
449449 end
450 xlabel('x');
451 ylabel('t');
450 xlabel('t');
451 ylabel('x');
452452 title('c_B(x,t)');
453453 compat.colorbar('SouthOutside');
454454 axis tight;
467467
468468% Yhtälön y_t = f(y, t) oikea puoli suurevektorille y.
469469%
470% Vain alueen sisäosa; tässä funktiossa ei käsitellä reunaehtoja.
471%
470472% s = tehtäväparametrit (struct)
471% y = kentän arvo hetkellä t
473% y = kentän y arvo hetkellä t
472474% t = aikamuuttujan arvo hetkellä t
473475%
474476function y_t = f(s, y, t)
494494 % Lasketaan ensin apusuureita (kaikki, mitä tarvitaan enemmän kuin kerran).
495495 p_x = ddx(y(s.I_p),s.h);
496496
497 % Sisäosa:
497 % Yhtälöiden oikea puoli on siis:
498498 y_t = [ g.*( p_x.^2 ) + 1/s.beta * ( ddx(g,s.h).*p_x + g.*d2dx2(y(s.I_p),s.h) ) ; ... % p_t
499499 (s.kappa/s.mu) * p_x .* ddx(y(s.I_cA),s.h) - s.k2 * y(s.I_cB) .* y(s.I_cA) ; ... % cA_t
500500 -s.k1*y(s.I_cA).*y(s.I_cB) ... % cB_t
501501 ];
502end
502503
503 % Käsitellään paineen reunaehto alueen loppureunalla x = L.
504 %
505 % Tiedetään p_x(L,t) = 0. Siispä yhtälö (p_t) tulee muotoon
506 %
507 % p_t(L,t) = g/beta * p_xx(L,t)
508 %
509 % missä p_xx(L,t) = ddx(p_x(L,t)).
510 %
511 p_x(end) = 0; % reunaehdosta!
512 % (Tarvittavat pari sisempää pistettä seuraavat derivaattakentästä alueen sisällä, joka laskettiin jo yllä.)
513 p_xx_at_L_t = ddx_bwd(p_x((end-2):end),s.h); % tarvitaan vain loppupisteen arvo; ei lasketa turhaan alueen sisällä
514 p_xx_at_L_t = p_xx_at_L_t(end); % poimitaan laskun tulos vektorista
504% Pakottaa reunaehdot toteutumaan suurevektorissa y.
505%
506% Huomaa, että tämän tehtävän reunaehdot eivät riipu ajasta,
507% joten parametri t on käyttämätön.
508%
509function y = bc_callback(s, y, t)
510 % Alkureuna x = 0.
511 %
512 % Dirichlet-ehdot ovat helppoja:
513 y(s.I_p(1)) = s.p_infty;
514 y(s.I_cA(1)) = s.cA0;
515515
516 % Siispä reunaehto toteutuu, kun määräämme, että
517 y_t(s.I_p(end)) = (g(end)/s.beta) * p_xx_at_L_t;
516 % Loppureuna x = L.
517 %
518 % Tässä helpointa on käyttää (esim. 2. kertaluvun) backward-differenssiä
519 % p_x:lle, josta voidaan ratkaista viimeisen solmun arvo, kun reunaehto
520 % sanoo, että p_x = 0:
521 %
522 % p_x(L,t) = ( 3*p(L,t) - 4*p(L-h,t) + p(L-2h,t) ) / (2h)
523 %
524 % mistä saadaan:
525 %
526 y(s.I_p(end)) = 2*s.h * ( 4*y(s.I_p(end-1)) - y(s.I_p(end-2)) ) / 3;
518527end
519528
520529%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
521530% Yleiskäyttöiset apufunktiot
522531%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
523532
524% Standardi Runge-Kutta 4. Integroi yhtälöä
533% Runge-Kutta 4. Tämä toteutus on suunniteltu alku-reuna-arvotehtäville (initial boundary value problem, IBVP).
525534%
535% Tämä integroi yhtälöä
536%
526537% y_t = f(y, t)
527538%
528539% yhden aika-askeleen. Suure y saa olla vektori, pysty- tai vaakaformaatilla
529540% ei ole merkitystä.
530541%
542% Poiketen pelkästä RK4:stä (joka integroi pelkästään tavallisten differentiaaliyhtälöiden systeemiä),
543% käyttäjän antamat reunaehdot pakotetaan toteutumaan jokaisella väliaskeleella.
544%
545% Huomaa, että tästä huolimatta vapausasteiden tulkinta alku-reuna-arvotehtävässä on täysin käyttäjän määräämä.
546% Tämä funktio "näkee" tehtävän ainoastaan systeeminä tavallisia differentiaaliyhtälöitä, kuten tavallinenkin RK4.
547% Ainoa ero on, että käyttäjän määräämä funktio pääsee muokkaamaan RK4:n välituloksia sellaisiksi, että ne
548% toteuttavat tehtävän reunaehdot. Tämä tarkentaa reuna-alkuarvotehtävän integrointia, estämällä välituloksia
549% ajautumasta omille teilleen reunaehtojen suhteen. Tämä puolestaan mahdollistaa tehtävän ratkaisun pienemmälläkin
550% määrällä vapausasteita.
551%
531552% Parametrit:
532553%
533554% y0 = suureen y arvo aika-askeleen alussa
534555% t0 = aika-askeleen alkuhetki
535556% dt = aika-askeleen pituus
536557% f = f(y, t) = yhtälön oikea puoli, function handle.
558% bc_callback = bc_callback(y,t) = funktio, joka pakottaa käyttäjän määräämät reunaehdot
559% toteutumaan datassa y hetkellä t.
537560%
538function y = rk4(y0, t0, dt, f)
561function y = rk4_ibvp(y0, t0, dt, f, bc_callback)
539562 k1 = f(y0, t0);
540 k2 = f(y0 + 1/2 * dt * k1, t0 + 1/2 * dt);
541 k3 = f(y0 + 1/2 * dt * k2, t0 + 1/2 * dt);
542 k4 = f(y0 + dt * k3, t0 + dt);
543 y = y0 + 1/6 * dt * (k1 + 2*k2 + 2*k3 + k4);
563
564 y1 = y0 + 1/2 * dt * k1;
565 t1 = t0 + 1/2 * dt;
566 y1 = bc_callback(y1, t1);
567
568 k2 = f(y1, t1);
569
570 y2 = y0 + 1/2 * dt * k2;
571 t2 = t0 + 1/2 * dt;
572 y2 = bc_callback(y2, t2);
573
574 k3 = f(y2, t2);
575
576 y3 = y0 + dt * k3;
577 t3 = t0 + dt;
578 y3 = bc_callback(y3, t3);
579
580 k4 = f(y3, t3);
581
582 y4 = y0 + dt * k4;
583 t4 = t0 + dt;
584 y4 = bc_callback(y4, t4);
585
586 % Standardi RK4:
587 % y = y0 + 1/6 * dt * (k1 + 2*k2 + 2*k3 + k4);
588
589 % Ekvivalentti ylläolevan kanssa, mutta toteuttaa reunaehdot paremmin:
590 % (Toivotaan, ettei y0:n vähentäminen aiheuta floating point -pyöristysongelmia.)
591 %
592 y = y0 + (1/3) * (y1-y0) + (2/3) * (y2-y0) + (1/3) * (y3-y0) + (1/6) * (y4-y0);
593
594 % Ja lopuksi varmistetaan vielä, että...
595 tf = t0 + dt;
596 y = bc_callback(y, tf);
544597end
545598
546599% Laskee skalaarikentän ensimmäisen paikkaderivaatan