Chalmers


Datorstudioövning 3: Optimering av kod

Målsättningen med dagens datorstudioövning är att snabba upp vår Matlab-kod. Vi kommer att använda oss av följande testproblem:


        du/dt - div(grad(u)) = x*(1-x)*y*(1-y)*cos(t) +

                         2*(x*(1-x) + y*(1-y))*sin(t),

  (1)     i Omega = [0,1] x [0,1], för 0 < t < pi/2,

        u(x,y,t) = 0, på randen till Omega, för 0 < t < pi/2,

        u(x,y,0) = 0, i Omega.

Den exakta lösningen är (verifiera detta)

  u(x,y,t) = x*(1-x)*y*(1-y)*sin(t)

För att kunna jämföra resultaten med varandra föreslås ni alla att använda följande diskretisering i rum och tid:

Rumsdiskretisering. Den s.k. standardtrianguleringen av enhetskvadraten i 128 stycken deltrianglar. Du skapar den så här i pdetool: Rita först Omega = [0,1] x [0,1] (enhetskvadraten). Innan du triangulerar Omega väljer du Parameters... från menyn Mesh. I rutan Maximum edge size skriver du inf och klickar på OK. Klicka nu på ikonen med en triangel för att triangulera. Du får nu en mesh som består av två trianglar! (Vi har ju just angivit att vi tillåter oändligt stora triangelsidor!). Förfina tre gånger. Du har nu skapat en så kallad standardtriangulering (med, i detta fall, 128 stycken trianglar). Exportera nu som vanligt matriserna p, e, och t till Matlab:s workspace.

Tidsdiskretisering. Använd en likformig indelning av tidsintervallet [0, pi/2] i femton stycken tidssteg med längd dt = pi/30.

Definiera nu samtliga funktionsfiler som behövs för att specificera data till problemet (1), d.v.s., för att specificera koefficienter till differentialekvationen, randvillkor samt begynnelsevillkor. Lös nu (1) med din tidsberoende lösare från DTA. Spara lösningen vid sluttiden pi/2 i kolonnvektorn U_slut. Tips: Flera av er har nog lösningen sparad kolonnvis (d.v.s. varje kolonn innehåller lösningen vid en viss tidpunkt) i en stor matris, med begynnelsevärdet i kolonn 1, o.s.v. Lösningen vid sluttiden ligger därmed i kolonn 16 (efter första tidssteget skriver vi ju i kolonn 2, o.s.v.) och vi kan, om vi antar att matrisen heter U (den kanske heter xi i ditt program), beräkna: >> U_slut = U(:, 16) (där ":" svarar mot samtliga rader).

Jämför nu U_slut med den exakta lösningen vid tiden pi/2:


  >> tid = pi/2;                       % Använd inte namnet t... :)
  >> x = p(1,:)';                      % Kolonnvektor.
  >> y = p(2,:)';
  >> u = x.*(1-x).*y.*(1-y).*sin(tid); % Blir också kolonnvektor...
  >> fel = u - U_slut;                 % ... så dimensionerna
                                       % stämmer överens!
  >> pdesurf(p, t, u)                  % plottar den exakta lösningen
  >> pdesurf(p, t, U_slut)             % plottar den approximativa lösningen
  >> pdesurf(p, t, fel)                % plottar felet
  >> max(abs(fel))                     % approximerar max-normen av felet

Om allt verkar vara i sin ordning är det nu dags att börja med optimeringen. VIKTIGT: Varje gång du gjort en förändring av ditt program bör du jämföra lösningen (vid tiden pi/2) som du får med den "uppsnabbade" koden, med U_slut för att kontrollera att de stämmer överens. Det viktigaste är ju ändå att programmet räknar rätt!

Vi skall börja med att "ta tid" på ditt nuvarande program. Lägg till raderna


  starttid = cputime;
  tic

först i din lösare, samt raderna

  sluttid = cputime;
  cputid = sluttid - starttid
  toc

sist i din lösare. Kommandona tic och toc ger dig helt enkelt en utskrift av elapsed_time, d.v.s, den tid (i sekunder) som ditt program tagit att köra. (Se >> help tic). Kommandot cputime returnerar CPU-tiden (i sekunder) som Matlab använt sedan Matlab startades, så cputid ger alltså den CPU-tid (i sekunder) som ditt program utnyttjat. (Se >> help cputime). Det kan vara intressant att jämföra dessa: om man är många användare som utnyttjar samma dator, eller om en användare har många processer igång samtidigt, är den senare förmodligen ett mer rättvist mått på ditt programs effektivitet.

Uppdelning av assemblerare. Eftersom det enda tidsberoendet hos data i ekvation (1) ovan återfinns i källtermen (övriga koefficienter samt randvillkor beror ej av tiden) vore det önskvärt att "en gång för alla", d.v.s., innan vi träder in i tidsloopen i vår lösare, assemblera samtliga matriser utom lastvektorn, vilken vi förstås är tvungna att assemblera om i varje tidssteg.

Mer allmänt kan vi tänka oss en situation där ett antal av matriserna kan beräknas en gång för alla, medan övriga måste räknas om i varje tidssteg. Hur kan vi åstadkomma detta?

Ett sätt är att vi skapar ett antal "små" assemblerare, där var och en har ansvar för "sin" matris. T.ex. skulle vi kunna skapa en funktion, minStyvhetsAssemblerare.m, som assemblerar styvhetsmatrisen:


  function A = minStyvhetsAssemblerare(p, t, tid);

  nnodes = size(p,2);         % Antal noder. (= antal kolonner i p)
  ntri = size(t,2);           % Antal trianglar. (= antal kolonner i t)

  A = zeros(nnodes, nnodes);  % Initiering av styvhetsmatris.

  for n = 1:ntri              % loopa över samtliga trianglar

  % Ta reda på hörnens nodnummer:  
  
    i = t(1,n);          % Nodnummer för hörn nummer 1 i triangel nummer n.
    j = t(2,n);          % Nodnummer för hörn nummer 2 i triangel nummer n.
    k = t(3,n);          % Nodnummer för hörn nummer 3 i triangel nummer n.
                         %
                         % 1 -> 2 -> 3 moturs riktning.
  
  % Ta reda på hörnens koordinater:		       
		       
    xnodi = p(1,i);      % x-koordinat för nod nummer i.
    xnodj = p(1,j);      % x-koordinat för nod nummer j.
    xnodk = p(1,k);      % x-koordinat för nod nummer k.
  
    ynodi = p(2,i);      % y-koordinat för nod nummer i.
    ynodj = p(2,j);      % y-koordinat för nod nummer j.
    ynodk = p(2,k);      % y-koordinat för nod nummer k.
		       
  % Räkna ut triangelns area:
  
    triangelarea = ((xnodj - xnodi)*(ynodk - ynodi) - ...
                    (xnodk - xnodi)*(ynodj - ynodi))/2;


    delomradesnummer = t(4,n);  % Rad 4 i t innehåller numret på delområdet
                                % i vilket triangel nummer n ligger.

  % Räkna ut diffusionskoefficienten i triangelns hörn vid tiden "tid":
  % (Vi skickar med numret på det delområde i vilket triangeln ligger
  % som en parameter, för att på ett enkelt sätt kunna specificera olika
  % uttryck för a i skilda delområden.)

    anodi = a_t_2D(xnodi, ynodi, tid, delomradesnummer);  % Värdet i nod i.
    anodj = a_t_2D(xnodj, ynodj, tid, delomradesnummer);  % Värdet i nod j.
    anodk = a_t_2D(xnodk, ynodk, tid, delomradesnummer);  % Värdet i nod k.

  % Räkna ut (de konstanta) gradienterna av phi_i, phi_j, och phi_k:
  
    Dphi_iDx = (ynodj - ynodk)/(2*triangelarea);
    Dphi_iDy = (xnodk - xnodj)/(2*triangelarea);
    gradphi_i = [Dphi_iDx; Dphi_iDy];
  
    Dphi_jDx = (ynodk - ynodi)/(2*triangelarea);
    Dphi_jDy = (xnodi - xnodk)/(2*triangelarea);
    gradphi_j = [Dphi_jDx; Dphi_jDy];
  
    Dphi_kDx = (ynodi - ynodj)/(2*triangelarea);
    Dphi_kDy = (xnodj - xnodi)/(2*triangelarea);
    gradphi_k = [Dphi_kDx; Dphi_kDy];
  
% Hörnkvadratur:

    A(i,i) = A(i,i) + ...
        (anodi*dot(gradphi_i,gradphi_i) + ...
         anodj*dot(gradphi_i,gradphi_i) + ...
         anodk*dot(gradphi_i,gradphi_i))*triangelarea/3;

    % här skall övriga åtta termer stå!

  end

Notera hur vi helt enkelt "klippt ut" just de delar i den "stora" assembleraren som behövs för att assemblera A, samt "klistrat in" dem i minStyvhetsAssemblerare.m. För att kontrollera att rutinen fungerar (viktigt!!!) kan vi t.ex. använda huset från förra studioövningen:

  >> p_hus = ...; % "point matrix"; den du skapade på förra övningen
  >> t_hus = ...; % "triangle matrix"
  >> tid = 0;
  >> A = minStyvhetsAssemblerare(p_hus, t_hus, tid)

Notera: Innan du anropar minStyvhetsAssemblerare måste du ha skapat en funktionsfil, a_t_2D.m, som definierar den (potentiellt tidsberoende) diffusionskoefficienten (du kan förstås kalla filen för vad du vill, och anropa funktionen med det namn du valt):

  function diffusion = a_t_2D(x, y, tid, delomradesnummer)

  % Funktionen  a(x,y,t)
  %
  % Exempel:  diffusion = 1;

  diffusion = 1;

Övning: Gör nu samma sak själv för övriga matriser och vektorer: "Tidsmassmatris", "Reaktiv massmatris", "Randbidragsmatris", "lastvektor", "randbidragsvektor". Kontrollera samtliga fem funktioner genom att göra "hustestet" ovan!

Tips: 1. Inga av matriserna i "husproblemet" beror av tiden. Det kan dock ändå vara bra att motsvarande assemblerare, liksom minStyvhetsAssemblerare.m, är "förberedda" för ett (potentiellt) tidsberoende. Om, liksom i "husproblemet" samt i problem (1), diffusionskoefficienten inte beror på tiden kan vi skicka in vad som helst i tidsargumentet (det används ju aldrig), t.ex. (som ovan) en 0:a. 2. Diffusionskoefficienten är densamma, a(x,y,t) = 1, i "husproblemet" och i problem (1). Tänk dock på att så i allmänhet ej är fallet: du måste alltså komma ihåg att ändra dina funktionsfiler från problem (1) till "husproblemet" när du gör dessa tester (och ändra tillbaka när du är klar; lite jobbigt kanske men bra träning! :). 3. Detsamma gäller matriserna p, t, och e. Ett smidigt sätt att komma runt detta är dock att (som ovan) kalla dem något annat, t.ex. p_hus, t_hus, och e_hus, under testet med "husproblemet". 4. Vi hade ingen "reaktiv massmatris" i "husproblemet" (och vi har det ej heller i problem (1)) men du kan testa "dess" assemblerare så här: Sätt "reaktionshastighetskoefficienten" (c(x,y,t) i termen c*u) till c = 1. Då skall ju den "reaktiva massmatrisen" och "tidsmassmatrisen" sammanfalla, och du kan alltså jämföra med "tidsmassmatrisen".

OBS! Detta betyder inte att vi "slänger ut" vår gamla assemblerare ur kursen! Det kan förekomma problem för vilka det är bäst att ha all assemblering samlad i en funktion.

Prova nu slutligen att använda dig av dina nya assembleringsrutiner för att lösa problem (1) ovan. Skriv alltså om din "Solver" så att funktionen minStyvhetsAssemblerare.m, samt assemblerarna för övriga tidsoberoende matriser, anropas innan du går in i tidsloopen, medan enbart minLastAssemblerare.m (eller vad du kallat den) anropas i varje tidssteg. Lös sedan (1) och ta tiden liksom ovan. Hur mycket snabbare blev programmet? (OBS! Glöm inte att kontrollera att din nya lösare räknar rätt! Jämför med U_slut.)

Tips: Skapa en kopia av din gamla "Solver" under ett nytt namn och gör ändringarna i denna. Dels kan det vara bra att ha en version av varje sort, dels skall man alltid se till att man har en fungerande version sparad när man sätter igång att skriva om sin kod, eftersom det finns en viss risk att man gör misstag... :) Det är då bra att ha en fungerande kod att gå tillbaka till, samt att jämföra resultaten med.

LU. Eftersom, för problem (1), samtliga matriser i vänsterledet är tidsoberoende och tidssteget är konstant, kan vi använda LU-faktorisering av "totalmatrisen" i vänsterledet. Om vi antar att tidsmassmatrisen kallas för M, styvhetsmatrisen kallas för A, den reaktiva massmatrisen kallas för Mc, och randbidragsmatrisen kallas för K, kan du beräkna LU-faktoriseringen av "totalmatrisen" genom att göra följande anrop efter att du assemblerat de tidsoberoende matriserna, men innan du går in i tidsloopen:


  dt = pi/30;
  A_tot = M + dt*(A + Mc + K);
  [Lower, Upper] = lu(A_tot);   % A_tot = Lower*Upper

(Det är möjligt att du tidigare kallat "totalmatrisen" för något annat, t.ex. VL.)

(Kommentar: Om du tittar på matriserna kanske du märker att Lower inte alltid är nedåt triangulär; det beror på att Matlab av numeriska skäl ibland väljer att byta plats på rader. Men det går, om man vill, att enkelt "byta tillbaka" så att Lower blir nedåt triangulär "på riktigt". Se: >> help lu)

Nu kan du inuti tidsloopen enkelt lösa det linjära ekvationssystemet i varje tidssteg (om vi vidare antar att lastvektorn kallas för F, randbidragsvektorn för G, aktuell lösning kallas för U_old och den sökta lösningen, när vi tagit ett tidssteg dt, kallas för U_new):


  b = M*U_old + dt*(F + G);
  y = Lower\b;
  U_new = Upper\y;

(Det är möjligt att du tidigare kallat "totalhögerledet" för något annat, t.ex. hl.)

Kommentar: Det kan kanske verka konstigt att det skulle vara "enklare" att göra på detta sätt, istället för att direkt skriva U_new = A_tot\b. Nu måste vi ju lösa två linjära ekvationssystem (istället för ett) i varje tidssteg. Poängen är att matriserna Lower och Upper är nedåt respektive uppåt triangulära. Därmed blir systemen mycket enklare att lösa för Matlab som använder Gauss-elimination: Tänk efter hur Gauss-elimination fungerar när matrisen är triangulär! (Möjligen blir det i detta fall med relativt få, 15, tidssteg och litet antal, 128, trianglar, inte en så stor tidsvinst. Testa i sådana fall, om du vill, med en finare mesh, t.ex. med 512 eller 2048 trianglar, och eventuellt ett mindre dt, så får du se!)

Formatet "Sparse". Eftersom matriserna är glesa (d.v.s. innehåller mest 0:or) kan du prova att deklarera dem som "sparse". Det gör du genom att t.ex. byta ut raden


  A = zeros(nnodes, nnodes);  % Initiering av styvhetsmatris.  

mot

  A = sparse(nnodes, nnodes);  % Initiering av styvhetsmatris.  

i funktionen minStyvhetsAssemblerare.m. (Se >> help sparse). Analogt för övriga tre matriser: M (tidsmassmatris), Mc (reaktiv massmatris), och K (randbidragsmatris).

Eftersom matriserna som uppkommer i en Finit Element-metod (FEM) typiskt är glesa är detta ett naturligt format: det lagrar och räknar endast med nollskilda element, vilket borde ge en tidsvinst. Men allting har ett pris: det åtgår också arbete till att hålla reda på vilka element som är skilda från noll och denna "book-keeping" tar också tid! Förhoppningsvis blir det en netto-vinst! Om inte så kan du också i detta fall prova med en finare triangulering där 0:orna blir ännu mer "dominanta".

Också bortsett från en (eventuell) tidsvinst har detta format en stor fördel: eftersom (som sagt) endast nollskilda element lagras fordras mindre minne vilket kan vara mycket betydelsefullt för att kunna räkna på stora problem, för vilka ens tillgängliga minne i datorn annars inte räckt till.

Kompakt (matris)notation. För de som vill: Studera sidan 24 i föreläsningsanteckningarna om Kodoptimering. Försök att realisera detta i någon eller några av dina assemblerare. Man kan vinna mycket tid på att utnyttja Matlab:s möjligheter att effektivt räkna med matriser. (Kommentar: Att skapa en 1 x ntri-vektor, area, innehållande alla triangelareor en gång för alla i början av sitt program, kan också vara bra av andra skäl: nu utför vi ju en hel del dubbelarbete genom att beräkna dessa varje gång en assemblerare anropas. Istället kan du då skicka med area som ett argument. Du kanske kan hitta (och eliminera!) fler liknande beräkningar som utförs "dubbelt"?!)

Kölista för studiopassen finns här: KD1,KD2.

Last modified: Thu Oct 24 14:28:22 MET DST 2002