Triwaq

VORtech Computing

Experts in Technisch Rekenwerk Postbus 260 2600 AG DELFT

MEMO Datum Auteur(s) Onderwerp

EV/M08.084 4 december 2008 Edwin Vollebregt

tel. 015-285 0125 fax. 015-285 0126 [email protected]

Onderzoek van ux-limiters voor het transportgedeelte van Waqua/Triwaq

Documentinformatie Versie Auteur

Datum

Opmerkingen

0.1 EV 0.2 EV Bestandslokatie:

24-11-2008 Eerste test-resultaten 04-12-2008 Implementatie in moederversie /v3/E05q bo simona/c87532-neg-transp-conc/report

Review

1 Inleiding Er zijn recent twee verschillende problemen opgetreden met negatieve concentraties in het transport-gedeelte van Waqua/Triwaq: 1. Het testmodel van Bas van 't Hof voor de OpenMI-koppeling (change c81249) levert negatieve zout-concentraties bij het einde van de Nieuwe Waterweg. 2. Bij het inspelen van het nieuwe CSM-model treden negatieve concentraties op. Deze negatieve concentraties blijken overshoots van het advectie-schema van het transportgedeelte te zijn, die optreden bij steile gradienten in het concentratie-veld. De vraag is of er verbeteringen aan de numerieke schema's mogelijk zijn met betrekking tot deze under/overshoots. Met name wordt hiervoor aan zogenaamde \ ux-limiters" gedacht. In dit memo beschrijven we drie eenvoudige test-cases voor het onderzoeken van discretisatieschema's voor het transportgedeelte van Waqua/Triwaq. Vervolgens laten we de eerste resultaten zien die met een \min-mod limiter" zijn bereikt.

2 Test-cases voor onderzoek van transport-schema's Numerieke schema's voor de advectieve termen van het transportgedeelte van Waqua/Triwaq kunnen eenvoudig gesoleerd worden onderzocht. Hiervoor gebruiken we een 1D kanaal met de volgende eigenschappen:

VORtech Computing Memo EV/M08.084

versie 0.2, 4 december 2008

2

schema formule voor concentratie op cell-interface e 2 orde centraal cm+1=2 = (cm + cm+1 )=2 1e orde upwind, qx > 0 cm+1=2 = cm 2e orde upwind, qx > 0 cm+1=2 = (3cm cm 1 )=2 3e orde upwind, qx > 0 cm+1=2 = (10cm 5cm 1 + cm 2 )=6 Tabel 1: Ruimtelijke discretisatie-schema's voor advectieve ux in transportgedeelte Waqua/Triwaq impliciet, trsadx expliciet, trsady Waqua 2e orde upwind 2e orde centraal Triwaq 3e orde upwind 2e orde centraal Tabel 2: Combinaties gebruikt in cyclische schema's in Waqua/Triwaq

Waterdiepte 1 m, een constante waarde valt weg uit de transportvergelijking;

roosterafstand x = y = 60 m, het Courant-getal voor het transportgedeelte is hiermee gelijk aan de tijdstap in minuten;

geen diusie in het transportgedeelte, pure advectie;

simulatie van een blok-pro el op m = 20 : 60 (siminp.bak1): direct onderzoeken van schokken, zowel positief als negatief. De maximale waarde is c = 1 kg=m3 , procentuele afwijkingen kunnen direct worden afgelezen.

simulatie van een glad cos-vormig pro el c = 0:5 (1 + cos(2 [0 : n]=n)), per periode, n = 60 in siminp.bak2 en n = 20 in siminp.bak3.

stroomsnelheid 1 m=s; praktisch geen bodemwrijving (Chezy-waarde 20:000 m1=2 =s), levert bij vlakke bodem een vlakke waterstand op;

beginconcentratie op roosterpunten 10 t/m 70, simulatie van twee uur (120 roosterpunten), rooster van 200 punten lang: geen eecten van randafhandeling.

n

punten

3 Initiele testresultaten Figuur 1 toont de eerste resultaten voor testcase 1, simulatie van een blok-pro el, bij een tijdstap (Courant-getal) van 0.5. In alle gevallen wordt het ADI-tijdsintegratieschema van Waqua/Triwaq gebruikt. De discretisaties voor de impliciete stap worden opgesteld in subroutine trsadx, voor de expliciete stap in trsady. De gebruikte schema's worden beschreven in Tabellen 1 en 2.



Pure advectie van blok−profiel, Dt=0.5


Exact centraal 1e orde upwind

1.2

1

1

0.8

0.8 Concentratie [kg/m3]

Concentratie [kg/m3]

1.2

0.6

0.4

0.4

0.2

0

0

−0.2

Exact 2e orde upwind 3e orde upwind

0.6

0.2

−0.2 80

100

120

140 Gridpunt [−]

160

180

200

90

100

110

120


1

1

140 150 Gridpunt [−]

160

170

180

190

200

170

180

190

200

Exact Waqua met limiter centraal met limiter 2e orde upwind met limiter Triwaq met limiter

0.8



130


Exact Waqua Triwaq

0.8

0.6

0.4

0.6

0.4

0.2

0.2

0

0

90

3

100

110

120

130


160

170

180

190

200

90

100

110

120

130


160

Figuur 1: Advectie van een blok-pro el, Courant-getal 0.5. De resultaten in Figuur 1 (links-boven) laten zien dat een gewone centrale discretisatie bij schokken sterke wiggles produceert en dat een eerste orde upwindschema veel arti ciele diusie oplevert. De tweede en derde orde upwind-schema's (rechts-boven) geven een overshoot die vooruit loopt op de schok. De schema's van Waqua en Triwaq (links-onder) geven veel kleinere overshoots, maar nog steeds bijna 10% van de grootte van de schok. Via min-mod ux-limiters (rechts-onder) kunnen deze overshoots worden voorkomen. De gebruikte ux-limiters zijn gebaseerd op de volgende slope-ratio's: c+ rm +1=2

= c rm +1=2 = u1+ rm+1=2 = u1 rm +1=2 =

( cm (cm+2 (cm+1 (cm+1

1 )=(cm+1 cm ) cm+1 )=(cm+1 cm ) cm )=(cm cm 1 ) cm )=(cm+2 cm+1 )

cm

(1) (2) (3) (4)



u2+ rm +1=2

= (cm+1 u2 rm +1=2 = (cm+1

cm )=(cm

2) cm+1 )

(5) (6)

cm

cm )=(cm+3

4

De de nitie van min-mod limiterfunctie is als volgt: (r )

= max(0; min(1; r))

(7)

Flux-limiter voor tweede orde centrale-schema: 0: u<0: u>

c+ = (rm +1=2 ); c 1 = (rm +1=2 );

1

cm+1=2

= ((2 1 ) cm + 1 cm+1 ) =2 cm+1=2 = (1 cm + (2 1 ) cm+1 ) =2

(8) (9)

Flux-limiter voor tweede orde upwind-schema: 0: u<0:

u>

u1+ = (rm +1=2 ); cm+1=2 = ((2 + 1 ) cm u1 1 = (rm +1=2 ); cm+1=2 = ((2 + 1 ) cm+1

1

1 ) =2 1 cm+2 ) =2

1 cm

(10) (11)

Flux-limiter voor derde orde upwind-schema: u2+ u1+ = (rm +1=2 ); 2 = (rm+1=2 ); cm+1=2 = ((6 + 31 + 2 ) cm + ( 31 22 ) cm 1 + 2 cm 2 ) =6 u2 u1 u<0: 1 = (rm +1=2 ); 2 = (rm+1=2 ); cm+1=2 = ((6 + 31 + 2 ) cm+1 + ( 31 22 ) cm+2 + 2 cm+3 ) =6

u>

0:

1

(12) (13) (14) (15)

De limiter-functie gebruikt het hoge orde schema bij = 1 en het lage orde schema bij = 0. Als de twee slopes die vergeleken worden verschillend teken hebben dan wordt het lage orde schema gebruikt. Als de slope verder weg van het discretisatiepunt kleiner is dan rond het discretisatiepunt zelf dan wordt het hoge orde schema gebruikt. Als de slope toeneemt met afstand dan wordt de verhouding gebruikt om te voorkomen dat er over/undershoot plaatsvindt. Het is niet geheel duidelijk waarom (cm+1 cm ) in de slope-ratio's voor het centrale schema in de teller staat en bij de upwind-schema's in de noemer moet staan. Dat heeft waarschijnlijk te maken met de neiging van de schema's om te over- of te onderschatten. De gebruikte limiter-functies zijn min of meer experimenteel bepaald, er is nog niet gekeken naar een nette/rigoreuze a eiding, bijvoorbeeld in overeenstemming met [1]. Er wordt een beveiliging gebruikt om deling door nul te voorkomen. Als deling door nul dreigt op te treden dan kan steeds het lagere orde schema worden gebruikt (?). De formules voor het derde orde upwind-schema zijn door onszelf geconstrueerd. Het is onduidelijk of hier nog betere benaderingen mogelijk zijn. Op Wikipedia staat een tiental verschillende limiter-functies. Deze leiden tot verschillende schema's. Het is onduidelijk welke limiter-functies in welke situaties toepasbaar zijn, welke



eisen ze stellen aan hoge- en lage-resolutie benaderingen, welke slope-ratio's er moeten worden ingevuld, etc. Hiervoor zouden we de literatuur moeten raadplegen of bij kennissen moeten langsgaan, met name bij het CWI. De schema's lijken niet sterk afhankelijk te zijn van de tijdstap die wordt gebruikt. Dit moet nog verder worden geveri eerd, met name voor de schema's met limiter waar de limiter-functie expliciet wordt geevalueerd.

4 Resultaten voor gladde pro elen In de vorige paragraaf is aangetoond dat ux-limiters negatieve concentraties kunnen voorkomen, terwijl de arti ciele diusie binnen de perken blijft. In deze paragraaf kijken we naar het eect voor gladdere pro elen. Hiervoor worden de test-cases siminp.bak2 en siminp.bak3 gebruikt, zie Figuur 2. De guur links-boven laat de berekende oplossingen zien voor het centrale schema en het eerste orde upwind schema voor het pro el met 20 gridpunten per periode. Deze leiden aan wiggles of enorme arti ciele diusie. De guur rechts-boven is voor de Waqua- en Triwaqschema's. Die behouden veel van het oorspronkelijke pro el en geven een beperkte over-shoot. Links-midden toont de resultaten voor schema's met limiters. Hierin zijn de overshoots weg maar komt daarvoor aardig wat inzakken van het pro el in de plaats. De mate van inzakken wordt getoond in de guur rechts-midden, waar het verschil met de analytische oplossing wordt getoond. De onderste twee guren zijn tenslotte voor het testgeval met 60 gridpunten per periode. Hierin is ook een kleinere tijdstap gebruikt (Courant-getal 0.5 in plaats van 1.25). Links-onder worden de fouten voor de Waqua- en Triwaq-schema's getoond, rechts-onder voor de schema's met ux-limiters. De ux-limiters doen het hier al een stuk beter dan bij 20 gridpunten per periode (rechts-midden), maar nog lang niet zo goed als de oorspronkelijke schema's.

5 Implementatie in de moederversie van SIMONA De hierboven beschreven experimentele schema's zijn in de moederversie van SIMONA opgenomen.

In de siminp- les zijn onder TRANSPORT - PROBLEM - METHODVARIABLES nieuwe keywords toegevoegd: ADVEC SCHEME, PARAM A, PARAM B en PARAM C.

Het eerste nieuwe keyword is het nummer van het schema dat moet worden gebruikt. De default is -1 voor de oorspronkelijke schema's. Waardes 0 t/m 5 betreen respectievelijk het Waqua-schema, 2e orde centraal, 1e orde upwind, 2e orde upwind, 3e orde upwind en het Triwaq-schema. Hierbij worden waardes 10, 20 of 30 opgeteld voor het aanzetten van de ux-limiters in respectievelijk de impliciete stap, de expliciete stap of beide stappen.

5



Pure advectie van cos−profiel met 20 gp, Dt=1.25


Exact centraal 1e orde upwind

1

0.8

0.8

0.6

0.6



1

6

0.4

0.2

Exact Waqua Triwaq

0.4

0.2

0 0

−0.2 90

100

110

120

130


160

170

180

190

−0.2 90

200

100

110


0.9

130


160

170

180

190

200

180

190

200

180

190

200

Pure advectie van cos−profiel met 20 gp, Dt=1.25 0.4

Exact Waqua met limiter centraal met limiter 2e orde upwind met limiter Triwaq met limiter

1

120

0.3

Waqua met limiter centraal met limiter 2e orde upwind met limiter Triwaq met limiter

0.8 0.2



0.7 0.6 0.5 0.4

0.1

0

−0.1

0.3 −0.2

0.2 0.1

−0.3

0 90

100

110

120

130


160

170

180

190

−0.4 90

200

100

110


0.04

0.04

0.03

0.03

0.02

0.02 Concentratie [kg/m3]

3


160

170

0.05 Waqua Triwaq

Concentratie [kg/m ]

130


0.05

0.01

0

−0.01

0

−0.01

−0.02

−0.03

−0.03

−0.04

−0.04

100

110

120

130


160

170

180

190

200

Waqua met limiter centraal met limiter 2e orde upwind met limiter Triwaq met limiter

0.01

−0.02

−0.05 90

120

−0.05 90

100

110

120

130


160

170

Figuur 2: Resultaten voor advectie van een cos-pro el met 20 gridpunten per periode en Courant-getal 1.25 (boven, midden) of 60 gridpunten per periode en Courant-getal 0.5 (onder). De eerste drie guren geven berekende oplossingen weer, de andere drie guren het verschil met de analytische oplossing.



De andere keywords zijn parameters die in de nieuwe schema's kunnen worden gebruikt. Op dit moment worden ze nog niet gebruikt.

De parameters worden ingelezen in Waqpre en op de SDS- le gezet. In Waqpro worden ze uitgelezen en doorgegeven naar de relevante subroutines trsadx en trsady.

De nieuwe schema's zijn alleen in TRIWAQ-berekeningen toegestaan, vanwege het gebruik van een zeven-breed stencil in de solver die in TRIWAQ wordt gebruikt.

De nieuwe schema's zijn volledig naast de oude schema's gemplementeerd. Ze hebben daardoor geen invloed op performance en resultaten van bestaande schematisaties. Wel leiden de aanpassingen tot andere compiler-optimalisaties en daarmee tot eecten van gevoeligheid.

De performance van de nieuwe schema's is nog niet geoptimaliseerd. Zo zijn alle schema's gemplementeerd in een enkele loop, en worden constante deel-expressies in iedere loopdoorgang opnieuw geevalueerd.

In de nieuwe schema's is nog geen aandacht besteed aan randafhandeling, noch voor open noch voor gesloten randen. Ze zijn vooralsnog beperkt tot schematische kanaaltjes zoals waarover hiervoor is gerapporteerd.

Referenties [1] W. Hundsdorfer, B. Koren, M. van Loon, and J.G. Verwer. A positive nite-dierence advection scheme. J.Comput.Phys, 117:35{46, 1995.

7

Triwaq

Recommend Documents