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.
|