Pokročilé programovací techniky v Matlabu T. Kozubek
Katedra aplikované matematiky Fakulta elektrotechniky a informatiky VŠB-Technická univerzita Ostrava
4.11. 2010
MI21
1
Obsah 1.
Matlab 1. 2. 3. 4. 5. 6. 7.
Funkce, řízení toku Homogenní pole (číselné matice, řetězce) Struktury Pole buněk OOP GUI Import, export, MEX
2.
Fourierova řada (spektrální analýza)
3.
Diskrétní Fourierova transformace (spektrální analýza)
4.
Konvoluce (vyhlazování signálu)
5.
Laplaceova transformace (řešení DR)
6.
Úlohy mechaniky – MatSol, paralelní programování
4.11. 2010
MI21
2
Matlab – funkce function [mean_x,std_x,var_x]=statistics(x) % STATISTICS Spocte zakladni statistiky vstupniho vektoru x % UZITI: [mean,std,var]=statistics(x) % Test na pocet vstupnich argumentu if nargin~=1; error('Spatny pocet vstupnich argumentu!'); end; % Test na korektnost vstupnich argumentu (x je ciselny vektor) if ~isvector(x); error('Vstupni argument x neni vektor!'); end; if ~isnumeric(x); error('Vstupni argument x neni ciselny vektor!'); end; % Spocti zakladni statistiky mean_x=mean(x); std_x=std(x); var_x=var(x); end % Podfunkce MEAN (upraveny vypocet stredni hodnoty) function f=mean(x) % Vypocet stredni hodnoty f=sum(x)/(length(x)+1); end 4.11. 2010
MI21
3
Matlab – homogenní pole Vytvoření matice
Přístup k prvkům, rušení
A=[1 2 3 234 5 6 7];
A(1:2,3)=[1;3]
B=[0,0,0; 2,3,4; 5,6,7];
A(2,:)=[]
A(A<7)=0
C=[A,B; B,A]; S=[‘abc’;’def’]
Řídké matice e=ones(5,1); A = spdiags([-e 2*e -e], -1:1, 5, 5) B=sprand(5,4,0.1) I=[1 1 2 3]; J=[1 3 2 4]; V=[1 1.5 2 3.7]; S=sparse(I,J,V,5,6) spy(S,'o',8); imagesc(S); colorbar; 4.11. 2010
MI21
4
Matlab - struktury • Struktura je heterogenní pole, jehož prvky jsou libovolné objekty – instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím názvů položek užitých v definici struktury. >>student=struct('jmeno','Jan Madlo', 'predmet',‘LAM','znamky',[1 2 1]); • Paměťové požadavky: není třeba souvislý blok paměti pro celou strukturu, ale pouze pro položky struktury. • Využití: snížení počtu vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu. • Přístup k položkám >>student.jmeno >>student.jmeno=‘Miro Vilik' 4.11. 2010
MI21
5
Matlab – pole buněk • Pole buněk je heterogenní pole, jehož prvky jsou libovolné objekty – instance tříd (řetězce, číselné pole, buněčné pole, struktury, atd). • Přístup k položkám je prostřednictvím indexů buněk. >>c={‘Hello’,rand(3),student}; • Paměťové požadavky: není třeba souvislý blok paměti pro celé pole, ale pouze pro jeho položky. • Využití: uložení řetězců různé velikosti, uložení matic různého řádu u metody rozložení oblastí, nahrazení vstupních a výstupních argumentů funkce, zvýšení čitelnosti kódu. • Přístup k položkám >>c{1}, c{1}(2:end), c{3}.jmeno >>c{1}=‘Ahoj’, c{3}.jmeno=‘Ivan Leden’ 4.11. 2010
MI21
6
Matlab – OOP
• • • •
MATLAB má zabudovánu kompletní podporu OO přístupu Jak systém roste, je čím dál složitější sledovat předávání dat mezi funkcemi. Opakování kódu copy&paste – redundance – nekonzistence. Deklarace třídních proměnných – předem se ví, jaké má proměnné. Srozumitelnost a přirozenost návrhu.
4.11. 2010
MI21
7
Matlab – OOP classdef MyClass < MySuperClass properties Property1 Property2 = sin(pi/12); % default value end properties (SetAccess = private, GetAccess = private) Stress Strain end methods function value = get.Property1(obj) %optionally implement get value end function obj = set.Property1(obj,value) %optionally implement set value end end end Dědičnost, přetěžování operátorů, konstruktory, proměnné s omezeným přístupem. 4.11. 2010
MI21
8
Matlab – GUI 1. Pomocí LAYOUT editoru (.fig, .m)
2. Pomocí Matlabovských příkazů (figure, uicontrol, uibuttongroup, guidata, guihandles, uimenu, …)
4.11. 2010
MI21
9
Matlab – import, export, MEX 1. Podpora práce s textovými i binárními soubory 2. Načtení obrázků, zvuků, videa, MS Excell, … Načtení a zobrazení obrázku rgb = imread('hokej.jpg'); image(rgb);
Načtení a přehrání zvukového souboru [y,fs,nbits]=wavread('babycry.wav'); sound(y,fs); plot(y);
MEX – přilinkování k Matlabu samostatných funkcí i celých knihoven v C, C++, Fortranu, Javě
4.11. 2010
MI21
10
Lineární zobrazení
f (z ) = a z, z ∈C, a ∈C f ( z) = a z = a e
i arg(a)
4.11. 2010
i arg(z )
ze
MI21
i (arg(a)+arg(z ))
=a ze
11
Fourierova řada (spektrální analýza)
{e
i ω nt
ω =
/
T
}
n∈ Z
... ortonorm.
2π ... úhlová T
f ( t ) ∈ L 2 ( 0 , T ),
systém
rychlost, f (t ) =
fcí L 2 ( 0 ,T )
T ... perioda +∞
∑
−∞
c n e i ω nt , t ∈ R
1 T − i ω nt ( ) cn = f t e dt , n ∈ Z , ∫ T 0
• Dvoustranné spektrum : {cn,ϕn}
4.11. 2010
MI21
12
Diskrétní Fourierova transformace (spektr. analýza) Fourierova řada +∞
f (t ) =
∑
−∞
N −1
N −1
∑c n =0
n
e
−i
k
k =0
i
2π nk ∆ t N∆t
1T cn = ∫ f (t ) e − iωnt dt , n ∈ Z , T 0
, t∈ R,
∑ f (t ) e
1 cn = N∆t
f ( tk ) =
cne
i ω nt
=
2π nk ∆ t N∆t
1 ∆t = N
N −1
∑c n =0
n
e
i
2π nk N
N −1
∑
k =0
fk e
−i
2π nk N
, n = 0 ,..., N − 1
, k = 0 ,..., N − 1
DFT fk =
N −1
∑c n =0
n
e
i
4.11. 2010
2π nk N
, k = 0 ,..., N − 1,
1 cn = N MI21
N −1
∑
k =0
fk e
−i
2π nk N
, n = 0 ,..., N − 1 13
Konvoluce posloupností (vyhlazování signálu)
{a }
∞ k k =0
{a }
,
M −1 k k =0
,
{b }
∞ l l =0
{b }
a0 a 1 a2 c = Ab = 0 ⋮ 0 0 4.11. 2010
{cn } = {ak }∗ {bl }, {c }
∞ n n =0
{c }
N −1 l l =0
N + M −1 n n =0
, cn = ∑ an −k bk k =0
, {c n } = {a k }* {bl }
0
0
⋯
0
0
a0
0
⋯
0
0
a1 0
a0 0
⋯ ⋯
0 a M −1
0 ⋯
⋮
⋮
⋱
⋮
⋱
0 0
0 0
⋯ ⋯
⋯ ⋯
a M −1 0
MI21
n
0 0 a0 ⋮ a M −2 a M − 1 0
b0 b 1 b2 ⋮ ⋮ b N −2 b N − 1 14
Laplaceova transformace (řešení DR)
an y
(n )
+ a n −1 y
( n −1 )
+ ⋯ + a1 y ′ + a 0 y = f (t ) ,
( n−1) ′ y( t0) =c0 , y ( t0 ) =c1 , ⋯,y ( t0 ) =cn−1 . ∞
LT { f } = F (s ) = ∫ f (t ).e − st dt , 0
DR
4.11. 2010
L(DR)
Y(s)
MI21
y(t)
15
Specifikace problému a generování sítě.
modelování geometrie
(typ úlohy, materiálové vlastnosti, počáteční a okrajové podmínky, …)
ANSA
Importer STL, IGIS, VRML, …
Exporter CDB format
CAD System Integration •Autodesk Inventor •CATIA •Pro/ENGINEER •Solid Edge •SolidWorks •Unigraphics
4.11. 2010
knihovna se škálovatelnými algoritmy založena na metodách rozložení oblasti FETI a BETI
T. Kozubek, T. Brzobohatý, A. Markopoulos Z. Dostál, V. Vondrák, R. Kučera, M. Sadowská MI21
16
problem
1.
F
Metody rozložení oblasti F
F
FETI
2.
FETI-DP (částečné rozdělení, regulární matice tuhosti)
F
FETIDP
3.
F
F
všechny podoblasti defekt matic K je znám
TFETI
4.11. 2010
podoblasti jsou volné nebo uchycené (předem není znám defekt matic K)
MI21
17
Introduction to Total FETI/BETI problem
1.
BETI
F
F
F
subdomains are fixed or free but with different defects
FETI
2.
FETI-DP (partial splitting, nonsingular)
F
FETI-DP
3.
F
TBETI
4.11. 2010
MI21
všechny podoblasti jsou volné se shodnou dimenzí jádra
18
Total FETI – primární proměnné min u Ku − u f 1 2
u K B c
T
T
funkcionál energie
pole posuvů matice tuhosti matice omezení vektor omezení
BI u ≤ cI BG u = cG BB u = cB
Total FETI
F
Dirichletovy okrajové podmínky jsou také vynuceny pomoci LM, BBu= cB
F
obstacle 4.11. 2010
MI21
19
Vše a mnohem více můžete najít Table of Contents - Preface. Part I. Background 1. Linear Algebra.- 2. Optimization. Part II. Algorithms 3. CG for Unconstrained Minimization 4. Equality Constrained Minimization 5. Bound Constrained Minimization 6. Bound and Equality Constrained Minimization Part III. Applications to Variational Inequalities 7. Solution of a Coercive Variational Inequality by FETI-DP method 8. Solution to a Semicoercive Variational Inequality by TFETI Method.- References.Index.
Seminář výpočetní mechaniky, 1-3. prosince 2010, VŠB-TU & Ústav geoniky AV ČR, v.v.i. 4.11. 2010
MI21
20
Škálovatelnost: kostka nad tuhou překážkou
4.11. 2010
MI21
21
4.11. 2010
MI21
22
23
1. ingredience: paralelní implementace blokově diagonální struktura matice tuhosti => vhodná pro paralelní implementaci
K1 O K = ⋯ O
O ⋮ O 2 K ⋮ O ⋯ ⋱ ⋯ O ⋮ K n
Mohou být řešeny koercivní a semikoercivní úlohy. 4.11. 2010
MI21
23
Paralelní programování Cluster HP model BLc7000 (c-class).
Typické případy • Paralelně smyčky for
Konfigurace výpočetních uzlů: 2x dual core CPU AMD Opteron 2210 HE 8GB ECC DDR2 667MHz RAM
– mnoho cyklů – dlouhé iterační cykly
• Offloading Work • Objemné datové sady IT4Innovations, Centrum excelence, Ostrava, Czech Republic, cluster with more than 30,000 processors. http://www.it4i.eu 4.11. 2010
MI21
24
Matlab Distributed Computing Server Worker Distributed Computing Server
Client Parallel Computing Toolbox
Scheduler or Job manager
• Local job manager • MathWorks job manager • Third-party schedulers (Windows CCS, PBS Pro, LSF) 4.11. 2010
MI21
Worker Distributed Computing Server
Worker Distributed Computing Server
25
Life cycle of a job Scheduler
createJob
Client
Queued
Running
Job
Job
Pending
Job
Job
Job
Job
Job
Job
submit
Job
Job
Job
Job Job
getAllOutputArguments
4.11. 2010
Finished
Worker Worker
Worker Worker Worker
Job
MI21
26
Distributed computations 1. Find scheduler or job manager 2. Create a job 3. Create tasks and associate them with the job 4. Send job to the front 5. Wait until job finishes 6. Gather results 7. Destroy job 4.11. 2010
>>sched=findResource(‘scheduler‘, ’type’,’local’); >>job=createJob(sched); >>createTask(job,@rand,1,{3,3}); >>createTask(job,@eye,1,{4}); >>createTask(job,@ones,1,{{4},{3}});
>>submit(job); >>waitForState(job); >>results=getAllOutputArguments(job) >>destroy(job); MI21
27
Parallel computations 1. 2.
3. 4. 5. 6. 7.
Find scheduler or job manager Create a parallel job and sets the numer of tasks Create parallel tasks Send job to the front Wait until job finishes Gather results Destroy job
>>sched=findResource(‘scheduler‘,… ’type’,’local’); >>pjob=createParallelJob(sched); set(pjob,‘MaximumNumberOfWorkers,4); set(pjob,‘MinimumNumberOfWorkers,4); >>createTask(pjob,@MyParTask,1,{}); >>submit(pjob); >>waitForState(pjob); >>results=getAllOutputArguments(pjob) >>destroy(pjob);
MPT functions: numlabs, labindex, labBarrierlabBroadcast, labProbe, labReceive, labSend, labSendReceive 28 MI21 4.11. 2010
Parallel computation scenario Send FEM models to workers and parallelize assembling of K, K+, R, fv, stresses, searching contact pairs, multiplication procedures. Slave 1 (subdom. 1)
Master
Slave 2 (subdom. 2)
Slave 3 (subdom. 3)
Send FEM Assemble K,K+,R,f Build contact conditions
Gather f Multiply by K+
Iterate
Gather result Compute stress Gather stress
4.11. 2010 X
MI21
29
Důlní průmysl: svěrný spoj důlní výztuže
tlak vyvolaný okolním materiálem R1
překryv
ocelová výztuž
R2
4.11. 2010
MI21
30
Důlní průmysl: svěrný spoj důlní výztuže Coulombovské tření A
f = 0.1
A
F
B
B
µ = 0.3
t2
C
tA = tD = 10 mm
t1
D
107mm
E = 2.1e5 MPa,
C
D
tB = tC = 40 mm
F
počet: binárních proměnných 65562 duálních proměnných
4.11. 2010
rovina symetrie
3112
MI21
31
Důlní průmysl: svěrný spoj důlní výztuže
8 těles 7 volných 4.11. 2010
MI21
32
Důlní průmysl: svěrný spoj důlní výztuže redukované napětí HMH
Posunutí
2D
3D
4.11. 2010
MI21
33
bez korekcí kontaktních směru
3 korekce
σN = 985.4 MPa
σN = 1020.6 MPa
1 NODAL SOLUTION STEP=2 SUB =34 TIME=.5 CONTSTOT (AVG) DMX =.100117 SMN =.012349 SMX =948.043
DEC 9 2008 11:06:26 PLOT NO. 8
MN
MX
.012349 252.821 126.416 379.225 PR03.db FKOP= -0.04 f=0.0 PLANE182
ANSYS (paralel) dynamický výpočet: 27 hodin 4.11. 2010
505.629
758.437 632.033
948.043
ANSYS (paralel) statický výpočet: 1 hodina MI21
MATSOL statická úloha, 3 korekce: sekvenční verze 15 min. paralelní 7 min. 34
Srovnání TBETI a TFETI TFETI
TBETI
1.548 1.611 0.774
Celková posunutí [mm]
0.387 0.017 TFETI
TBETI
1.548 1.611 0.774 0.387 0.017
4.11. 2010
MI21
35
Srovnání TBETI a TFETI 250 podoblastí f = 0.1 Čas výpočtu
2.41 h.
Čas výpočtu
2.54 h.
Celkový čas
2.55 h.
Celkový čas
3.33 h.
násobení matice vektor
2 126
násobení matice vektor
1 882
primární proměnné
1 592 853
duální proměnné
261 553
primární proměnné
713 751
duální proměnné
261 553
max(abs(u))
1.303 mm
max(abs(u))
1.378 mm
max(totalU)
1.548 mm
max(totalU)
1.616 mm
FETI 4.11. 2010
BETI MI21
36
Automobilový průmysl: kuličkové ložisko vnější kroužek klec vnitřní kroužek hřídel zatížení
uchycení
4.11. 2010
MI21
37
3
Válečkové ložisko větrné elektrárny
Statistics počet těles
73
počet podoblastí
700
primární proměnné
2,73 M
duální proměnné
459,8 k
počet iterací
4.11. 2010
4270
MI21
38
Automobilový průmysl: pneumatika 3.9 mil stupňů volnosti, 4000 iterací
celková posunutí 4.11. 2010
MI21
redukované napětí HMH 39
Světlomet – Visteon (Autopal) 12 těles
4.11. 2010
VŠB – TU Ostrava MI21
40
Světlomet – Visteon 64 Oblastí
4.11. 2010
VŠB – TU Ostrava MI21
41
Světlomet – Visteon (Test solution-Total Displacement)
4.11. 2010
VŠB – TU Ostrava MI21
42
Dveře auta Sestava složená z 29 těles
4.11. 2010
VŠB – TU Ostrava MI21
43
Dveře auta
4.11. 2010
VŠB – TU Ostrava MI21
44
Dveře auta 128 Oblastí
4.11. 2010
VŠB – TU Ostrava MI21
45
Dveře auta (Test solution-Total Displacement)
4.11. 2010
VŠB – TU Ostrava MI21
46
Auto – Sestava obsahuje okolo 700 těles
4.11. 2010
VŠB – TU Ostrava MI21
47
3000 subdomains
4.11. 2010
VŠB – TU Ostrava MI21
48
700 bodies
3000 subdomains 4.11. 2010
VŠB – TU Ostrava MI21
49
Rozšíření na dynamické úlohy
4.11. 2010
MI21
50
Děkuji za pozornost
4.11. 2010
MI21
51
