TRIWAQ Cursusboek Module 3

Programmeren aan Parallel WAQUA/TRIWAQ Cursusboek Module 3

Programmeren aan Parallel WAQUA/TRIWAQ Cursusboek Module 3

Version Maintenance Copyright

: : :

1.0, oktober 2002 see www.helpdeskwater.nl/waqua Rijkswaterstaat

Cursus Domein Decompositie en Parallel Rekenen in SIMONA

Cursusboek Module 3

Inhoud PROGRAMMEREN AAN PARALLEL WAQUA/TRIWAQ................................................................I CURSUSBOEK MODULE 3.................................................I 1

INTRODUKTIE PARALLEL PROGRAMMEREN .. 4 1.1 1.2 1.3 1.4 1.5 1.6

ONDERDELEN VAN PARALLEL PROGRAMMEREN ......... 4 PARALLELLE COMPUTER ARCHITECTUREN ................. 5 PARALLELLE PROGRAMMEERPARADIGMAS ................ 7 PROBLEEM DECOMPOSITIE STRATEGIËN ................... 11 CLASSIFICATIE VAN REKENPROBLEMEN................... 11 ANALYSEREN VAN PARALLELLISME EN DATAAFHANKELIJKHEDEN .......................................................... 12 2 GLOBAAL ONTWERP VAN PARALLEL WAQUA/TRIWAQ............................................................. 18 2.1 UITGANGSPUNTEN BIJ DE PARALLELLISATIE ............ 18 2.2 OVERZICHT VAN DE SYSTEEM-STRUCTUUR .............. 20 2.3 UITWERKING VAN DE ROOSTER-OPDELING............... 22 2.4 IMPLEMENTATIE VAN NUMERIEKE BEREKENINGEN (“LOOPS”) IN WAQUA/TRIWAQ ..................................... 24 2.5 PARALLELLISATIE VAN NUMERIEKE BEREKENINGEN (“LOOPS”) IN WAQUA/TRIWAQ ..................................... 26 3

NUMERIEKE ASPECTEN ........................................ 33 3.1 3.2 3.3

AANPASSINGEN AAN SOLVERS I .............................. 33 AANPASSINGEN AAN SOLVERS II ............................. 42 VEREISTE HERSTRUCTURERING VAN DE SEQUENTIËLE PROGRAMMATUUR ............................................................. 48 3.4 PARALLELLISATIE VAN CONVERGENTIE-CRITERIA ... 52 3.5 VOORTGEZETTE DATA-ANALYSE ............................. 55 4 THEORIE EN PRAKTIJK VAN COMMUNICATIE MET COCLIB .................................................................... 64 4.1 4.2 4.3 4.4 4.5

OVERZICHT VAN COCLIB...................................... 64 COCLIB EN COEXEC........................................... 65 INDEXSETS EN STENCILS ......................................... 66 DE COCLIB ROUTINES........................................... 71 GEBRUIK VAN GLOBALE COMMUNICATIES ............... 77

5 UITBREIDINGEN VOOR DOMEIN DECOMPOSITIE ............................................................... 80 5.1 5.2 5.3 6

INTERNE WERKING VAN COEXEC ..................... 94 6.1 6.2

Versie 1.0, oktober 2002

BASISFILOSIFIE ....................................................... 80 INTERPOLATIE IN COCLIB ..................................... 83 INTERFACES BIJ HORIZONTALE VERFIJNING.............. 91 DE FUNCTIES VAN COEXEC .................................. 94 IN- EN UITVOER VAN COEXEC............................... 94

2


6.3 6.4 6.5

Cursusboek Module 3

UITGEVOERDE CHECKS ........................................... 96 HET OPSTARTEN VAN PROCESSEN............................ 97 BOODSCHAPPEN TUSSEN WAQPRO EN COEXEC... 97

7 INTERNE WERKING VAN DE PARTITIONER COPPRE ........................................................................... 100 7.1 7.2 7.3 7.4

DE FUNCTIES VAN COPPRE ................................. 100 INVOER EN UITVOER VAN COPPRE....................... 101 DE WERKING VAN COPPRE.................................. 106 TOEVOEGEN VAN NIEUWE DATA STRUCTUREN ....... 112

8 INTERNE WERKING VAN DE COLLECTOR COPPOS............................................................................ 113 8.1 8.2 8.3

FUNCTIES VAN COPPOS ...................................... 113 HET COLLECTIE ALGORITME ................................. 113 SPECIALE GEVALLEN ............................................ 115

9 CASE-STUDY: PARALLELISATIE VAN HET HORIZONTAAL K- MODEL........................................ 117 REFERENTIES ................................................................ 118


3


1

Cursusboek Module 3

Introduktie Parallel Programmeren Voorwoord: het doel van module 3 is om inzicht te geven in de manier waarop parallel rekenen en domein decompositie in WAQUA/TRIWAQ zijn geïmplementeerd. Na afloop van de cursus zouden de deelnemers de gehanteerde aanpak globaal moeten kunnen plaatsen ten opzichte van mogelijke andere aanpakken, de consequenties van parallel rekenen en domein decompositie kunnen beoordelen voor wijzigingen aan WAQUA/TRIWAQ, en voor een aantal gevallen de parallellisatie van nieuwe functionaliteit zelf kunnen uitvoeren. In dit eerste hoofdstuk geven we een overzicht van wat het in het algemeen inhoudt om programmatuur te ontwikkelen voor parallelle computers. Daarbij ligt de nadruk op de gehanteerde aanpak voor parallellisatie van WAQUA/TRIWAQ, maar wordt deze aanpak ook in een bredere context geplaatst.

1.1

Onderdelen van parallel programmeren Hoewel parallel rekenen op grote schaal succesvol wordt toegepast bestaat er niet een enkele manier voor parallel programmeren. In plaats daarvan zijn er vele verschillende programmeertalen die kunnen worden gebruikt (honderden), ieder met hun eigen sterke en zwakke kanten. De verschillende alternatieven zijn moeilijk met elkaar vergelijkbaar doordat ze zijn gebaseerd op sterk verschillende uitgangspunten, die vaak voortkomen uit verschillende toepassingsgebieden waarop men zich richt, of verschillende soorten parallelle computers die men wil kunnen gebruiken. In deze cursus richten we ons op numerieke simulatieproblemen. In deze problemen staat de te gebruiken oplosmethode vaak centraal: simulatieprogramma’s zijn gericht op het uitvoeren van een bepaalde taak, en die taak wordt beschreven door de stappen die er moeten worden uitgevoerd. Ook binnen deze klasse van problemen blijkt er een veelheid aan methoden te zijn voor parallel rekenen, en blijken methoden vanuit allerlei verschillende invalshoeken te kunnen worden bekeken. Er is dan ook geen enkel recept te geven voor het parallelliseren van een applicatie. Wel kunnen we een vijftal belangrijke aspecten van parallel rekenen identificeren: 1) Bepalen van het parallellisme in de gestelde rekentaak, d.w.z. welke berekeningen in principe gelijktijdig kunnen worden uitgevoerd, en eventueel vergroten van het parallellisme door veranderen van het gebruikte algoritme. 2) Classificeren van het soort parallellisme: afhankelijk van de structuur en regelmaat van de berekeningen hebben bepaalde methoden voor parallellisatie


4


Cursusboek Module 3

duidelijke voordelen of kunnen ze juist direct worden uitgesloten. 3) Kiezen van een mechanisme voor het verdelen van de berekeningen, bijvoorbeeld op basis van de verschillende soorten berekeningen (stroming vs. stoftransport) of op basis van de data-objecten van het probleem. 4) Kiezen van een programmeerparadigma. Dit is min of meer de software architectuur: de programmeertaal, omgeving, en de structuur van het programma of de programma’s. 5) Kiezen van de parallelle computer(s) waarop de berekening moet kunnen worden uitgevoerd. Bij dit laatste punt aangekomen blijkt direct dat de punten niet in de gegeven volgorde kunnen uitgewerkt: de te gebruiken hardware kan een grote invloed hebben op de te gebruiken programmeertaal, probleemopdeling en zelfs de numerieke algoritmes die men wil gebruiken. Daarom geeft bovenstaande geen recept maar een soort checklist van onderdelen die een goede parallellisatiestrategie bevat. De vijf aspecten moeten goed van elkaar worden onderscheiden en voor alle aspecten moeten duidelijke keuzes worden gemaakt. In de volgende paragrafen worden de onderdelen een voor een uitgewerkt, te beginnen met de laatste. In het volgende hoofdstuk wordt dit stramien gebruikt voor het beschrijven van de parallellisatiestrategie die voor WAQUA/TRIWAQ is gebruikt.

1.2

Parallelle computer architecturen In de afgelopen jaren is er een sterke convergentie geweest in de architecturen van commercieel succesvolle parallelle computers. Eigenlijk kan men stellen dat alle parallelle computers tegenwoordig zijn opgebouwd uit een aantal onafhankelijke volledig functionele processoren die met elkaar zijn verbonden door een krachtig interconnectie-netwerk en die zijn voorzien van een grote hoeveelheid geheugen en schijfruimte. Dit in tegenstelling tot begin jaren negentig van de vorige eeuw, toen er een groot verschil was tussen (parallelle) vectorcomputers (Cray YMP, C90), clusters van volledig functionele processoren (Ncube, Intel Paragon, Transputers), en machines waarin er een veel groter aantal zeer eenvoudige processing elements werden gebruikt (Maspar-I, CM-5). De toenmalig gangbare classificatie van computers als SIMD of MIMD (single of multiple instruction, multiple data) is voor parallel rekenen van minder belang geworden, maar is nog wel relevant voor speciale processoren zoals voor grafische bewerkingen. Ook de structuur van de gebruikte interconnectie-


5


Cursusboek Module 3

netwerken (ring, mesh, torus, hypercubus) is tegenwoordig een stuk minder relevant door de sterke verbetering van routingchips, operating systemen (bijv. virtueel shared memory) en communicatie-bibliotheken (PVM, MPI). Toch zijn er nog de nodige verschillen tussen parallelle computers die van belang zijn voor het parallel programmeren. Deze verschillen zitten vooral in de hierarchische opbouw van het geheugen en in de snelheid van het interconnectie-netwerk. Computergeheugens zijn meestal hierarchisch georganiseerd omdat er geheugen van verschillende snelheden wordt toegepast: een klein beetje zeer snel en kostbaar geheugen en een grote hoeveelheid langzaam en goedkoper geheugen. Het snelste geheugen bestaat uit de registers van de processoren zelf. Daarna komen het “cache”-geheugen en het centrale geheugen, en tenslotte worden ook harde schijven als geheugen gebruikt (swappen, out of core rekenen). Cache geheugen wordt gebruikt om gegevens en programma-instructies die mogelijk veel keer achter elkaar nodig zijn “dicht bij” de processor te bewaren. Op basis van de geheugenstructuur kunnen parallelle computers in de volgende drie categoriën worden onderverdeeld: Symmetric multiprocessors (SMP). Dit zijn computers waarin alle processoren even snel bij alle stukken van het centrale geheugen kunnen komen. Bijvoorbeeld dual-Pentium PC’s bevatten een enkele systeembus waar aan de ene kant de twee processoren op zijn aangesloten en aan de andere kant het centrale geheugen ligt. Vergelijkbare principes worden toegepast in multiprocessor werkstations zoals de HP K460, SUN HPC450 en SGI Onyx. NUMA multiprocessors, ofwel Virtually (distributed) shared memory computers (VSM, DSM). Dit zijn computers die zijn opgebouwd uit meerdere processoren die wel allemaal bij alle stukken van het centrale geheugen kunnen komen, maar waarbij de toegangstijden variëren voor verschillende stukken (NUMA=non-uniform memory access). Bijvoorbeeld doordat de machine fysiek is op gebouwd uit processormemory combinaties die met elkaar zijn verbonden via een interconnectie-netwerk, en waar de systeemhardware of software ervoor zorgt dat de geheugens van andere processoren toegankelijk zijn. Typische voorbeelden van deze classe zijn de SGI Origin2000/3800 series (Unite, Teras) en de Cray T3E. Multicomputers – parallelle computers die bestaan uit afzonderlijke processor-memory combinaties die zijn verbonden via een of ander netwerk, en waarbij iedere


6


Cursusboek Module 3

processor alleen zijn eigen geheugen kan benaderen. De processor-memory combinaties kunnen zijn geïntegreerd in een enkele machine, zoals bij de IBM SP2, maar vaker worden parallelle computers van dit type gevormd door het met elkaar verbinden van afzonderlijke sequentiële of parallelle computers. Clusters van Linux PC’s vallen in deze categorie. Deze worden ook wel netwerken van werkstations genoemd (NOW), en als er voornamelijk standaard hardware wordt gebruikt dan wordt ook de benaming Beowulf cluster toegepast. De verschillen tussen deze architecturen hebben enkele belangrijke consequenties voor de manier waarop men ze kan programmeren. Op multicomputers is geen gezamelijk geheugen beschikbaar, dus moet men altijd enige vorm van “message passing” gebruiken (zie volgende paragraaf), terwijl op SMP’s ook een “shared memory” benadering mogelijk is. Op NUMA multiprocessors zijn in theorie beide aanpakken mogelijk. Maar in de praktijk eisen de verschillende geheugensnelheden vaak toch dat de programmeur zelf in de gaten houdt wanneer welke gegevens worden uitgewisseld. Hoewel de communicatie misschien niet hoeft te worden geprogrammeerd heeft dit toch veel weg van de “message passing” benadering. De snelheid van het interconnectie-netwerk van de parallelle computer variëert van supersnel (bus of andere structuur in een SMP-machine) tot erg langzaam (samen gebruikte supercomputers van verschillende rekencentra). In veel gevallen is de communicatietijd voor een boodschap goed te benaderen met een opstarttijd (latency) en een bandbreedte (1/tijd per byte). Dit onafhankelijk van de werkelijke afstand tussen processoren in het netwerk. In een mesh-structuur zou je verwachten dat het aantal “hops” zou meewegen in de communicatietijd, maar tegenwoordig blijken dit soort effecten van ondergeschikt belang. Gangbare getallen voor latency en bandbreedte zijn 800 s en 5MB/sec voor een Linux cluster met fast-Ethernet (100Mbit/sec), en 15 s, 150MB/sec voor de TERAS machine van SARA.

1.3

Parallelle programmeerparadigmas Een programmeerparadigma is een stel van regels dat men kiest voor het structureren van een programmeertaak. Bijvoorbeeld de programmeertaal die men kiest en de conventies die men daarbij afspreekt voor het vormgeven van modules. Ook de structurering van de programmatuur hoort erbij: het executiemodel, de onafhankelijke eenheden van het programma (processen, threads, objecten). Veel van deze aspecten horen ook bij de gebruikte software architectuur.


7


Cursusboek Module 3

Het is belangrijk om te beseffen dat er verschillende programmeerparadigmas mogelijk zijn en om het gebruikte paradigma goed te beschrijven. De keuze voor een paradigma heeft grote invloed op verschillende kwaliteitsaspecten van de uiteindelijke programmatuur (onderhoudbaarheid, uitbreidbaarheid, performance). De beschrijving van het paradigma geeft veel houvast bij de implementatie van de programmatuur. Het eerste onderscheid tussen verschillende programmeerparadigmas dat van belang is voor parallel rekenen is het onderscheid tussen impliciet en expliciet parallel programmeren. Met impliciet parallel programmeren wordt bedoeld dat de programmeur zich op een vrij abstract niveau bezig houdt met het programma zonder concreet te denken in termen van de rekenprocessen die het uiteindelijke werk gaan doen. Een eerste vorm hiervan is dat de programmeur een sequentiëel programma maakt en dat de parallellisatie automatisch gebeurt door een parallelliserende compiler. Dit gaat moeizaam als de programmeur geen rekening heeft gehouden met parallellisatie (“legacy” Fortran77-code) en levert dan vaak slechts een beperkte versnelling op. Het gaat beter als de programmeur het parallellisme in zijn programma meer zichtbaar maakt voor de compiler. Dit kan door toevoegen van compiler-directives, of bijvoorbeeld door het gebruik van “data parallelle” statements: array-syntax, “forall”-loops, etc. Dit is in het bijzonder sterk uitgewerkt in de programmeertaal High Performance Fortran (HPF), welke tegenwoordig is opgenomen in Fortran95. De data parallelle benadering is vooral voor regelmatige problemen goed geschikt. Een andere vorm van impliciet parallel programmeren is bijvoorbeeld dat de programmeur de code schrijft voor een heleboel objecten die met elkaar interacteren zonder erbij te bepalen hoe deze over de beschikbare processoren moeten worden verdeeld. Nog een andere vorm is de benadering waarin de programmeur programma-code schrijft voor aparte taken en hun volgorde-afhankelijkheden geeft (taakgraaf), maar het daadwerkelijk uitvoeren van de taken overlaat aan het systeem. Expliciet parallel programmeren is dat de programmeur zich sterk bezig houdt met de daadwerkelijke rekenprocessen die uiteindelijk worden uitgevoerd, inclusief hun onderlinge interacties. Aangezien het managen van meerdere processen extra complicaties met zich meebrengt heeft dit natuurlijk niet de voorkeur. Maar in de praktijk blijkt de extra complexiteit redelijk te kunnen worden beheerst, en blijken impliciete benaderingen veelal niet afdoende voor het bereiken van tevredenstellende performance.


8


Cursusboek Module 3

Het belangrijkste verschil tussen verschillende methoden voor expliciet parallel programmeren is hoe de verschillende processen met elkaar interacteren. De twee extremen hierbij zijn als volgt: Shared data paradigma – alle rekenprocessen werken op een enkele globale data-ruimte. Alle processen kunnen in principe alle variabelen uit die data-ruimte opvragen zonder tussenkomst van andere processen. Binnen dit paradigma speelt synchronisatie van de berekeningen van verschillende processen een belangrijke rol. Process/channel paradigma – ieder proces werkt op een eigen lokale data-ruimte die niet direct toegankelijk is voor de andere processen. De uitwisseling van gegevens tussen de processen is expliciet zichtbaar in hun programmas. Hiervoor worden “channels” gebruikt: deze zorgen voor ontkoppeling van welke data er gecommuniceerd wordt en met wie er gecommuniceerd wordt. Deze twee methoden kunnen ook “shared memory programming” en “message passing programming” worden genoemd. Maar dat is enigszins misleidend omdat het process/channel paradigma net zo goed op een shared memory computer kan worden toegepast, terwijl het shared data model ook op een (distributed memory) multicomputer kan worden ondersteund. Verder kleven aan de termen shared memory en message passing connotaties als “low level” en moeilijk, terwijl dit bij gebruik van een geschikte software-omgeving allerminst het geval hoeft te zijn. De voor grootschalige rekenproblemen meest relevante programmeeromgevingen die bij bovenstaande paradigmas horen zijn respectievelijk OpenMP en MPI. Dit zijn beiden standaarden die door een breed forum zijn ontwikkeld en worden ondersteund. OpenMP is een set van compiler-directives en library-functies waarmee een sequentiëel Fortran of C-programma kan worden geparallelliseerd binnen het shared data paradigma. Het executiemodel dat hierbij wordt gebruikt is dat er in eerste instantie een enkele “thread” wordt opgestart voor het programma, maar dat er voor verschillende secties van de code meerdere threads kunnen worden gebruikt. De programmeur geeft aan welke secties dat zijn, en hoe de threads aan het begin en einde van deze secties moeten synchroniseren. De variabelen van het programma zijn in principe globaal beschikbaar voor alle verschillende threads, maar kunnen door de programmeur als thread-local worden gedeclareerd. Communicatie gebeurt impliciet via de variabelen die de verschillende processen met elkaar delen. OpenMP is primair gericht op de “single program multiple data” (SPMD) manier


9


Cursusboek Module 3

van programmeren: een enkel programma dat meerdere keren wordt opgestart voor verschillende sets van invoer-data. MPI is een bibliotheek van communicatiefuncties die in Fortran en C-programmas kunnen worden gebruikt. Hier worden er meerdere instanties van een enkele executable tegelijk opgestart en deze gaan dan al communicerend samen een probleem oplossen. Maar ook kunnen er onderweg nog nieuwe processen worden opgestart, men is dus niet beperkt tot het SPMD-model. Er is een groot aantal communicatie-routines waarmee een proces data kan versturen naar en ontvangen van andere processen, evenals routines waarmee kan worden gecommuniceerd met groepen van andere processen. De synchronisatie van de processen gebeurt impliciet op de momenten dat ze wachten op boodschappen van de anderen. Het grootste verschil tussen de shared data en process/channel paradigmas zit in de manier waarop data is gekoppeld aan processen. In het process/channel paradigma is er een sterke koppeling, zijn de gegevens op ieder moment het exclusieve bezit van een enkel proces. Hierdoor zijn de gegevens in de geheugenhierarchie steeds “dicht bij” de processor die ze nodig heeft, is er een goede “data localiteit”. Dit is belangrijk voor het bereiken van goede performance. Verder zijn zulke programmas gemakkelijker te analyseren en debuggen, interacties tussen verschillende processen kunnen niet gemakkelijk over het hoofd worden gezien (de programmeur heeft niet te maken met “anti-dependencies” en “race conditions”, zie paragraaf 1.6). In het shared data paradigma is de koppeling tussen data en processen minder strak. De programmeur hoeft daardoor niet precies te vertellen welke gegevens er moeten worden gecommuniceerd, waardoor dit model in het algemeen als gemakkelijker te gebruiken wordt ervaren. Verder is dit voordelig als de werkverdeling dynamisch moet worden aangepast of als het probleem minder regelmatig is. Dit echter ten koste van een stuk controle over de performance, en ten koste van de overdraagbaarheid naar computers zonder een (virtueel) shared memory systeem. Tenslotte noemen we paar bekende structuren van samenwerkende processen die ook als onderdeel van het paradigma worden beschouwd. Dit zijn de “farmer-worker” organisatie, ook wel “master-slave” en “process-farm” genoemd, waarbij een coördinerend proces steeds taken uitdeelt aan de werker-processen, de “pipeline” organisatie, waarbij ieder proces zijn eigen bewerking uitvoert op een stel gegevens en deze dan doorgeeft aan het volgende proces in de rij, en de “agenda-parallelle” manier van werken, waarbij alle rekenprocessen dezelfde agenda van activiteiten uitvoeren en daarin steeds hun eigen gedeelte uitvoeren. Dit laatste is min of meer gelijk aan de SPMD benadering.


10


1.4

Cursusboek Module 3

Probleem decompositie strategiën Het volgende aspect van parallel programmeren dat we bespreken is de strategie die wordt gevolgd voor het verdelen van het rekenwerk in afzonderlijke taken. Drie generieke mechanismes die hiervoor kunnen worden gebruikt zijn het verdelen van de data van het probleem, het verdelen van de operaties van het probleem, of het verdelen van taken, operaties plus data, van het probleem. Het verdelen van de data van een probleem is een aantrekkelijke methode omdat er in het algemeen bij grote rekenproblemen ook grote aantallen verschillende data-items worden gebruikt, en omdat er vaak sprake is van een enkele “agenda” van “super-stappen”, berekeningen voor een groot aantal data-items tegelijkertijd. Het principe van datadecompositie wordt ook wel data parallellisme genoemd, of geometrische decompositie (als de data een ruimtelijke betekenis hebben), rooster-opdeling (als niet de data-items zelf maar een onderliggend abstract rooster wordt verdeeld) of domein decompositie. Relevante aspecten van een datadecompositie zijn of deze statisch is of dynamisch moet variëren, en of alle delen gelijk of onderling verschillend zijn. Het verdelen van de operaties van een probleem wordt veelal functioneel parallellisme genoemd. Dit heeft vooral toepassingen in pipelines – produktielijnen waar steeds dezelfde serie operaties moet worden uitgevoerd op verschillende hoeveelheden data. Het verdelen van alle berekeningen in een stel taken met volgorde-restricties (taakgraaf) is praktisch omdat het nauw aansluit bij het analyseren van het parallellisme in een berekening (zie paragraaf 1.6). Met name als de hoeveelheid werk per taak moeilijk van te voren is in te schatten. Denk bijvoorbeeld aan een directe solver voor een in blokken gepartitioneerde ijle matrix: het is gemakkelijk om taken aan te wijzen (vegen blokken) en hun volgorde-relaties te specificeren, en het is een stuk moeilijker om een geschikte data-decompositie van het probleem te maken (verdeling van de matrix of van de onbekenden) die het rekenwerk gelijkmatig over de processoren verdeelt.

1.5

Classificatie van rekenproblemen De keuze tussen verschillende decompositie-strategiën en programmeerparadigmas wordt sterk bepaald door een paar algemene karakteristieken van het probleem dat moet worden opgelost. Volgens (Fox, 1994) zijn hierbij vooral de ruimtelijke en temporele structuur van het algoritme van belang en moeten er vijf categoriën worden onderscheiden:


11


Cursusboek Module 3

Natuurlijk parallelle algoritmes (embarrassingly parallel) – algoritmes die bestaan uit een groot aantal onafhankelijke berekeningen. Bijvoorbeeld het evalueren van zetten in een schaakprogramma. Synchrone en los-synchrone algoritmes – dit zijn algoritmes waarin steeds gelijktijdige updates worden gedaan aan een groot aantal data-items. Het onderscheid tussen synchroon en los-synchroon is of de updates van alle data-items identiek zijn of niet, en bepaalt of een algoritme op een SIMD parallelle computer kan worden uitgevoerd. Dit is is tegenwoordig minder belangrijk. Los-synchrone problemen kunnen goed met een data parallelle benadering worden aangepakt. Veel tijdsintegratie-procedures en iteratieve solvers vallen in deze klasse. A-synchrone algoritmes – dit zijn algoritmes met verschillende data-items (verschillende updateberekeningen) en met in de tijd variërende interactiepatronen. Voorbeelden hiervan zijn discrete-event simulaties met een groot aantal verschillende actoren en parallelle directe solvers voor ijle matrices. Samengestelde problemen – dit zijn algoritmes die zijn opgebouwd als samenstelling van meerdere verschillende algoritmes die ieder voor zich in een van de bovenstaande categoriën kunnen vallen. Denk bijvoorbeeld aan een gekoppelde simulatie van waterbeweging en sedimenttransport.

1.6

Analyseren van parallellisme en dataafhankelijkheden Nu dat we weten waar we naartoe willen, wat er mogelijk is voor het implementeren van parallelle algoritmes, gaan we tenslotte kijken naar het parallelliseren van een algoritme zelf. In de praktijk zal dit juist een van de eerste stappen zijn, is de toepassing belangrijker dan de computer waarmee wordt gewerkt, maar bovenstaande uitleg is een nuttige achtergrond voor het huidige onderwerp. Het probleem van het parallelliseren van een algoritme is dat je niet zomaar opeenvolgende berekeningen tegelijkertijd mag uitvoeren: als de volgende berekening het resultaat van een vorige gebruikt dan is er een data-afhankelijkheid die ook bij parallelle verwerking gerespecteerd moet worden, zoals bijvoorbeeld in het volgende code fragment:

100


nfac = 1 do 100 i=2,n nfac = nfac * i continue

12


Cursusboek Module 3

De data-afhankelijkheden van een algoritme kunnen worden weergegeven in een gerichte graaf, de dataafhankelijkheidsgraaf. De knopen in de graaf zijn de berekeningen van het algoritme, de pijlen geven een informatie-stroom weer: van de berekening die een resultaat produceert naar de berekening die dit resultaat gebruikt. De data-afhankelijkheidsgraaf geeft al het parallellisme weer dat er bestaat bij het gekozen niveau van de te onderscheiden berekeningen. Die berekeningen kunnen zo klein zijn als een enkele optelling, maar kunnen ook wat groter gekozen worden zoals bijvoorbeeld de sommatie van de elementen van een vector, de berekening van een inprodukt e.d. Wanneer de berekeningen verder worden samengevoegd tot grotere eenheden dan wordt de graaf ook wel een taakgraaf genoemd, zie paragraaf 1.4. Twee berekeningen in een dataafhankelijkheidsgraaf of twee taken in een taakgraaf kunnen tegelijkertijd worden uitgevoerd als er geen pad bestaat in de graaf van een van beiden naar de andere. We bekijken als voorbeeld verschillende iteratieve oplosmethoden voor voor de Poisson-vergelijking u xx u yy f op een vierkant domein [0,1] [0,1] met (Dirichlet) randvoorwaarde u 0 . Bij gebruik van centrale differenties en een constante roosterafstand h 1 (n 2) levert dit het bekende stelsel vergelijkingen op: ui , j

1

ui

1, j

4ui , j

ui

1, j

ui , j

1

h 2 fi , j

(1.1)

De resulterende matrix is ijl, heeft vijf niet-nul elementen per rij en heeft bandbreedte 2n 1 . De matrix is positief definiet (alle eigenwaarden groter dan nul) en net aan diagonaal dominant: eigenschappen die van belang zijn voor de convergentie-snelheid van iteratieve solvers. De eerste solver die we bekijken is de Jacobi-methode. In iedere iteratie q 1, 2, wordt de volgende berekening uitgevoerd: for j 1 to n

(1.2)

for i 1 to n uiq, j

( h 2 fi , j

uiq, j 11 uiq 1,1 j

uiq 1,1 j

uiq, j 11 ) / 4

Het is duidelijk dat alle berekeningen binnen een iteratie volledig onafhankelijk van elkaar kunnen worden uitgevoerd. Maar ook kunnen berekeningen van verschillende iteraties tot op zekere hoogte tegelijkertijd worden uitgevoerd. Bijvoorbeeld kan worden nagegaan dat de berekening van uiq, j niet afhankelijk is, ook niet indirect, van de berekening van uiq s, j , voor s 1.


13


Cursusboek Module 3

Met de data-afhankelijkheidsgraaf kunnen in principe allerlei analyses worden uitgevoerd. Bijvoorbeeld kan de lengte van het langste pad worden bepaald. Hiervoor wordt aan iedere berekening bepaalde kosten toegevoegd, en de kosten van het langste pad geven dan de kleinst mogelijke rekentijd bij parallelle verwerking. Daarbij wordt er overigens wel vanuit gegaan dat er een onbeperkt aantal processoren beschikbaar is en dat communicatie oneindig snel verloopt. Verder kunnen er verschillende toewijzingen van de berekeningen aan abstracte processoren worden gemaakt en op dezelfde manier worden beoordeeld. Op verschillende manieren kan hiermee invulling worden gegeven aan de hoeveelheid parallellisme in het algoritme (maximale breedte, gemiddelde breedte, langste pad t.o.v. totale hoeveelheid werk, e.d.). Een relevant getal voor de Jacobi-methode is de breedte: er kunnen steeds n n processoren worden benut. Een tweede iteratie-methode die kan worden gebruikt voor het Poisson-probleem is de Gauss-Seidel methode: for j 1 to n

(1.3)

for i 1 to n u

q i, j

2

( h fi , j

q i, j 1

u

q i 1, j

u

q 1 i 1, j

u

q 1 i, j 1

u

)/4

Deze versie van Gauss-Seidel gebruikt de “lexicografische” nummering waarbij, in een sequentiële implementatie, van links naar rechts alle punten binnen een rij, en daaromheen van onder naar boven alle rijen een voor een worden afgelopen. Steeds worden zo veel mogelijk waardes op het nieuwe iteratieniveau q gebruikt, n.l. voor de punten links en onder het punt waarvoor de berekening wordt uitgevoerd. Als er alleen binnen een iteratie wordt geparallelliseerd, bijvoorbeeld omdat er na iedere iteratie op basis van alle nieuwe waardes wordt bekeken of er wel of geen volgende iteratie nodig is, dan is de lengte van het langste pad in de data-afhankelijkheidsgraaf 2n 1 stappen, en zijn er n processoren nodig om deze minimale executietijd te bereiken. Een mogelijke mapping en scheduling van de berekeningen waarmee dit wordt bereikt is rijgewijze verdeling van de roosterpunten over de processoren, en gelijktijdige berekening van de waardes uiq, j voor i j s (diagonalen van het rekenrooster). Een belangrijke eigenschap van een opdeling van de graaf in stukken is de granulariteit (omvang) van de stukken. Iedere afhankelijkheid tussen berekeningen van verschillende processoren vereist communicatie en/of synchronisatie, en hieraan is steeds een systeemafhankelijke opstarttijd verbonden. Een fijnkorrelige opdeling van een algoritme kan hierdoor alleen op een computer met een kleine opstarttijd efficient worden uitgevoerd, een grofkorrelige opdeling is minder kritisch ten aanzien van de gebruikte computer. Deze


14


Cursusboek Module 3

granulariteit verschilt sterk tussen de Jacobi en Gauss-Seidel methodes: bij gebruik van n processoren kunnen deze bij de Jacobi-methode steeds n roosterpunten achter elkaar updaten, in de Gauss-Seidel methode moet er (bij n processoren) na iedere update-stap afzonderlijk worden gecommuniceerd. Bovenstaande laat zien dat de data-afhankelijkheidsgraaf een nuttig hulpmiddel kan zijn bij het bepalen van het potentieel parallellisme in een algoritme en bij het zoeken van geschikte manieren om dit parallellisme kan worden uitgebuit. Het nut van deze methode moet overigens niet worden overschat: in de praktijk blijkt er vooral op relatief kleine aantallen processoren te worden gewerkt, en de communicatietijd wordt in de analyse grotendeels genegeerd. Verder is er hierboven geen aandacht besteed aan de numerieke efficiëntie van de verschillende algoritmen, terwijl de verschillen hierin enorm groot kunnen zijn. Zo zijn er met de Gauss-Seidel methode voor het Poissonprobleem half zoveel iteraties nodig als met de Jacobi-methode, en kan er met over-relaxatie in de Gauss-Seidel methode nog een enorme winst worden behaald. Dit is snel belangrijker de verschillen in speedup, de efficiëntie van een parallelle implementatie, zeker wanneer er niet al te veel processoren worden gebruikt. In de praktijk van parallel WAQUA/TRIWAQ wordt dataafhankelijkheidsanalyse vooral op de volgende manier gebruikt. Aan de hand van deze analyse kan worden bepaald wat er moet worden gecommuniceerd tussen de verschillende processoren. Dit wordt geïllustreerd aan de hand van de redblack Gauss-Seidel methode voor het Poisson-probleem.

for j 1 to n for i 2 to n step 2 (odd(i q i, j

u

2

q 1 i, j 1

( h fi , j u

u

for j 1 to n for i 1 to n step 2 (even(i uiq, j

( h 2 fi , j uiq, j

1

j )) q 1 i 1, j

uiq 1,1 j

uiq, j 11 ) / 4

(1.4)

j ))

uiq 1, j uiq 1, j uiq, j 1 ) / 4

Hoe deze methode is afgeleid (aangepaste nummering van roosterpunten volgens schaakbord-patroon, matrix-splitting L D vs. U) of hoe de convergentie zich verhoudt tot die van de Gauss-Seidel methode met lexicografische nummering (asymptotisch hetzelfde als de matrix symmetrisch positief definiet is) is voor de volgende analyse niet van belang. Belangrijk is om te herkennen dat alle berekeningen in de eerste loop onafhankelijk van elkaar zijn, evenals alle berekeningen in de tweede loop. De tweede loop gebruikt echter wel de resultaten van de eerste loop, en de eerste loop gebruikt resultaten van de tweede loop uit de voorgaande iteratie (q-1).


15


Cursusboek Module 3

Een geschikte manier om bovenstaande berekening te parallelliseren binnen het process/channel-paradigma (expliciet communicatie aangeven) is dus om een statische roosteropdeling te maken in het gewenste aantal deelroosters, om per processor de voor het eigen deelrooster benodigde fi , j en ui0, j te bepalen, en om dan voorafgaand aan iedere loop een communicatie-stap in te lassen om de benodigde waardes van buurpunten te verkrijgen. Welke waardes zijn dit? Het zijn steeds de nieuwste iteranden voor de roosterpunten buiten het eigen deelrooster. Welke roosterpunten buiten het eigen deelrooster? Dit is analytisch uit te schrijven indien we rechthoekige deelroosters veronderstellen. Maar dat is een zware beperking voor WAQUA/TRIWAQ, omdat er veel landpunten zijn waarvoor geen berekeningen hoeven te worden gedaan. Daarom gebruiken we de volgende techniek. Initialiseer een array met 1-en in de roosterpunten van het eigen deelrooster, en met 0-en daarbuiten. Initialiseer een tweede array op 0. Voer een loop uit over alle (i, j ) van het eigen deelrooster waarvoor i j oneven is, en markeer daarvoor steeds de vier naburige punten met 1 in het tweede array. De punten waarvoor er moet worden gecommuniceerd voorafgaand aan de eerste van de twee loops van de red-black Gauss-Seidel methode zijn nu die punten waarvoor het tweede array een 1 bevat en het eerste array een 0 bevat: nodig in de berekening èn niet van het eigen deelrooster. Dit is het principe van data-analyse met behulp van index-sets en stencils, wat verder wordt uitgewerkt in hoofdstuk 3 van deze syllabus. Tenslotte wordt opgemerkt dat afhankelijkheidsanalyse ook een belangrijk hulpmiddel is voor het doorgronden van stukken programma-code. Een verschil ten opzichte van het analyseren van algoritmes is dat in een algoritme een waarde xiq steeds uniek is gedefiniëerd, terwijl variabelen in een programma steeds nieuwe waardes kunnen krijgen. Hierdoor zijn er meerdere soorten afhankelijkheden van belang: Data-afhankelijkheid: een berekening gebruikt het resultaat van een vorige berekening. Eerstgenoemde berekening mag de waarde pas lezen als de ander klaar is met schrijven. Anti-afhankelijkheid: een berekening overschrijft de waarde van een variabele die nog nodig is voor een andere berekening. De overschrijvende berekening mag pas worden uitgevoerd als de andere berekening klaar is met lezen. Output-afhankelijkheid: twee berekeningen wijzigen dezelfde variabele. De berekeningen moeten in de goede volgorde worden uitgevoerd om hetzelfde eindresultaat op te leveren. Versie 1.0, oktober 2002

16


Cursusboek Module 3

In een process/channel programmeermodel moet de programmeur alleen de eerste soort afhankelijkheden in het programma detecteren. In een shared data programmeermodel zijn ook de andere afhankelijkheden van belang. Indien deze niet goed worden onderkend dan kunnen er race-condities ontstaan: het resultaat van het programma is afhankelijk van de precieze timing en snelheid van de verschillende rekenprocessen.


17


2

Cursusboek Module 3

Globaal ontwerp van parallel WAQUA/TRIWAQ In het vorige hoofdstuk is een introductie gegeven in de verschillende aspecten die er komen kijken bij parallel programmeren, en zijn er per onderdeel verschillende benaderingen gepresenteerd. Nu beschrijven we hoe die theorie is toegepast voor parallel WAQUA/TRIWAQ. Enerzijds geeft dit een globaal beeld van hoe parallel WAQUA/TRIWAQ in elkaar zit, anderzijds is het een voorbeeld van hoe men te werk kan gaan bij parallellisatie van een willekeurig ander probleem.

2.1

Uitgangspunten bij de parallellisatie Het globaal ontwerp van de huidige versie van parallel WAQUA/TRIWAQ is grotendeels gemaakt in 1995. Er was toen ervaring met een prototype parallelle versie van de eerste versie van TRIWAQ met vaste lagen, en de eerste versie van WAQUA/TRIWAQ in SIMONA was zojuist gereed gekomen. Er was nog weinig zicht op de huidige convergentie van parallelle computer-architecturen, en daarom was overdraagbaarheid (portabiliteit) een belangrijk doel. Netwerken van werkstations waren sterk in opmars, en zouden als het eventueel mogelijk was moeten kunnen worden gebruikt. De eerste optie die is onderzocht en is afgeschoten is die van automatische parallellisatie met behulp van een parallelliserende compiler. Automatische parallellisatie lijkt niet principiëel onmogelijk, maar er is wel een echte SMP machine nodig (geen NUMA multiprocessor) om een goede speedup te behalen op meer dan een handvol processoren. Bovendien is er een heel goede compiler nodig om het parallellisme te herkennen en het werk slim te verdelen of moet de code fors worden geherstructureerd. Herstructurering was sterk ongewenst vanwege de omvang van de programmatuur en daaraan verbonden onderhoudskosten, ook de gevectoriseerde versie van TRIWAQ is zeer moeilijk up-to-date te houden geweest met de standaardversie. De twee grootste problemen voor automatische parallellisatie zijn de “double sweeps” voor de tridiagonale stelsels voor waterstanden (weinig parallellisme, steeds verschillende verdelingen van roosterpunten over de processoren), en alle verschillende soorten arrays en loops voor verschillende soorten punten uit het grid (randpunten, barriers, source-punten e.d.). De handmatige parallellisatie van WAQUA/TRIWAQ wordt beschreven aan de hand van de vijf onderwerpen die in paragraaf 1.1 zijn onderscheiden.


18


2.1.1

Cursusboek Module 3

Bepalen van het parallellisme Als eerste is natuurlijk gekeken naar de numerieke methoden die in de nieuwe versie van TRIWAQ in SIMONA werden gebruikt. Deze leken sterk op die van de vaste-lagen versie van TRIWAQ, behalve dat het aantal lagen per roosterpunt in het horizontale vlak constant is, en dat de voor parallellisatie erg vervelende “CUE”-methode van de eerste versie van TRIWAQ in de nieuwe versie was vervangen door de “red-black Jacobi” methode.

2.1.2

Classificeren soort parallellisme De gebruikte algoritmen kunnen worden geclassificeerd als lossynchroon. Veel berekeningen in TRIWAQ zijn van de vorm “voor alle roosterpunten doe...”. Eerst worden voor alle roosterpunten de coëfficiënten van de op te lossen stelsels vergelijkingen bepaald, daarna wordt voor alle roosterpunten een iteratie uitgevoerd, enz. De berekeningen verschillen tussen verschillende roosterpunten, bijvoorbeeld in de gebruikte differentie-benadering voor de advectieve termen, of t.g.v. aanpassingen aan de discretisaties nabij open/gesloten randen. Een aantal “kernels” kan ook als a-synchroon worden getypeerd: de double-sweep solver voor tri-diagonale stelsels, en de Gauss-Seidel gebaseerde solvers voor de turbulente grootheden en voor de impulsvergelijking in WAQUA. De oorspronkelijke implementatie van WAQUA is ook meer asynchroon dan de berekeningen van TRIWAQ. De verschillende stadia van de berekeningen waren sterk met elkaar vervlochten (bijv. opstellen/oplossen vergelijkingen). Om parallellisatie met de hieronder beschreven aanpak mogelijk te maken is de implementatie geherstructureerd, zodanig dat deze grote overeenkomsten heeft met die van TRIWAQ.

2.1.3

Parallelle computer-architecturen Zoals hierboven al is aangegeven was portabiliteit een belangrijke doelstelling voor de parallellisatie. Verder waren er clusters van werkstations beschikbaar bij Rijkswaterstaat en bij de TU Delft, en werden deze als een relatief goedkope en krachtige rekenfaciliteit gezien. Daarom werd ervoor gekozen om de parallelle versie geschikt te maken voor multi-computers (zonder shared memory), tenzij dat echt door performance of complexiteit onmogelijk zou zijn.

2.1.4

Parallel programmeerparadigma Vanwege de wens om netwerken van werkstations te kunnen gebruiken en vanwege de onmogelijkheid van automatische parallellisatie lag de keuze voor de een of andere variant van


19


Cursusboek Module 3

“process/channel” paradigma voor de hand. Ieder rekenproces krijgt zijn eigen subdomein toegewezen, en doet daar alle berekeningen voor. De processen communiceren met elkaar ten aanzien van de waardes van roosterfuncties nabij subdomeinranden. Daarbij wordt niet gezegd of die communicatie via “message passing” gebeurt of via bufferarrays in een shared memory. Een ander aspect van het programmeerparadigma betreft het opstarten van de rekenprocessen en verdelen van de data van het probleem. Hierbij is ervoor gekozen om ieder rekenproces zijn eigen versie van de SIMONA-omgeving te geven (geheugenadministratie). Dit om te voorkomen dat er een compleet gedistribueerde implementatie van die omgeving moet worden gemaakt voor bepaalde parallelle computers. Daarom is de gevolgde aanpak goed te karakteriseren als “het koppelen van verschillende instanties van het oorspronkelijke programma”.

2.1.5

Probleemdecompositie Binnen de beschreven aanpak zijn er nog verschillende mechanismen mogelijk voor het verdelen van het rekenwerk. Er is gekozen voor een aanpak waarbij het rekenrooster wordt verdeeld in deelroosters, en waarbij een rekenproces volledig verantwoordelijk is voor de berekeningen van zijn deelrooster. Dit omdat er weinig verschillende functies in het programma zijn die onafhankelijk van elkaar, gelijktijdig kunnen worden uitgevoerd. De roosteropdeling wordt alleen bepaald voor het horizontale vlak omdat het aantal lagen in TRIWAQ-simulaties relatief klein is en omdat een aantal berekeningen in vertikale richting sterk is gekoppeld. De opdeling wordt eenmalig bepaald en is dan vast gedurende de hele berekening (statische i.p.v. dynamische roosteropdeling). Het dynamisch kunnen variëren van de roosteropdeling zou voordelig kunnen zijn in situaties waarin veel droogvallen en onderlopen optreedt (bijv. hoogwatergolven in rivieren), maar is moeilijk te implementeren binnen de gekozen strategie van het koppelen van instanties van het oorspronkelijke programma. Twee manieren om ook bij veel droogvallen en onderlopen een redelijke load-balance te handhaven zijn om bij het verdelen rekening te houden met de bodemligging, en om meer subdomeinen/rekenprocessen te gebruiken dan er processoren zijn.

2.2

Overzicht van de systeem-structuur Het ontwerp dat gaandeweg de parallellisatie tot stand is gekomen wordt gepresenteerd in Figuur 1.


20


Cursusboek Module 3

Title: opratovw4.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 1 – Globale WAQUA/TRIWAQ

systeem-structuur

van

parallel

Een aantal belangrijke aspecten van dit ontwerp zijn als volgt. 1) Het externe gedrag van de parallelle versie is hetzelfde als van de sequentiële versie, in die zin dat een enkele globale invoer (SDS-)file wordt geconverteerd in een enkele globale uitvoerfile. Hierdoor kan parallel rekenen worden gebruikt in combinatie met alle standaard preen postprocessing programma’s, en zijn de details van de parallellisatie grotendeels voor de eindgebruiker. Twee aspecten waarvoor dit niet geldt zijn de message-output en de report-file (per subdomein apart). 2) De partitioner is geïntroduceerd om het rekenrooster te verdelen in min of meer evengrote stukken, en om de gegevens van de initiële globale SDS-file te distribueren over de deeldomeinen. Hij zorgt ervoor dat ieder rekenproces een SDS-file krijgt toegewijzen die vrijwel volledig aan het format van sequentiëel WAQUA/TRIWAQ voldoet. De collector doet het tegenovergestelde, voegt de uitvoerdata per deeldomein samen tot een globale uitvoerfile. Het distribueren en collecteren van de gegevens is strikt gescheiden van het rekenen door de rekenprocessen. Dit zorgt ervoor dat er binnen een rekenproces geen verwarring kan ontstaan


21


Cursusboek Module 3

tussen de globale en lokale dimensies van allerlei soorten arrays. 3) Het master-proces (executive) heeft slechts een zeer beperkte functie. Hij zorgt voor het opstarten van de rekenprocessen, voor het verzamelen van de diagnostische uitvoer in een enkele file, en voor het signaleren en afhandelen van problemen met de rekenprocessen. 4) De WAQUA/TRIWAQ rekenprocessen zijn ieder apart instantiaties van het volledige programma, inclusief file I/O, initialisaties e.d. Het gezichtspunt dat hierbij wordt gebruikt is dat van het koppelen van verschillende simulaties, in plaats van dat een enkele simulatie wordt opgesplitst. Ten behoeve van deze koppeling is het oorspronkelijke programma uitgebreid, onder andere met een nieuw type randvoorwaarde en met communicatie van waardes bij subdomeinranden. Verder wordt de oorspronkelijke programma-code volledig hergebruikt. Het gezichtspunt van gekoppelde zelfstandige simulaties wordt in praktijk gebruikt voor het realiseren van domein decompositie en voor on-line koppelingen. Voor zulke koppelingen is het nuttig als de rekenprocessen zo min mogelijk aannames maken over de omgeving waarbinnen ze draaien. Mede hierom hebben de rekenprocessen in een parallelle run slechts beperkte informatie over de overall-structuur van de berekening. 5) De communicatie tussen de verschillende rekenprocessen is strikt gescheiden van de berekeningen in het programma. Aan de ene kant is dit gedaan vanwege portabiliteit, opdat het precieze communicatiesysteem van de parallelle computer kan veranderen zonder grote gevolgen voor het parallelle programma. Een belangrijkere reden voor het introduceren van een aparte communicatie-bibliotheek is echter het verstoppen van details voor applicatie-programmeurs. Bijvoorbeeld de boekhouding van de precieze vorm van subdomein-randen en buurprocessen voor applicatieprogrammeurs. Met de communicatie-bibliotheek wordt volledige vrijheid bij het opdelen van grillige rekenroosters geïmplementeerd.

2.3

Uitwerking van de rooster-opdeling In een parallelle simulatie wordt ieder deeldomein door een apart rekenproces afgehandeld alsof het een zelfstandig domein is. Hiervoor is het oorspronkelijke programma uitgebreid om te kunnen samenwerken met andere instanties van zichzelf.


22


Cursusboek Module 3

Het basis-principe voor parallellisatie van de berekeningen van WAQUA/TRIWAQ is data-parallellisme. Ieder rekenproces voert in principe dezelfde reeks van stappen uit (de gezamelijke agenda), maar dan alleen op het eigen gedeelte van de data. Het eigen gedeelte van de data wordt gegeven door de opdeling van het globale rooster. Deze opdeling is geïmplementeerd als een array HOWNER(1:nmax,1:mmax) waarin per roosterpunt (array-index) van het globale domein wordt aangegeven welk rekenproces hiervan de eigenaar is. De rekenprocessen werken echter allemaal op een eigen lokale geheugenruimte, en het zou zonde zijn om ze allemaal arrays met de afmetingen van het globale domein te laten bewaren. In plaats daarvan wordt er per deeldomein een uitsnede gemaakt uit het globale domein die voldoende is voor de berekeningen voor het deeldomein. Deze uitsnede van het globale domein bevat alle roosterpunten van het subdomein, plus aan alle kanten drie extra rijen en kolommen van punten van de buurdomeinen. Deze extra punten van de buurdomeinen worden gebruikt voor de opslag van gegevens van de buurdomeinen die nodig zijn in de berekeningen, zoals waterstanden en zout-concentraties, en worden de guard band van het deeldomein genoemd. De guardbandbreedte is instelbaar, en staat standaard op 3 ingesteld ten behoeve van het stoftransport-gedeelte van TRIWAQ. Beschouw als voorbeeld een globaal rooster van 10x30 punten dat stripsgewijs wordt verdeeld in drie deeldomeinen. Stel dat deeldomein 1 loopt van n=1 tot en met 12, deeldomein 2 loopt van n=13 t/m 18 en dat deeldomein 3 loopt van n=19 tot en met 30. Dan worden de uitsnedes van het globale domein respectievelijk [1..10] [1..15] , [1..10] [10..21] en [1..10] [15..30] . Merk op dat er in parallelle berekeningen geen guard band wordt toegevoegd buiten het oorspronkelijke globale domein, dit is anders bij gebruik van domein decompositie met horizontale verfijning. Vervolgens worden de uitsnedes per deeldomein verschoven zodanig dat het eerste punt uitkomt op (1,1). Dit is nodig omdat de ondergrens van allerlei arrays in WAQUA/TRIWAQ is vastgelegd op 1. Tenslotte worden de afmetingen mmax en nmax voor de subdomeinen berekend. In het gegeven voorbeeld zijn de verschuivingen van de drie deeldomeinen de vectoren (0,0), (0,9) en (0,14), en zijn de afmetingen nmax respectievelijk 15, 12 en 21. Deze getallen worden berekend door de partitioner COPPRE bij het aanmaken van de SDS-files per deeldomein. Op deze SDS-files worden ook de bijbehorende uitsnedes van het HOWNER-array gezet: per SDS-file is dit een array POWNER(1:nmax,1:mmax), met de voor dat subdomein bepaalde afmetingen nmax en mmax.


23


2.4

Cursusboek Module 3

Implementatie van numerieke berekeningen (“Loops”) in WAQUA/TRIWAQ De berekeningen in WAQUA/TRIWAQ bestaan grotendeels uit “loops” over alle roosterpunten. Bijvoorbeeld een loop voor het berekenen van de horizontale termen in de impuls-vergelijking, een loop voor de dichtheidsterm, enz. Deze “data-parallelle stappen” worden in WAQUA/TRIWAQ op een specifieke manier geïmplementeerd, namelijk via de “irogeo-tabel”. Deze implementatie wordt toegelicht voordat we de parallellisatie van zulke berekeningen bespreken aan de hand van de gekozen roosteropdeling. De irogeo-tabel geeft een beschrijving van de rijen en kolommen van het rekenrooster. Het rekenrooster is oorspronkelijk als volgt tot stand gekomen: 1) De gebruiker specificeert de afmetingen mmax en nmax van de rechthoek waarbinnen het rooster ligt. 2) De gebruiker geeft de grid-enclosures op. Dit zijn polygonen met horizontale, vertikale of diagonale lijnstukken die het actieve gedeelte van het rooster omsluiten of juist eilanden eruit wegsnijden. De enclosure-punten zelf horen niet tot het echte rekenrooster. Alle punten (m,n) die op deze manier worden geselecteerd worden “cellen” van het rekenrooster. Op alle hoekpunten van de cellen worden kromlijnige coördinaten en dieptecijfers opgegeven, de hoekpunten zijn de diepte-punten van een WAQUA-rooster. In de middens van alle zijvlakken van de cellen worden snelheidspunten gedefinieerd, en de middens van de cellen worden waterstandspunten genoemd. 3) De gebruiker kan op segmenten van de grid-enclosures openingen definieren. Ten behoeve van de randafhandeling wordt er een ring van extra rand-cellen rondom de echte cellen van het rekenrooster gelegd. Verschillende stukken van deze ring mogen niet met elkaar samenvallen; eilanden moeten tenminste twee cellen breed zijn, en er mogen geen scherpe hoeken in enclosures zitten. 4) De programmatuur bepaalt voor welke roosterpunten (m,n) er echte cellen zijn gedefiniëerd en waar er randcellen bestaan. Dit wordt bijvoorbeeld opgeslagen in het kenmerk-array KCS: de codes 0, 1 en 2 worden gebruikt om aan te geven of een punt een landpunt is, of dat het een intern punt of een open randpunt is. Voor de interne punten (KCS=1) wordt ook de irogeo-tabel opgesteld. Figuur 2 geeft een voorbeeld van een irogeo-tabel.


24


Cursusboek Module 3

Title: irogeo.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 2 - Voorbeeld van een rekenrooster inclusief randpunten, met bijbehorende irogeo-tabel (rijen en kolommen)

Aan de hand van de irogeo-tabel kunnen snel alle actieve roosterpunten worden gevonden, terwijl de landpunten worden genegeerd. Bij de implementatie van loops moet wel worden uitgekeken of de loop loopt over waterstandspunten of snelheidspunten, en of de randpunten wel of niet moeten worden meegenomen. Daarom zijn er verschillende soorten loops: Interne waterstandspunten: m=mfu..ml Alle waterstandspunten, inclusief m=mf..mlu, mf=mfu-1, mlu=ml+1.

randpunten:

U-snelheidspunten: m=mf..ml. De variabele-namen mf, mfu, ml en mlu worden altijd op deze manier gebruikt in het gehele WAQUA/TRIWAQ-programma. Let op dat snelheidspunt (m,n) aan de rechterkant van waterstandspunt (m,n) ligt. Dit is conform de WAQUA staggered-grid conventie, de manier waarop de roosters van waterstands, snelheids en dieptepunten op arrays met integer indices worden afgebeeld. Als voorbeeld van een berekening volgens de irogeo-tabel beschouwen we een loop uit de niet-geparallelliseerde versie van TRIWAQ (medio 1999), voor het toevoegen van de bijdrage van de atmosferische druk-term aan de gediscretiseerde impulsvergelijking (subroutine trscue.f, versie 3.5). c c c

--- Calculation of atmospheric pressure gradient term (i.c. isvwp=1). --- Remark: no atmospheric pressure on boundaries if (isvwp.eq.1) then do 670 k=1,kmax do 660 irk=1,norows n=irogeo(1,irk)


25


Cursusboek Module 3

mf=irogeo(2,irk)-1 ml=irogeo(3,irk) nmf=iadres(n,mf) nml=iadres(n,ml) do 650 nm=nmf,nml,minc nmu=nm+minc if (kfu(nm).eq.1 .and. kcs(nm).eq.1 .and. + kcs(nmu).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 650 continue 660 continue 670 continue end if

De buitenste loop (670) is voor alle lagen van het rooster, in alle roosterpunten in het horizontale vlak zijn kmax lagen gedefinieerd. De tweede loop (660) gaat over alle rijen van de irogeo-tabel. Dan worden de coördinaten van de rij gelezen uit de irogeo-tabel: de n-coördinaat van de hele rij, en de mcoördinaten van het eerste en laatste punt van de rij. Merk op dat niet het getal “mfu” maar het getal “mf” wordt bepaald. De loop gaat dus over snelheidspunten en niet over waterstandspunten. Dat is logisch, omdat de impulsvergelijking voor snelheidspunten wordt gediscretiseerd (coëfficiënt ddk). Dan volgen de 1D-adresberekeningen waarmee dezelfde rekenroutine gespiegeld kan worden voor de tweede halve tijdstap, en in loop 650 worden dan alle snelheidspunten van de rij afgelopen. Vervolgens wordt een controle uitgevoerd: de bijdrage aan de impulsvergelijking moet alleen worden gemaakt als het snelheidspunt nat is (geen schotje, kfu(nm)=1), als de twee naastgelegen waterstandspunten interne punten zijn (kcs(nm)=kcs(nmu)=1), en als het snelheidspunt niet is afgesloten door een barrier in laag k (brlsu(nm,k)<=1e6). Tenslotte volgt de berekening zelf: een bijdrage aan het rechterlid van de impulsvergelijking, het verschil van de luchtdrukwaardes in de twee naasgelegen waterstandspunten gedeeld door de afstand tussen die punten en vermenigvuldigd met hdtdw, de halve tijdstap gedeeld door de dichtheid van water.

2.5

Parallellisatie van numerieke berekeningen (“Loops”) in WAQUA/TRIWAQ Na bovenstaande introductie van de irogeo-tabel beschrijven we nu de parallellisatie van berekeningen in WAQUA/TRIWAQ. Daarbij kijken we met name naar dataparallelle berekeningen zoals hierboven geïntroduceerd. Speciale dingen zoals de parallellisatie van solvers en bijvoorbeeld berekeningen voor barriers worden later behandeld in hoofdstuk 3. De parallellisatie is primair gebaseerd op de verdeling van het globale rekenrooster, zie paragraaf 2.3. Als voorbeeld nemen


26


Cursusboek Module 3

we de berekening van de luchtdruk-term. Iedere processor voert de berekeningen uit voor alle roosterpunten van het eigen deeldomein. Daarbij gebruikt hij de irogeo-tabel voor zijn eigen deeldomein, zie Figuur 3. Title: irogeo_sub1.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Title: irogeo_sub2.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 3 - Deelroosters van twee subdomeinen met bijbehorende irogeo-tabellen per subdomein.

De afspraak die wordt gehanteerd bij de parallellisatie van subroutine trscue.f is dat ieder subdomein het stelsel discretiseert en de oplossing bepaalt voor alle Usnelheidspunten van het eigen subdomein. Hoe krijgt een subdomein dit voor elkaar? Hij voert loops uit over alle rijen van zijn subdomein, van punt “mf” tot en met “ml”, en als het punt “mf” van een ander subdomein is dan doet hij hier geen berekening. Het bepalen of er een subdomein-rand is aan het begin van een rij wordt gedaan via het array van rand-codes irobou. Dit gebruikt een van de belangrijkste uitbreidingen aan het WAQUA/TRIWAQ programma voor het samenwerken van de rekenprocessen: de introductie van een speciale rand-code voor het aangeven van een subdomeinrand. De rand-codes per rij/kolom geven normaal gesproken aan hoe het rooster is


27


Cursusboek Module 3

gekoppeld aan de buitenwereld. De gewone rand-codes zijn “gesloten rand”, “waterstandsrand”, “snelheidsrand”, “debietrand” en “Riemann-rand”, in de parallelle versie is hier de “subdomein-interface” bijgekomen. Een waterstandsrand vereist soms andere berekeningen dan een gesloten rand of een snelheidsrand. En een subdomein-interface moet vaak weer anders worden afgehandeld. Het is namelijk geen echt randpunt maar een intern punt van het globale rekendomein waar de differentiaalvergelijkingen moeten worden opgelost, en tegelijkertijd is het ook een punt waarvoor de berekeningen primair door een ander rekenproces zullen worden uitgevoerd. Met behulp van de subdomein-randcodes wordt het stukje code uit trscue.f m.b.t. de luchtdrukterm als volgt geparallelliseerd: c c c

--- Calculation of atmospheric pressure gradient term (i.c. isvwp=1). --- Remark: no atmospheric pressure on boundaries

if (isvwp.eq.1) then do 670 k=1,kmax do 660 irk=1,norows n =irogeo(1,irk) msta=irogeo(2,irk)-1 ml =irogeo(3,irk) if (irobou(1,irk).ge.ibitfc) msta=msta+1 nmsta=iadres(n,msta) nml =iadres(n,ml) do 650 nm=nmsta,nml,minc nmu=nm+minc if (kfu(nm).eq.1 .and. kcs(nm).eq.1 .and. + kcs(nmu).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 650 continue 660 continue 670 continue end if

Het enige verschil in dit stukje code is dat de variabele-naam “mf” is vervangen door “msta”, en dat “msta” met een wordt opgehoogd ten opzichte van zijn default-waarde “mf” in het geval dat de randcode voor het begin van de rij aangeeft dat het een subdomeinrand is. Om de getoonde berekening goed te laten werken in een parallelle run moet echter wel aan een extra voorwaarde (preconditie) worden voldaan. Namelijk, de arrays kcs en pres moeten de goede waardes bevatten in alle elementen die in de berekening worden gelezen. Dat is geen probleem voor de array-elementen die bij het eigen subdomein horen, maar er wordt ook verwezen naar roosterpunten van naburige subdomeinen. Namelijk, voor subdomein 1 van Figuur 3, voor punt m=ml=4 van alle rijen, wordt verwezen naar “mu”=m+1=5, een waterstandspunt van het naburige subdomein. De verwijzingen zitten in de array-referenties


28


Cursusboek Module 3

kcs(nmu) en pres(nmu). De verwijzingen zijn een soort van

data-afhankelijkheden zoals in paragraaf 1.6 zijn besproken. Het zoeken en volgen van deze afhankelijkheden wordt dataanalyse genoemd. Hier wordt uitvoerig op ingegaan in paragraaf 3.5. De data-afhankelijkheden worden gerespecteerd als aan de preconditie wordt voldaan. Dat kan in het algemeen worden gedaan door het invoegen van communicaties voorafgaand aan de getoonde berekening. Hiervoor wordt de “update”-operatie gebruikt. Deze kopiëert de waardes uit de naburige subdomeinen naar de guard band van het eigen subdomein. Ofwel: ieder subdomein verstuurt waardes uit zijn eigen binnengebied naar de andere subdomeinen die hierin zijn geïnteresseerd, en ontvangt daarna de waardes van de buursubdomeinen voor zijn eigen guard band. In het getoonde voorbeeld is communicatie wordt een statisch array gebruikt en een invoer-array. Array kcs wordt eenmalig gevuld gedurende de initialisatie-fase van WAQUA/TRIWAQ. Aan het einde van de initialisatie-fase wordt het eenmalig gecommuniceerd (in subroutine waspiu.f), en voorafgaand aan de getoonde loop is er dan geen communicatie meer nodig. Array pres wordt direct afgeleid uit gegevens van de windSDS-file per subdomein. Bij het aanmaken van deze deel-SDSfiles voegt de partitioner COPPRE al de guard band toe, en zorgt die er ook voor dat de guard band is gevuld met de bijbehorende waardes van naburige subdomeinen. Omdat de waardes in array pres door programma WAQPRO niet worden veranderd, is communicatie van pres voorafgaand aan de bestudeerde berekening niet nodig. In het bovenstaande voorbeeld is beschreven hoe een dataparallelle loop kan worden geparallelliseerd. Daarbij is de “owner computes”-regel gebruikt om te bepalen welk subdomein welke berekeningen van het globale algoritme uitvoert. De benodigde communicaties zijn bepaald via dataanalyse. Nu laten we eerst zien hoe hetzelfde principe uitpakt voor andere soorten loops, en daarna beschrijven we het principe van “redundante berekening” om communicatie te kunnen uitsparen. Een tweede te parallelliseren loop wordt verkregen door de loop van het eerdere voorbeeld te herstructureren. Dit kan door te herkennen dat de checks op kcs(nm) en kcs(nmu) er precies toe leiden dat punten mf en ml worden uitgesloten van de loop. Equivalent met de getoonde loop voor de parallelle versie is daarom: if (isvwp.eq.1) then do 670 k=1,kmax do 660 irk=1,norows n =irogeo(1,irk) mfu =irogeo(2,irk)


29


Cursusboek Module 3

mend=irogeo(3,irk)-1 if (irobou(2,irk).eq.ibitfc) mend=mend+1 nmfu =iadres(n,mfu) nmend=iadres(n,mend) do 650 nm=nmfu,nmend,minc nmu=nm+minc if (kfu(nm).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 650 continue 660 continue 670 continue end if

De loop gaat nu over alle inwendige U-snelheidspunten van het globale domein: m=mfu g..mldg (=mlg-1). Ieder subdomein neemt de doorsnede van deze range met zijn eigen subdomein. Het bepaalt de inwendige snelheidspunten van zijn eigen range mfl..mll. De subscripts zijn hier toegevoegd om onderscheid te maken tussen globale en per-subdomein-waardes. De punten mfl..mldl van het subdomein zijn altijd inwendige snelheidspunten van het globale domein, punt mll is het niet per sé, maar kan het wel zijn. Namelijk als er een subdomeinrand is aan het einde van de lokale rij, en als die subdomeinrand niet samenvalt met een echte rand dan is het punt mll een inwendig punt van de globale rij. In parallel WAQUA/TRIWAQ is het begrip van ownership strikt doorgevoerd, ook voor de echte randpunten van het globale domein. Iedere (m,n)-coördinaat heeft een unieke eigenaar p. Dit betekent dat rekenprocessen soms alleen eigenaar zijn van het randpunt van een rekenrij van het globale domein, terwijl de rest van de rij toebehoort aan een andere processor. Dit is met name niet te vermijden bij diagonale randen (waterstandsranden en gesloten randen). Om het rekenproces te kunnen laten bepalen welke randafhandeling er precies nodig is, is ervoor gekozen om extra rand-codes te introduceren voor samenvallende fysieke randen en subdomeininterfaces. Hiervoor wordt de rand-code voor de gebruikte fysieke rand opgeteld bij de rand-code voor een interface. Dus de rand-code 12 signaleert een interface (code ibitfc=10) die samenvalt met een waterstandsrand (code 2). In het voorbeeld is dit gebruikt door alleen bij nietsamenvallende subdomeinranden (.eq.ibitfc) het randpunt ml mee te nemen in de loop. Merk op dat in het eerdere voorbeeld de conditie .ge.ibitfc is gebruikt. In Figuur 3 is een andere complicatie te zien die ontstaat bij het samenvallen van subdomeinen fysieke randen. In de getoonde roosteropdeling is punt (4,3) van subdomein 2 (globale coördinaten (5,3)) het eindpunt van een rij van het globale rekenrooster. Dit punt zorgt er nu voor dat subdomein 2 een rekenrij bevat van 0 inwendige punten: mfu=4, ml=3! In subdomein 2 heeft deze rekenrij de randcode ibitfc aan het


30


Cursusboek Module 3

begin van het eigen stuk (subdomein interface, punt mf is een inwendig punt van het globale domein), en bijvoorbeeld randcode 1 aan het einde van de rij (gesloten rand). De randcode aan het einde van rij 1 van subdomein 1 is dan ibitfc+1. In de bovenstaande loop levert de speciale rij van subdomein 2 geen problemen op. Loop 650 gaat van m=4 tot m=2, dus wordt geen enkele keer doorlopen. In andere situaties moet men er wel verdacht op zijn dat er rijen kunnen voorkomen zonder interne punten. Een andere manier van parallelliseren dan de “owner computes” regel toe te passen is via “redundante berekening”. Je zorgt er dan voor dat loops niet alleen voor de roosterpunten van het eigen subdomein worden uitgevoerd, maar dat ze ook punten uit de guard band meenemen. Hiermee bereik je dat de uitvoer-variabele van een loop ook in de punten van de guard band wordt ingevuld (post-conditie), en mogelijk kun je hiermee communicatie uitsparen voor een volgende berekening. Dit wordt onder andere gebruikt bij het discretiseren van de impulsvergelijking in subroutine trsumo.f, waarbij ook de luchtdrukterm weer aan de orde is. Bij het bepalen van tridiagonale stelsels voor de waterstanden in subroutine trssuw.f zijn per roostercel de coëfficiënten AA, CC en DD van de impulsvergelijking nodig van de twee naastgelegen Usnelheidspunten. Bij subdomeinranden aan het begin van een rij wordt hierbij verwezen naar roosterpunten van naburige subdomeinen, en moet dus aan een pre-conditie worden voldaan van de vorm “coëfficiënten in punten van de guard band ingevuld”. In plaats van het gebruiken van communicatie voor het invullen van waardes in de guard band (de updateroutine), wordt hier aan de pre-conditie voldaan door de benodigde coëfficiënten in subroutine trsumo.f redundant te berekenen. De post-conditie die wordt opgesteld voor subroutine trsumo.f luidt nu “het berekenen van AAK, CCK en DDK voor alle punten mf..ml van alle rijen van het eigen subdomein” (AAK is per laag, AA is de vertikaal geïntegreerde coëfficiënt). Ten behoeve van deze post-conditie wordt de luchtdrukterm als volgt geparallelliseerd: if (isvwp.eq.1) then do 1960 k=1,kmax do 1955 irk=1,norows n =irogeo(1,irk) msta=irogeo(2,irk) mend=irogeo(3,irk)-1 if (irobou(1,irk).eq.ibitfc) msta=msta-1 if (irobou(2,irk).eq.ibitfc) mend=mend+1 nmsta=iadres(n,msta) nmend=iadres(n,mend) do 1950 nm=nmsta,nmend,minc


31


Cursusboek Module 3

nmu=nm+minc if (kfu(nm).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 1950 continue 1955 continue 1960 continue end if

De meeste loops in trsumo.f gaan nu over de range mf..ml. De getoonde loop moet echter de echte randpunten uitsluiten, dus gaat in een sequentiële run van mfu tot en met mld. In parallelle runs worden de punten mf en ml hier weer bijgetrokken, àls dat inwendige snelheidspunten van het globale domein zijn. De pre-conditie van deze loop verandert door het gebruik van redundante berekening. Doordat de loop kan worden uitgevoerd voor punten “mf” van andere subdomeinen, moeten ook de volgende variabelen in de guard band beschikbaar zijn: kfu(nmf), brlsu(nmf,1:kmax), pres(nmf), gvu(nmf). Het lijkt alsof er nu door gebruik van redundante berekening juist meer wordt gecommuniceerd dan voorheen! Dat is uiteindelijk echter niet het geval. Array pres hoeft helemaal niet te worden gecommuniceerd, gvu is statisch dus slechts eenmalig, en kfu en brlsu moeten toch worden gecommuniceerd, ook als er geen gebruik wordt gemaakt van redundante berekening. Daarentegen moeten AA, CC en DD juist in iedere iteratie van de oplosprocedure voor de continuïteitsvergelijking opnieuw worden gecommuniceerd. Redundante berekening kan vrij gemakkelijk op een vrij kleine schaal worden gebruikt: berekening van laagdiktes uit de laagposities, berekenen van Chezy-waardes uit Manningcoëfficiënten e.d. Voor omvangrijkere toepassing van redundante berekening wordt de benodigde data analyse echter snel zeer complex. Daarom wordt dit in het algemeen niet aangeraden, en zal ook voor parallel WAQUA/TRIWAQ worden onderzocht hoe nodig de redundante berekening van de impulsvergelijking is voor het behalen van goede parallelle performance.


32


3

Cursusboek Module 3

Numerieke aspecten In het vorige hoofdstuk is een overzicht gegeven van de overall strategie die wordt gevolgd in parallel WAQUA/TRIWAQ. Deze bestaat kort gezegd uit het koppelen van min of meer WAQUA/TRIWAQ rekenprocessen,

onafhankelijke

het gebruik van een roosteropdeling voor verdeling van de verantwoordelijkheden over de verschillende rekenprocessen, het gebruik van een guard band voor het opslaan van benodigde waardes voor de implementatie van subdomein-randvoorwaarden, het gebruik van abstracte communicatie-opdrachten voor het up-to-date maken van waardes in de guard band rond subdomeinen, en uit het gebruik van aparte programma’s voor het maken van aparte SDS-files per deeldomein en weer samenvoegen hiervan tot een globale uitvoerfile. In dit hoofdstuk worden de numerieke aspecten hiervan verder uitgewerkt. Het gaat daarbij met name om speciale algoritmes die niet een-twee-drie met data-parallellisme zijn te verdelen, en om het bepalen van welke communicatie op welke punten in het algoritme nodig is.

3.1

Aanpassingen aan solvers I Als eerste bekijken we de Jacobi- en Gauss-Seidel-achtige solvers die in WAQUA/TRIWAQ worden gebruikt: De zgn. “red-black Jacobi” methode voor de impulsvergelijking in TRIWAQ (subroutine trsjam.f). De “red-black Jacobi” methode voor de transportvergelijking in TRIWAQ (subroutine trsjac.f). De “flow-dependent Gauss-Seidel methode” voor het horizontale transport van de turbulente grootheden in TRIWAQ (subroutine trsdgs.f). De “flow-dependent Gauss-Seidel methode” voor de impulsvergelijking in WAQUA (subroutine wasuxc.f). Wat deze solvers met elkaar gemeen hebben is dat ze zonder aanpassingen (red-black Jacobi) of met beperkte aanpassingen (flow-dependent Gauss-Seidel) op een data-parallelle manier kunnen worden geparallelliseerd.


33


3.1.1

Cursusboek Module 3

Red-black Jacobi methode voor impulsvergelijking TRIWAQ Bij de impliciete behandeling van de impulsvergelijking in subroutine trscue.f van TRIWAQ wordt voor de horizontale advectieve termen een tweede orde upwind-schema gebruikt. Dit leidt tot een koppeling van negen onbekenden in het horizontale vlak en een tridiagonale structuur in de vertikale richting:

a0um, n ,k b 2 xum

1 2,n ,k

b 2 y um, n

2, k

b0um, n,k

c0um, n,k

b 1 xum

b 1x um

1, n , k

b 1 y um , n

1 1, n , k

b 1 y um ,n

1, k

1,k

b 2 x um

(3.1)

2,n , k

b 2 y um ,n

2, k

d0

De coëfficiënten a-d zijn verschillend voor alle roosterpunten. Van de coëfficiënten b 2 x , b 2 y zijn er twee nul afhankelijk van de stroomrichting, maar omwille van de eenvoud van de solver wordt dit niet gebruikt. De “red-black Jacobi methode” is een iteratie-methode die is gebaseerd op een schaakbord-verdeling in rode en zwarte punten (m,n) van het horizontale vlak. Iedere iteratie is verdeeld in twee stappen: een update voor alle rode punten, en een update voor alle zwarte punten. De update bestaat uit het oplossen van tridiagonale stelsels voor de vertikale richting, waarbij de horizontale termen “expliciet” op het laatst bekende iteratie-niveau worden gebruikt. Voor iteratie q voor de rode punten wordt dit: e0

d 0 b 2 xumq 12, n , k b 2 y umq , n1

solve a0 , b0 , c0

b 1 xumq 11, n , k

b 1 xumq 11, n , k

b 2 xumq 12, n , k

b 1 y umq , n1 1,k

b 1 y umq ,n1 1,k

b 2 y umq ,n1

2, k

k

umq , n, k

k

e0

2,k

k

De accolades geven hierin aan dat er coëfficiënten en onbekenden zijn voor alle lagen van het model, en dat deze samen een tridiagonale matrix-vergelijking vormen. Voor de zwarte punten, als de waardes uq voor rode punten al berekend zijn wordt de volgende berekening gedaan:

e0

d 0 b 2 x umq 12, n ,k b 2 y umq ,n1

solve a0 , b0 , c0

k

2,k

b 1x umq

b 1 y umq ,n

umq , n ,k

k

b 1xumq

1,n , k

b 1 y umq , n

1, k

e0

1,n , k 1, k

b 2 x umq 12, n , k b 2 y umq ,n1

2,k

k

De red-black Jacobi methode is geïmplementeerd als een sequentie van data-parallelle stappen. Voorafgaand aan de eerste iteratie wordt voor alle roosterpunten een gedeelte van de (constante) tridiagonale matrices geveegd en wordt de Versie 1.0, oktober 2002

34


Cursusboek Module 3

beginschatting ingevuld ( u 0 u oud ). Per halve iteratie wordt eerst voor alle betrokken punten het rechterlid berekend, in de volgende loop wordt het rechterlid voorwaarts geveegd, daarna wordt de terugsubstitutie gedaan, en in een laatste loop wordt de update ten opzichte van de vorige iteratie bepaald ( || u q u q 1 || )ten behoeve van het stop-criterium. In de verschillende loops zijn er afhankelijkheden tussen de verschillende lagen, maar per loop zijn de berekeningen voor verschillende punten van het horizontale rooster volledig onafhankelijk van elkaar. Dit is een terugkerende structuur in TRIWAQ: de koppelingen in de vertikale richting zijn in het algemeen een stuk sterker dan in de horizontale richtingen. Mede hierom is er in parallel TRIWAQ gekozen voor een roosteropdeling van alleen het horizontale vlak, om alle lagen van een bepaald punt steeds aan hetzelfde deeldomein toe te wijzen. De koppeling tussen de verschillende lagen wordt geïllustreerd aan de hand van een code-fragment voor het vegen van de tridiagonale matrices per roosterpunt: c c

First a LU-decomposition of the tri-diagonal structure in vertical direction is made.

1000 1010 1020

do 1020 k=2,kmax do 1010 irk=1,norows n =irogeo(1,irk) mf =irogeo(2,irk)-1 ml =irogeo(3,irk) nmf =iadres(n,mf) nml =iadres(n,ml) do 1000 nm=nmf,nml,minc if (kfu(nm).eq.1) then aak(nm,k)=aak(nm,k)/bbk(nm,k-1) bbk(nm,k)=bbk(nm,k)-cck(nm,k-1)*aak(nm,k) endif continue continue continue

In dit fragment wordt een berekening voor alle Usnelheidspunten uitgevoerd met behulp van de irogeo-tabel (loops 1010 en 1000, zie paragraaf 2.4). De controle kfu.eq.1 selecteert alle natte punten, waar schotjes staan hoeft geen berekening te worden uitgevoerd. Loop 1020 gaat over de lagen van het rooster. Deze loop bevat twee afhankelijkheden. Ten eerste is de berekening van bbk(nm,k) in het tweede statement van de body afhankelijk van het eerste door het gebruik van aak(nm,k). Ten tweede is er een “loop-carried dependency” van het tweede statement terug naar het eerste. De berekening van aak voor laag k=3


35


Cursusboek Module 3

hangt af van de berekende waarde van bbk voor k=2. Loop 1020 is dus puur sequentiëel. In de horizontale richting zijn er daarentegen juist geen afhankelijkheden. Dit komt doordat alle rijen irk die in loop 1010 worden afgelopen tuples (n,mf,ml) opleveren die nooit naar dezelfde nm in loop 1000 verwijzen. Verder gebruikt iedere doorgang van loop 1000 alleen gegevens met arrayindex nm, wordt er niet naar omliggende roosterpunten verwezen. Met betrekking tot het horizontale rooster is de berekening dus volledig data-parallel, en kan ze eenvoudig worden geparallelliseerd door aanpassing van mf bij subdomeinranden (zie paragraaf 2.5). Een vergelijkbaar verhaal kan worden gehouden voor alle andere loops van de red-black Jacobi solver. In veel gevallel zijn er daarbij in het geheel geen referenties naar naburige roosterpunten (zoals in het code-fragment hierboven), zodat er geen pre-condities m.b.t. waardes in de guard band van toepassing zijn. Alleen voor de berekening van het rechterlid voor de tridiagonale stelsels voor de vertikaal zijn waardes van naburige roosterpunten nodig. Daarbij hoeft steeds alleen voor rode of voor zwarte punten te worden gecommuniceerd, en kunnen twee strategiën worden gevolgd. Ofwel communiceren net voordat de waardes nodig zijn (aan het begin van een halve iteratie), ofwel direct nadat ze berekend zijn (aan het einde van een halve iteratie). Dit levert alleen kleine verschillen in de preen post-condities aan het begin en einde van het totale iteratieproces. De guard band moet worden ingevuld voor “twee rijen en kolommen van roosterpunten” rondom het eigen subdomein. Dit wordt geformaliseerd in paragraaf 3.5 (dataanalyse). Het programmeren van de communicatie, inclusief het initialiseren van communicatie-tabellen voor rode en zwarte punten wordt in paragraaf 4.3 beschreven. Tenslotte moet er nog iets met het stop-criterium worden gedaan, dit wordt uitgewerkt in paragraaf 3.4. Een tricky punt bij de implementatie van de red-black Jacobi methode in parallelle runs is dat de onderverdeling rood/zwart wordt bepaald op basis van roostercoördinaten (m,n), terwijl er in parallelle runs met een lokaal verschoven coördinatenstelsel wordt gewerkt. Bij de rood/zwart-verdeling moet de verschuiving dus ongedaan worden gemaakt, moet er met globale (m,n)-coördinaten worden gerekend, en ook de communicatie-tabellen voor alleen rode of alleen zwarte punten moeten consistent hiermee worden gedefiniëerd.

3.1.2

Red-black Jacobi methode voor transportvergelijking TRIWAQ In de transportvergelijking in TRIWAQ (subroutine trsdif.f) wordt een ADI-benadering gebruikt die leidt tot impliciete


36


Cursusboek Module 3

koppeling langs rekenrijen en in de vertikale richting, en tot ontkoppeling van de vergelijkingen in de dwarsrichting. Voor de horizontale advectieve termen wordt een hogere orde upwind-schema gebruikt. Dit leidt tot een stelsel van vergelijkingen met de volgende structuur:

a 1 x c 'm

1, n , k 1

b 3 x c 'm

3, n ,k

b 1 x c 'm

1,n , k

c 1 x c 'm

1, n , k 1

a0 c 'm, n,k

1

a 1x c 'm

1,n , k 1

b 2 x c 'm

2, n , k

b 1x c 'm

1,n , k

b 2 x c 'm

2, n ,k

b 3 x c 'm

3,n ,k

c0 c 'm, n, k

1

c 1 x c 'm

1,n ,k 1

b0c 'm ,n, k

(3.2)

d0

De notatie c ' voor de onbekende concentraties is gebruikt om onderscheid te maken met de coëfficiënten c van het stelsel. De “prime” verwijst naar het nieuwe tijdsniveau, maar dit is bij de bespreking van de impulsvergelijking hierboven weggelaten. De coëfficiënten a-d zijn wederom verschillend voor alle roosterpunten. Van de coëfficiënten b 2 x , b 3 x zijn er twee nul afhankelijk van de stroomrichting, maar omwille van de eenvoud van de solver wordt dit niet gebruikt. De kruistermen a 1 x , c 1x zijn afkomstig van de anti-creeptermen. De red-black Jacobi methode voor de transportvergelijking (subroutine trsjac.f) volgt precies hetzelfde stramien als de solver voor de impulsvergelijking uit de vorige paragraaf. Er wordt een ondervereling gemaakt in rood en zwart, alle horizontale termen worden steeds naar het rechterlid gebracht, en daarna worden er tridiagonale stelsels opgelost voor de vertikaal. Ook wordt de matrix van het tridiagonale stelsel wederom vooraf geveegd (op een iets andere manier dan in trsjam.f). Een verschil in de parallellisatie van deze solver ten opzichte van die voor de impulsvergelijking is dat er nu geen koppeling is tussen verschillende rijen van het rekenrooster. Dit is te zien in vergelijking (3.2), doordat er geen verwijzing wordt gemaakt naar onbekenden met andere n-coördinaten. Dit zou kunnen worden uitgebuit, als de roosterpunten voorafgaand aan de berekening worden herverdeeld dan is er in het iteratieproces geen communicatie meer nodig. Dat herverdelen leidt echter tot een complexe implementatie en vergt ook zelf veel communicatie. Daarom wordt de gewone roosteropdeling gebruikt voor de parallellisatie van deze berekeningen. Een ander verschil met de vorige solver is dat er nu communicatie nodig is voor drie kolommen van roosterpunten links en rechts van het eigen subdomein, en dat er voor het iteratieproces geen communicatie nodig is voor punten boven/onder het eigen subdomein.


37


3.1.3

Cursusboek Module 3

Flow-following Gauss-Seidel methode voor turbulentie-vergelijking TRIWAQ In het k- turbulentie model in TRIWAQ wordt een tijdstapsplitsingsmethode gebruikt voor de tijdsintegratie. In een eerste stap worden tussenwaardes berekend door alleen de horizontale advectieve termen in rekening te brengen, daarna worden de waares op het nieuwe tijdstip berekend door de overige termen in de tijd te integreren (productie, dissipatie, buoyancy, vertikaal transport). De tweede stap leidt tot een aantal data-parallelle stappen voor het berekenen van coëfficiënten, en daarna tot het oplossen van tridiagonale stelsels voor de vertikale richting. Deze berekeningen zijn eenvoudig te parallelliseren op basis van de gegeven roosteropdeling. In de eerste stap is er wel een horizontale koppeling. De horizontale advectieve termen worden met een eerste orde upwind-schema benaderd. Dit leidt tot vergelijkingen van de volgende vorm: b 1 y km* ,n

1, k

b 1 x km*

1, n , k

b0 km* ,n ,k

b 1 x km*

1,n ,k

b 1 y km* ,n

1,k

d0

Hierin is k * de onbekende turbulente energie op de fractionele tijdstap, de vergelijking voor de energie-dissipatie * heeft dezelfde vorm. Door de upwind-benaderingen en door de afwezigheid van een diffusie-term is steeds een van de coëfficiënten b 1x en een van b 1y gelijk aan nul. Welke dat is is afhankelijk van de stroomrichting. Merk verder op dat in deze vergelijking slechts een enkel laagnummer k voorkomt, de vergelijkingen voor verschillende laag-interfaces zijn volledig onafhankelijk van elkaar. Door de speciale structuur met nullen op twee buitendiagonalen is de Gauss-Seidel methode uitermate effectief voor deze stelsels. Wanneer de stroming namelijk overal dezelfde richting heeft, bijvoorbeeld um, n ,k 0, vm ,n ,k 0 voor alle punten binnen een bepaalde laag, dan kan de oplossing van het stelsel in een keer worden berekend door linksonder te beginnen en punt voor punt naar rechtsboven toe te werken. De vergelijking reduceert dan namelijk in ieder punt tot b 1 y km* , n

1, k

b 1x k m*

1, n ,k

b0 k m* , n ,k

d0

en de volgorde van rekenen zorgt er dan voor dat er alleen waardes op het nieuwe iteratie-niveau gebruikt worden, uitgaande van de exacte oplossing op de linkerrand en de benedenrand. Zou je in dit geval in de rechterbovenhoek van het domein beginnen dan worden de buitendiagonalen steeds juist op het oude iteratie-niveau genomen:


38


b0 kmq ,n ,k

d 0 b 1 y kmq ,n1 1, k

Cursusboek Module 3

b 1 x kmq

1 1, n , k

Dit levert effectief de Jacobi-methode op. De exacte oplossing propageert dan slechts met 1 roosterpunt per iteratie naar 1 rechts/boven. Op de rand is k0 de exacte oplossing, dus ook k1,1 2 2 is goed, dan k1,2 , k2,1 enz. Dit laat zien dat de nummering van roosterpunten in de Gauss-Seidel methode cruciaal is voor advectie-gedomineerde problemen met upwind discretisaties.

De manier waarop bovenstaande analyse wordt gebruikt in de solver voor de vergelijkingen van het turbulentiemodel is als volgt. Iedere iteratie wordt onderverdeeld in vier stappen. In de eerste stap worden alle roosterpunten berekend waarvoor u 0, v 0 , in de tweede stap de punten met u 0, v 0 , dan de punten met u 0, v 0 , en tenslotte de punten met u 0, v 0 . In paragraaf 1.6 is besproken dat de Gauss-Seidel methode zonder red-black strategie slechts weinig parallellisme bevat. Dat geldt evengoed voor de hier gebruikte flow-following Gauss-Seidel methode. Bovendien is moeilijk te voorspellen welke roosterpunten tegelijk kunnen worden berekend, en strookt dit niet met de gehanteerde roosteropdeling. In plaats van de methode van het sequentiële programma exact te parallelliseren, precies dezelfde numerieke eigenschappen te behouden, kiezen we er daarom voor om in parallelle sommen een iets andere methode te hanteren. De methode voor parallel rekenen wordt geconstrueerd vanuit het stramien dat wordt gehanteerd binnen parallelle berekeningen: rekenstappen voor hele subdomeinen afgewisseld met updates van de guard band. We kiezen ervoor om aan het begin van iedere iteratie de guard band te updaten voor een ring van een rij/kolom rondom het hele subdomein, en om hiermee per subdomein het sequentiële algoritme toe te passen. Deze wijziging leidt tot een mengeling van de Gauss-Seidel en Jacobi-methodes. Bijvoorbeeld wordt voor een roosterpunt met u 0, v 0 bij aan de onderrand van een subdomein de volgende update-vergelijking toegepast: b0 kmq ,n ,k

d 0 b 1 y kmq ,n1 1, k

b 1 x kmq

1, n , k

De waarde van het roosterpunt (m-1,n) links, van het eigen subdomein, is al berekend en wordt op het nieuwe iteratieniveau genomen, de waarde van het roosterpunt (m,n-1) onder van het buursubdomein wordt op het oude iteratieniveau genomen. Deze aanpassing leidt tot een iets slechtere convergentie in parallelle berekeningen dan in sequentiële berekeningen. Maar Versie 1.0, oktober 2002

39


Cursusboek Module 3

dit valt in het algemeen sterk mee: doordat de advectieve termen van ondergeschikt belang zijn in het turbulentiemodel reikt de invloed van een roosterpunt niet zo heel ver, terwijl de informatie per iteratie wel een heel subdomein door wordt getransporteerd. Vanwege de eventuele verslechtering van de convergentie is er bij de parallellisatie van TRIWAQ een stop-criterium ingevoerd in subroutine trsdgs.f, voorheen werden altijd precies twee iteraties uitgevoerd. Dit stop-criterium moet ervoor zorgen dat een parallelle berekening met veel subdomeinen een kwalitatief net zo goede oplossing oplevert als een sequentiële run. De implementatie hiervan wordt in paragraaf 3.4 behandeld.

3.1.4

Flow-following Gauss-Seidel methode voor impuls-vergelijking WAQUA In de impliciete impulsvergelijking in WAQUA worden de advectieve termen gediscretiseerd met centrale differenties in de rekenrichting en tweede orde upwind-differenties in de dwarsrichting. Dit leidt tot vergelijkingen van de volgende vorm: b 1 xum

1, n , k

b 2 y um, n

2, k

b0um, n ,k b 1 y um, n

b 1x um 1, k

1, n ,k

b 1 y um ,n

1,k

b 2 y um ,n

2, k

d0

(3.3)

Hierin is een van de coëfficiënten b 2 y nul en een van de coëfficiënten b 1y klein door het gebruik van upwinddifferenties. Daarom is de verwachting dat het volgen van de stroomrichting in y-richting voordelig is voor de convergentie van de solver. Binnen een rij van het rooster is er echter een sterke koppeling door het gebruik van centrale differenties, en kunnen er tridiagonale stelsels worden opgelost. En verder moet er rekening mee worden gehouden dat de stroomsnelheid V meestal niet uniform gericht is. De solver die nu in de sequentiële versie van WAQUA gebruikt wordt bestaat uit het bepalen van de dominante richting van de V-snelheid (sommeren van Vm,n), het een voor een updaten van de rijen van het rekenrooster, en per rij oplossen van een tridiagonaal stelsel. De rijen worden in oneven iteraties behandeld in de richting van de overheersende V-snelheid, en in even iteraties in tegengestelde richting afgelopen. Het parallellisme in deze methode is zeer beperkt. In het oplossen van tridiagonale stelsels zit vrijwel geen parallellisme (zie volgende paragraaf), en met iedere volgende rij kan pas worden begonnen als de vorige is afgehandeld. Behalve als er meerdere rijen zijn met gelijke n-coördinaat, maar dit parallellisme is vrijwel niet te exploiteren (te weinig, te onvoorspelbaar, te ongelijkmatig).


40


Cursusboek Module 3

We moeten dus wederom een andere oplosmethode construeren die wel voldoende parallellisme bevat. Dit gebeurt op een vergelijkbare manier als voor het turbulentie-model van TRIWAQ: een keer per iteratie communiceren, en binnen een iteratie alle waardes van buurdomeinen op het oude iteratieniveau gebruiken (Jacobi-koppeling). Deze keuze wordt verder gemotiveerd in (Vollebregt en Roest, 1998). De berekeningen die per subdomein moeten worden uitgevoerd zijn nu niet meer “data-parallel” te noemen. De gevolgde benadering voor de parallellisatie kan nog wel “agenda parallel” worden genoemd. Alle processoren hanteren dezelfde agenda van stappen. Dit wordt ook vaak “SPMD” genoemd (single program, multiple data), maar die term is minder scherp afgebakend. De benodigde communicatie is in dit geval 1 kolom en twee rijen om het eigen subdomein heen, zie vergelijking (3.3). Een speciaal aandachtspunt bij de implementatie van deze solver betreft de start/eindwaardes van de m-coördinaten voor de tridiagonale stelsels per rij. In een sequentiële berekening lopen deze stelsels per rij van mf tot ml, is er een vergelijking voor alle snelheidspunten van de rij. Voor punten mf en ml zijn dit randvoorwaarden, voor mfu t/m mld wordt de volledige impulsvergelijking gediscretiseerd. In geval van subdomeinranden aan het begin en einde van de rij dan wordt dit anders. Het punt mf is dan van een buurdomein, hier moet een Jacobi-vergelijking worden opgesteld: q 1 umf

umq *,1neighbour

Deze vergelijking wordt ook gebruikt als punt mf ook een echt randpunt is van het globale domein (randcodes ibitfc+1..6). Het punt ml is van het subdomein zelf. Als dit een intern punt is van het globale domein (randcode ibitfc), dan moet hier de volledige impulsvergelijking worden gediscretiseerd en is er een Jacobi-vergelijking nodig voor punt mlu. Anders is het een randpunt van het globale domein (samenvallende subdomein/echte rand, randcodes ibitfc+1..6). Dan volstaat het om de echte randvoorwaarde op te leggen voor punt ml, en is er geen vergelijking voor punt mlu. Tenslotte moet er worden opgelet met rijen die bestaan uit een enkel randpunt van het globale domein. In dat geval is mf gelijk aan ml, en dreigen de randvoorwaardes voor het begin en einde van de rij elkaar te overschrijven. Als een rij van een subdomein bestaat uit alleen het startpunt van een rij van het globale domein dan is de echte randvoorwaarde voor mf nodig en bestaat de randafhandeling voor ml uit niets anders dan de Jacobi-vergelijking voor mlu. Voor het eindpunt van een rij wordt de Jacobi-vergelijking voor mf worden gebruikt;


41


Cursusboek Module 3

toepassen van de echte randvoorwaarde voor punt ml stuit op problemen met de vereiste snelheid voor punt mld voor waterstandsranden.

3.2

Aanpassingen aan solvers II Vervolgens bekijken we de solvers voor de continuïteitsvergelijking in WAQUA/TRIWAQ en de solver voor de transportvergelijking in WAQUA. Deze solvers zijn allemaal gebaseerd op gekoppelde vergelijkingen voor rijen van het rekenrooster. Voor de continuïteitsvergelijking zijn dit tridiagonale stelsels, voor de transportvergelijking in WAQUA pentadiagonale stelsels. In sequentiële runs kunnen deze stelsels worden opgelost met de zogenaamde “double sweep” methode. Dit is in feite GaussEliminatie zonder gebruik van pivoting. Deze methode bevat weinig parallellisme. Daarom zijn speciale aanpassingen nodig voor de parallellisatie van deze rekenkernen. Deze aanpassingen worden achtereenvolgens besproken voor de continuïteitsvergelijking in TRIWAQ (subroutine trssuw.f), de continuïteitsvergelijking in WAQUA (subroutine wassuc.f), en de transportvergelijking inWAQUA (subroutine wasdfc.f).

3.2.1

Solver voor de continuïteitsvergelijking in TRIWAQ De algemene strategie die wordt gevolgd voor het oplossen van de gekoppelde impuls- en continuïteitsvergelijkingen staat beschreven in de technische documentatie van TRIWAQ (Zijlema, 1998). Deze bestaat globaal uit het discretiseren van de impuls-vergelijking tot aakm, k

um , k

m

cckm, k

m 1

ddkm ,k , m

mf ..ml

Deze vergelijking wordt vertikaal geïntegreerd met de nieuwste schatting voor de laagdiktes, op basis van de nieuwste iterand voor de waterstanden. De geïntegreerde vergelijking wordt in de continuïteitsvergelijking gesubstitueerd en dit levert een lineair tridiagonaal stelsel. am

q m 1

bm

q m

cm

q m 1

dm , m

mf ..mlu

Dit stelsel wordt opgelost met de double sweep, met de nieuwe waterstand worden droogvalcontroles uitgevoerd en nieuwe laagdiktes berekend. Hiermee wordt opnieuw vertikaal geïntegreerd en opgelost totdat de berekende waterstanden of de resulterende stroomsnelheden nauwelijks meer veranderen (stop-criterium). Het discretiseren van de impulsvergelijking, het berekenen van laagdiktes, de droogvalcontrole en dergelijken zijn allemaal data parallelle berekeningen die op de gewone manier kunnen worden geparallelliseerd. De moeilijkheid zit in het oplossen Versie 1.0, oktober 2002

42


Cursusboek Module 3

van de tridiagonale stelsels met de double sweep methode, een puur sequentiële berekening. Verschillende mogelijkheden zijn overwogen als alternatief voor het oorspronkelijk gehanteerde algoritme. Zo zijn er bijvoorbeeld parallelle algoritmes die in een keer tridiagonale stelsels kunnen kraken: Wang’s partitiemethode, cyclische reductie. Deze hebben als nadeel dat ze twee keer zo veel rekenwerk vergen als de double sweep, en dat ze meerdere communicatie-stappen na elkaar vereisen. De verwachting is daarom dat ze geen goede performance opleveren op computers waar de opstarttijd van communicaties relatief groot is (netwerken van werkstations). Een tweede strategie is om de coëfficiënten van de stelsels te herdistribueren zodanig dat iedere processor een aantal complete stelsels kan oplossen zonder interactie met de anderen. Dit wordt een vrij complexe strategie als je ook het rekenwerk van de droogvalcontrole e.d. wilt herverdelen, en als je dat niet doet dan moet je in iedere iterate veel communiceren. Daarom is in eerste instantie gekozen voor een strategie à la de impulsvergelijking in WAQUA: gebruiken van Jacobi-vergelijkingen op subdomeinranden. In een later stadium is daar een verbetering op gemaakt die de convergentie in parallelle runs sterk verbetert: combinatie met de TWGE-methode. TWGE staat voor two-way Gaussian-Elimination. Deze methode is erop gebaseerd dat het in principe niet uitmaakt of je een double-sweep aan de ene kant van een rij begint of aan de andere kant begint. Je kunt dan ook met twee processoren tegelijk beginnen te vegen aan het stelsel. Bij de subdomeinrand wisselen de twee processoren de geveegde vergelijkingen met elkaar uit, en kunnen ze de terugsubstitutie uitvoeren. Dit staat verder uitgewerkt in de technische documentatie van TRIWAQ en in (Vollebregt, 1997). Wanneer een rij van het globale rekenrooster is verdeeld over meer dan twee processoren dan gaat deze truc weer niet op. In dat geval worden de sub-rijen in paren van twee verdeeld: (1,2), (3,4), (5,6). Per paar wordt TWGE toegepast, en op de overgangen tussen twee paren worden Jacobi-vergelijkingen gebruikt. In de volgende iteratie wordt een andere verdeling in paren gebruikt: (1), (2,3), (4,5), (6). Dit zorgt ervoor dat er in de ene iteratie een sterke koppeling/informatieuitwisseling is op de ene helft van de subdomeinranden, en in de andere iteratie op de andere helft van de subdomeinranden. Numeriek gezien kan dit BJ-TWGE algoritme worden beschouwd als een domein decompositie-algoritme met sterk overlappende domeinen. De overlap zorgt voor een grote verbetering van de convergentieeigenschappen ten opzichte van het blok-Jacobi algoritme dat equivalent is met een domein decompositie-algoritme zonder overlap van de domeinen.


43


Cursusboek Module 3

De BJ-TWGE methode is geïmplementeerd met behulp van een rood-zwart kleuring van de deel-rijen per subdomein. In die kleuring wordt onderscheid gemaakt tussen globale rijen die in 1, 2 of meer stukken zijn verdeeld. In het eerste geval is de rij in de parallelle som “altijd rood”, d.w.z. in alle iteraties wordt de forward/backward sweep uitgevoerd. In het tweede geval is de eerste deelrij “altijd rood” en de tweede deelrij “altijd zwart”. Voor zwarte deelrijen wordt de richting van de double sweep veranderd, wordt een backward/forward sweep uitgevoerd. Als een rij in meer dan twee stukken is verdeeld dan wordt de eerste bijbehorende deelrij “afwisselend rood/zwart”, de tweede “afwisselend zwart/rood”, en zo om en om tot de laatste deelrij. Een deelrij die “afwisselend rood/zwart” is is “rood” in oneven iteraties (forward/backward sweep) en “zwart” in even iteraties (backward/forward sweep). Merk op dat het TWGE-algoritme leidt tot exacte oplossing van lineaire tridiagonale stelsels, terwijl de Jacobi-koppeling leidt tot een benaderende oplossing. Een effect van de iteratie-fout is dat de overall iteratiemethode voor de gekoppelde impuls- en continuïteitsvergelijkingen niet langer massabehoudend is op iteratie-niveau. Om toch massabehoud te bereiken in parallelle berekeningen wordt een correctie-stap uitgevoerd na afloop van het iteratieproces. Deze massacorrectie zou niet nodig zijn als er voor een directe parallelle solver was gekozen zoals Wang’s methode. De analyse van de benodigde communicatie is voor deze solver een stuk lastiger dan voor de solvers uit paragraaf 3.1. Naast de waterstanden (binnen iteratieproces voor 1 kolom aan weerszijden naast het subdomein, na afloop ook voor 1 rij boven/onder het subdomein) moeten ook de coëfficiënten voor het TWGE-algoritme worden gecommuniceerd. Tenslotte is er communicatie nodig voor de schotjes en stroomsnelheden, maar deze communicatie wordt uitgespaard door het gebruik van redundante berekening. De gedetailleerde data-analyse voor subroutine TRSSUW is te vinden in de technische documentatie van TRIWAQ.

3.2.2

Solver voor de continuïteitsvergelijking in WAQUA De oplosmethode voor de gekoppelde impuls- en continuïteitsvergelijkingen in subroutine WASSUC van WAQUA is grotendeels hetzelfde als de overeenkomstige solver in TRIWAQ. Wederom wordt eerst de impulsvergelijking gediscretiseerd, wordt deze per iteratie in de continuïteitsvergelijking gesubstitueerd, worden de resulterende tridiagonale stelsels per rij opgelost, en worden droogvalcontroles uitgevoerd en stroomsnelheden berekend.


44


Cursusboek Module 3

Het grootste verschil tussen deze solver en die van TRIWAQ is dat het iteratieproces in WAQUA per rij wordt uitgevoerd in plaats van voor alle rijen tegelijk. Dit gold met name voor de niet-geparallelliseerde versie van WAQUA, en is nog slechts gedeeltelijk waar voor de parallelle versie. De methode werkte vroeger voor iedere rij van het rekenrooster afzonderlijk:

voor rij = 1, .., norows doe - voer droogval/onderloop-controles uit - bepaal coefficienten impuls/continuiteitsvgl. - los stelsel op met iteratieve procedure en met droogval-controles In dit schema zijn de berekeningen van verschillende rijen van elkaar afhankelijk door de droogvalcontroles. Namelijk, de coëfficiënten van de impulsvergelijking hangen af van de actuele geometrie in naastgelegen rijen. Dus als er in rij 1 iets droogvalt dan heeft dat effect op rij 2, en dus kunnen de berekeningen voor rijen 1 en 2 niet tegelijkertijd worden uitgevoerd. Vanwege deze afhankelijkheid had de oorspronkelijke methode te weinig parallellisme: alle processoren zouden steeds samen aan een enkele rij moeten werken, en per rij is er niet zo heel veel gelijktijdig te doen. Wanneer de processoren toch samen aan een rij gaan werken, dan treedt meteen het tweede probleem op. Ze moeten dan ontzettend vaak met elkaar gaan communiceren, voor iedere iteratie voor iedere rij afzonderlijk. Vanwege de opstarttijd van communicaties in parallelle computers (latency) zou er vrijwel geen enkele versnelling mogelijk zijn. Om dit te verbeteren wordt momenteel het volgende schema gebruikt:

voor rij = 1, .., norows doe - voer droogval/onderloop-controles uit voor rij = 1, .., norows doe - bepaal coefficienten impuls/continuiteitsvgl. voor iteratie = 1, 2, ... doe voor rij = 1, .., norows doe - bepaal coefficienten tridiagonale stelsels voor rij = 1, .., norows doe - los tridiagonale stelsels op ..... droogval-controles, stroomsnelheden, ....

(3.4)

Dit schema voldoet aan hetzelfde stramien als de andere berekeningen in WAQUA/TRIWAQ: min of meer dataparallelle berekeningen (agenda parallel) voor het hele veld.


45


Cursusboek Module 3

Het iteratieproces wordt toch nog min of meer per rij apart uitgevoerd. Hiervoor wordt er per rij apart bepaald hoeveel iteraties er nodig zijn, op basis van het stop-criterium max m um[ q,]n um[ q, n1] vel . Dit is ingepast in het schema (3.4) via een array ianthq(1:norows) dat per rij aangeeft of de rij al geconvergeerd is of niet: c----------------------------------------------------------c calculate new depth at u-points c----------------------------------------------------------do 3000 ic=1,norows if (ianthq(ic).ne.0) then n =irogeo(1,ic) mf =irogeo(2,ic)-1 ml =irogeo(3,ic)

+ 2900 3000

do 2900 m=mf,ml ... hu(nm)=0.5*(dp(nm)+dp(ndm))+ tetau(nm)*sq(nm)+(1.0-tetau(nm))*sq(nmu) continue end if continue

De waardes van ianthq worden per iteratie bepaald door middel van communicatie tussen de verschillende subdomeinen (zie paragraaf 3.4 over convergentie-criteria). Het stop-criterium per rij zorgt voor onvoorspelbare verschillen tussen de werklast van de verschillende subdomeinen. Het is namelijk vooraf niet te zeggen hoeveel iteraties er per rij nodig zullen zijn, dit variëert per rij en gedurende de simulatie. Wanneer er in een iteratie al een aantal rijen geconvergeerd zijn, en de ene processor heeft een groter deel dan de andere processor, dan wordt de initiëel gelijkmatige werkverdeling verstoord. Daarnaast wordt er bij het communiceren van waterstanden, schotjes en coëfficiënten voor het TWGE-algoritme geen rekening gehouden met de rijen die al geconvergeerd zijn; de hele guard band wordt bijgewerkt. Tenslotte kost de communicatie voor het stop-criterium per rij in een aantal situaties relatief veel tijd. Door deze factoren is het in parallelle runs relatief minder voordelig om per rij te itereren dan in sequentiële runs. Een stuk van het snelheidsverschil met TRIWAQ wordt erdoor opgeheven, ofwel de speedup die met TRIWAQ kan worden behaald is in het algemeen hoger dan wanneer hetzelfde model met WAQUA wordt doorgerekend. De data-analyse ten behoeve van de benodigde communicatie in subroutine WASSUC van WAQUA komt sterk overeen met die van subroutine TRSSUW van TRIWAQ. Wederom wordt er redundante berekening gebruikt om communicatie uit te


46


Cursusboek Module 3

sparen. De complete analyse documentatie van WAQUA.

3.2.3

staat

in

de

technische

Solver voor de transportvergelijking in WAQUA De transportvergelijking wordt in WAQUA op een vergelijkbare manier gediscretiseerd als in TRIWAQ, namelijk met een ADI-benadering. Dit leidt tot gekoppelde stelsels vergelijkingen per rij van het globale rekenrooster en tot ontkoppeling in de dwarsrichting. Omdat de transport-vergelijking lineair is, omdat er slechts een laag wordt gerekend, en vanwege het tweede orde upwindschema voor de advectieve termen is dit een lineair pentadiagonaal stelsel: b1

c1

f1

x1

d1

a2

b2

c2

f2

x2

d2

e3

a3

b3

c3

f3

x3

d3

e4

a4

b4

c4

f4

x4

d4

e5

a5

b5

c5

x5

d5

e6

a6

b6

x6

d6

(3.5)

Hierin is de notatie voor de coëfficiënten afgestemd op de variabele-namen in de programmatuur, en staat x voor de onbekende concentraties. Deze pentadiagonale stelsels worden in sequentiële berekeningen opgelost met een double sweep algoritme (Gauss Eliminatie). Dat is uitermate praktisch omdat dat algoritme in een keer de exacte oplossing oplevert. Merk op dat bij de continuïteitsvergelijking alsnog moet worden geïtereerd vanwege de niet-lineariteiten in het daarbij op te lossen stelsel. In parallelle runs wordt de double sweep op dezelfde manier geparallelliseerd als voor de continuïteitsvergelijking, namelijk via de BJ-TWGE methode. Er wordt weer onderscheid gemaakt tussen rijen van het globale rekenrooster die in 1, 2 of meer stukken zijn verdeeld. Rijen die niet in stukken zijn verdeeld worden “altijd rood” genoemd en hiervoor wordt het sequentiële double sweep algoritme gebruikt. Rijen die in twee stukken zijn verdeeld worden met het TWGE-algoritme aangepakt: het eerste stuk heet “altijd rood” en hiervoor wordt een forward/backward sweep uitgevoerd, het tweede stuk is “altijd zwart” en gebruikt een backward/forward sweep. Merk op dat het TWGE-algoritme hiervoor in een keer de exacte oplossing bepaalt. Rijen die in drie of meer stukken worden verdeeld worden opgelost met Jacobi-vergelijkingen erbij. De deelrijen worden afwisselend “afwisselend rood/zwart” en “afwisselend zwart/rood” gelabeld. Door het geruik van Jacobi-


47


Cursusboek Module 3

vergelijkingen is de oplossing na 1 sweep niet exact, en moet er verder worden geïtereerd. Per iteratie wordt een verdeling in paren bepaald voor toepassing van TWGE. De introductie van een iteratieproces ten behoeve van parallelle berekeningen vergt ook de introductie van een stop-criterium. Dit is op dezelfde manier uitgewerkt als voor de continuïteitsvergelijking in WAQUA: het benodigde aantal iteraties wordt per rij afzonderlijk bepaald. In dit stop-criterium wordt rekening gehouden met rijen van het globale rooster die niet of slechts in twee stukken worden onderverdeeld: voor deze rijen wordt slechts 1 iteratie uitgevoerd. In het geval dat alle rijen van het globale rooster in ten hoogste twee stukken zijn verdeeld dan is iteratie overbodig, en dan wordt ook het stop-criterium niet ge-evalueerd. Per iteratie moeten de concentraties worden gecommuniceerd voor twee kolommen aan weerszijden van ieder subdomein. Waar er TWGE is gebruikt zou deze communicatie kunnen worden uitgespaard, behalve dat dit verschillen oplevert op het niveau van afrondfouten. In iedere iteratie moeten verder de TWGE-coëfficiënten worden gecommuniceerd. Na afloop van het iteratieproces wordt er gecommuniceerd voor de dwarsrichting, voor twee rijen boven/onder het eigen subdomein.

3.3

Vereiste herstructurering van de sequentiële programmatuur In de vorige paragrafen is de strategie beschreven voor de parallellisatie van alle iteratieve solvers in WAQUA/TRIWAQ. Daarbij is aan de orde gekomen dat de numerieke methoden uit de sequentiële versie soms moeten worden aangepast, en dat ook vaak de implementatie van de numerieke technieken moet worden geherstructureerd. In deze paragraaf geven we een overzicht van aanpassingen aan sequentiëel WAQUA/TRIWAQ.

3.3.1

Solvers en convergentie-criteria De gemaakte aanpassingen aan numerieke technieken zijn vooral het gebruik van Jacobi-vergelijkingen op subdomeinranden en de introductie/wijziging van stop-criteria. Daarnaast zou men ten behoeve van de parallellisatie andere differentiebenaderingen kunnen gebruiken op subdomeinranden. Dit wordt niet gedaan, omdat een belangrijk uitgangspunt bij de parallellisatie is geweest dat gebruikers niets van de parallellisatie zouden moeten merken behalve de hogere rekensnelheid. De nauwkeurigheid van een parallelle run moet daarom hetzelfde zijn als die van een sequentiële run. De Jacobi-vergelijkingen zorgen voor een ander verloop van de convergentie van iteratieve procedures, maar leiden in principe


48


Cursusboek Module 3

niet tot een ander eindresultaat indien er wordt doorgeïtereerd tot op machine-nauwkeurigheid. Een uitzondering hierop komt door het algoritme voor droogvallen en onderopen: per iteratie wordt er gekeken of er roosterpunten moeten worden drooggezet, en dit kan in parallelle runs nu net iets anders uitpakken dan in sequentiële runs. Dit is een beperking van het huidige algoritme voor droogvallen en onderlopen die in de toekomst zou kunnen worden opgelost. Omdat het verloop van de convergentie in parallelle runs anders kan zijn dan in sequentiële runs zijn ook extra stopcriteria ingevoerd: in het turbulentie-gedeelte van TRIWAQ, in de impulsvergelijking in WAQUA, en in de transportvergelijking in WAQUA. Het convergentie-gedrag is verder afhankelijk van hoeveel en waar er Jacobi-vergelijkingen worden gebruikt, dus van de partitionnering van het globale rekenrooster en van het gebruikte aantal subdomeinen. In plaats van een vast aantal iteraties (turbulentie: 2, impulsvgl. WAQUA: ITERMOM, transport WAQUA: 1) wordt er nu daarom doorgeïtereerd totdat de aanpassingen aan de schatting van de oplossing kleiner zijn dan een vooraf ingestelde iteratienauwkeurigheid. In de transporten turbulentievergelijkingen zijn deze toleranties hard ingecodeerd, voor de impulsvergelijking kan de gebruiker de nauwkeurigheid opgeven in de simulatie-invoerfile. Bij de parallellisatie is ook een nieuw stopcriterium toegevoegd voor de continuïteitsvergelijkingen in WAQUA en TRIWAQ. Het nieuwe criterium kijkt naar de updates in de waterstanden in plaats van naar de updates in de stroomsnelheden. Deze wijziging heeft niet direct iets te maken met de parallellisatie. Het criterium is ingevoerd omdat hiermee de berekening van stroomsnelheden in iedere iteratie kan worden uitgespaard, terwijl de uiteindelijke oplossing na hetzelfde aantal iteraties er niet door verandert. Het aantal iteraties verandert wel licht, afhankelijk van de nauwkeurigheid die wordt ingesteld, maar deze nauwkeurigheid wordt toch min of meer arbitrair gekozen.

3.3.2

Aanpassingen routines voor QH-openingen Naast de beschreven aanpassingen aan de solvers zijn er op verschillende andere plaatsen in WAQUA/TRIWAQ stukken berekeningen geherstructureerd. Daarbij ging het vooral om wijzigingen van de implementatie, zonder gevolgen voor de berekende resultaten: Herstructurering van subroutine WASSUC van WAQUA (continuïteitsvergelijking), loop over alle rijen tegelijk in plaats van per rij afzonderlijk. Deze aanpassing heeft wel invloed op het rekenresultaten, en is besproken in paragraaf 3.2.2.


49


Cursusboek Module 3

Herstructurering van subroutine WASCBA voor de dynamische barriersturing (openen/sluiten van barriers op basis van een “regime”, afhankelijk van voorspelde waterstanden, snelheden, concentraties e.d. in plaats van vooraf opgegeven via een tijdserie). Herstructurering van subroutines WASFQH en WASFQA voor QH-openingen (randvoorwaarde dynamisch bepaald uit een QH-relatie) en QADopeningen (debiet-randvoorwaarde met automatische verdeling langs een opening). In al deze gevallen is het doel van de herstructurering geweest om dezelfde berekeningen te kunnen uitvoeren met een kleiner aantal communicaties in parallelle berekeningen. Dit vanwege de latency, de wachttijd die in rekening wordt gebracht bij iedere communicatie in de meeste parallelle computers. Ter illustratie van de gemaakte aanpassingen bekijken we subroutine WASFQH voor QH-openingen. Het oorspronkelijke algoritme kan als volgt worden samengevat (voor een vertikale opening): voor alle QH-openingen i doe: sommeer Qi

m (i )

m

H mU U m langs de opening

bepaal H i uit QH-tabel voor opening i zet H i in array van randvoorwaarden m(i )

Het probleem met dit algoritme voor de parallellisatie is dat openingen over meerdere subdomeinen kunnen worden verdeeld, en dat er voor iedere opening afzonderlijk moet worden gecommuniceerd om het totale debiet te bepalen. De strategie die in de parallelle versie wordt gevolgd is om eerst alle debieten te bepalen (in een werk-array in plaats van een werk-variabele), inclusief communicatie, en daarna voor alle openingen de bijbehorende H op te zoeken en de randvoorwaarde in te vullen. Dit vergt slechts 1 communicatie in plaats van per opening afzonderlijk. De berekening van het debiet per opening is als volgt geparallelliseerd: ieder subdomein initialiseert het werk-array op 0, telt voor iedere QH-opening op wat de bijdrage aan het debiet is voor de eigen snelheidspunten, en de som over alle processoren is dan het totale debiet per opening van het globale domein:


50


Cursusboek Module 3

voor alle QH-openingen i doe: sommeer Qi

m (i )

m

H mU U m in eigen subdomein

sommeer Qi over alle subdomeinen (communicatie) voor alle QH-openingen i doe: bepaal H i uit QH-tabel voor opening i zet H i in array van randvoorwaarden m(i) eigen subdom.

Het sommeren van debieten voor snelheidspunten voor het eigen subdomein is automatisch geïmplementeerd door de manier waarop globale SDS-files worden gepartitionneerd. Een subdomein-SDS-file bevat hierdoor alleen informatie over het eigen aantal openingen, en de coördinaten-arrays (kbm,kbn) voor openingspunten bevatten alleen de coördinaten van de punten per opening van het eigen subdomein. Hierin is het uitgangspunt te zien dat subdomeinen zoveel mogelijk als complete zelfstandige rekenprocessen worden beschouwd. Voor het implementeren van de sommatie over alle subdomeinen is wel informatie over het globale domein nodig. Ten eerste moet het totale aantal QH-openingen bekend zijn voor de lengte van het gebruikte werk-array, ten tweede moet er een mapping worden gegeven van de eigen QH-openingen naar de bijbehorende globale QH-openingsnummers. Er kunnen bijvoorbeeld 5 QH-openingen in het globale domein zijn terwijl ntoqh=2 voor een subdomein; dat subdomein moet dan weten van welke van de 5 globale openingen hij een gedeelte heeft. Voor de lengte van het werk-array wordt het totaal aantal openingen ntoglb van het globale domein gebruikt, dit wordt door partitioner COPPRE aan de SDS-file toegevoegd (in array PCONST). De mapping van lokaal naar globaal wordt voor alle openingen opgeslagen in array iglopn. Dit array wordt in WAQPRE aangemaakt met de identiteit “globaal openingsnummer = eigen openingsnummer”, en wordt dan automatisch door COPPRE per subdomein aangemaakt door de waardes van de deel-openingen per subdomein te selecteren. De sommatie over de subdomeinen zelf wordt uitgevoerd met een zogenaamde “globale communicatie-operatie”, de globale som. Dit is een subroutine die door alle betrokken subdomeinen tegelijk wordt aangeroepen met een array van eigen bijdragen, en die bij return de globale waardes retourneert: elementsgewijs de som over alle subdomeinen. In het communicatie-schema van globale communicaties spelen alle betrokken processoren een actieve rol. Pas recentelijk is de communicatie-routine uitgebreid zodat de betrokken processoren kunnen worden gekozen door het definiëren van groepen van processen. In de huidige operationele versie van WAQUA/TRIWAQ moeten nog steeds alle rekenprocessen meedoen in globale communicaties. Vanwege deze beperking


51


Cursusboek Module 3

wordt subroutine WASFQH momenteel door alle rekenprocessen aangeroepen als er ergens in het globale domein een QH-opening bestaat, ook door rekenprocessen waarvoor het eigen subdomein geen openingen of administratie van QH-tabellen bevat. Deze beschrijving van de parallellisatie van QH-openingen is typerend voor een bepaald soort van herstructureringen van de niet-parallelle versie van WAQUA/TRIWAQ. Bijvoorbeeld gaat de herstructurering van de routines voor dynamische barrier-sturing vrijwel analoog: de oorspronkelijke niet-parallelle berekening bestaat uit een enkele loop (over alle “condities”) die communicatie vereist tussen verschillende subdomeinen (gemiddelde waterstand over een lijn, waardes uit verschillende punten van het globale domein), de parallellisatie-strategie is eerst een loop voor het bepalen van lokale waardes, dan globale communicatie, en dan het gebruik van de waardes, en de implementatie vereist een extra werk-array, een stukje extra administratie (bijv. globale nummers), en het splitsen van een enkele loop in twee of meer aparte loops.

3.4

Parallellisatie van convergentie-criteria In de diverse iteratieve procedures in WAQUA/TRIWAQ worden verschillende convergentie-criteria gebruikt. Deze criteria zijn steeds van de vorm “als maximum update tussen opeenvolgende iteraties groter dan grenswaarde dan nog een iteratie”:

zolang( x [ q ]

x [q

1]

max m, n xm[ q,]n

xm[ q,n1]

abs

) doe ...

Hierin wordt steeds het maximum genomen over alle roosterpunten (m,n) van het globale rekenrooster. In parallelle berekeningen is het praktisch om in alle subdomeinen precies evenveel iteraties uit te voeren. Dit vereenvoudigt de communicatie van waardes nabij subdomeinranden, omdat er “blind” waardes naar alle buurdomeinen kunnen worden verstuurd zonder er rekening mee te houden of die buurdomeinen wel op waardes zitten te wachten. Verder is het aantal iteraties in sequentiële berekeningen ook voor alle roosterpunten hetzelfde. Het in alle subdomeinen precies evenveel iteraties uitvoeren is in lijn met de “agenda-parallelle” manier van parallel rekenen, waarbij alle rekenprocessen dezelfde agenda van activiteiten uitvoeren, en steeds min of meer in dezelfde fase van de berekeningen zitten, met elkaar in de pas blijven lopen. Dit vereenvoudigt de analyse van gebeurtenissen in een parallelle run.


52


Cursusboek Module 3

Het basisprincipe voor de parallellisatie van iteratieve procedures en stop-criteria in parallel WAQUA/TRIWAQ is dus om na afloop van iedere iteratie een globaal stop-criterium te evalueren, en afhankelijk daarvan doen ofwel alle of geen van alle processoren een volgende iteratie. Het evalueren van de globale stop-criteria kan op een paar verschillende manieren worden gedaan met behulp van globale communicatie. Er kunnen residuen of vlaggen worden gecommuniceerd, en de laatsten kunnen met plus of met max met elkaar worden gecombineerd. Tenslotte is het voor de analyse van het convergentie-gedrag van solvers en meer in het algemeen voor het onderzoeken van model-instabiliteiten prettig om de lokatie ((m,n)-coördinaten) van het grootste residu te printen, en dit moet worden meegenomen in de globale communicatie. Voor deze verschillende opties kunnen diverse globale communicatie-routines worden gebruikt, zie paragraaf 4.5. In parallel WAQUA/TRIWAQ zijn twee subroutines geprogrammeerd voor alle verschillende convergentie-checks in de diverse iteratieve procedures (waschk.f, wasck2.f). Een aantal speciale features van deze routines zijn als volgt: In de transporten turbulentiegedeeltes van WAQUA/TRIWAQ kunnen de absolute groottes van de onbekenden sterk variëren, zodat het moeilijk is vooraf een goede waarde voor de iteratie-nauwkeurigheid in te stellen. Daarom is het hier beter om een relatief stopcriterium zoals

zolang( x [ q ]

x [q

1] rel

x [ q ] ) doe ...

of

zolang( max m, n ( xm[ q,]n

xm[ q,n1] ) / xm[ q,]n

rel

) doe ...

te gebruiken. De eerste van deze twee varianten vereist aparte combinatie van de maximale updates x [ q ] x [ q 1] en de maximale waardes x [ q ] per subdomein, pas daarna kan het criterium worden geevalueerd. Hiervoor zijn de convergentie-check routines uitgebreid met extra argumenten voor het activeren van het relatieve stop-criterium en het doorgeven van de maximale waarde per subdomein. De tweede variant is een stuk strenger ten aanzien van roosterpunten met relatief kleine waardes van xm[ q,]n en wordt daarom in praktijk niet gebruikt. In de solvers voor de continuïteitsvergelijking in WAQUA/TRIWAQ kan het communiceren voor de stop-criteria worden gecombineerd met de


53


Cursusboek Module 3

communicatie over al of niet droogvallen in het domein en met de controle op instabiliteit van de simulatie (onrealistisch grote waterstanden). Het is efficiënt om deze verschillende globale controles te combineren omdat globale communicatie-operaties relatief duur zijn: er is een aantal opeenvolgende communicatiestappen nodig van orde 2 log( p) , met p het aantal rekenprocessen/subdomeinen. De convergentie-controle per iteratie vergt relatief veel communicatie-tijd, en deze tijd neemt steeds verder toe naarmate er met meer processoren wordt gewerkt (slecht schaalbaar). Daarentegen is het benodigde aantal iteraties van te voren al aardig in te schatten: ongeveer gelijk aan het vereiste aantal iteraties van de vorige tijdstap. Deze twee gegevens zijn gecombineerd in het zogenaamde minit-schema. Om de kosten van de convergentie-controle te beperken doen we in parallelle runs steeds tenminste evenveel iteraties als er in de vorige tijdstap nodig waren (het minimum aantal iteraties minit), en slaan we de convergentie-check in de eerste minit-1 iteraties over. In de eerste tijdstap is minit 1 en wordt er direct in de eerste iteratie gecommuniceerd. In een latere tijdstap is minit bijvoorbeeld 4. In dat geval worden in de vierde iteratie alle residuen van iteraties 1 t/m 4 gecommuniceerd. Op basis daarvan kan worden gekeken of er minder dan vier, vier of meer dan vier iteraties nodig zijn/waren. Het precieze schema staat in de technische documentatie van TRIWAQ uitgewerkt. In WAQUA wordt het benodigde aantal iteraties per rij afzonderlijk bepaald. Hiervoor wordt een array van vlaggen over alle subdomeinen gecombineerd: per rij een vlag of er wel of niet aan het stop-criterium wordt voldaan. Pas zodra er convergentie is bereikt voor alle rijen dan worden de grootste residuen over het hele veld gecommuniceerd ten behoeve van print-uitvoer. De bijbehorende subroutine wasck2.f ondersteunt ook droogvallen per rij en een variant van het minit-schema per rij, en is daardoor behoorlijk complex geworden. Bovendien neemt de communicatie-tijd toe naarmate er meer rijen in het globale rooster zijn en meer subdomeinen worden gebruikt (complexiteit 2 norows log( p) ), dus is deze routine slecht schaalbaar. Hiervoor moet eens een beter alternatief worden uitgewerkt.


54


3.5

Cursusboek Module 3

Voortgezette data-analyse In hoofdstuk 1 is als algemene strategie voor het parallelliseren van een numerieke applicatie gegeven het op elkaar afstemmen van parallelle numerieke algoritmes, applicatie-klasse, probleem decompositie strategie, programmeer-paradigma en te gebruiken parallel rekenplatform. In hoofdstuk 2 is beredeneerd dat voor WAQUA/TRIWAQ een strategie op basis van roosteropdeling (=probleem decompositie), en met onafhankelijke rekenprocessen met lokaal geheugen en communicatie van guard band waardes (=programmeer paradigma) het meest geschikt is. De implementatie van de gekozen strategie is gedeeltelijk in hoofdstuk 2 gepresenteerd (implementatie van de roosteropdeling, parallellisatie van berekeningen m.b.v. de irogeo-tabel). Daarnaast is er in de voorgaande paragrafen aangegeven welke aanpassingen aan numerieke algoritmes er zijn gebruikt: ter verbetering van het parallellisme en om beter te passen binnen het gebruikte programmeer paradigma (“agenda parallel”). In deze paragraaf bespreken we tenslotte hoe de benodigde communicatie kan worden bepaald met behulp van data-analyse. Het doel van data-analyse is het bepalen van welke communicatie en synchronisatie er nodig is om te voldoen aan de data-afhankelijkheden van het gekozen parallelle algoritme. Binnen het gekozen stramien van min of meer data-parallelle “super-stappen” (agenda parallelle benadering) bestaat dataanalyse uit het bepalen van de preen post-condities van alle superstappen, en ervoor zorgen dat aan deze preen postcondities wordt voldaan. In de gehanteerde “process/channelbenadering” wordt de benodigde synchronisatie automatisch uitgevoerd wanneer processen wachten op de te ontvangen data van buur-processen. De preen post-condities bestaan dus alleen uit de gebruikte en gewijzigde data per super-stap, en het voldoen aan de preen postcondities bestaat uit het invoegen van communicatie-operaties. Als voorbeeld nemen we eerst weer de discretisatie van de luchtdrukterm uit subroutine trscue.f (zie paragraaf 2.4 voor de uitleg van berekeningen m.b.v. de irogeo-tabel, en paragraaf 2.5 voor de parallellisatie van dit stukje code m.b.v. de “ownercomputes” regel).

c c c

--- Calculation of atmospheric pressure gradient term (i.c. isvwp=1). --- Remark: no atmospheric pressure on boundaries if (isvwp.eq.1) then do 670 k=1,kmax do 660 irk=1,norows n =irogeo(1,irk) msta=irogeo(2,irk)-1


55


Cursusboek Module 3

ml =irogeo(3,irk) if (irobou(1,irk).ge.ibitfc) msta=msta+1 nmsta=iadres(n,msta) nml =iadres(n,ml) do 650 nm=nmsta,nml,minc nmu=nm+minc if (kfu(nm).eq.1 .and. kcs(nm).eq.1 .and. + kcs(nmu).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 650 continue 660 continue 670 continue end if

In dit stukje code wordt een loop uitgevoerd over alle Usnelheidspunten van het eigen gedeelte van het globale rekenrooster (n.l., m=mf..ml in sequentiële runs, mf is van een ander subdomein i.g.v. randcodes .ge.ibitfc). De variabelen kfu, brlsu, ddk en gvu worden alleen gelezen voor arrayindices i gelijk aan waarden van de loop-counter nm. Voor deze arrays geldt dus de pre-conditie in

: kfu, brlsu, ddk en gvu voor eigen U-snelheidspunten De arrays kcs en pres worden gelezen voor waterstandspunten (m,n) en (m+1,n), met (m,n) de logische coördinaten van de eigen U-snelheidspunten nm van de loop. Voor deze arrays kan nu de volgende pre-conditie worden opgeschreven:

in

: kcs, pres voor waterstandspunten die links of rechts naast een eigen U-snelheidspunt liggen De post-conditie voor de loop is dat array ddk is gewijzigd:

out

: ddk voor eigen U-snelheidspunten Bovenstaande preen post-condities geven per loop en per array aan welk gedeelte van het array er wordt gelezen (in) of geschreven (out). Voordat we verder gaan om de notatie te verbeteren en uit te breiden eerst een paar opmerkingen: 1. preen post-condities zijn in het algemeen krachtige hulpmiddelen om berekeningen van elkaar te ontkoppelen. Ze kunnen worden beschouwd als een interface-beschrijving van een stuk code, beschrijven het externe gedrag van de code zonder in te gaan op de interne werking ervan. 2. Net als interface-beschrijvingen mogen preen postcondities extra input/output-argumenten beschrijven die niet daadwerkelijk in de code worden gebruikt, ze mogen strengere eisen stellen dan feitelijk nodig is. In het bovenstaande voorbeeld is dit daadwerkelijk gebruikt door niet te letten op de nat/droog-status van roosterpunten (ddk wordt niet gelezen of geschreven


56


Cursusboek Module 3

voor droge punten) of op de intern/rand-status van waterstandspunten (pres wordt niet gelezen voor randpunten). De verbeteringen in de notatie die hieronder worden uitgewerkt zijn als volgt: De namen van de verschillende relevante sets van punten worden geformaliseerd. De verwijzingen geformaliseerd.

naar

links/rechts

worden

Er wordt een tussenniveau ingevoerd voor beschrijving van loops waaruit de preen post-condities kunnen worden afgeleid. De verschillende relevante verzamelingen van punten worden geformaliseerd aan de hand van “index sets”. Een index set is een algemene verzameling van elementen waarbij de elementen een uniek label hebben (de index). De index hoeft geen scalaire waarde te zijn, binnen WAQUA/TRIWAQ worden meestal de logische (m,n)-coördinaten van roosterpunten gebruikt. Voor de data-analyse zijn er een paar verzamelingen van punten met speciale relevantie, zie onderstaande tabel. “eigen U-snelheidspunten”

irogeo-tabel, msta..ml

“alle U-snelheidspunten”

irogeo-tabel, mf..ml

“eigen V-snelheidspunten”

icolge-tabel, nsta..nl

“alle V-snelheidspunten”

irocolge-tabel, nf..nl

“eigen waterstandspunten in X-richting”

irogeo-tabel, msta..mend

“alle waterstandspunten in X-richting”

irogeo-tabel, mf..mlu

“eigen waterstandspunten in Y-richting”

icolge-tabel, nsta..nend

“alle waterstandspunten in Y-richting”

icolge-tabel, nf..nlu

“alle interne waterstandspunten”

irogeo-tabel, mfu..ml

In de tabel is te zien dat punten worden geselecteerd op basis van verschillende kenmerken: eigen/alle, U/V/S, X/Y-richting. Daarnaast zijn er ook andere kenmerken van belang. Voor snelheidspunten is dit de nat/droog-status, waterstandspunten kunnen actief zijn of niet (kcs=1,2 of 0). Deze verschillende selectie-criteria worden gebruikt om symbolische namen voor de verschillende verzamelingen punten te vormen. Zo wordt de verzameling van eigen U-snelheidspunten verder als “own-up” aangeduid, en is de verzameling van alle actieve waterstandspunten in X-richting “all-actv-wlp-xdir”. Iedere verzameling correspondeert met een bepaalde loop en selectie met kenmerk-arrays, en voor iedere verzameling kan hiermee


57


Cursusboek Module 3

ook een masker-array imask(1:nmax,-2:mmax+3) worden gevuld. De verwijzingen naar punten in de omgeving van een ander punt worden geformaliseerd via “stencils”. Binnen de dataanalyse is alleen het patroon van interacties van belang, en worden er geen coëfficiënten (-1, -1, 4, -1, -1) aan de verschillende punten van het stencil gekoppeld. Het standaard vijfpuntsstencil wordt daarom geschreven als de verzameling van offsets stenc 1

( 1, 0), (0, 0), (1, 0), (0, 1), (0,1)

Door het gebruik van staggering in WAQUA/TRIWAQ is er een groot aantal verschillende stencils van belang. Bijvoorbeeld gebruikt de berekening van de luchtdrukgradiënt in Usnelheidspunten de waardes van array pres van de twee naastgelegen waterstandspunten (links/rechts), en door de gebruikte nummering van waterstands- en snelheidspunten zijn dit offsets links=(0,0), rechts=(1,0). De volgende tabel geeft een overzicht van een aantal veelgebruikte stencils. stenc1

( 1, 0), (0, 0), (1, 0), (0, 1), (0,1)

…

stenc1

( 1, 0), (0, 0), (1, 0)

nmd, nm, nmu

stenc+1

stencsu

(0, 0), (1, 0)

nm, nmu

stenc-1

stencus

( 1, 0), (0, 0)

nmd, nm

stenc1

(0, 1), (0, 0), (0,1)

ndm, nm, num

stencvu

(0, 0), (1, 0), (0, 1), (1, 1)

nm, nmu, ndm, ndmu

stencuv

(0, 0), ( 1, 0), (0,1), ( 1,1)

nm, nmd, num, numd

( 1, 0), (0, 0), (1, 0), (0, 2), (0, 1), (0,1), (0, 2)

…

stenc1

2

De tabel illustreert tegelijkertijd de gebruikte systematiek voor de naamgeving van stencils. De subscript verwijst naar het aantal punten dan wordt meegenomen in iedere richting (1) of in specifieke richtingen (+1 ), of naar de soort conversie die wordt uitgevoerd: van S- naar U-punten (su). De laatste kolom laat zien hoe de stencils vaak in de code zijn te herkennen, welke array-indices er meestal voor de offsets worden gebruikt. Het meest opvallende hierin is dat in de programmatuur de volgorde (m,n) is omgedraaid naar (n,m), zodat logische coördinaat (m+1,n) (offset (1,0)) wordt aangeduid met nmu en niet met mun. De derde wijziging in de notatie is de invoering van een aparte beschrijving voor loops. In die beschrijving worden


58


Cursusboek Module 3

achtereenvolgens de loop-range, de gewijzigde variabelen met offsets ten opzichte van de loop-index, en de gebruikte variabelen met offsets ten opzichte van de loop-index gespecificeerd. Voor de berekening van de luchtdrukterm is dit: own-up

loop modf.

ddk for stenc0

using

kfu, brlsu, ddk, gvu for stenc0, kcs, pres for stenc+1

Hierin staat stenc0 voor het stencil dat alleen offset (0,0) bevat, voor variabelen die alleen worden gelezen of gewijzigd voor de verschillende waarden van de loop-index. Uit deze beschrijving kunnen de preen post-condities voor de loop direct worden afgeleid. Dit is gebaseerd op het “vermenigvuldigen” van een index set met een stencil. Namelijk, array pres wordt gelezen voor alle array-indices die worden bereikt door een offset uit stencil stenc+1 op te tellen bij een loop-index uit index-set own-up. De nieuwe verzameling die met deze operatie wordt gecreëerd hangt sterk samen met het Cartesisch produkt van elementen uit de index set en het stencil, vandaar de naam “vermenigvuldigen”. Eigenlijk is het beter om te spreken over de “dilatatie” ofwel de “verschuiving” van een index set met een stencil. De resulterende index set wordt aangeduid met het symbool : “own-up stenc+1 ”. Tenslotte worden een aantal afkortingen ingevoerd ten behoeve van de data-analyse. Het gebruik van stenc0 wordt weggelaten bij gewijzigde variabelen (modf.). Voor uitvoervariabelen wordt de waarde stenc0 veel gebruikt, en als deze er niet bij wordt vermeld dan is dit de enige logische keuze. Variabelen die “nooit veranderen” worden meestal niet meegesleept in de hele data-analyse: kcs, guu, gvu, etc. De niet-eigen roosterpunten die worden bereikt door toepassing van een stencil in alle punten van een indexset worden aangeduid met de guard band ten opzichte van dat stencil:

guard1 ( p) = own - up

stenc 1

own - up

Hiermee zijn de notaties voor de grondige data-analyse allemaal beschreven. De toepassing ervan vergt enige oefening aan de hand van de beschikbare voorbeelden. Hiervoor kan bijvoorbeeld een routine van parallel TRIWAQ worden genomen en naast de complete data-analyse uit de technische documentatie van TRIWAQ worden gehouden. Ook wordt later in dit document een compleet voorbeeld uitgewerkt, namelijk de parallellisatie van een stuk nieuwe functionaliteit, betreffende het horizontaal k- turbulentie-model. Tenslotte Versie 1.0, oktober 2002

59


Cursusboek Module 3

geven we in deze paragraaf nog twee voorbeelden van de toepassing van de notatie voor twee moeilijke loops: een loop die variabelen wijzigt m.b.v. een stencil, en een loop die wordt geparallelliseerd m.b.v. redundante berekening. In de meeste loops in de rekenroutines van WAQUA/TRIWAQ is de loop-index ook de index die wordt gebruikt voor de resultaat-arrays. Dit soort loops is van de vorm “voor alle indices i bereken x(i)”. Er zijn een paar situaties wanneer er van dit schema wordt afgeweken. Met name bij de berekening van fluxen voor het discretiseren van de transportvergelijkingen. De fluxen leven in snelheidspunten, terwijl de coëfficiënten van de transportvergelijking met waterstandspunten worden geassocieerd. In dat geval is het volgende schema praktisch:

voor alle U-snelheidspunten doe bereken hulpvariabelen: flux tel op bij coefficienten voor waterstandspunten links/rechts Een loop waarin dit schema wordt gebruikt is de berekening van de advectieve en diffusieve transporten in de y-richting (Vpunten, onder/boven i.p.v. links/rechts) in subroutine wasdfc.f van WAQUA. c-----------------------------------------------------------------c --- for all active points is computed : c - explicit advective and diffusive transport (contut and vipat) c - contributions are added to wl-points top and bottom c transports are set to zero i.c. of dry velocity points, c such that temporarily dry points keep their old c concentration c-----------------------------------------------------------------c do 3000 irk=1,nocols m =icolge(1,irk) nf =icolge(2,irk)-1 nl =icolge(3,irk) c c --- loop over all v-points (using computational columns) c do 2900 n=nf,nl nm =iadres(n,m) num =nm+ninc nm =lgrid(nm) num =lgrid(num) c c --- compute dp/dy and d2p/dy2 in v-point (m,n), add c contributions to wl-points (m,n) and (m,n+1) c (coefficient dl) c if ( khv(nm).gt.0 ) then qyt =hdt*qy(nm) contut=hdt*qy(nm)*.5*(rp(nm,l)+rp(num,l)) chezy =sqrt(chzfac/czv(nm)) vipat =( cdag*abs(qyt)/(r*chezy*gvv(nm)) + + .5*(difco(nm)+difco(num)) )*rdx*hv(nm)* + gvv(nm)/gev(nm) * (rp(num,l)-rp(nm,l)) dl(nm) =dl(nm)-contut+vipat


60


Cursusboek Module 3

dl(num)=dl(num)+contut-vipat endif 2900 continue 3000 continue

De beschrijving van deze loop ten behoeve van de data-analyse is als volgt: all-wet-vp

loop modf.

dl for stenc+1(

using

khv, qy, czv, hv for stenc0, rp, difco for stenc+1(

Deze beschrijving is eenvoudig te verifiëren. Het is iets moeilijker om de bijbehorende preen post-condities te destilleren. Dit hangt op de interpretatie van de verschuiving “all-wet-vp ( stenc+1(”. Wanneer echter de betekenissen van all-wet-vp en stenc+1( erbij worden gehaald (resp. icolge, nf..nl, khv>0, en {(0,0), (0,1)}), dan is duidelijk dat het resultaat van de verschuiving kan worden ge-enumereerd met een loop over icolge, punten nf..nlu. Dit is index set all-wlpydir, behalve dat hierin de nat/droog-status van V-punten onder tafel is geschoven. Ten behoeve van de communicatie in parallelle runs is het praktisch om details van droogvallen en onderlopen niet weer te geven in de preen post-condities. Dit is omdat het in de communicaties ondoenlijk is om onderscheid te maken tussen natte en droge punten: dat leidt tot tijdsafhankelijke verzamelingen van te communiceren punten, en tot iedere tijdstap herberekenen van de vereiste communicatie-tabellen. Dit kost een stuk meer tijd dan het in iedere tijdstap communiceren van waardes voor droge punten. De preen post-condities van de loop worden nu bepaald door de loop-indexset en de gebruikte stencils met elkaar te vermenigvuldigen en door het afsplitsen van de roosterpunten van het eigen subdomein: in

: khv, qy, czv, hv for own-vp ( guard-1((p), rp, difco for own-wlp-ydir

out

guard1((p)

: dl for own-wlp-ydir ( guard1((p) De pre-conditie (in) laat zien voor welke delen van de guard band de verschillende arrays moeten zijn gecommuniceerd of consistent zijn berekend. De post-conditie (out) laat zien voor welke punten er gerekend wordt. Hier valt op dat er ook voor punten van andere subdomeinen wordt gerekend, en verder is ook te zien dat er voor open randpunten van het domein aanpassingen aan dl worden gemaakt. Dit klopt niet met de randvoorwaarde die voor de transportvergelijking wordt gebruikt. Eigenlijk zou er steeds moeten worden gecontroleerd of indices nm en num van array dl niet verwijzen naar


61


Cursusboek Module 3

randpunten nf en nlu. Deze controles worden in subroutine wasdfc.f achterwege gelaten omdat de randvoorwaarden pas later worden opgesteld, waarbij bovenstaande waardes van dl worden overschreven. Het tweede voorbeeld betreft een stukje code waarin redundante berekening op een zodanige manier wordt toegepast dat dit leidt tot samengestelde stencils. Het betreft de bekende berekening van de luchtdruk-term voor de impulsvergelijking, uit subroutine trsumo.f van TRIWAQ. if (isvwp.eq.1) then do 1960 k=1,kmax do 1955 irk=1,norows n =irogeo(1,irk) msta=irogeo(2,irk) mend=irogeo(3,irk)-1 if (irobou(1,irk).eq.ibitfc) msta=msta-1 if (irobou(2,irk).eq.ibitfc) mend=mend+1 nmsta=iadres(n,msta) nmend=iadres(n,mend) do 1950 nm=nmsta,nmend,minc nmu=nm+minc if (kfu(nm).eq.1 .and. brlsu(nm,k).le.1e6) + ddk(nm,k)=ddk(nm,k) + hdtdw*(pres(nmu)-pres(nm))/gvu(nm) 1950 continue 1955 continue 1960 continue end if

De post-conditie voor heel subroutine trsumo.f luidt “de waardes van de coëfficiënten AAK, CCK en DDK van de gediscretiseerde impulsvergelijking per laag zijn bekend voor alle punten van index set all-up x all-layers. Aan deze postconditie wordt voldaan door alle termen van de impulsvergelijking redundant te berekenen voor de punten van all-up die niet tot het eigen subdomein behoren (guard-1 ). De luchtdrukterm wordt in het globale domein alleen berekend voor interne snelheidspunten. In parallelle runs moet deze term dus worden berekend voor index set all-int-up. In sequentiële runs is dit range mfu..mld, in parallelle runs kunnen hier de punten mf en ml bijkomen, àls dat inwendige snelheidspunten van het globale domein zijn. De beschrijving van de loop wordt aldus: all-int-up x all-layers

loop modf.

ddk

using

kfu, brlsu, ddk for stenc0, pres for stenc+1

De standaard manier om de pre-conditie hieruit af te leiden levert voor array pres de index-set “all-int-up stenc+1 ”. Voor welk gedeelte van de guard band moet pres nu worden gecommuniceerd? Dit valt te beantwoorden met enige manipulaties met index sets. De verzameling all-int-up is te schrijven als Versie 1.0, oktober 2002

62


(own-int-up

Cursusboek Module 3

guard-1 ) = (own-int-up

stenc-1 )

Hiermee kan het produkt als volgt worden geschreven: all-int-up

stenc+1 = (own-int-up

stenc-1 )

= own-int-up

(stenc-1

= own-int-up

stenc1 = own-int-up

stenc+1 =

stenc+1 ) = guard1 (p)

De pre-conditie van de getoonde loop wordt hiermee in

: pres for own-int-up

guard1 (p),

kfu, brlsu for own-int-up

guard-1 (p)

De strekking van dit voorbeeld is dat redundante berekening leidt tot producten van stencils, en dat die producten met manipulaties kunnen worden herleid tot nieuwe enkelvoudige stencils.


63


4

Cursusboek Module 3

Theorie en praktijk van communicatie met COCLIB De WAQUA/TRIWAQ rekenprocessen moeten met elkaar gegevens uitwisselen om bijvoorbeeld de waarden uit te wisselen die ze voor hun deel van het gebied berekend hebben en om met elkaar te controleren of convergentie bereikt is. Deze communicaties worden uitgevoerd met COCLIB, een bibliotheek van communicatie-operaties. Deze bibliotheek is zodanig geschreven dat de gebruiker niet te maken krijgt met het daadwerkelijk verzenden en ontvangen van gegevens: de gebruiker specificeert alleen welke gegevens hij beschikbaar heeft voor andere processen en welke gegevens hij wil ontvangen. De communicatiebibliotheek zorgt ervoor dat de benodigde gegevensuitwisseling plaatsvindt. Om te specificeren welke gegevens er beschikbaar zijn voor andere processen en welke gegevens er ontvangen moeten worden, zijn een aantal krachtige mechanismen ontwikkeld, waardoor de programmeur met weinig moeite heel nauwkeurig kan aangeven wat er gecommuniceerd moet worden. In dit hoofdstuk worden die mechanismen behandeld en wordt getoond hoe COCLIB gebruikt dient te worden.

4.1

Overzicht van COCLIB Voordat COCLIB gebruikt kan worden, dient ze eerst geïnitialiseerd te worden door het aanroepen van de routine cocini. Deze routine initialiseert de interne datastructuren van de bibliotheek maar geeft ook de informatie door die het proces normaal gesproken vanaf de command-line zou ontvangen. De tegenhanger van cocini is cocend, waarmee het gebruik van de communicatie bibliotheek wordt beëindigd. Naast cocini en cocend kent COCLIB in essentie een drietal soorten routines: 1. configuratieroutines, waarmee men een aantal zaken specifeert die bij de communicatie een rol spelen. 2. nabuur communicatie routines, waarmee een proces informatie uitwisselt met zijn directe buren. 3. globale communicatie routines, waarmee een proces informatie uitwisselt met alle andere processen. Met de configuratieroutines definieert men bijvoorbeeld indexsets en stencils, met behulp waarvan men later aan de nabuur-communicatieroutines nauwkeurig kan opgeven wat er gecommuniceerd moet worden. Verderop in dit hoofdstuk wordt hier nog uitgebreid bij stilgestaan. Bij WAQUA/TRIWAQ is het configureren van COCLIB een flink stuk programmatuur. Dit gebeurt echter eenmalig bij het opstarten


64


Cursusboek Module 3

van WAQUA/TRIWAQ en heeft als voordeel dat er tijdens het rekenproces weinig extra uitgezocht hoeft te worden, zodat de communicaties vlot kunnen verlopen. In de categorie nabuur-communicatieroutines valt momenteel slechts één routine: cocupd. Deze routine wordt gebruikt om een te specificeren stuk van de guardband up-to-date te brengen met de meest recente waarden die buurprocessen daarvoor uitgerekend hebben. Globale-communicatieroutines worden vooral gebruikt om te bepalen of er convergentie is bereikt. Alle processen bepalen apart van elkaar of er binnen hun subdomein convergentie is bereikt, en door middel van een globale communicatie wordt bepaald of er subdomeinen zijn waarvoor convergentie nog niet is bereikt (zie paragraaf 3.4).

4.2

COCLIB en COEXEC De communicatie-bibliotheek COCLIB is nauw verbonden met het master/executive-proces COEXEC en de twee kunnen dan ook niet los van elkaar gezien worden. Op het moment dat COEXEC een proces opstart, stuurt het meteen een opstartboodschap erachteraan. De opstartboodschap bevat algemene informatie over het aantal processen dat opgestart wordt, maar bevat daarnaast ook de informatie die het proces in een niet-parallelle run via de command-line binnen zou krijgen (naam van SDS-file, experiment, print uitvoerfile). Het proces dient te beginnen met het aanroepen van de routine cocini, die deze opstartboodschap opvangt. Behalve het opvangen van de opstartboodschap en het terugleveren van de input-argumenten, dient cocini ook om een aantal interne datastructuren van COCLIB te initialiseren. Als tegenhanger van cocini dient op het eind van een proces cocend aangeroepen te worden. Deze routine stuurt een boodschap naar COEXEC waarin gemeld wordt dat het proces gaat termineren. COEXEC weet dan dat het proces op een ordelijke manier tot een einde komt. Als COEXEC op zeker moment constateert dat een proces verdwenen is, waarvan hij geen afmelding gekregen heeft, dan zal hij besluiten dat dat proces blijkbaar gecrasht is, wat een reden kan zijn om de hele parallelle run af te breken. In de communicatieroutines is ook een stuk functionaliteit ingebouwd om eventuele problemen die bij het communiceren ontstaan op te vangen. Op het moment dat een COCLIB merkt dat een communicatie niet vlot genoeg verloopt, wordt er contact gezocht met COEXEC om na te gaan wat de oorzaak van de problemen kan zijn. Zoals gezegd houdt COEXEC bij welke processen er


65


Cursusboek Module 3

actief zijn. Als de oorzaak van het communicatieprobleem gelegen is in het feit dat het proces waarmee gecommuniceerd moet worden er niet (meer) is, dan kan COEXEC dat constateren en terugmelden aan het proces dat de communicatie initieerde. Als een COCLIB routine een probleem gemeld heeft bij COEXEC dan blijft de routine in principe proberen om de communicatie alsnog tot stand te brengen, tenzij COEXEC meldt dat dat waarschijnlijk niet gaat lukken. Als de communicatie uiteindelijk toch lukt, dan meldt de COCLIB routine dat weer aan COEXEC, zodat die weet dat het probleem is opgelost. COEXEC houdt ook bij hoe vaak bepaalde communicatie problemen gemeld worden en kan op een zeker moment besluiten dat de hoeveelheid problemen dermate groot is geworden dat het niet zinvol lijkt om de run door te zetten. Op dat moment stuurt COEXEC een panic-boodschap naar alle processen. Zodra een proces constateerd dat een communicatie traag verloopt, dan wordt als eerste gekeken of COEXEC wellicht al een panic-boodschap heeft gestuurd. Als dat inderdaad het geval is dan wordt het proces binnen COCLIB afgebroken.

4.3

Indexsets en Stencils In COCLIB wordt heel nauwkeurig gespecificeerd wat er gecommuniceerd moet worden, zodat er nooit teveel informatie wordt uitgewisseld en de performance dus optimaal blijft. Dit specificeren van wat er gecommuniceerd moet worden, gebeurt met behulp van indexsets en stencils. In deze paragraaf zullen deze begrippen eerst nader toegelicht worden.

4.3.1

Indexsets Elk data-array in een programma heeft een zogeheten indexset, waarmee op het meest elementaire niveau de verzameling van zijn toegestane indices wordt aangeduid. Een array van 100 elementen heeft als indexset [1,100]. Een dergelijke indexset kan een bepaalde interpretatie kennen. De elementen van de indexset [1,mnmaxk] bijvoorbeeld zijn elk gekoppeld aan één van de rekenpunten in het WAQUA/TRIWAQ rooster. De elementen van de indexset [1,nsrc] zijn elk gekoppeld aan één van de bronpunten. Het is met name deze interpretatie die de ene indexset van de andere onderscheidt. Stel dat men een array met lengte 10 heeft waarin de waterstanden in de 10 meetpunten liggen opgeslagen en een array met lengte 10 waarin de waterstanden bij de 10 barrierpunten liggen opgeslagen, dan hebben deze twee arrays toch verschillende indexsets, hoewel ze dezelfde lengte hebben.


66


Cursusboek Module 3

Anderzijds, als men een ander array met lengte 10 heeft, waarin de U-snelheid bij elke barrier wordt opgeslagen, dan heeft dit array dezelfde indexset als het array waarin de waterstanden bij de barriers zijn opgeslagen. De indexsets die tot nu toe aan de orde gekomen zijn, hebben allemaal een ruimtelijke interpretatie. Dat wil zeggen: de elementen van de indexsets zijn steeds gekoppeld aan bepaalde locaties in het rekenrooster. Daarnaast kent men uiteraard ook indexsets die geen ruimtelijke interpretatie kennen, bijvoorbeeld de indexset [1,LMAX], waarvan elk element op een bepaalde opgeloste stof betrekking heeft. Tenslotte zijn er nog arrays waarvan de indexset geen enkele zinvolle interpretatie heeft, waarbij bijvoorbeeld gedacht kan worden aan het array idimen. Bij een indexset behoren coördinaten en een eigenaarschap. Deze twee attributen van een indexset komen nu eerst aan de orde.

4.3.2

Coördinaten van indices Voor het huidige betoog zijn uitsluitend die indexsets van belang die een ruimtelijke interpretatie hebben. Gezien de ruimtelijke interpretatie van dergelijke indexsets is het mogelijk om aan elk van de elementen van de indexset logische coördinaten toe te kennen die de plek van het element in het rekenrooster weergeven (zie ook Figuur 4). Voor de indexset {1..mmax} {1..nmax} is dit heel eenvoudig: de index van het element zelf (d.i. (m,n)) is zelf de coördinaat. Een dergelijke indexset, waarvan de indices zelf direct ook de coördinaten zijn, wordt een regelmatige indexset genoemd. Voor de indexset {1..mnmaxk} zijn de bijbehorende (m,n)coördinaten te vinden via de lgrid-tabel (eigenlijk is de inverse nodig, het opzoeken van (m,n) bij gegeven nm is nu veel werk). Voor de indexset {1..NSRC} zijn de coördinaten te vinden in de arrays mdis en ndis. Dergelijke indexsets, waarbij de coördinaten middels een tabel opgezocht moeten worden, worden onregelmatige indexsets genoemd. De reden waarom het nuttig is om deze coördinaten toe te kennen komt later aan de orde. Nu volstaat het om op te merken dat het te maken heeft met het feit dat WAQUA/TRIWAQ gebaseerd is op eindige-differentieberekeningen, waarin bijvoorbeeld de waarde in een roosterpunt berekend wordt uit de vier omliggende roosterpunten. De informatie in arrays wordt dus vaak gebruikt op basis van de met de elementen geassocieerde coördinaten.


67


Cursusboek Module 3

Title: ixset.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 4 Betekenis van coördinaten en eigenaarschap van indices

4.3.3

Het eigenaarschap van indices In parallel WAQUA/TRIWAQ kent het rekenrooster interne punten en guardband punten. Guardbandpunten zijn punten van het rekenrooster waarvan het deeldomein geen eigenaar is. Normaalgesproken zal het er niet rekenen maar het kan er wel variabelen lezen. Via zijn coördinaten kan een element van een indexset geassocieerd zijn met een intern punt of met een guardband punt. Daarmee wordt dan tevens het eigenaarschap van dat element bepaald (zie ook Figuur 4). Het deeldomein is ‘eigenaar’ van die indices die geassocieerd zijn met interne punten. De indices die geassocieerd zijn met punten in de guardband hebben, zijn eigendom van de eigenaar van dat guardbandpunt. Hierbij is het van belang om rekening te houden met het feit dat in WAQUA/TRIWAQ een gestaggerd rooster wordt gebruikt, waardoor het eigenaarschap van de U-punten van het rekenrooster bijvoorbeeld kan verschillen van het eigenaarschap van dieptepunten. De indexset van een array met lengte mnmaxk waarin waarden voor U-punten staan (bijvoorbeeld de U-snelheden) kan dus een ander eigenaarschap hebben dan de indexset van een array waarin waarden voor dieptepunten staan, ook al heeft dat laatste array ook een lengte mnmaxk. In dit geval is er sprake van twee aparte indexsets met


68


Cursusboek Module 3

Title: stencil.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 5: De array elementen die ververst moet worden voor het gegeven stencil

elk mnmaxk elementen. Dit detail is niet van belang voor parallel WAQUA/TRIWAQ maar speelt wel een rol bij domein decompositie.

4.3.4

Stencils In numerieke berekeningen kent men het begrip stencil. Stel dat men voor de berekening van een nieuwe waterstandswaarde in punt (m,n) de waterstanden nodig heeft uit punten (m-1,n), (m+1,n), (m,n), (m,n-1) en (m,n+1). In dat geval wordt door COCLIB onder het begrip stencil de verzameling offsets {(1,0),(0,0),(1,0),(0,-1),(0,1)} verstaan. In woorden: het stencil is de verzameling offsets (of offsetvectoren) waarmee men vanuit een punt (m,n) alle punten kan bereiken die in een berekening nodig zijn. Het stencil dat hier als voorbeeld wordt gebruikt, wordt gewoonlijk een vijfpunts stencil genoemd. In het vervolg zal ook gesproken worden over het toepassen van een stencil op een punt (m,n). Dit levert alle punten op die door middel van een offset uit het stencil bereikt worden vanuit het punt (m,n). Toepassen van het bovengenoemde vijfpunts stencil op het roosterpunt (6,3) levert dus de punten (5,3), (6,3), (7,3), (6,2) en (6,4) op, voorzover die bestaan.


69


Cursusboek Module 3

In het betoog tot nu toe is het stencil toegepast op de punten van het rekenrooster. Zoals boven beschreven is, hebben veel indexsets een relatie met het rekenrooster middels de coördinaten die met de elementen van de indexset zijn geassocieerd. Daardoor kan men ook een stencil toepassen op een willekeurige indexset (mits die een ruimtelijke interpretatie kent). Het toepassen van een stencil op een element van een indexset levert al die punten van dezelfde indexset op waarvan de coördinaten bereikt worden door het stencil toe te passen op de coördinaten van het betreffende element. Overigens kan men ook een stencil toepassen op een indexset en het resultaat daarvan vertalen naar een andere indexset. Stel dat men een stencil toepast op ixset1 en het resultaat vertaalt naar ixset2, dan is het resultaat daarvan alle punten in ixset2 waarvan de coördinaten worden bereikt door het stencil toe te passen op de coördinaten van de elementen van ixset1.

4.3.5

Interfaces Als men nu wil weten welke punten (m,n) er allemaal gebruikt worden in een berekening met het vijfpuntsstencil waarbij nieuwe waarden berekend worden voor alle interne punten van subdomein 1, dan kan men het stencil toepassen op alle interne punten van subdomein 1. Dat levert een set van punten op waarin ook punten voorkomen die geen interne punten van subdomein 1 zijn (zie Figuur 5). Dit zijn punten die nodig zijn voor de berekening van nieuwe waarden vlak bij de grens tussen twee subdomeinen. Deze punten die geen eigendom zijn van subdomein 1 worden aangeduid als de externe interface van subdomein 1 voor de berekening met het vijfpuntsstencil. Er bestaat ook een interne interface van subdomein 1; daarmee wordt de set van punten in subdomein 1 bedoeld die subdomein 2 nodig heeft voor de berekening met het vijfpuntsstencil. In het voorbeeld van Figuur 5 is de externe interface van subdomein 1 gelijk aan de interne interface van subdomein 2 en omgekeerd. De interne en externe interface worden samen aangeduid als de interface voor de berekening. In Figuur 5 wordt ook getoond dat men ook kan spreken over een interface van een bepaalde indexset. Die interface wordt weer verkregen als alle niet-eigen indices die worden bereikt door het stencil toe te passen op elk van de eigen indices. Merk op dat hetzelfde stencil, toegepast op twee verschillende indexsets ook twee verschillende interfaces oplevert.

4.3.6

Interfaces verversen met COCLIB Voorafgaand aan de berekening met het vijfpuntsstencil, moet subdomein 1 ervoor zorgen dat het de meest recente waarden (bijv. waterstanden) verkrijgt in zijn externe interface. Dat wil zeggen: de waarden in de externe interface punten (in de guard


70


Cursusboek Module 3

band) moeten ververst worden (tenzij men al zeker weet dat de meest recente waarden al aanwezig zijn). Dit verversen wordt in COCLIB aangeduid met de Engelse term update. COCLIB biedt voor dit verversen één enkele, zeer krachtige operatie: cocupd. De invoer voor deze operatie bestaat uit een array waarvan waarden ververst moeten worden, een specificatie van de indexset van die array en een specificatie van de interface waarvoor het verversen moet gebeuren. Deze operatie stuurt alle waarden uit de interne interface naar de buurdomeinen die ze nodig hebben en ontvangt van alle buurdomeinen nieuwe waarden voor de punten in de externe interface. Het argument waarmee men in cocupd de indexset van het array specificeert is een handle, die verkregen wordt uit één van de operaties waarmee men een indexset kan definiëren. Elke indexset moet dus van te voren gedefinieerd worden. Dat gebeurt tijdens de initialisatie van WAQUA/TRIWAQ op het moment dat alle daarvoor benodigde informatie is ingelezen. Verderop in dit hoofdstuk wordt daar nog nader op ingegaan. Zo ook is het argument waarmee men in de cocupd aanroep de interface specificeert een handle naar een vooraf gedefiniëerde interface. Het definiëren van de interfaces gebeurt ook weer bij het initialiseren van WAQUA/TRIWAQ. Een interface wordt doorgaans aangeduid met een naam die verwijst naar het stencil waarmee de interface bepaald is. In slordig spraakgebruik wordt de interface soms ook wel aangeduid als het stencil. Ook het definiëren van interfaces met behulp van stencils komt verderop in dit hoofdstuk nog aan de orde.

4.4

De COCLIB routines Een volledig overzicht van de routines in COCLIB is te vinden in de “Programmer’s Guide COCLIB”. Deze paragraaf geeft een overzicht van de voornaamste routines plus wat voorbeelden van het gebruik ervan, toegespitst op parallel rekenen. De voorzieningen voor domein decompositie komen in een later hoofdstuk aan de orde.

4.4.1

Starten, stoppen en command-line argumenten De eerste routine van COCLIB die in een programma moet worden aangeroepen is cocini. Deze initialiseert een aantal interne datastructuren van COCLIB en ontvangt bovendien de opstartboodschap van COEXEC, waarin ondermeer de argumenten zijn opgenomen die het programma in een nietparallelle omgeving via de commandline binnen zou krijgen. Deze argumenten worden teruggeleverd in een array van strings. De elementen in dit array kunnen worden opgevraagd


71


Cursusboek Module 3

met de routine cocarg, die daarbij ook direct het type van het element teruggeeft (keyword, string, integer of real). Deze routine correspondeert met de routine SIRDEL die in een nietparallelle omgeving voor hetzelfde doel gebruikt wordt. Het aanroepen van cocend beëindigt het gebruik van COCLIB in een programma. Deze routine stuurt een boodschap naar COEXEC om te melden dat het proces verder geen communicatie zal doen. Deze routine moet helemaal aan het eind van het programma (of vlak voor een exit()) aangeroepen worden omdat er op sommige platforms vanuit gegaan wordt dat het programma direct na het aanroepen van cocend zal termineren. Op die platforms is het ook zo dat cocend pas eindigt op het moment dat alle parallelle processen cocend hebben aangeroepen.

4.4.2

Opslagruimte voor interne tabellen COCLIB maakt op grote schaal gebruik van allerhande interne tabellen: indexsets en interfaces wordt opgeslagen als een tabel, en ten behoeve van domein decompositie worden er ook diverse tabellen opgeslagen die gebruikt worden bij de interpolaties. Hiervoor is flink wat opslagruimte nodig, waarvan de grootte uiteraard samenhangt met de grootte van het model dat doorgerekend wordt. Om te voorkomen dat ook kleine modellen met een heel grote interne opslagruimte moeten werken, is het wenselijk om deze opslagruimte dynamisch te laten creëren. COCLIB maakt daarom gebruik van een groot werkarray dat de applicatie die COCLIB gebruikt zelf moet aanmaken en die steeds meegegeven moet worden bij de aanroep van COCLIB routines. Dit werkarray staat bekend onder de naam icocad. COCLIB biedt wel een routine waarmee uitgerekend kan worden hoelang dit array moet zijn voor WAQUA/TRIWAQ op basis van een aantal relevante afmetingen van het model (mmax, nmax, kmax). Na het aanroepen van cocini zal dus vrij snel eerst cocmem aangeroepen worden om na te gaan hoelang het array icocad moet zijn, waarna het array aangemaakt kan worden. Pas daarna kunnen de meeste overige COCLIB routines gebruikt worden.

4.4.3

Configuratie en gebruik van nabuurcommunicaties met COCLIB Voor het verversen van waarden in de guardband biedt COCLIB de routine cocupd. Met deze routine kan men de guardbandwaarden in maximaal vier (integer of real) arrays tegelijkertijd laten verversen. Het is voordelig om zoveel mogelijk arrays tegelijk te laten verversen omdat het altijd voordeliger is om één keer een lange boodschap te sturen dan


72


Cursusboek Module 3

viermaal een boodschap met een kwart van de lengte (vanwege de latency). cocupd heeft de volgende heading: subroutine cocupd(numint, numrea, ifld1, setnm1, itfnm1, ifld2, setnm2, itfnm2, ifld3, setnm3, itfnm3, ifld4, setnm4, itfnm4, icocad, name)

intpl1, intpl2, intpl3, intpl4,

De eerste twee argumenten geven aan hoeveel integer arrays er ververst moeten worden en hoeveel real arrays. Daarna volgen de arrays zelf met hun bijbehorende argumenten. Als eerste komen de numint integer arrays, daarna de numrea real arrays. Voor elk van de arrays moet men opgeven wat de indexset van het array is, voor welke interface de guardband ververst moet worden en welke interpolatiemethode er gebruikt moet worden. Dat laatste, de interpolatiemethode, komt nu nog niet aan de orde. Wel zal worden besproken hoe de indexset en de interface worden gespecificeerd. Beiden worden aan cocupd opgegeven door middel van hun naam. Een aanroep van cocupd kan er dus uitzien als: Call cocupd(1,0, kcs,’fullbox’,’stcmax’,’ 0 ,’ ’ ,’ ’ ,’ 0 ,’ ’ ,’ ’ ,’ 0 ,’ ’ ,’ ’ ,’ icocad,name1(1:ilen))

’, ’, ’, ’,

waarmee dus de guardbandwaarden van het array kcs worden vervangen. De indexset van dit array is de indexset met de naam ‘fullbox’ en de guardband moet vervangen worden voorzover de punten vallen binnen de interface ‘stcmax’. Zowel de indexset als de interface moet vooraf gedefinieerd zijn, anders volgt een foutmelding.

4.4.4

Het definiëren van indexsets Indexsets kunnen gedefinieerd worden met de routines cocidr, cocidi en cocprd. De eerste twee routines dienen om respectievelijk een regelmatige en een onregelmatige indexset te definiëren. De laatste routine (cocprd) dient om het cartesisch produkt van twee eerder gedefiniëerde indexsets te definiëren. Zoals eerder al werd uitgelegd, is een regelmatige indexset een indexset waarvan de elementen zelf direct hun coördinaten vormen. Een bekend voorbeeld is de indexset “fullbox”, met grootte [1..nmax]x[-2..mmax+3]. Bij het definiëren van een regelmatige indexset hoeft men alleen de twee afmetingen mee te geven en een ownership array, d.i. een array met de afmetingen van de indexset waarin voor elk indexpaar is


73


Cursusboek Module 3

aangegeven welk proces de eigenaar van het element is. De punten waarvan de eigenaar het proces zelf is, zijn de interne punten. Punten waarvan andere processen eigenaar zijn, zijn guardband punten. Voor de indexset 'fullbox', is dit in principe het array POWNER dat door COPPRE op de deeldomein SDS-files wordt gezet. Het kan echter zijn dat men onderscheid wil maken tussen de fullbox voor v-punten en voor u-punten. In dat geval kan het ownership array voor de fullbox van v-punten iets verschillen van dat voor de fullbox van upunten. Het volgende code fragment toont hoe een proces met procesnummer 1 een zeer eenvoudige regelmatige indexset van 3 bij 3 elementen wordt gedefiniëerd:

+ +

+

integer idmmin(2), idmmax(2), iowner(3,3) data owner/2, 2, 2, 3, 1, 4, 3, 3, 3/ data idmmin/1,1/ data idmmax/3,3/ call cocidr(‘ix_exmpl’,2,idmmin,idmmax,iowner, icocad,name(1:ilen))

Merk op dat het eigenaarschaparray aangeeft dat alleen het middenpunt daadwerkelijk toebehoort aan het proces en dat deoverige punten vooral toebehoren aan processen 2 en 3 (en een enkel punt nog aan proces 4). Bovendien zijn de indices van deze indexset eenvoudigweg de indices 1 tot en met 3. Als de indices in de eerste dimensie van –2 tot 0 zouden lopen, dan zou idmmin(1)=-2 zijn en idmmax(1)=0. Een onregelmatige indexset is een indexset waarbij de coördinaten die met elke index geassocieerd zijn, apart moeten worden opgegeven. Een dergelijke indexset is in principe 1dimensionaal (d.i. de indexset heeft de vorm {1..card}). De coördinaten die met elke index geassocieerd zijn, kunnen echter meerdimensionaal zijn. Meestal zijn de coördinaten tweedimensionaal (ze beelden af op de set {1..mmax} {1..nmax}). De routine cocidi, waarmee een onregelmatige indexset gedefinieerd wordt, vraagt behalve de lengte van de indexset ook de dimensie van de coördinaten en de coördinaten zelf. Daarnaast moet ook nu weer het eigenaarschap array meegegeven worden. Stel dat men de regelmatige indexset die zojuist gebruikt is als voorbeeld van het definiëren van een regelmatige indexset nu als onregelmatige indexset wil definiëren (wat overigens niet erg zinvol zou zijn). Dan zou het code fragment er als volgt uit komen te zien.

+ +


integer icoor(2,9), iowner(9) data icoor/1, 1, 1, 2, 1, 3,

74


+ + + + + +

+

Cursusboek Module 3

2,1, 2,2, 2,3, 3,1, 3,2, 3,3/ data ionwer/2,2,2,3,1,4,2,2,2/ call cocidi(‘ix_exmpl’,9,2,icoor,iowner, icocad,name(1:ilen))

Met cocprd kan men het Cartesisch produkt van twee eerder gedefinieerde indexsets definiëren. Men zou kunnen denken dat de indexset ‘fullbox’ ook als Cartesisch produkt gedefinieerd kan worden. Dat kan evenwel niet omdat bij cocprd geen nieuw ownership array meegegeven kan worden. Bovendien wordt in cocprd gecontroleerd of ten hoogste één van de indexsets punten bevat die geen eigendom zijn van het aanroepende proces zelf (omdat anders het eigenaarschap van de elementen van het produkt dubbelzinnig zou worden). COCPRD wordt vooral gebruikt voor het definiëren van 3D indexsets op basis van 2D indexsets; bijvoorbeeld de indexset ‘fullbox*kmax’ op basis van ‘fullbox’ en ‘kmax’.

4.4.5

Het definiëren van interfaces Nadat de indexsets gedefiniëerd zijn, kunnen de interfaces gedefiniëerd worden. Voor elk stencil waarvoor men in het programma een nabuur-communicatie wil doen, moet een interface gedefiniëerd worden. Het definiëren van een interface gebeurt met de routine cocitf. De routine is zeer krachtig: hij kan gebruikt worden om zeer complexe interfaces te definiëren. In de meest eenvoudige vorm hoeft men slechts de naam van een indexset, een beschrijving van het stencil en de naam van de te creëren interface op te geven. Bijvoorbeeld:

+ + + + + +

integer mask(3,3), ioffs(2,5) data mask/9*1/ data ioffs/-1,0, 0,-1, 1,0, 0,1, 1,1/ call cocitf(‘ix_exmpl’,mask,5,ioffs,2, ‘ix_exmpl’,mask,’ ‘,’ ‘,’ ‘, ‘stc_5’,icocad,name(1:ilen))

In dit voorbeeld wordt een interface aangemaakt voor een vijfpunts stencil. Merk overigens op dat de naam van de interface (‘stc_5’) refereert aan het stencil dat gebruikt is. Dit is conventie in de huidige versie van parallel WAQUA/TRIWAQ en is ook prettig voor het begrijpen van de code, hoewel een interface feitelijk iets anders is dan een stencil.


75


Cursusboek Module 3

Het array mask en de argumenten in de middelste regel van de call naar cocitf lijken op het eerste gezicht overbodig, maar dienen om meer complexe situaties aan te kunnen. Daarop zal nu verder worden ingegaan. De externe interface die met cocitf bepaald wordt, wordt gevonden door het volgede algoritme: for alle elementen i1 in set1 if owner1(i1)=me and mask1(i1)==1 for elke offsetvector v van het stencil einde = i1 + v if einde in set2 and owner2(einde)= not me and if mask2(einde)==1 einde behoort tot externe interface endif endif endfor endif endfor

waarin set1 en mask1 de eerste twee argumenten van cocitf zijn, en set2 en mask2 verderop in de argumentenlijst voorkomen, na de specificatie van het stencil. De owner arrays horen bij de indexset en worden bij het definiëren ervan meegegeven (zie boven). Voor de interne interface geldt m.m. hetzelfde algoritme. In het zojuist gegeven voorbeeld waren set1 en set2 aan elkaar gelijk en was mask in beide gevallen overal 1. Dit zal echter in het algemeen niet zo zijn. Het array mask1 heeft dezelfde afmetingen als indexset set1 waar hij bij hoort en dient om ervoor te zorgen dat het stencil niet in elk punt van de indexset wordt toegepast, maar slechts in een beperkt aantal punten (te weten: die punten waarvoor mask de waarde 1 heeft). Een bekend voorbeeld is bijvoorbeeld een Red-Black Jacobiiteratie methode. Daarbij worden ombeurten nieuwe waarden berekend in rode en in zwarte punten, waarbij de zwarte en rode punten een schaakbordpatroon vormen. Als men een nieuwe waarde wil berekenen in een rood punt, heeft men daarvoor alleen de waarden nodig die rondom rode punten liggen. In dit geval zal men voorafgaand aan het herberekenen van de rode punten de guardband willen verversen voor een interface die gevonden wordt door het stencil toe te passen op alle rode punten. Dit kan men doen door bij het definiëren van deze interface een masker mee te geven waarin alleen de rode elementen de waarde 1 hebben. De tweede indexset, set2, die men aan COCITF mee kan geven, zal in veel gevallen dezelfde zijn als set1. In bepaalde gevallen kan het echter nodig zijn om hiervan af te wijken. Stel dat men bijvoorbeeld een nieuwe waarde wil berekenen voor


76


Cursusboek Module 3

een attribuut van een barrier en dat men daarbij waterstanden nodig heeft die om de barrier heen liggen. In geval wil men voorafgaand aan de berekening van waterstandsveld alleen die guardband punten verversen dicht bij een barrier liggen.

de dat het die

Dit kan met cocitf bereikt worden door voor set1 de indexset van alle barriers te nemen en voor set2 de set van alle waterstandspunten. De interface die dan bepaald wordt bestaat uit alle waterstandspunten in de guardband die bereikt worden door het stencil toe te passen op een barrierpunt. Deze interface kan gebruikt worden om het waterstandsveld te updaten. De overige argumenten op de tweede regel van de aanroep van cocitf in bovenstaand voorbeeld (de lege strings) zijn nodig in het geval van horizontale verfijning en zullen in een volgend hoofdstuk aan de orde komen. Voor parallel rekenen en domein decompositie met vertikale verfijning zijn ze niet nodig. In het bovenstaande is aangegeven hoe een interface voor een 2D indexset aangemaakt kan worden. In veel gevallen (in TRIWAQ) zal men echter een interface nodig hebben voor een 3D indexset (het produkt van een 2D indexset met bijvoorbeeld de indexset kmax). Uitgaande van een interface voor een 2D indexset kan een interface voor een 3D indexset aangemaakt worden met de routine coctpr. In feite wordt hiermee geen nieuwe interface aangemaakt, maar wordt aangegeven dat de interface ook gebruikt mag worden voor de 3D indexset.

4.5

Gebruik van globale communicaties Naast de voorzieningen voor nabuurcommunicaties biedt COCLIB ook voorzieningen voor zogeheten globale communicaties. Dit zijn communicaties waarbij alle processen (of een bepaalde subgroep van alle processen) gezamenlijk informatie uitwisselen. Een globale communicatie wordt bijvoorbeeld gedaan als een alle processen een iteratie van een iteratief proces hebben volbracht en voor hun eigen deeldomein hebben vastgesteld of convergentie al dan niet bereikt is. Alle processen roepen dan een globale communicatieroutine aan, waarbij ze hun eigen convergentievlag als argument meegeven en als resultaat het minimum van de vlaggen van alle deelnemende processen krijgen. Als dit minimum aangeeft dat er een proces is waarvoor convergentie nog niet bereikt is, dan doen alle processen nog een nieuwe iteratie. In pseudo-code: if convergentie bereikt cnvflg = 1 else cnvflg = 0 endif call cocigl(cnvflg,1,COCMIN,’all’,info,


77


Cursusboek Module 3

name(1:ilen)) if(info==0) if(cnvflg == 1) doe nog een iteratie else error: communicatie mislukt endif

In bovenstaand voorbeeld geeft de constante COCMIN aan dat het minimum van de waarden van alle processen moet worden teruggegeven. Zo zijn er ook constanten COCMAX en COCSUM, met voor de hand liggende betekenissen. Met deze globale communicatieroutines kan men meer dan één enkele integer (of real) tegelijk communiceren. Als men in bovenstaand voorbeeld in plaats van cnvflg een array van waarden ingestop had, dan zou ieder element van dat array na afloop gevuld zijn met het minimimum van de overeenkomstige array-elementen van alle deelnemende processen. Dit werd in oudere versies van COCLIB ondermeer gebruikt om expliciet informatie te krijgen van elk van de deelnemende processen. Stel dat men weet dat er nproc processen deelnemen aan de globale communicatie. Dan kan het proces dat van zichzelf weet dat het proces iproc is een array toeleveren ter lengte van nproc waarin element iproc gevuld is met de waarde die het proces wil toeleveren en verder overal nul. Als alle processen ditzelfde doen, dan zal na een globale communicatie met COCSUM dit array gevuld zijn met een waarde voor elk proces. Overigens gaat dit gebruik enigzins in tegen de algemene filosofie van COCLIB, waarbij een proces zo min mogelijk moet afweten van de context waarin het draait (dus eigenlijk ook niets moet afweten van het aantal processen die er aktief zijn). Voor deze vorm van communicitie is de nieuwere operatie COCNWS eigenlijk meer geschikt. In bovenstaand voorbeeld is het argument ‘all’ opgenomen, wat betekent dat alle aktieve processen geacht worden aan de communicatie mee te doen. Dit betekent dan ook dat de communicatie pas eindigt als alle processen hun bijdrage hebben geleverd. In sommige gevallen is dat niet handig. Stel dat men bijvoorbeeld een globale communicatie wil doen om het debiet door alle stukken van een opening te verkrijgen. Dan wil men deze communicatie alleen uitvoeren met processen die ook daadwerkelijk een stuk van de opening hebben. Om dit mogelijk te maken is het concept van groepen geïntroduceerd. Met de routine COCDGR kan men een procesgroep definiëren die bestaat uit een subgroep van een andere groep. COCDGR is zelf ook een globale communicatie, waarbij alle processen die eraan deelnemen aan elkaar


78


Cursusboek Module 3

doorgeven of ze wel of niet tot een subgroep behoren. Bijvoorbeeld: if barriers in dit subdomein ihavbr = 1 else ihavbr = 0 endif call cocdgr(‘have_barriers’,ihavbr,’all’, name(1:ilen))

waarin alle aktieve processen met elkaar communiceren om te bepalen welke van hen barriers bevatten. Deze nieuwe groep (de groep ‘have_barriers’) kan vervolgens gebruikt worden om binnen die groep een globale communicatie te doen. In bovenstaand voorbeeld werd de routine COCIGL gebruikt voor een globale communicatie van integer waarden. Naast deze routine zijn er nog drie andere globale communicatierouitnes: COCRGL (analoog aan COCIGL maar dan voor real waarden), COCRGS en COCNWS. COCRGS lijkt in veel opzichten op COCRGL (en dus op COCIGL) met als enige uitzondering dat elk element van het array dat gecommuniceerd moet worden een tuple is van met elkaar samenhangende waarden. COCRGS ondersteunt slechts één reductie-operatie die aangegeven wordt met de constante COCMAX1. Deze operatie zorgt ervoor dat elk deelnemend proces de tuple terugkrijgt waarvan het eerste element het grootste is van alle toegeleverde tuples. Een voorbeeld waarbij deze operatie gebruikt wordt is een van de convergentiekriteria in WAQUA/TRIWAQ. Elk proces levert een tuple waarin achtereenvolgens het maximale residu, het eigen procesnummer en de (m,n)-locatie van het maximale residu zijn opgenomen. Na de globale communicatie met COCRGS bezit elk proces het tuple waarin het grootste residu voorkwam, met daarbij het procesnummer dat het grootste residu gevonden heeft en de locatie van het grootste residu. Tenslotte is er nog de operatie COCNWS. Deze operatie ondersteunt een ‘newspaper’ vorm van globale communicatie waarbij elk proces een aantal waarden (vgl. een artikel in een krant) toelevert en het resultaat bestaat uit een array waarin alle waarden die door alle processen geleverd zijn achterelkaar zijn opgenomen (vgl.de krant die uit artikelen bestaat). Als resultaat komt bovendien een array terug waarin te zien is welke waarden door welk proces zijn bijgedragen.


79


5

Uitbreidingen voor Domein Decompositie

5.1

Basisfilosifie

Cursusboek Module 3

Tot nu toe is uitsluitend aandacht besteed aan de voorzieningen voor parallel rekenen. Er is uitgebreid stilgestaan bij indexsets, stencils en interfaces. Met de introduktie van vertikale en horizontale verfijningen zijn deze concepten weliswaar gehandhaafd, maar ze hebben een uitgebreidere betekenis gekregen. Daarnaast zijn er enkele nieuwe concepten geïntroduceerd die vooral verband houden met de benodigde interpolaties. Deze paragraaf zal ingaan op de uitgebreidere betekenis die het begrip guardband gekregen heeft. De volgende paragraaf bespreekt de gevolgen daarvan voor het begrip van interfaces.

5.1.1

De noodzaak voor een guardband Bij de behandeling van dit concept in het vorige hoofdstuk is aangegeven dat een indexset de verzameling van indices is waarmee elementen van een array kunnen worden geadresseerd. De indexsets die voor COCLIB van belang zijn, hebben een ruimtelijke interpretatie: elke index is gekoppeld aan een logische coördinaat in het rekenrooster. Bij parallel rekenen hebben alle deeldomeinen een stuk van een enkelvoudig rekenrooster. Afgezien van de guardband overlappen de delen van de verschillende processen elkaar niet. Er zou ook geen guardband nodig zijn, als men de waarden die een proces nodig heeft uit een ander deeldomein rechtstreeks zou kunnen lezen uit de datastructuren van dat deeldomein. Er zijn echter twee complicaties waarom er toch een guardband nodig is. Ten eerste is de communicatie tussen twee processen te traag om iedere waarde apart over te halen, zodat men in één keer alle benodigde waarden wil communiceren. Daarvoor is tijdelijke opslag nodig binnen het gebruikende proces. Hiervoor dient de guardband. Ten tweede wil men bij een berekening alle gebruikte punten gewoon adresseren alsof ze in het eigen rooster liggen en niet per punt uitzoeken waar dat punt precies te vinden is. Bij parallel rekenen is de guardband een eenvoudige uitbreiding van het eigen rekenrooster, zodat men normale adressering kan gebruiken. De guardband dient dus zowel als tijdelijke opslag van opgehaalde waarden alsook om het adresseren van die waarden gemakkelijk te maken.


80


Cursusboek Module 3

Title: overlap.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 6: De overlapset bij de communicatie van grof naar fijn

5.1.2

De guardband bij verfijning Zodra roosters onderling verfijnd zijn, wordt de situatie wat complexer. De waarden die ontvangen worden uit een ander proces passen niet meer op het eigen rekenrooster en dus ook niet in de eigen arrays. Stel dat een proces waarden ontvangt van een buurproces met 6 lagen, terwijl het proces zelf 3 lagen heeft, dan zijn de arrays van dit proces zelfs met guardband niet toereikend om de gegevens van het buurproces op te slaan. Bovendien wil men niet steeds als er een waarde wordt gebruikt die eigenlijk bij een buurproces thuishoort, ter plekke de interpolatie naar het eigen rekenrooster uitvoeren.

5.1.3

De overlapset Dit probleem is opgelost door elk proces naast zijn eigen rekenrooster, uitgebreid met een guardband, nog een extra rooster (het overlaprooster) te geven dat bestaan uit die delen van het eigen rooster die overlappen met nabuurroosters (inclusief guardband) tezamen met die delen van nabuurroosters die overlappen met het eigen rooster (inclusief guardband). Er is een indexset geassocieerd met dit extra rooster en deze indexset heet de overlapset. Elk element van de overlapset is geassocieerd met een roosterpunt in het overlaprooster, wat op zijn beurt weer geassocieerd kan zijn met ofwel een punt van het eigen rekenrooster of een punt van een nabuurrekenrooster. Met de overlapset is ook een stuk geheugenruimte (een werkarray) verbonden, dat gebruikt wordt voor de communicaties tussen processen waarbij sprake is van verfijning.


81


Cursusboek Module 3

Communicaties waarbij interpolaties plaatsvinden, verlopen via de overlapset (zie Figuur 6). Het proces dat gegevens wil versturen, kopiëert die eerst naar de overlapset, vervolgens worden de gegevens van de overlapset gecommuniceerd naar het ontvangende proces, waar ze weer op de overlapset (van het ontvangende proces) geplaatst worden. Op de overlapset van het ontvangende proces vindt vervolgens de interpolatie plaats. Merk op dat hierbij ook waarden van het ontvangende proces gebruikt worden, die ook op de overlapset geplaatst zijn. Bovendien wordt de interpolatie pas uitgevoerd nadat van alle toeleverende processen de bijdragen zijn ontvangen. Het is dus mogelijk dat voor de berekening van een waarde in de guardband van een proces waarden worden gebruikt van het proces zelf en van meer dan één ander proces. Dit maakt het mogelijk om op geometrisch zeer complexe manieren te koppelen. Na de interpolatie worden de gegevens gekopiëerd van de overlapset naar (de guardband van) het eigen rekenrooster. De filosofie van COCLIB gaat er van uit dat een proces geen weet heeft van zijn context. Door middel van de update operatie ontvangt het nieuwe waarden in de guardband, maar het proces zou geen bemoeienis moeten hebben met de wijze waarop die waarden tot stand komen. Bij vertikale verfijning is deze filosofie onverkort gehandhaafd: de overlapset is geheel een interne aangelegenheid en is nergens zichtbaar voor de applicatie programmeur. Ook intern in COCLIB is de overlapset niet expliciet aanwezig (maar wel conceptueel). Bij horizontale verfijning is dit principe enigszins losgelaten, omdat het correct bepalen van de overlapset heel veel inhoudelijke kennis van WAQUA/TRIWAQ en van de gebruikte interpolatiemethoden vergt. Daarom wordt de overlapset in WAQUA/TRIWAQ met horizontale verfijning expliciet aangemaakt. Ook komt men de overlapset tegen in een aantal operaties (bijvoorbeeld bij het definiëren van interfaces en bij het definiëren van interpolatiemethoden).

5.1.4

Interpolatiepunten en discretisatiepunten Het begrip van eigenaarschap van punten in de overlapset en van punten in gewone indexsets is ook gewijzigd. Bij parallel rekenen kon men van ieder punt in een indexset eenduidig vaststellen wie de eigenaar was. Bij horizontale verfijning is dit niet meer zo eenvoudig. Er zijn immers punten in de indexset van het ene proces die eenvoudigweg niet voorkomen in de overeenkomstige indexset bij een ander proces. Welk proces is hiervan dan de eigenaar. De oplossing voor dit probleem is dat er binnen een indexset (behalve de overlapset) niet meer gesproken wordt over een eigenaarschap maar over interpolatiepunten en discretisatie-


82


Cursusboek Module 3

punten. Een interpolatiepunt is een punt waar waarden niet door het proces zelf verkregen worden maar worden vastgesteld door interpolatie van waarden van een ander proces. Een discretisatiepunt is een punt waarin het proces zelf daadwerkelijk rekent. Een indexset is conceptueel in zijn geheel eigendom van het proces (behoudens stukken waar het proces met een parallelle buurman gekoppeld is of met een buurman waarbij vertikale verfijning optreedt). Op de overlapset ligt dit net even complexer, omdat daarin punten van meerdere processen voorkomen, waaronder dus interpolatiepunten en discretisatiepunten van elk van die processen. Daarom speelt het eigenaarschap op de overlapset nog wel een rol. Een element van een overlapset is eigendom van een van de processen en is ofwel een interpolatiepunt ofwel en discretisatiepunt. Overigens biedt COCLIB de mogelijkheid om per indexset de grens tussen interpolatiepunten en discretisatiepunten per gebruik vast te leggen. Bij het definiëren van horizontale interpolaties (zie onder) moet men expliciet aangeven wat de interpolatiepunten zijn en wat de discretisatiepunten zijn.

5.2

Interpolatie in COCLIB Bij een update operatie kan men specificeren welke interpolatie er moet worden gebruikt om binnenkomende gegevens naar het eigen rooster te interpoleren. De interpolatiemethode die daarbij wordt gespecificeerd is een samengestelde interpolatiemethode, d.i. een interpolatiemethode die bestaat uit een aantal onderdelen. Alvorens in te gaan op de betekenis en het gebruik van die onderdelen, zal eerst aangegeven worden wat een samengestelde interpolatiemethode inhoudt.

5.2.1

Interpolaties in de update operatie De update operatie bestaat intern uit vier fasen: Fase 1: Voor alle te communiceren velden Indien nodig Kopieer eigen gegevens naar de overlapset Fase 2: Voor alle buurprocessen Voor alle te communiceren velden Bepaal wat verstuurd moet worden en verstuur Fase 3: Voor alle binnenkomende boodschappen Pak de gegevens uit Plaats de gegevens op de overlapset en/of de eigen arrays eventueel met interpolatie Fase 4: Voor alle gegevens op de overlap set Interpoleer de gegevens Kopieer de gegevens naar de eigen arrays


83


Cursusboek Module 3

Een samengestelde interpolatiemethode is opgebouwd uit interpolatiemethoden voor elk van de fasen 1, 3 en 4. Voor fase 1 is de interpolatiemethode altijd een betrekkelijk simpele, want het betreft hier feitelijk geen interpolatie maar een kopiëeraktie. Dit wordt toch een interpolatie genoemd omdat het door COCLIB op vrijwel dezelfde manier afgehandeld wordt. Een dergelijke interpolatiemethode wordt een reshape methode genoemd (alleen de structuur van de data verandert, niet de inhoud). Bij fase 3 is wel echt sprake van interpolatie, want hier worden gegevens met vertikale verfijning direct geïnterpoleerd naar de juiste laagverdeling. Als er geen sprake is van vertikale verfijning dan worden de gegevens direct op de eigen arrays of op de overlapset geplaatst. Bij fase 4 vindt zowel horizontale interpolatie (binnen de overlapset) als een reshape (van overlap naar eigen arrays) plaats. De reshape operatie is (de inverse van) de methode die in fase 1 gebruikt wordt. Een samengestelde interpolatiemethode wordt met de routine cocmth opgebouwd uit vooraf gedefiniëerde enkelvoudige interpolatiemethoden. Daarbij moet men een reshape operatie opgeven en interpolatiemethoden voor fasen 3 en 4 (of de string ‘none’ als er geen interpolatie nodig is). Een aanroep van cocmth kan er als volgt uitzien: call cocmth(‘bilin’,’fullbox’,’ddh-refine’, ‘from fullbox to overlap_d’, ‘none’,’bilin’,name1(1:ilen))

De argumenten op de eerste regel zijn respektievelijk de naam van de samengestelde interpolatie methode, de naam van de indexset waarop de methode werkt en de koppelingscategorie (waarover later meer). Op de tweede regel staat de naam van de reshape methode die in fase 1 en fase 4 van de update gebruikt wordt . Op de derde regel staan de interpolatiemethoden voor vertikale verfijning (‘none’) en voor horizontale verfijning (‘bilin’). De reshape methode en de interpolatiemethoden zelf zijn van te voren gedefiniëerd. Daarop wordt verderop in dit hoofdstuk teruggekomen.

5.2.2

Koppelingscategoriën In het algemene geval kan een deeldomein gekoppeld zijn aan verschillende deeldomeinen waarbij meerdere vormen van verfijning voorkomen. Een domein kan bijvoorbeeld gekoppeld zijn aan drie andere domeinen, waarbij het ene buurdomein horizontaal verfijnd is, het tweede deeldomein vertikaal verfijnd is en het derde deeldomein niet verfijnd is (gewoonlijk aangeduid als een parallel buurproces).


84


Cursusboek Module 3

In sommige gevallen wil men van geval tot geval andere vormen van interpolatie gebruiken. Bij een parallel buurproces hoeft er helemaal geen interpolatie te gebeuren en het zou onnodige tijdverspilling zijn om dat wel te doen (waarbij dan alle gewichten gelijk aan 1.0 zouden zijn). Bij een vertikaal verfijnd buurproces wil men niet gehouden zijn om ook een (zinloze) horizontale interpolatie uit te voeren. Daarom wordt een interpolatiemethode altijd gedefinieerd voor een bepaalde koppelingscategorie, waarbij een koppelingscategorie het soort verfijning aangeeft dat het proces heeft met een buurproces. Bekende koppelingscategorieën zijn bijvoorbeeld: parallel (geen verfijning), ddvert (vertikale verfijning), ddh-refine (horizontale verfijning). Omgekeerd moet men elke samengestelde interpolatiemethode definiëren voor elke koppelingscategorie. Als men een samengestelde interpolatiemethode voor een bepaalde koppelingscategorie niet definiëert, dan wordt er aangenomen dat buurprocessen die tot die categorie behoren, moeten worden behandeld als parallelle buurprocessen. Uiteraard doet zich de vraag voor hoe COCLIB weet tot welke koppelingscategorie een buurproces behoort. Dit wordt bepaald door het aanroepen van de routine coctyp, waaraan het aanroepende proces een aantal eigenschappen van het eigen rooster en model meegeeft (bijvoorbeeld de eigen laagverdeling). De routine coctyp wisselt deze eigenschappen uit met alle buurprocessen en roept dan per buurproces een routine aan waarin op basis van de eigenschappen van het eigen domein en dat van het buurproces wordt bepaald wat de koppelingscategorie voor het buurproces is. De routine die daarvoor aangeroepen wordt is de routine coccat en deze moet in principe door de gebruiker van COCLIB geïmplementeerd worden (vgl de gebruikersroutines zoals die WAQUA voorkomen). Men is dus vrij om op basis van de uitgewisselde domeineigenschappen zelf te bepalen tot welke categorie een buurproces behoort.

5.2.3

Het definiëren van reshape methoden Een reshape methode is een methode om gegevens van de ene indexset (lees: van de ene arraystructuur) over te zetten naar een andere indexset (arraystructuur). Met name komt dit soort van methode aan de orde in fase 1 en 4 van de update operatie, waar gegevens uit de eigen arrays overgezet moeten worden van en naar een werkarray op de overlapset. Technisch gesproken is een reshape methode niets anders dan een tabel waarin staat welk element uit het bron-array moet worden gekopiëerd naar welk element van het doel-array. Dit soort tabel (die intern in COCLIB veel gebruikt wordt), wordt een relatietabel genoemd. COCLIB biedt de gebruiker enkele


85


Cursusboek Module 3

routines om dergelijke relatietabellen (of reshape methodes) aan te maken. De reshape methode die in samengestelde interpolatiemethoden wordt gebruikt, wordt gewoonlijk aangemaakt met de routine cocrss. Aan deze routine geeft men de namen mee van de bron-indexset en de bijbehorende overlapset, waarna op basis van coördinaten die met de indices zijn geassocieerd de gewenste tabel wordt opgesteld (en intern opgeslagen onder de naam die men aan de tabel verbindt). Er is één relatietabel die niet apart gedefiniëerd hoeft te worden (en ook niet mag worden) en dat is de relatietabel ‘unity’. Deze relatietabel legt tussen twee willekeurige indexsets van dezelfde lengte een 1:1 relatie tussen de elementen. Er zijn nog enkele andere routines in COCLIB voor het definiëren van relatietabellen. Deze routines spelen vooral een rol tijdens het construeren van een overlapset, en zullen hier nu niet verder behandeld worden.

5.2.4

Basismethoden en coëfficiëntensets COCLIB beschouwt een interpolatiemethode als een combinatie van een interpolatieformule, de verzameling van coëfficienten die daarin gebruikt worden en de indexset waarop de methode werkt. De interpolatieformule wordt aangeduid met de term basismethode, de gebruikte coëfficienten worden coëfficientenset genoemd. De basismethoden zijn hard ingeprogrammeerd in COCLIB en daarmee ligt ook vast welke coëfficienten er gebruikt moeten worden. De waarden van die coëfficienten kan men echter zelf bepalen en, indien nodig, ook veranderen zo vaak men maar wil. Eenmalig, aan het begin van het programma definiëert men een bepaalde coëfficientenset en vervolgens kan men die set herhaaldelijk vullen met nieuwe waarden. Dit laatste is bijvoorbeeld nodig voor vertikale interpolatie voor het geval dat een sigma laag gekoppeld is aan tenminste één vaste laag. In dat geval hangen de interpolatiegewichten af van de actuele waterstand (die de dikte van de sigma-laag bepaalt). De interpolatiegewichten zijn uitgedrukt in coëfficiënten die de actuele waterstand voorstellen. In WAQPRO wordt met regelmaat de actuele waterstand ingevuld in de betreffende coëfficiëntenset, zodat de interpolatie altijd de juiste waarden gebruikt. COCLIB zorgt ervoor dat als er nieuwe waarden voor een coëfficiëntenset worden gegeven, deze automatisch doorgegeven worden aan alle buurprocessen die daar belang bij hebben. Voor de zojuist beschreven sigma-vast interpolatie heeft een proces bijvoorbeeld de waterstanden nodig in de eigen punten en in de punten van elk van de buurprocessen.


86


5.2.5

Cursusboek Module 3

Het definiëren en vullen van een coëfficiëntenset Een coëffientenset wordt gedefinieerd met de routine cocico. Een aanroep van deze routine kan er als volgt uitzien: call cocico(‘zku_full’,’full_u*kmax0’,’stcmax’, ‘full_u*kmax0’,’unity’,’REAL’,0, icocad,name1(1:ilen))

Het eerste argument is de naam van de coëfficiëntenset, die later weer gebruikt wordt bij het definiëren van interpolatiemethoden. Daarna komt de naam van de indexset waarop de coëfficienten worden toegeleverd (d.i. de indexset van het array waarmee de coëfficiënten aan COCLIB worden doorgegeven. Vervolgens komt de naam van een interface (in dit geval ‘stcmax’). Deze interface is de interface waarop buren hun coëfficiënten uitwisselen. Als deze te klein gekozen wordt, dan kan het voorkomen dat een proces nieuwe waarden voor zijn coëfficiënten specificeert (met coccoe, waarover later meer), deze nieuwe coëfficienten niet gebruikt worden in de interpolatie die in buurprocessen worden uitgevoerd. Vervolgens komt en nogmaals de naam van een indexset (in dit geval dezelfde als de indexset waarop de waarden toegeleverd worden). Deze tweede indexset is de indexset waar de coëfficiënten op gedefiniëerd zijn. Hierbij is het voor coëfficiëntensets die gebruikt worden voor vertikale interpolatie van belang of dit een één-dimensionale indexset is of een meerdimensionale indexset. Als een interpolatiemethode met een één-dimensionale indexset wordt toegepast op een array met een meerdimensionale indexset, dan wordt de coëfficientenset in ieder horizontaal punt opnieuw gebruikt. Met andere woorden de coëfficiënten worden horizontaal constant verondersteld. Dit speelt weer met name een rol bij een koppeling van een sigma laag aan een vaste laag. In dit geval zijn de coëfficiënten niet horizontaal constant (de waterstand is immers niet horizontaal constant) en dus moet men een meer-dimensionale indexset gebruiken om de coëfficiënten of te definiëren. In de bovenstaande aanroep van cocico staat vervolgens de naam van een reshape methode waarmee gegevens van de eerstegenoemde indexset (die waarop de coëfficiënten worden toegeleverd) worden overgezet op de tweede (die waarop de coëfficiëntenset feitelijk is gedefiniëerd). In dit geval is dat de voorgedefiniëerde reshape methode ‘unity’ omdat beide indexsets aan elkaar gelijk zijn. Het volgende argument in de aanroep van cocico geeft het datatype van de coëfficiënten aan. Tenslotte komt er nog een getal, waarmee aangegegeven wordt of er nog speciale eigenschappen met de coëfficiënten verbonden zijn. Het kan


87


Cursusboek Module 3

bijvoorbeeld zijn dat men bij het gebruik van een coëfficiëntenset er van uitgaat dat de coëfficiënten niet meer wijzigen nadat ze eenmaal gedefiniëerd zijn, zodat evenuteel afgeleide grootheden slechts eenmalig berekend hoeven te worden. In dat geval moet men bij de aanroep van cocico aangeven dat de coëfficiëntenset niet ververst mag worden. Nadat een coëfficiëntenset gedefiniëerd is door het aanroepen van cocico moeten tenminste eenmaal de waarden van de coëfficienten gevuld worden met de routine coccoe. Het gebruik van deze routine spreekt voor zich en zal hier niet verder worden toegelicht.

5.2.6

Interpolatiemethoden voor vertikale verfijning Voor het definiëren van interpolatiemethoden voor vertikale verfijning is er slechts één routine nodig in COCLIB en dat is de routine cocint. Een voorbeeld van een aanroep van deze routine is: call cocint(‘meth_w(sigfix)’,’full_s*kmax’, ‘ddv-sigfix’, ‘linear for layer-interface samples’, ‘zks_full’, icocad, name1(1:ilen))

Het eerste argument is de naam van de nieuw te definiëren methode. Daarna komen de naam van de indexset waarop de methode werkt en de naam van de koppelingscategorie waarvoor de methode te gebruiken is. Het feit dat men de koppelingscategorie hier ook moet opgeven (en nogmaals bij het definiëren van een samengestelde methode) is een overblijfsel uit de COCLIB versie waarin samengestelde methoden nog niet voorkwamen. In elk geval moet de koppelingscategorie die hier opgegeven wordt dezelfde zijn als de koppelingscategorie van een samengestelde methode waarin deze vertikale interpolatie gebruikt wordt. Vervolgens wordt de naam van één van de basismethoden gegeven. Zoals boven al beschreven specificeert dit feitelijk de formule van de interpolatie. Er worden drie methoden ondersteund door COCLIB: ‘constant for layer-averages’ Deze methode wordt gebruikt om variabelen te interpoleren die een laaggemiddelde voorstellen zoals bijvoorbeeld de horizontale stroomsnelheden u en v. ‘constant for layer-integrals’ Deze methode wordt gebruikt om variabelen te interpoleren die een laagtotaal voorstellen zoals bijvoorbeeld debieten. ‘linear for layer-interface samples’ Deze methode wordt gebruikt om variabelen te interpoleren die op een bepaald punt van de laag gedefinieerd zijn, zoals bijvoorbeeld de verticale stroomsnelheid w.


88


Cursusboek Module 3

Indien er ooit andere methoden nodig zouden zijn, dan moeten die in COCLIB bijgeprogrammeerd worden. Tenslotte wordt, voorafgaand aan de normale COCLIB argumenten icocad en name, de naam van de te gebruiken coëfficiëntenset opgegeven. In dit geval is dat de coëfficiëntenset van laaginterface-locaties. Dit is in dit geval een coëfficiëntenset die gedefinieerd is op de 3D indexset ‘fullbox*kmax0’, zodat de coëfficiënten per 3D index kunnen verschillen. Dit is nodig omdat deze interpolatie methode gebruik wordt voor de zogenoemde sigma-vast koppeling, waarbij een vaste laag aan een sigma-laag gekoppeld wordt. In dat geval moeten verschillen de interpolatiegewichten zowel per laag als ook in de horizontaal van elkaar (omdat ze immers direct afhangen van de waterstand ter plaatse). In het geval dat er alleen sigma lagen aan sigma lagen en/of vaste lagen aan vaste lagen gekoppeld worden, kan worden volstaan met een 1D coëfficiëntenset gedefinieerd op ‘kmax0’, omdat de interpolatiegewichten horizontaal niet van elkaar zullen verschillen. COCLIB ziet aan de afmeting van de coëfficiëntenset of deze voor elk horizontaal punt ongewijzigd mag worden gebruikt (1D coëfficiëntenset) of dat er per horizontaal punt andere gewichten nodig zijn (3D coëfficiëntenset). Merk op dat de coëfficiëntenset ‘zks_full’ die in dit geval gebruikt wordt ook iedere tijdstap ververst moet worden om ervoor te zorgen dat alle interpolatiegewichten correct zijn.

5.2.7

Interpolatiemethoden voor horizontale verfijning Voor horizontale verfijning is er niet, zoals bij verticale verfijning, één enkele routine om interpolatiemethoden te definiëren. Dit is omdat de argumentenlijsten per basismethode nogal van elkaar verschillen, waardoor het niet goed mogelijk is om daar één enkele routine voor te gebruiken. Waar men dus bij verticale verfijning de naam van de basismethode als argument moet doorgeven, moet men bij horizontale verfijning de routine aanroepen die bij de gewenste basismethode hoort. Er worden vijf basismethoden voor horizontale verfijning ondersteund: Lokaal constant waarbij de waarde in de gehele horizontale roostercel constant wordt verondersteld. Bij interpolatie van grof naar fijn wordt de waarde van de grove cel in elke fijne cel gezet. Bij interpolatie van fijn naar grof worden de waarden van de fijne cellen gemiddeld, gewogen naar het oppervlakte van de cel. Interpolaties gebaseerd op deze basismethode worden aangemaakt met routine coccel.


89


Cursusboek Module 3

Lin(weigh) die gebruikt wordt voor interpolatie van grootheden die op cellfaces gedefiniëerd zijn (zoals horizontale stroomsnelheden). Deze methode interpoleert lokaal constant in de richting van de cellface en lineair in de richting loodrecht daarop. Interpolaties op basis van deze basismethode worden aangemaakt met routine cocclf. Bilineair hetgeen voor zich spreekt. Interpolaties met deze basismethode worden aangemaakt met de routine cocbil. Sum die gebruikt wordt voor debieten. Bij interpolatie van fijn naar grof worden de waarden van de fijne cellen bij elkaar opgeteld. Bij de interpolatie van grof naar fijn worden de waarden op een andere manier geïnterpoleerd (bijvoorbeeld lin(weigh)). Bij het definiëren van interpolatiemethoden op basis van de Sum methode moet men dan ook altijd de naam van een andere interpolatiemethode opgeven. Max voor het interpoleren van integer waarden. Bij de interpolatie van fijn naar grof wordt de grootste waarde van alle fijne cellen binnen een grove cel genomen. Bij de interpolatie van grof naar fijn wordt de waarde van de grove cel verdeeld over de fijne cellen volgens één van de andere interpolatiemethoden (zoals ook gebeurt bij de Sum methode). Deze methode wordt gebruikt om masker-arrays te interpoleren. Een uitvoerige beschrijving van elk van de interpolatiemethoden, inclusief gebruikte formules, is opgenomen in (Vollebregt et al. 2001), dat desgewenst opgevraagd kan worden bij de projectleider van SIMONA*SYS. De routines coccel, cocclf en cocbil vergen veel meer argumenten dan wat bij vertikale verfijning het geval is. Daarom is het illustratief om een voorbeeld nader te bekijken. Een aanroep van de routine coccfl kan er als volgt uitzien: Call coccfl(‘lin(weigh)’,’overlap_u’, mdscru,mskitp,‘mnmaxk_u’, ’from overlap_u to mnmaxk_u’, ‘maximum interface’, ‘coordmd_u’,’coordm_u’, ’coordndu’,’coordn_u’, ’cellface_u’, icocad, name1(1:ilen))

Als eerste argument wordt de naam van de nieuw te definiëren methode opgegeven (‘lin(weight)’) gevolgd door de naam van de overlapset die gebruikt wordt voor de interpolatie. Daarna volgen twee masker arrays, waarin respectievelijk de discretisatiepunten en de interpolatiepunten met de waarde 1 zijn aangegeven. Deze masker arrays leven op de indexset ‘mnmaxk_u’ die als derde argument op de tweede regel wordt opgegeven, gevolgd door de reshape methode waarmee waarden van deze indexset naar de overlapset overgezet


90


Cursusboek Module 3

worden. Vervolgens wordt er een argument (‘maximum interface’) opgegeven dat aangeeft welke interface er gebruikt moet worden om masker-arrays op de overlapset te communiceren. Het begrip van interfaces bij horizontale verfijning komt verderop in dit hoofdstuk nog uitgebreid aan bod. Tenslotte komen de vijf coëfficiëntensets die voor deze interpolatiemethode nodig zijn en de gebruikelijke COCLIB argumenten icocad en name. De coëfficiëntensets moeten allemaal ‘read-only’ zijn (zie het definiëren van coëfficiëntensets) omdat ze gebruikt worden om een veel efficiëntere opslag van gewichten te berekenen en deze berekening niet bij elke communicatie herhaald wordt. De coëfficiëntwaarden mogen dus niet veranderen tijdens een run.

5.3

Interfaces bij horizontale verfijning In het vorige hoofdstuk is het begrip interfaces al ter sprake gekomen omdat het een belangrijke rol speelt bij de parallelle communicaties. Voor domein decompositie met horizontale verfijning is dit begrip iets uitgebreid, hetgeen ook betekent dat er bij het definiëren van interfaces een aantal extra argumenten meegegeven moeten worden. Het belangrijkste verschil met interfaces voor parallelle communicaties (en ook voor verticale verfijning) bestaat eruit dat bij communicaties met horizontale verfijning gebruik wordt gemaakt van de overlapset omdat de roosters van naburige domeinen in horizontale zin gewoonlijk niet dezelfde celgrootte hebben. Bij parallel rekenen kent ieder proces het (logische) rooster van zijn buurproces (omdat dit immers zijn eigen guardband is). Het kent dus ook de wijze waarop het buurproces het opgegeven masker (waarmee de toepassing van het stencil beperkt wordt tot de punten waarvoor echt nieuwe waarde ontvangen moeten worden) op dat rooster toepast. Bij horizontale verfijning weet een proces niet precies welk masker een buurproces toepast. Het zou zijn eigen masker kunnen overzetten naar de overlapset en dan interpoleren naar de roosterpunten van de betreffende buurman binnen de overlapset, maar het is nooit zeker of dit geïnterpoleerde masker precies gelijk is aan het masker dat het buurproces zelf specificeert. Dat betekent dat processen elkaar moeten laten weten welke maskers ze toepassen, zodat daarna elk proces precies kan bepalen welke waarden er naar het buurproces toegestuurd moeten worden. Binnen het definiëren van een interface is dus een communicatie geïntroduceerd. Deze communicatie is feitelijk een gewone cocupd en die heeft informatie nodig over


91


Cursusboek Module 3

de te gebruiken overlapset, de reshape operatie waarmee gegevens van de indexset zelf naar de overlapset (en terug) kunnen worden gezet en tenslotte nog een interface. Merk op dat men dus voor het definiëren van een interface al een interface nodig heeft. Hierover zal zodadelijk nog iets meer gezegd worden. Voor het communiceren van maskers bij het definiëren van interfaces wordt een interface gebruikt die in elk geval groot genoeg is: het masker moet immers gecommuniceerd worden voor alle punten die gebruikt worden binnen de interface die gedefiniëerd wordt. Daarom wordt altijd een zogeheten ‘maximale interface’ gebruikt. Merk op dat deze interface alleen gebruikt wordt voor het communiceren op de overlapset zelf; het is dus een interface op de overlapset. Een voorbeeld van het definiëren van een interface voor horizontale verfijning is de volgende: Call cocitf(‘mnmaxk_s’, mskall, noffs, istcof,idimof, ‘mnmaxk_s’, mskall, + ‘overlap_s’,‘from mnmaxk_s to overlap_s’, + ‘maximum interface’, ‘stcmax’, icocad, name1(1:ilen))

Op de eerste regel staan de argumenten die al bekend zijn van de interfaces voor parallel rekenen. Op de tweede regel staat als eerste de naam van de overlapset die gebruikt wordt, gevolgd door de naam van de reshape methode die gebruikt wordt om gegevens van de indexset ‘mnmaxk_s’ naar de overlapset ‘overlap_s’ over te zetten (en terug). Vervolgens, staat op de derde regel de naam van interface die gebruikt wordt om maskers te communiceren (‘maximum interface’), waarna weer de bekende overige argumenten volgen. Zoals zojuist al opgemerkt werd, lijkt er een probleem te ontstaan met het feit dat men een interface nodig heeft om een andere interface te kunnen maken. Dit roept de vraag op hoe de ‘eerste’ interface dan gemaakt wordt. Daarvoor heeft COCLIB een aparte routine (cocifd). Met deze routine creëert men een interface direct, dus niet via het gebruik van een stencil. Bij de set waarop de interface gedefinieerd moet worden, geeft men dan direct lijsten op waarin per buurproces aangegeven wordt voor welke indices er waarden ontvangen moeten worden en voor welke indices er waarden naar dat buurproces toe gestuurd moeten worden. De lijsten zelf worden aangemaakt met de routine cocrdm. Een voorbeeld van het definiëren van een interface met cocifd is het volgende, waarbij alleen de relevante code fragmenten weergegeven zijn: do 100 iprc=1,numprc c c

fill iwork such that indices that are to be sent to process iprc are marked with ‘1’ write(tabnam,*) ‘internal itf overlap_s’,


92


+ + +

Cursusboek Module 3

iprc call cocrdm(tabnam,’overlap_s’, iwork, 1, itftab(iprc,1), icocad, name1(1:ilen)) :

100

continue call cocifd(‘overlap_s’, ‘maximum interface’, + itftab, numprc, icocad, + name1(1:ilen))

Als eerste wordt in een loop over alle buurprocessen steeds een tabel aangemaakt waarin staat welke indices van de overlapset ‘overlap_s’ verzonden moeten worden naar de betreffende buurman. Een handle naar deze tabel wordt opgeslagen in het array itftab(numprc,2), waarin allereerst handles naar de interne interfaces per buurproces zijn opgeslagen en vervolgens handles naar de externe interface per buurproces. Als voor elke buurman tabellen zijn aangemaakt voor de interne en externe interfaces, wordt tenslotte de interface zelf aangemaakt met cocifd. Overigens is het zeer waarschijnlijk dat de meeste programmeurs zelf nooit de routine cocifd zullen gebruiken.


93


Cursusboek Module 3

6

Interne werking van COEXEC

6.1

De functies van COEXEC COEXEC (vroeger ook wel master genoemd) is het programma dat de parallelle WAQUA/TRIWAQ processen (verder rekenprocessen genoemd) opstart en beheert. Het programma is bedoeld als een generiek programma dat geen kennis heeft van WAQUA/TRIWAQ en dus toepasbaar is voor elke applicatie die ontwikkeld is op basis van COCLIB. Met de introductie van DDHOR heeft COEXEC echter ook een beperkte inhoudelijke rol gekregen, doordat het programma een aantal checks uitvoert op de consistentie van de aangeboden domeinen. Daarbij heeft het bovendien de mogelijkheid gekregen om alleen die checks uit te voeren en niet daadwerkelijk over te gaan tot het starten van de rekenprocessen, zodat COEXEC ook gebruikt kan worden om de invoer te controleren. Het opstarten van de processen omvat niet alleen het daadwerkelijk starten van de processen, maar ook het toesturen van een opstartboodschap, waarin de informatie is opgenomen die een rekenproces normaalgesproken via de commandline binnen zou krijgen. In principe is COEXEC erop ingericht om een goede plaatsing (mapping) van processen op processoren te bepalen, maar deze functionaliteit is niet uitgewerkt. Het beheer van processen houdt een aantal zaken in. Ten eerste handelt COEXEC situaties af waarbij één van de rekenprocessen in de problemen komt. Daarnaast verzamelt het de uitvoer die de processen naar standaard uitvoer en standaard error printen, voegt aan die uitvoer een label toe waaraan te zien is welk proces de uitvoer genereerde en schrijft die uitvoer vervolgens naar zijn eigen standaard uitvoer.

6.2

In- en uitvoer van COEXEC Figuur 7 toont de in- en uitvoer van COEXEC. De informatie over welke processen er opgestart moeten worden, leest COEXEC uit een configuratiefile. Voor parallel rekenen en domein decompositie met vertikale verfijning bevat deze configuratiefile een lijst van op te starten processen met daarbij vermeld onder welke naam ze bekend moeten zijn, welke executable er moet worden uitgevoerd, hoeveel exemplaren van het proces er moeten worden opgestart en welke invoerargumenten ze moeten hebben.


94


Cursusboek Module 3

Title: coexec.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 7: overzicht van de in- en uitvoer van COEXEC

Het aantal exemplaren dat van een proces moet worden opgestart is bedoeld om een deeldomein parallel door te rekenen. Deze mogelijkheid is wel voorzien maar wordt verder niet ondersteund. De mogelijkheid om per proces een andere executable te laten uitvoeren, wordt gebruikt om processen met verschillende buffersizes te gebruiken, zodat niet alle processen met de maximaal benodigde buffersize hoeven te draaien. Dit is vooral van belang als een domein decompositie som op een enkele processor wordt uitgevoerd. De configuratiefile voor domein decompositie met horizontale verfijning bevat vergelijkbare maar net iets andere informatie. Dit is vooral gedaan omdat de configuratiefile voor domein decompositie met horizontale verfijning normaalgesproken expliciet door de gebruiker geschreven wordt, zodat de inhoud ervan zoveel mogelijk moet aansluiten bij de beleving van de gebruiker.


95


Cursusboek Module 3

Om het onderscheid tussen de twee soorten files te kunnen maken begint een configuratiefile voor parallel rekenen en domein decompositie met vertikale verfijning met het keyword PROCESSES, terwijl een configuratiefile voor domein decompositie met horizontale verfijning begint met het keyword DOMAINS. Dit geeft direct ook het verschil in benadering aan. Bij domein decompositie met horizontale verfijning volgt na het hoofdkeyword DOMAINS een lijst met domeinen. De domeinen komen hierbij overeen met de aparte siminpfiles/roosters die aan elkaar gekoppeld worden (merk op dat er in parallelle sommen en in DDVERT-berekeningen een enkel globaal domein/siminp-file is). Per domein worden opgegeven: de logische naam van het domein, het runid van de SDS-file voor het domein, de naam van het experiment op die SDS-file en de naam van de executable die uitgevoerd moet worden. Bovendien kan er nog een parameter opgegeven worden aan COEXEC, waarmee ingesteld wordt hoe nauwkeurig COEXEC moet kijken of de aangeboden roosters op elkaar aansluiten. Behalve de configuratiefile gebruikt COEXEC in het geval van domein decompositie met horizontale verfijning ook de 000file(s) van de domeinen. Daaruit wordt alleen het aantal deeldomein gelezen dat uit het betreffende domein gehaald is door COPPRE. Op basis daarvan bepaalt COEXEC hoeveel processen er voor het domein moeten worden opgestart en welke invoerargumenten deze processen moeten hebben.

6.3

Uitgevoerde checks Zoals gezegd, voert COEXEC een aantal controles uit in het geval van domein decompositie met horizontale verfijning. Deze controles zijn nodig omdat de gebruiker de SDS-files voor de verschillende domeinen onafhankelijk van elkaar aanmaakt. COEXEC is het eerste programma dat de invoer als geheel kan overzien. Dit in tegenstelling tot de situatie bij parallel rekenen en domein decompositie met vertikale verfijning, waar COPPRE in principe de gehele configurtie kent en alle benodigde controles kan uitvoeren. COEXEC kan ook in een ‘checkonly’ mode gebruikt worden, waarbij uitsluitend de controles uitgevoerd worden. Op die manier kan eerst gekeken worden of de invoer valide is voordat er daadwerkelijk een run opgestart wordt. COEXEC controleert: of alle domeinen hetzelfde procesmodel gebruiken (dus allemaal met transport of juist allemaal zonder, allemaal


96


Cursusboek Module 3

hetzelfde turbulentiemodel, allemaal ofwel WAQUA ofwel TRIWAQ, etc.) of de roostertypes gecombineerd mogen worden (rechlijnig of kromlijnig mogen gecombineerd worden, bolcoördinaten mogen niet met andere roostertypes samen gebruikt worden) of, in geval van transportberekening, het aantal opgeloste stoffen hetzelfde is of de tijdstapinformatie (begin- en eindtijd en tijdstap) en referentiedatum hetzelfde is of dezelfde iteratieconstanten gebruikt worden of aangrenzende domeinen die horizontaal verfijnd zijn dezelfde laagverdeling hebben (dus niet ook vertikaal verfijnd zijn of een andere verdeling van de lagen hebben) of aangrenzende domeinen die horizontaal verfijnd zijn een geldige verfijning hebben. Deze laatste controle is een zeer complexe waarin erg veel situaties worden afgevangen die tot problemen bij de berekening zou kunnen leiden.

6.4

Het opstarten van processen In COEXEC zijn voorzieningen getroffen om expliciet te bepalen welk proces op welke processor moet draaien. Deze voorzieningen zijn echter niet uitgewerkt en ook nooit nodig gebleken. Op dit moment worden processen eenvoudigweg round-robin op de beschikbare hosts geplaatst. Voor het opstarten van rekenprocessen wordt de PVM operatie pvmfspawn gebruikt. Als deze operatie om enige reden niet lukt, dan wordt een foutcode teruggegeven. De betekenis daarvan is terug te vinden in de PVM documentatie. De foutcode –7 komt relatief veel voor en verdient daarom enige speciale aandacht. Deze foutcode houdt in dat de executable die opgestart moet worden niet gevonden is. Dit kan erop duiden dat het zoekpad waarin COEXEC de executables zoekt niet goed is ingesteld (dit moet in de PVM hostfile gebeuren, waarvan de default versie in de SIMONA bindirectory te vinden is). Andere mogelijke oorzaken zijn dat de naam van de executable niet goed is ingevuld, of dat er iets is misgegaan bij het hercompileren van WAQPRO met een afwijkende buffersize.

6.5

Boodschappen tussen WAQPRO en COEXEC Zodra een proces door COEXEC is opgestart, ontvangt het van COEXEC een opstartboodschap. Deze boodschap wordt opgevangen in de COCLIB routine COCINI. De boodschap


97


Cursusboek Module 3

bevat in de eerste plaats de argumenten die de rekenprocessen normaalgesproken via de commandline binnen zouden krijgen. Deze argumenten leest COEXEC uit de configuratiefile. COCINI geeft deze argumenten terug als argument zodat het aanroepende proces ze kan gebruiken. Behalve de invoerargumenten bevat de opstartboodschap ook informatie over de configuratie (aantal deeldomeinen, de runids van de deeldomeinen) en informatie over het de plaats van het proces binnen die hele configuratie (domeinnummer, nummer van het proces in de lijst van alle processen). Tenslotte bevat de opstartboodschap ook nog de waarde van de parameter waarmee in het geval van horizontale verfijning de nauwkeurigheid wordt aangegeven waarmee roosters aan elkaar moeten passen. Nadat de processen opgestart zijn, gaat COEXEC over tot het afhandelen van boodschappen van de rekenprocessen. Dit is een passieve rol, waarbij COEXEC vrijwel uitsluitend actief is op het moment dat een rekenproces een boodschap gestuurd heeft. Als COEXEC gedurende langere tijd geen boodschappen ontvangt, dan zal het een controle uitvoeren om na te gaan of alle processen nog actief zijn. Als daarbij geconstateerd wordt dat alle processen onverwacht verdwenen zijn, dan wordt de run afgebroken met een foutmelding. Een dergelijke afloop is zeer ongebruikelijk, maar kan een enkele keer voorkomen. Meestal is daarbij sprake van een programmeerfout of een probleem met het operating systeem. De boodschappen die uitgewisseld worden tussen de rekenprocessen en COEXEC vallen uiteen in twee typen. Enerzijds zijn er boodschappen die door PVM (de elementaire communicatie bibliotheek waarop COEXEC/COCLIB gebaseerd is) worden gegenereerd. Anderzijds zijn er boodschappen die door COEXEC/COCLIB zelf gegenereerd worden. De PVM-boodschappen die de rekenprocessen naar COEXEC sturen omvatten boodschappen waarmee het rekenproces meldt dat het actief geworden is dat het termineert dat het iets naar standaard uitvoer of standaard error geschreven heeft. De laatste boodschap bevat ook de string die naar standaard uitvoer of –error gestuurd is. Deze string wordt voorzien van een label met een identificatie van het versturende proces en vervolgens naar COEXEC’s eigen standaard uitvoer geschreven.


98


Cursusboek Module 3

De COCLIB boodschappen die de rekenprocessen naar COEXEC sturen, omvatten boodschappen waarmee het rekenproces meldt dat het verdacht lang moet wachten in een communicatiestap en aan COEXEC vraagt om na te gaan of de communicatiepartner nog gewoon functioneert dat een voorheen gemeld communicatieprobleem is opgelost dat het proces COCLIB heeft afgesloten Als een proces problemen meldt met een communicatie, dan zal COEXEC allereerst nagaan of de communicatiepartner nog actief is. Als dat het geval is, dan meldt COEXEC aan de probleemmelder terug dat die het nogmaals moet proberen. Het rekenproces zal deze procedure een aantal keer doorlopen totdat ofwel de communicatie alsnog slaagt (waarop dit aan COEXEC wordt gemeld) ofwel COEXEC dermate vaak een probleemmelding over deze communicatie heeft ontvangen dat hij besluit dat het geen zin meer heeft om verder te proberen. In dat geval ontvangen alle processen een zogeheten panic message, op grond waarvan ze onmiddelijk zullen besluiten om te termineren zodra een communicatie niet lukt. Op sommige platforms is het nodig om een proces te laten wachten met termineren totdat alle processen klaar zijn. Dit omdat anders het operating system tegelijk met het eerst terminerende proces alle processen beëindigt. Daarom heeft COCLIB een voorziening om in COCEND te blijven wachten totdat COEXEC meldt dat alle processen COCEND hebben aangeroepen en de hele run kan worden beëindigd. De boodschappen die COEXEC aan de rekenprocessen stuurt zijn hierboven al aan terloops aan de orde geweest. Deze boodschappen omvatten boodschappen waarmee COEXEC de opstartinformatie doorgeeft aan elk rekenproces meldt dat functioneert

een

beoogde

communicatiepartner

nog

meldt dat er een serieus probleem is ontstaan in de run (bijvoorbeeld doordat een communicatie langdurig niet lukt), de zogenaamde panic message meldt dat de run beëindigd kan worden.


99


Cursusboek Module 3

7

Interne werking van de partitioner COPPRE

7.1

De functies van COPPRE COPPRE (ook wel de de partitioner genoemd) is het programma dat gebruikt wordt om uit een SDS-file die door WAQPRE wordt gegegeneerd een SDS-file voor een deeldomein te maken. Voor parallel rekenen heeft COPPRE een drieledige functie: 1. het bepalen van een geschikte opdeling van het rekenrooster, 2. het aanmaken van SDS-files voor elk van de stukken (deeldomeinen) waarin het rooster is opgedeeld 3. en het toevoegen van enkele extra gegevens aan de SDS-files voor de deeldomeinen, die betrekking hebben op het parallel rekenen. COPPRE kan werken met een extern toegeleverde opdeling van het rekenrooster (partitionering). Als de gebruiker zelf een partitionering toelevert, dan zal COPPRE daar niets aan toevoegen. COPPRE kan ook zelf de partitionering bepalen op basis van enkele ingebouwde heuristieken, waaruit de gebruiker moet kiezen (een beschrijving is te vinden in de “User’s guide for Parallel WAQUA/TRIWAQ and for Domain Decomposition”). De opdeling van het rekenrooster die COPPRE bepaalt, wordt weggeschreven in de resultaat file (coppre-r.runid). Men kan COPPRE dan ook gebruiken om uitsluitend een partitionering te maken (dus zonder deeldomein SDS-files aan te maken). De coppre-r.runid file kan men wijzigen en vervolgens weer aan COPPRE geven als extern toegeleverde partitie. Op die manier kan men zelf ingrijpen in de partitionering die COPPRE maakt. Voor het aanmaken van de SDS-files voor de deeldomeinen bouwt COPPRE eerst een aantal tabellen op waarin is vastgelegd hoe gegevens voor de deeldomein SDS-file te vinden zijn in de gegevens op de originele SDS-file. Deze gegevens worden bewaard in een aparte SDS-file (de zogenoemde 000-file, omdat de naam ervan de structuur SDSrunid-000 heeft). De collector COPPOS gebruikt deze zelfde tabellen om de rekenresultaten in de deeldomein SDS-files weer terug te zetten in de originele SDS-file. COPPOS kan dus niet werken zonder dat COPPRE eerst de tabellen heeft aangemaakt. In die zin kan men het aanmaken van tabellen als een aparte functie van COPPRE zien. Voor domein decompositie heeft COPPRE nog een tweetal extra functies gekregen.


100


Cursusboek Module 3

Figuur 8: het toevoegen van een guardband voor DDHOR

In de eerste plaats voert COPPRE voor DDVERT ook een aantal consistentie controles uit. Zo moet de laagverdeling van aan elkaar grenzende deeldomeinen op elkaar aansluiten, waarbij alleen een verfijning van 1:n per laag van het grove deeldomein is toegestaan. Bij het aanmaken van een deeldomein voor DDVERT zal COPPRE controleren of het aan te maken deeldomein inderdaad een goede laagverdeling heeft ten opzichte van aangrenzende deeldomeinen die eerder zijn gegenereerd. Ten tweede heeft COPPRE voor DDHOR als extra functie de taak gekregen om een guardband toe te voegen op plaatsen waar dat nodig is. Bij DDHOR wordt ieder deeldomein door een aparte siminp-file gemodelleerd. Dat betekent dat er aan de rand van het deeldomein in principe geen guardband bestaat. Als men zo’n deeldomein dan toch wil koppelen aan een ander domein, zal er een guardband aangemaakt moeten worden door het gehele deeldomein groter te maken en en stuk te verschuiven. Dit wordt geïllustreerd in Figuur 8.

7.2

Invoer en uitvoer van COPPRE Figuur 9 geeft een overzicht van de in- en uitvoer van COPPRE. De verschillende onderdelen daarvan zullen nu kort beschreven worden.

7.2.1

De te partitioneren SDS-file(s) Aan de invoerzijde is er allereerst uiteraard de SDS-file die gepartitioneerd moet worden. Deze file kan eventueel een referentie bevatten naar een SVWP file waarin ruimtelijk variërende wind en luchtdruk opgeslagen liggen. Als een dergelijke referentie voorkomt, dan wordt ook die SVWP file gepartitioneerd.


101


Cursusboek Module 3

Title: coppre.eps Creator: fig2dev Version 3.2 Patchlevel 3c Preview: This EPS picture was not saved with a preview included in it. Comment: This EPS picture will print to a PostScript printer, but not to other types of printers.

Figuur 9 overzicht van de in- en uitvoer van COPPRE

7.2.2

De configuratie file Verder is er een configuratiefile, die normaalgesproken door het runscript wordt aangemaakt. Hierin wordt gespecificeerd hoe de opdeling in deeldomeinen is. Het format van de configuratiefile is beschreven in de “Users Guide for Parallel WAQUA/TRIWAQ and for Domain Decomposition”. De configuratiefile kan ofwel als hoofdkeyword het woord PARTITIONING bevatten, waarna een partionering voor parallel rekenen volgt, ofwel als hoofdkeyword het woord DECOMPOSITION bevatten, waarna een partitionering voor domein decompositie volgt. In hoofdzaak zijn er drie manieren om een partitionering te specificeren: 1.


door het specificeren van één van de ingebouwde partitioneringsmethoden

102


Cursusboek Module 3

2.

door het toeleveren van een areas file of

3.

door het toeleveren van een partitioneringsfile

De eerste methode wordt niet ondersteund voor domein decompositie omdat het nooit zinvol is om COPPRE zelf een geschikte domein decompositie te laten bepalen. Naast een specificatie van de partitionering bevat de configuratiefile optioneel ook een stuk informatie over de parallelle computer waarop gerund gaat worden. Op basis van deze gegevens kan COPPRE eventueel een partitionering optimaliseren, maar deze functionaliteit wordt niet gebruikt.

7.2.3

De LDS-beschrijving(en) Tenslotte heeft COPPRE ook nog informatie nodig over de gegevens die op de SDS-file staan. Deze wordt toegeleverd via de file coplds.ldsnaam, waarbij op dit moment alleen de ldsnamen waqua en svwp ondersteund worden. Het format van deze file is beschreven in de “User’s guide for Parallel WAQUA/TRIWAQ and for Domain Decomposition”. De file bevat drie hoofdstukken: 1. parameters 2. indexsets 3. arraybeschrijvingen Onder het keyword parameters kan men de waarde zetten van een parameter die gebruikt wordt in de verdere beschrijving van de LDS structuur. Deze waarde kan men ofwel direct zetten, of laten inlezen uit de SDS-file. Zo zou men bijvoorbeeld een parameter KURFLG kunnen definiëren, waarvan COPPRE de waarde inleest vanuit de toegeleverde SDS-file. Onder het keyword indexsets wordt beschreven welke indexsets er gebruikt worden in de beschrijving van arrays verderop, en hoe die indexsets gedistribueerd moeten worden. Beschouw bijvoorbeeld het array LGRID, met afmeting mmax bij nmax. De indexset van LGRID is dan de set [1,mmax] [1,nmax]. De set [1,mmax] kan men een naam geven. In de huidige file coplds.waqua heeft deze indexset de naam mmax gekregen. Dit is enerzijds wat verwarrend omdat de indexset wordt aangeduid met de naam die doorgaans de grootte van de indexset aangeeft. Anderzijds is dit wel een comfortabele conventie bij het beschrijven van de structuur van arrays (zie onder). Bij het definiëren van een index set moet men zijn grootte aangeven (d.i. de cardinaliteit van de set), wat gewoonlijk ook weer gebeurt met een symbolische verwijzing naar een waarde op de SDS-file. Bijvoorbeeld: de indexset kmax heeft een


103


Cursusboek Module 3

grootte die gegeven wordt door ‘MESH_IDIMEN(18)’, waarbij COPPRE zelf opzoekt welke waarde op die plek op de SDS-file staat. Behalve de grootte van de indexset moet men ook definiëren hoe de indexset gedistribueerd moet worden. Daarbij kan men de indexset ofwel repliceren (d.i. niet distribueren maar integraal overnemen in alle deeldomein SDS-files, denk aan KMAX), ofwel distribueren. In dat laatste geval moet men de distributietabel (of ook MAP) aangeven. Dat begrip komt in de volgende paragraaf uitgebreider aan de orde. Bovendien kan men bij een gerepliceerde indexset aangeven welke interpolatiemethode er gebruikt moet worden bij het samenvoegen van resultaten in COPPOS. Dit komt uitgebreider aan de orde bij de beschrijving van van COPPOS. Tenslotte bevat de LDS-beschrijving de specificatie van de feitelijke datastructuren. Dat gaat in twee stappen: eerst worden alle compounds gespecficieerd en daarna alle leaf-arrays. Bij de compounds moet worden aangegeven of de compound tijdsafhankelijk is of niet door middel van het keyword NTIMES. De naam van het keyword lijkt te suggereren dat men het aantal tijdstippen moet opgeven, maar feitelijk wordt elke waarde groter dan nul opvat als een indicatie van tijdsafhankelijkheid. Als men voor NTIMES een negatieve waarde opgeeft, dan geeft men daarmee aan dat het array niet tijdsafhankelijk is maar wel door COPPOS moet worden behandeld (zie ook de beschrijving van COPPOS). Bij de leaf-arrays geeft men niet alleen aan op welke plek onder welke compound het array thuishoort, maar ook wat de structuur van het array is in termen van indexsets. Van het array LGRID zal men dus aangeven dat het thuishoort onder MESH en dat zijn indexset de structuur MMAX*NMAX heeft, waarbij de namen MMAX en NMAX de namen zijn van indexsets die eerder gedefiniëerd zijn. Het ‘*’-teken geeft het Carthesisch produkt van twee indexsets aan. Behalve dit teken worden ook de ‘+’ en de max() ondersteund. Met de plus-operator wordt de concatenatie van twee sets wordt bedoeld, de max-operator selecteert uit een aantal index-sets de set met de grootste afmeting. In de beschrijving van een array kan men naast expliciet gedefiniëerde indexsets ook impliciete indexsets gebruiken, die aangegeven worden met een getal dat de grootte van de set aangeeft: de indexset ‘3’ stelt de set [1..3] voor.

7.2.4

De COPPRE-R-file Na het bepalen van de partitionering exporteert COPPRE de gevonden opdeling van het rooster in een report file, copprer.runid. Deze file heeft het format van een partitioning-file, en kan dus desgewenst ook in het vervolg aangeboden worden aan COPPRE voor nieuwe runs met hetzelfde model. Dit is vooral


104


Cursusboek Module 3

zinvol omdat men de file zelf kan bewerken, hetzij met de hand, hetzij met VISIPART, zodat men verbeteringen kan aanbrengen aan de partitionering die COPPRE zelf bepaalt. COPPRE heeft dan ook de mogelijkheid om uitsluitend de partitionering te bepalen (optie PARTONLY in de COPPRE configuratie file), zodat men COPPRE kan gebruiken om een initiële partitionering te maken met een van de ingebouwde methoden en die vervolgens verder te optimaliseren. Deze manier van werken is enigszins achterhaald doordat VISIPART een enclosure uit een siminp-file kan lezen en ook automatische verdeling van (deel-)domeinen in stukken ondersteunt.

7.2.5

De 000-file Naast de SDS-files voor de deeldomeinen maakt COPPRE ook een file aan waarin allerhande informatie wordt opgeslagen die nodig is voor COPPOS. Deze informatie kan soms ook nuttig zijn voor debug doeleinden als bij het aanpassen van COPPRE blijkt dat iets niet goed (meer) werkt. De inhoud van de 000-file (evenals de structuur van een aantal arrays die COPPRE alleen op run-time gebruikt en niet opslaat) is beschreven in “Local data structure COUPLE”, die te vinden is in de sysdoc directory van de SIMONA installatie. Hier zal de inhoud alleen op hoofdlijnen geschetst worden. Op de 000-file is een viertal compounds aanwezig, waarvan twee (PARTITIONING en LDSDSCR) de interne opslag vormen van de input (respectievelijk de configuratie-file en de LDS-beschrijving). De andere twee compounds (OWNER en CONV) worden door COPPRE zelf gecreëerd. Het OWNER array bevat een specificatie van de opdeling van het rekenrooster (zoals die ook wordt geëxporteerd naar de COPPRE-R-file) en van de opdeling van de gedistribueerde indexsets. De compound CONV bevat de vertaaltabellen waarmee per deeldomein wordt aangegeven hoe de gegevens in de datastructuren gevonden kunnen worden in de structuren voor het domein als geheel.

7.2.6

De deeldomein SDS-files Tenslotte produceert COPPRE uiteraard de SDS-files voor de deeldomeinen. Dit zijn volwaardige SDS-files waarin de beschrijving van de deeldomeinen is vastgelegd alsof elk deeldomein met een eigen siminp-file was aangemaakt. Dit houdt ondermeer in dat de (m,n)-coördinaten per deeldomein tellen vanaf 1. De inhoud is dan ook vrijwel volledig conform de LDSbeschrijving van WAQUA is (en zoals die door WAQPRE wordt aangemaakt). Als enige uitzondering geldt dat er een tweetal arrays aan toegevoegd is: PCONST en POWNER. Het


105


Cursusboek Module 3

eerste array bevat een aantal constanten die soms nodig zijn bij het parallel programmeren (het totaal aantal deeldomeinen, het nummer van het eigen deeldomein, etc.). Het array POWNER is een array ter grootte van mmax nmax waarin per (m,n)coördinaat aangegeven is tot welk deeldomein het behoort.

7.3

De werking van COPPRE In deze paragraaf zal allereerst worden beschreven wat de principes zijn achter het aanmaken van arrays voor deeldomeinen. Daarna wordt getoond wat de structuur (op hoofdlijnen) van COPPRE is, waarna de verschillende elementen daaruit nader worden toegelicht.

7.3.1

Principes van COPPRE Het principe waarop COPPRE gebaseerd is, kan het best worden geïntroduceerd aan de hand van een voorbeeld. Dit voorbeeld zal tevens de belangrijkste termen en begrippen introduceren. Stel men heeft een array D met een lengte van 100 elementen. Dit array bevat de bodemdieptes op een lijn van 100 roosterpunten. De verzameling van indices [1,100] van dit array, de indexset, wordt aangeduid met een naam, ixset. Stel dat men het domein in twee delen verdeelt, zodanig dat de eerste 40 roosterpunten in het eerste domein vallen en de resterende 60 in het tweede domein. Men wil nu voor beide deeldomeinen een versie van het array D aanmaken waarvan de inhoud bestaat uit die elementen van het originele array die voor het betreffende deeldomein van belang zijn. Het originele array wordt meestal aangeduid als het globale array (d.i. globaal over alle deeldomeinen heen), en de arrays per deeldomein worden dan aangeduid als lokale arrays (dus lokale versies van het globale array). De verdeling van het domein specificeert in feite een partitionering van ixset. De partitionering is een functie op ixset, waarvan de waarde gegeven wordt door het nummer van het deeldomein waarin het element terecht komt. De partitionering is dus part(i )

1 voor i =[1,40] 2 voor i =[41,100]

Het lokale array D in deeldomein 1 zal een lengte van 40 krijgen en in deeldomein 2 een lengte van 60. De cardinaliteit van indexset ixset in deeldomein 1 is dus 40 en die in deeldomein 2 is 60.


106


Cursusboek Module 3

program coppre call coprin

Lees configuratie file

call copart

Bepaal partitionering

if(create SDS-files) call coprlm

Lees LDS-beschrijving

call copcnv

Bouw conversietabellen op

call copdlw

Creëer deeldomein SDS-files

call copwme

Creëer de 000-file

endif Figuur 10 Hoofdstructuur van COPPRE

Kijkend naar het array dat voor deeldomein 1 aangemaakt moet worden, dan bestaan de elementen daarvan uit de eerste 40 elementen van het originele array. Kortweg kan men voor deeldomein 1 schrijven do i=1,local_cardinality1(ixset) local_D(i)= global_D(globnr1(i)) enddo

waarbij local_cardinality1() de lokale cardinaliteit van ixset geeft en globnr1() voor i=1,40 de waarden 1 tot en met 40 teruggeeft. De zelfde methode kan men voor het lokale array D voor deeldomein 2 toepassen, als geven local_cardinality2() en globnr2() uiteraard andere waarden terug. Men kan dus per deeldomein voor dit array een globnr-tabel aanmaken waarin wordt vastgelegd hoe men de waarde van een element kan vinden in het globale array. De globnr-tabellen worden aangeduid met de naam conversietabellen. In feite is de tabel globnr niet zozeer gekoppeld aan het array dat gedistribueerd wordt maar aan de indexset van het array. Stel dat men nog een array heeft met dezelfde indexset, dan kan men dat met dezelfde tabel distribueren.

7.3.2

De structuur van COPPRE In Figuur 10 is de hoofdstructuur van COPPRE weergegeven. Aan het inlezen van de invoergegevens en aan het bepalen van de partitionering zal hier verder geen aandacht geschonken worden omdat het voor het begrip van de verdere werking van COPPRE niet van belang is. De overige onderdelen worden hieronder nader beschreven.


107


7.3.3

Cursusboek Module 3

Het aanmaken van conversietabellen In het zojuist behandelde voorbeeld is al aangegeven dat men de partitionering van de indexset moet kennen om een array te kunnen opdelen in arrays per deeldomein. COPPRE bepaalt feitelijk alleen de partitionering van de indexset mmax nmax en leidt de partitionering van alle andere indexsets daaruit af. In deze paragraaf zal aangegeven worden hoe dat gebeurt. De partitionering van de indexset mmax nmax is in feite een tabel HOWNER die aan elk element (m,n) van mmax nmax het nummer toevoegt van het deeldomein waartoe die (m,n) coordinaat behoort. Stel dat men nu de partitionering van de indexset mnmaxk wil weten. Zoals bekend is mnmnaxk een gecomprimeerde vorm van de indexset mmax nmax, waaruit alle landpunten verwijderd zijn. De LGRID tabel legt een relatie tussen de indexset mmax nmax en de indexset mnmaxk. Langs de omgekeerde weg geeft lgrid impliciet een relatie tussen de indexset mnmaxk en de indexset mmax nmax. Als men de omgekeerde lgrid-tabel aanduidt met invlgrid, dan kan de partitionering van mnmaxk gevonden worden door partmnmaxk(i) = howner(invlgrid(i))

Op dezelfde manier kan men de partitionering van bijvoorbeeld de bronpunten bepalen. De lokatie van bronpunten is vastgelegd in de arrays ndis en mdis: partsrc(i) = howner(mdis(i),ndis(i))

Kortom, de partitionering van elke willekeurige indexset kan worden afgeleid uit de partitionering van de indexset mmax nmax als de relatie van die indexset tot de indexset mmax nmax bekend is. In het vervolg zal deze relatie aangeduid worden met de naam relation. Dit is dus een relatie tussen de elementen van een indexset en de elementen van mmax nmax. Het bepalen van de conversietabel voor een indexset gaat nu als volgt in zijn werk. Bij elk element i van de indexset wordt via relation in partmn het nummer opgezocht van de eigenaar van dat element (d.i. het nummer van het deeldomein waarin het element thuishoort). Vervolgens wordt de lokale indexset voor die eigenaar met eentje uitgebreid en wordt er in de globnr tabel voor die eigenaar voor dat element de waarde i gezet, d.i. het nummer van het originele element in in de globale indexset. In pseudo code komt dat neer op het algoritme dat is weergegeven in Figuur 11.


108


Cursusboek Module 3

forall p: local_cardinality(p)=0 for i=1,global_cardinality owner = howner(relation(i)) local_cardinality(owner) = local_cardinaltiy(owner)+1 globnr(owner)(local_cardinality(owner))=i endfor Figuur 11 Algoritme voor het bepalen van conversietabellen

Op deze manier worden in COPPRE de conversietabellen opgebouwd en tegelijkertijd de lokale cardinaliteiten van de indexsets bepaald. Het algoritme als zodanig is generiek. Omdat de relatie relation echter in WAQUA/TRIWAQ niet op een uniforme manier geïmplementeerd is, is het stuk van COPPRE waarin deze tabellen worden aangemaakt toch nog enigzins afhankelijk van de structuur van de LDS. Er is overwogen om ook de in de LDS voorkomende relaties extern te specificeren, om zodanig COPPRE geheel generiek te maken, maar iets dergelijks is tot op heden niet geïmplementeerd. Op dit moment maakt COPPRE daarom een vast aantal conversietabellen aan op basis van ingeprogrammeerde kennis van de LDS. In de extern toegeleverde beschrijving van de LDS (in de file coplds.waqua) worden deze tabelnummers gekoppeld aan indexsetnamen door middel van het keyword MAP in de definitie van de indexsets. Als de LDS wordt uitgebreid met nieuwe indexsets dan bestaat er dus een kans dat er ook een nieuwe conversietabel moet worden toegevoegd. In dat geval moet er dus een stuk intern in COPPRE worden geprogrammeerd. Dit programmeerwerk kan beperkt blijven als de relatie is vastgelegd op een manier die vaker in WAQUA/TRIWAQ wordt gebruikt, maar kan ook lastiger zijn als de relatie in een heel ongebruikelijke vorm wordt gelegd.

7.3.4

Conversietabellen voor rijen en kolommen In het voorgaande is beschreven wat conversietabellen voorstellen en hoe ze gemaakt worden. Daarbij is er steeds van uit gegaan dat bij elk element van een indexset precies één element in mmax nmax hoort (bijvoorbeeld: een bronpunt heeft precies één (m,n)-coördinaat). Er zijn echter ook indexsets waarbij dat niet het geval is. Denk aan de indexset [1,norows], de indexset van de gridrijen. Elk element van deze indexset correspondeert met een rij van punten in mmax nmax. De relatie wordt gelegd via de irogeo tabel.


109


Cursusboek Module 3

In dit geval behoort één rij (een element van de indexset [1,norows]) niet meer noodzakelijk uniek tot één bepaald deeldomein: de rij zelf kan tot meerdere deeldomeinen behoren. In dit geval zijn er dus in meerdere deeldomeinen lokale rijen die elk een deel vormen van een bepaalde globale rij. In dit geval bevat de globnr tabel voor een deeldomein niet alleen het nummer van de globale rij waarbij een lokale rij behoort, maar ook het begin en einde van de lokale rij binnen de globale rij. Bijvoorbeeld: stel dat globale rij 3 begint in het stuk van mmax nmax dat aan deeldomein 1 is toegekend, vervolgens vanaf het vierde punt van de rij tot aan het 26e punt door deeldomein 3 loopt en tenslotte eindigt in deeldomein 2. Dan zullen de globnr tabellen van deeldomeinen 1, 2 en 3 elk een entry bevatten dat verwijst naar de globale rij 3. In de globnr tabel van deeldomein 1 staat dan aangegeven dat de lokale rij overeenkomt met punten 1 tot en met 3 van de globale rij. Zo ook staart in de globnr tabel van deeldomein 3 aangegeven dat de lokale rij overeenkomt met punten 4 tot en met 25 van globale rij 3. Deze informatie is vooral van belang tijdens de correctie-fase die hieronder beschreven wordt.

7.3.5

Het opbouwen van de deeldomein LDS Nadat COPPRE de invoer heeft gelezen en de conversietabellen heeft gecreëerd begint het programma met het opbouwen van de LDS voor elk van de deeldomeinen. Dit gebeurt in een loop over alle deeldomeinen waarin voor elk deeldomein een SDSfile wordt aangemaakt, gevuld en gesloten. Het vullen van de SDS-file doet COPPRE op basis van de beschrijving van de LDS die het uit de coplds.waqua leest. In een loop over alle compounds wordt steeds een hele compound aangemaakt en weggeschreven. Daarbij moeten ook lokale versies van de data-arrays aangemaakt worden. Dit doet COPPRE op basis van de beschrijving van de indexset van de data-array. Dat is een vaak complexe expressie van Cartesische produken en concatenaties van allerlei indexsets. Sommige van die indexsets zullen gedistribueerd zijn (d.i. gekoppeld zijn aan een conversietabel) en anderen zullen gerepliceerd zijn. In elk geval kan COPPRE op basis van de lokale cardinaliteiten van de indexsets met behulp van de beschrijving van de indexset van een data-array de lokale lengte van het array bepalen. Zodra die lengte bekend is, kan het array worden gealloceerd. Het vullen van het array gebeurt vervolgens door bij elk element ervan, middels de beschrijving van de indexset en de


110


Cursusboek Module 3

aanwezige conversietabellen het overeenkomstige element in de originele, globale array te vinden en de waarde daarvan te kopiëren.

7.3.6

Correcties Boven is beschreven hoe met behulp van conversietabellen en de beschrijving van de LDS automatisch een deeldomein versie van de LDS kan worden gegenereerd. Daarbij worden waarden uit de originele, globale arrays rechtstreeks gekopiëerd naar de deeldomein arrays. Voor de meeste gegevens in de LDS van WAQUA/TRIWAQ is dat voldoende. Maar het zal onmiddelijk duidelijk zijn dat er ook waardes aangepast moeten worden. Denk bijvoorbeeld aan de waarde van mnmaxk in MESH_IDIMEN. De waarde daarvan is in elk van de deeldomeinen anders dan in het volledige domein. Daarom wordt er helemaal op het einde van COPPRE nog een correctieslag uitgevoerd, waarin dit soort aanpassingen plaatsvinden. Voor een deel zijn deze aanpassingen onafhankelijk van WAQUA/TRIWAQ, voor een deel zijn ze ook heel specifiek. Dit is dus een tweede stuk (naast het aanmaken van de conversietabellen) van COPPRE dat minder generiek is dan het mogelijk zou kunnen zijn. Ook in dit deel kunnen aanpassingen nodig zijn als de LDS wordt uitgebreid met een nieuwe datastructuur, die speciale conversies vereist. Zoals gezegd is een deel van de correcties niet erg inhoudelijk van aard. Ten eerste is in de LDS-beschrijving (file coplds.waqua) beschreven waar de cardinaliteiten van de indexsets zijn opgeslagen. Bij de indexset mnmaxk is bijvoorbeeld vastgelegd dat deze de cardinaliteit ervan is opgeslagen in MESH_IDIMEN(5). Op basis van deze informatie kan COPPRE zonder verdere informatie zien dat MESH_IDIMEN(5) moet worden aangepast aan de lokale cardinaliteit van de indexset mnmaxk. Omdat de (m,n)-coördinaten van de deeldomeinen beginnen bij 1, moeten ook de (m,n)-coördinaten van alle bronpunten, barriers etc. aangepast worden. Deze coördinaatarrays zijn precies die arrays die ook gebruikt zijn bij het aanmaken van de conversietabellen. Deze correctie is dus sterk gerelateerd aan het WAQUA/TRIWAQ specifieke stuk van COPPRE bij het aanmaken van de conversietabellen. Bij de indexsets van rijen (en soortgelijke structuren zoals kolommen en crossecties) moet er bovendien nog een extra correctie uitgevoerd worden. De (m,n)-coördinaten in de irogeo tabel moeten niet alleen aangepast worden aan de verschuiving van de oorsprong van de (m,n)-coördinaten, maar


111


Cursusboek Module 3

bovendien aan het feit dat de lokale rij slechts een deel van de globale rij vormt. Het zou te ver voeren om alle aanpassingen die gedaan worden hier in detail te beschrijven. Hoewel het een relatief groot deel van de programmacode van COPPRE betreft, gaat het uiteindelijk slechts om een zeer klein deel van de LDS waarop correcties moeten worden uitgevoerd.

7.4

Toevoegen van nieuwe data structuren Als men aan WAQUA/TRIWAQ een nieuw array wil toevoegen en dat array ook in de parallelle versie wil gebruiken, zal men te maken krijgen met COPPRE. In het gunstigste geval (en waarschijnlijk ook de meeste gevallen) kan men volstaan met het uitbreiden van de beschrijving van de LDS in coplds.waqua. Als men daarbij echter een nieuwe indexset nodig heeft, en als deze indexset een ruimtelijke interpretatie heeft (dat wil zeggen: op de een of andere manier gekoppeld aan de indexset mmax nmax) dan zal men echter ook aan COPPRE zelf moeten programmeren. In dat geval zal men een nieuwe conversietabel moeten introduceren. De programmatuur is op dit punt zodanig opgezet dat het betrekkelijk gemakkelijk zou moeten zijn om een nieuwe conversietabel aan te maken. Bovendien zal men dan wellicht toevoegingen moeten doen aan de correctieroutines. Voor beide is het handig om de relatie tussen de nieuwe index set en de set mmax nmax op te slaan op een manier die vaker gebruikt wordt in WAQUA/TRIWAQ, omdat men dan gebruik kan maken van bestaande programmastructuren.


112


Cursusboek Module 3

8

Interne werking van de collector COPPOS

8.1

Functies van COPPOS Het programma COPPOS (ook wel de collector genoemd) is in vergelijking met COPPRE functioneel gezien een relatief eenvoudig programma. Het heeft slechts één functie en dat is het samenvoegen van de resultaten die elk van de rekenprocessen naar zijn eigen SDS-file geschreven heeft. De samengevoegde resultaten worden opgeslagen in de SDS-file waaruit COPPRE aanvankelijk de deeldomein SDS-files heeft gegenereerd. De methode waarop COPPOS gebaseerd is, is feitelijk het omgekeerde van de methode waarop COPPRE gebaseerd is. COPPOS maakt dan ook uitgebreid gebruik van de tabellen die COPPRE bepaald heeft en in de SDS-runid-000 file heeft opgeslagen. De enige functie die COPPOS extra heeft ten opzichte van COPPRE is gelegen in de interpolatie van gegevens. Dit speelt uitsluitend een rol bij domein decompositie met vertikale verfijning, waarbij deeldomeinen met verschillende aantallen lagen worden samengevoegd tot een deeldomein met het maximaal voorkomende aantal lagen. Stel dat men bijvoorbeeld een experiment heeft uitgevoerd waarin het ene deel van het domein met 5 lagen is doorgerekend en het andere deel van het domein met 10 lagen, dan wordt een SDS-file gecreëerd met 10 lagen, waarin gegevens uit het gebied waarin 5 lagen zijn gebruikt naar 10 lagen toe geïnterpoleerd worden. Dit geeft in het gebied waar 5 lagen gebruikt zijn een resolutie die hoger lijkt dan er feitelijk bereikt is. Anderzijds heeft dit wel als voordeel dat de normale postprocessing tools gebruikt kunnen blijven worden. Ten behoeve van deze interpolatie bevat de LDS-beschrijving die door COPPRE gelezen wordt bij enkele indexsets die aan KMAX zijn gerelateerd informatie over de wijze waarop die indexset geïnterpoleerd moet worden. Bij domein decompositie met horizontale verfijning wordt COPPOS per globaal domein gedraaid, om per domein een resultaat SDS-file te genereren met de gegevens uit alle deeldomeinen die uit dat domein gegenereerd zijn.

8.2

Het collectie algoritme Zoals al werd opgemerkt is de methode waarop COPPOS gebaseerd is min of meer het omgekeerde van de methode die COPPRE gebruikt. Op hoofdlijnen werkt COPPOS zoals wordt weergegeven in Figuur 12. Van binnenuit redenerend wordt steeds een data-


113


Cursusboek Module 3

vul de tijdsonafhankelijke interpolatietabellen voor alle compounds voor alle tijdstippen waarvoor de compound bestaat vul de tijdsafhankelijke interpolatietabellen lees de lokale compound in uit alle deeldomein SDS-files lees de globale compound uit de doel-SDS-file voor alle data-arrays van de compound zet alle elementen van het globale data-array op nul voor alle deeldomeinen als interpolatie nodig creëer een tijdelijke versie van het locale array met het aantal lagen waarheen geinterpoleerd wordt interpoleer deeldomeingegevens naar het tijdelijk array voeg het resultaat in in de globale data-array vernietig het tijdelijk array anders voeg direct in in de globale data-array eind als einde loop voor alle deeldomeinen einde loop voor alle data-arrays van de compound schrijf het globale array weg naar de doel-SDS-file einde loop voor alle tijdstippen waarvoor de compound bestaat einde loop voor alle compounds voer correcties uit Figuur 12: het samenvoegingsalgoritme van COPPOS

array van de doel-SDS-file gevuld met de resultaten van elk van de deeldomeinen totdat alle data-arrays van een compound voor een bepaald tijdstip gevuld zijn. Op dat moment wordt de compound voor het betreffende tijdstip naar de SDS-file geschreven en wordt eventueel dezelfde compound voor een volgend tijdstip gevuld met waarden uit de deeldomein SDSfiles. Merk op dat er dus steeds uitgegaan wordt van een array van de globale SDS-file, die op nul gezet wordt en vervolgens aangevuld met de waarden uit de deeldomein SDS-files. Deze procedure leidt ertoe dat delen van een globaal array die niet gevuld worden vanuit deeldomein SDS-files op waarde nul blijven staan. In Figuur 12 is ook aangegeven hoe de interpolatie in het algoritme geplaatst is. Er wordt gebruik gemaakt van interpolatietabellen. Deze tabellen geven in compacte vorm de uit te voeren interpolatie weer en zijn feitelijk dezelfde tabellen als die gebruikt worden in COCLIB voor het interpoleren van gegevens tijdens de communicatie tussen deeldomeinen. Interpolatietabllen kunnen tijdsonafhankelijk zijn, of Versie 1.0, oktober 2002

114


Cursusboek Module 3

tijdsafhankelijk. Bij tijdsafhankelijke interpolatietabellen hangen de interpolatiegewichten af van de waterstand op het moment waarvoor geïnterpoleerd moet worden. Dit speelt een rol in situaties waarbij een sigma-laag in het ene domein verbonden is met een vaste laag in een ander domein. Voor deeldomeinen waarvoor interpolatie nodig is, worden eerst de gegevens geïnterpoleerd naar het gewenste aantal lagen alvorens ze worden ingevoegd in de globale array. Op deze manier zijn het interpoleren en het samenvoegen gescheiden. Het invoegen van gegevens uit een deeldomein SDS-file in de globale data-array is een bewerking die sterkt lijkt op de bewerking die COPPRE uitvoert om lokale data-arrays te vullen. Kortweg komt deze operatie neer op: do i=1,local_cardinality1(ixset) global_array(globnr1(i))= local_array(i) enddo

waarbij weer dezelfde globnr (conversie-)tabellen gebruikt worden als in COPPRE. Net als COPPRE kent ook COPPOS een correctiefase na afloop van het samenvoegen, al is deze veel beperkter dan in COPPRE. Bij COPPOS wordt de correctiefase gebruikt om een kopie van het array HOWNER (zie COPPRE) aan de SDS-file toe te voegen, zodat men achteraf kan zien welk deel van het domein tot welk deeldomein behoorde en de arrays KHU en KHV te repareren, waarin na samenvoeging van de deeldomein arrays nog het aantal lagen is terug te vinden dat in elk van de deeldomeinen gebruikt is. Deze waardes worden vervangen door het aantal lagen van de resultaat-SDS-file.

8.3

Speciale gevallen In principe hoeft COPPRE alleen tijdsafhankelijke arrays samen te voegen. Tijdsonafhankelijke arrays worden immers aangemaakt door WAQPRE en blijven onveranderd tijdens de run. Er is echter één uitzondering: de compound-array HARMONIC_TIDE die op runtime wordt aangemaakt en die dus alsnog gecollecteerd moet worden. In dit geval maakt COPPOS eerst de compound(-structuur) aan met daaronder alle data-arrays op de goede plek en met de goede lengte, voordat met het daadwerkelijke collecteren van gegevens van de deeldomein-SDS-files wordt begonnen. Voor het overige kan COPPOS dit array dan vrijwel hetzelfde behandelen als tijdsafhankelijke arrays. COPPRE voegt de deeldomein SDS-files samen voor alle tijdstippen, met uitzondering van het eerste tijdstip. Dit omdat


115


Cursusboek Module 3

het eerste tijdstip (de initiële waarde) altijd door WAQPRE is aangemaakt en dus al goed op de globale SDS-file staat.


116


9

Case-study: parallelisatie horizontaal k- model

Cursusboek Module 3

van

het

In deze cursushandleiding is een introductie gegeven in parallel rekenen in het algemeen en in de manier waarop dit wordt gebruikt binnen WAQUA/TRIWAQ. Daarbij zijn allerlei aspecten aan de orde geweest die van belang zijn voor de verdere ontwikkeling van parallel WAQUA/TRIWAQ, zoals de strategie voor het parallelliseren van de verschillende numerieke solvers, de manier waarop data-analyse wordt toegepast, en de interne werking van de hulpprogrammatuur COCLIB, COEXEC, COPPRE en COPPOS. In dit hoofdstuk wordt een en ander toepasbaar gemaakt voor de programmerende gebruiker die nieuwe functionaliteit wil toevoegen aan parallel WAQUA/TRIWAQ. Hiervoor wordt de toevoeging van het horizontaal k- turbulentie-model als voorbeeld gebruikt. Deze functionaliteit is door WL/RIKZ ontwikkeld voor sequentiële versies van TRIWAQ en is door VORtech Computing geparallelliseerd. In de rapportage van dit laatste project (VORtech Computing TR01-10, “Parallellisatie van het horizontaal k- turbulentie-model”) wordt een gedetailleerde beschrijving gegeven van de stappen die zijn uitgevoerd voor de parallellisatie en van de manier waarop verschillende keuzes worden gemaakt. Het rapport TR01-10 is bijgevoegd bij deze cursushandleiding.


117


Cursusboek Module 3

Referenties (Zijlema, 1998)

M. Zijlema, TRIWAQ – three dimensional shallow water flow model, Technisch Rapport SEPRA in opdracht van RIKZ.

(Vollebregt, 1997)

E.A.H. Vollebregt, Parallel Software Development Techniques for Shallow Water Models, Proefschrift, Technische Universiteit Delft, 1997.

(Vollebregt en Roest, 1998)

E.A.H. Vollebregt en M.R.T. Roest, Prospects of the Parallellisation of WAQUA, Technisch Rapport TR98-02, VORtech Computing, Postbus 260, 2600 AG Delft

(Vollebregt et al. 2001)

E.A.H. Vollebregt, M.R.T. Roest en B. van ’t Hof, Detailontwerp domein decompositie met horizontale verfijning, Technisch Rapport TR01-06, VORtech Computing, Postbus 260, 2600 AG Delft


118

TRIWAQ Cursusboek Module 3

Recommend Documents