II.3.  Diszkretizáció

Az előrejelzési feladat analitikus megoldása nem ismert, ezért a hidro-termodinamikai egyenletrendszert numerikus módszerek segítségével oldjuk meg. A folytonos egyenleteket diszkretizáljuk, azaz egy háromdimenziós rács rácspontjaiban tekintjük a meteorológiai állapothatározókat, s az előrejelzést (a modellintegrálást) az időtáv időlépcsőkre osztásával, lépésenként készítjük el. A diszkretizációval kapcsolatban az egyik legfontosabb kérdés, hogy az egyenletekben szereplő térbeli és időbeli differenciálást milyen numerikus módszerekkel végezzük el. A differenciál-operátorok közelítésére két módszercsaládot alkalmazhatunk: (i) a véges különbséges módszereknél a deriváltakat az állapothatározók adott időpontbeli illetve rácspontbeli értékeinek segítségével állítjuk elő, míg (ii) a Galjorkin módszer-család (Haltiner és Williams, 1980) esetében a meteorológiai változókat olyan ún. bázis függvények segítségével írjuk fel, amelyek révén a differenciálás analitikusan elvégezhető. A Galjorkin-módszerek két típusa a véges elem és a spektrális módszer, és többnyire globális modellekben alkalmazzák őket, míg a véges különbséges módszert elsősorban korlátos tartományú modellekben használják. A jegyzet jelen fejezetében először a véges különbséges módszerek, majd a spektrális technika részleteibe nyújtunk betekintést.

II.3.1. Véges különbséges közelítés

Szépszó Gabriella

A vizsgált parciális differenciálegyenletek

Elsőként áttekintjük a lineáris parciális differenciálegyenletek típusait. Ugyan az előrejelzési feladat alapját bonyolultabb, nem-lineáris parciális differenciálegyenlet-rendszer képezi, mégis lényeges ismerni az egyszerűsített lineáris egyenletek megoldására alkalmazott módszereket, mert hasonlóakat használunk az összetett problémákon is. Kiindulásként tekintsük a vertikális irányban homogén, összenyomhatatlan, súrlódásmentes, forgó folyadék egyenleteit (Vreugdenhil, 1994):

(II.26.)

ahol u(xyt) és v(xyt) a horizontális áramlási sebesség komponensei, h(xyt) a hullám magassága, g a gravitációs gyorsulás,  f a Coriolis-paraméter és a horizontális divergencia. A fenti egyenletek olyan folyadék mozgását írják le, melynek vertikális kiterjedése elhanyagolható a horizontális méretéhez képest, ezért ezeket sekélyvíz-egyenleteknek nevezik. A sekélyvíz közelítés jól használható a nagytérségű légköri folyamatok leírásánál, pl. a Rossby-hullámok esetében. Példáinkban a fenti egyenletek egyszerűsített, lineáris változatain fogjuk bevezetni és vizsgálni a véges különbséges módszereket, mégpedig:

  • Az egydimenziós lineáris advekciós egyenleten, amelyben az advekciós sebesség konstans (c):

(II.27.)

  • Az egydimenziós lineáris gravitációs hullám-egyenleten, amelyben H a folyadék átlagos vastagsága (az óceán esetében 4 km, a légkör esetében 10 km):

(II.28.)

  • Az advekciót és a gravitációs hullám-tagokat egyaránt tartalmazó feladaton:

(II.29.)

Az egyenletekben kétfajta differenciál-operátor szerepel, az idő és a hely szerinti derivált, s ezek sajátosságai – mint az alábbiakban látni fogjuk – a diszkretizálásukra alkalmazott véges differencia sémák jellegében is megmutatkoznak.

Véges differencia sémák az időbeli és térbeli deriváltak közelítésére

A diszkretizációhoz először is definiáljunk egy rácsot, az egyszerűség kedvéért ekvidisztáns Δx rácsfelbontással, továbbá osszuk fel az időtávot Δt hosszúságú időlépésekre. A folytonos feladathoz konstruálunk egy olyan feladatot, melyet ezen a diszkrét téren értelmezünk (a továbbiakban ezt hívjuk diszkrét vagy véges differencia feladatnak). Ahogyan már említettük, a véges különbséges diszkretizáció során adott változó időbeli és térbeli deriváltjait a változó különböző szempontok szerint kiválasztott időlépcsőkben illetve rácspontokban felvett értékei alapján számítjuk. Az egyenletekben szereplő térbeli differenciál-operátoroknak az adott j rácspontra vonatkozó diszkretizációjára a legelterjedtebben a következő módszereket alkalmazzuk:

  • Bal oldali séma:

(II.30.)

  • Jobb oldali séma:

(II.31.)

  • Középponti vagy centrált séma:

(II.32.)

A sémák értelemszerűen arról kapták a nevüket, hogy a j-edik pontbeli derivált kiszámításához mely pozíciójú rácspontokat használjuk fel.

A numerikus megoldással kapcsolatban elvárjuk, hogy a diszkretizált operátorok elegendően nagy pontossággal közelítsék a folytonos operátort. A diszkrét és folytonos feladat közelségét a konzisztencia adja meg: a véges differencia feladat konzisztens a folytonos feladattal, ha mellett a két feladat különbségéből adódó maradéktagtart a 0-hoz (ezt csonkítási hibának is nevezik, de ne keverjük össze a spektrális módszernél alkalmazott csonkítással). A konzisztencia rendjét a csonkítási hiba vezető tagjának fokszáma adja meg, s ez minél magasabb, annál pontosabban közelíti az adott véges differencia séma a folytonos problémát. Ha elvégezzük a konzisztencia vizsgálatát (ld. az alfejezet végén az 1. feladatot), akkor megállapíthatjuk, hogy a középponti séma magasabb (másod-) rendű pontossággal közelíti a folytonos differenciál-operátort, mint a bal és jobb oldali sémák.

Az időbeli deriválás esetében tekintsük a következő feladatot:

(II.33.)

A fenti ún. Cauchy-probléma (azaz a differenciálegyenlet és a kezdeti feltétel együttese) esetében a feladat tulajdonképpen u(t) előrejelzése a kiindulási állapot ismeretében. Az egyenletben szereplő differenciál-operátor közelítésére explicit vagy implicit sémákat használhatunk. Az explicit sémák az adott időlépcsőbeli u meghatározásához csak ismert időlépcsőbeli értékeket használnak fel, míg az implicit sémák a még nem ismert időlépésekből is felhasználnak értékeket. Az alábbiakban néhány egyszerű véges differencia sémát mutatunk az időbeli derivált diszkretizációjára (ahol n az n-edik időlépést jelöli):

  • A legegyszerűbb az explicit Euler- vagy Euler-forward séma, aminél a forward elnevezés arra utal, hogy a séma a következő időlépcsőbeli értéket az ismert értékekből, tehát időben előre határozza meg:

(II.34.)

  • Implicit Euler-séma vagy Euler-backward séma:

(II.35.)

  • Leapfrog-séma:

(II.36.)

  • Másodrendű implicit séma:

(II.37.)

A felsorolt véges differencia módszerek közül a két Euler-séma elsőrendű, míg a leapfrog- és az utolsó, implicit módszer másodrendű (tehát nagyobb) pontosságot biztosít. Az Euler-forward és a leapfrog-módszer explicit sémák, (II.35.) és (II.37.) esetében viszont un+1 meghatározásához felhasználjuk az  f-nek az (+ 1)-edik időlépésben felvett értékét. Ahhoz, hogy un+1-et ilyen módon ki tudjuk számítani, operátor-invertálást vagy iterációt kell alkalmazni. Ez  f alakjától függően bonyolulttá is teheti az implicit sémák használatát, alkalmazásuk esetenként mégis előnyös lehet a számítási hatékonyság szempontjából (erre később még visszatérünk).

A leapfrog módszer az (+ 1)-edik időlépésbeli u meghatározásához nemcsak az n-edik, de a két lépéssel korábbi (– 1)-edik időlépésből is felhasznál információkat, azaz a módszer egy három-időszintes séma. Erre utal az elnevezése is: az aktuális u kiszámításához a kettővel korábbi időlépcsőből használja fel az u értékét és a megelőző időlépcsőbeli értéket csak az  f-en keresztül veszi figyelembe („átugorja” a megelőző időlépcsőbeli u-t), így a páros és páratlan időlépcsőbeli értékek gyengén vannak csatolva. Ebből következően a leapfrog-sémát csak a második időlépés után lehet teljes egészében alkalmazni, s az első időszint értékeit más numerikus séma segítségével kell meghatározni – erre használhatunk például egylépéses Euler-módszert (a leapfrog séma ezen tulajdonságából eredő következményeket később tárgyaljuk).

Konvergencia és stabilitás

A numerikus megoldással kapcsolatban nemcsak azt várjuk el, hogy a diszkretizált egyenletek elegendően nagy pontossággal közelítsék a folytonos feladatot, de azt is, hogy (U) megoldásuk konvergáljon a(z) (u) folytonos megoldáshoz. Konvergencia esetén mellett a diszkrét feladat megoldása tart a folytonos feladat megoldásához bármely j rácspontban és > 0 időpontban, azaz

(II.38.)

A konvergenciának tehát az előrejelzési időtáv minden időpontjában fenn kell állnia (azaz például a hatórás előrejelzések esetében csakúgy, mint a kétnapos prognózisoknál).

A konvergencia matematikai feltételének teljesülését nehéz belátni, ezért helyette legtöbbször a stabilitás teljesülését vizsgáljuk. Stabil feladatról akkor beszélünk, ha a feladat megoldása „folytonosan függ” a kiindulási feltételektől, azaz kis eltérés (hiba) a kezdeti feltételben nem vezet lényegesen eltérő megoldásra. Világos, hogy a stabilitás külön értelmezhető a folytonos és a diszkrét feladatban – előbbinél fizikai stabilitásról, utóbbinál számítási stabilitásról beszélünk. A meteorológiai problémák esetében, ahol az állapothatározók korlátos értékkészletű függvények, a stabilitást a hiba korlátosságán keresztül vizsgáljuk. Adott pontbeli és időlépésbeli Uj,n megoldás stabilitásához szükséges, hogy rögzített Δx rácsfelbontás mellett és esetén az

(II.39.)

alakban definiált hiba nem nő az idővel, azaz Az alkalmazott véges különbséges módszer pedig akkor stabil, ha bármely kezdeti feltételhez tartozó megoldás kielégíti a fenti feltételt (Mesinger és Arakawa, 1976). Lényeges, hogy ez csak a stabilitás szükséges feltétele, elégséges feltételt ugyanis nehéz megadni, s csupán néhány speciális esetben lehetséges.

A stabilitás és a konvergencia között a Lax-Richtmyer tétel (1956) teremt kapcsolatot, amely kimondja, hogy egy konzisztens véges differencia sémákkal megadott feladat akkor és csak akkor konvergens, ha stabil. Azaz a konzisztencia és a konvergencia együttes fennállása esetén a séma stabil, illetve a konzisztens és stabil séma egyben konvergens is. A tétel gyakorlati szempontból bír jelentőséggel, mert külön-külön a konzisztencia és a stabilitás fenti feltételének vizsgálata egyszerűbb, mint közvetlenül a konvergencia teljesülését ellenőrizni.

A stabilitás vizsgálata

A stabilitás feltételét egy gyakorlati példán keresztül szemléltetjük. Tekintsük a (II.27.) egydimenziós lineáris advekciós egyenletet. Ennek analitikus megoldása φ (– ct, 0), amely a kezdeti feltételben megadott hullám c sebességgel való haladását írja le (II.12. ábra).

II.12. ábra. Az egydimenziós lineáris advekciós egyenlet megoldásának sematikus rajza: a kezdeti feltétellel meghatározott hullám c > 0 sebességgel x-irányba való áthelyeződése egy időlépcső alatt.

Diszkretizáljuk a feladatot úgy, hogy az időbeli deriválásra explicit Euler-, a térbeli differenciálásra pedig bal oldali sémát alkalmazunk Δx rácsfelbontás és Δt időlépés mellett:

(II.40.)

Vezessük be az Courant-számnak nevezett mennyiséget, s rendezzük át az egyenletet az időlépcsők szerint:

(II.41.)

A fenti összefüggés alapján a séma a φ  következő időlépcsőbeli értékét adott rácspontban a φ előző lépésbeli értékeiből állítja elő, az aktuális (φ j,n) és a tőle balra lévő (φ j-1,n) rácspont felhasználásával.

Vizsgáljuk meg a Courant-szám értékét! Ha r 0 és 1 közé esik (), akkor a megoldás interpoláció eredményeként áll elő (II.13. ábra) és a séma pontosan úgy működik, ahogyan az advekció folyamata. Ellenben ha < 0 vagy > 1, akkor a séma extrapoláció segítségével állítja elő a megoldást φ j,n és φ j-1,n értékéből (tehát az említett két pont értékével ír le egy a tartomány másik részéről érkező információt; II.13. ábra). Az extrapolációval az a probléma, hogy az így nyert megoldás nem korlátos, ugyanis a háromszög-egyenlőtlenség alapján a (II.41.) egyenletre teljesül az alábbi:

(II.42.)

Abban az esetben, ha a φ j,n együtthatója 1-nél nagyobb, a megoldás minden időlépésben növekszik és nem marad stabil – a diszkrét feladat néhány időlépés után „felrobban”. Az együttható akkor nem nagyobb 1-nél, ha teljesül. Ha a koordináta-rendszert úgy választjuk, hogy benne a c sebesség pozitív irányú, akkor ez egyenértékű a

(II.43.)

feltétellel.

A feltétel általánosítható más feladatokra is, s általános alakjában Courant–Friedrichs–Lewy vagy CFL-kritériumnak (Courant et al., 1928) nevezik. Kimondja, hogy a diszkretizált feladat stabilitásához az időlépcsőt nem választhatjuk korlátlan hosszúságúra, ennek határt szab az alkalmazott térbeli felbontás és a feladat által leírt mozgásformák leggyorsabb terjedési sebessége. (Ismét kiemeljük, hogy a stabilitásnak ez szükséges, de nem elégséges feltétele.) Richardson 1920-as években végrehajtott kísérlete, az első kézzel végzett „numerikus” előrejelzés többek között azért végződött kudarccal, mert ez a kritérium akkor még nem volt ismert, s a számítások során Richardson túl hosszú időlépcsőt alkalmazott. A CFL-feltételnek a számítási hatékonyság szempontjából van jelentősége, a meteorológiai előrejelzések készítésénél a számítási műveletek és az adatok rendkívüli mennyisége miatt ugyanis az alkalmazott numerikus módszerek hatékonysága is lényeges szempont. Minél nagyobb időlépcsőt tudunk használni, az integrálási műveleteket annál kevesebb lépésben kell megismételni, s az előrejelzést annál gyorsabban tudjuk előállítani. Célunk tehát az, hogy az előrejelzési feladatot olyan numerikus sémák segítségével oldjuk meg, melyek adott felbontás mellett a lehető legnagyobb időlépcső használatát engedik meg.

II.13. ábra. Az egydimenziós lineáris advekciós egyenletben szereplő c advekciós sebesség, valamint a diszkretizációhoz használt Δx rácsfelbontás és Δt időlépés kapcsolatának sematikus rajza: az r Courant-szám értékétől függően az (n + 1)-edik időszinten a megoldás interpolációval (első panel) vagy extrapolációval (második és harmadik panel) áll elő (Kalnay, 2003 nyomán).

A numerikus stabilitást az előzőekben az ún. direkt vagy maximum-kritérium módszer segítségével vizsgáltuk, ami viszont szemléletessége ellenére csak korlátozottan használható. Helyette kiterjedten alkalmazzák a Neumann-módszert lineáris (vagy linearizált) problémák esetében. A módszert az advekciós egyenlet példáján keresztül mutatjuk be. Tekintsük tehát a (II.27.) feladatot a hozzá tartozó kezdeti feltétellel. A lineáris differenciálegyenletek megoldása kifejezhető függvénysorok segítségével, a Neumann-módszer alkalmazásánál a kezdeti feltételt és a megoldást Fourier-sor alakban keressük. Kihasználjuk még, hogy lineáris feladat esetében a Fourier-sor minden tagja megoldás, így elegendő egyetlen (jelen esetben a k-adik) Fourier-komponenst tekinteni. Ugyanezt a módszert követjük a folytonos feladathoz konstruált véges differencia feladat esetében is. Ez alapján a diszkrét feladat kezdeti feltétele tetszőleges j rácspontban, valamint megoldása az n-edik időlépésben az alábbi módon írható fel:

(II.44.)

ahol k és ck az adott hullám-módushoz tartozó hullámszám illetve együttható, λk pedig a k hullámszámhoz tartozó és a numerikus sémától függő komplexértékű, amplitúdó-jellegű mennyiség. A fenti összefüggésben tehát a numerikus megoldás úgy áll elő, hogy a kezdeti feltétel minden időlépésben egy amplitúdó-jellegű mennyiséggel szorzódik. Ahhoz, hogy a megoldás korlátos maradjon, szükséges, hogy ne legyen nagyobb 1-nél – ez a stabilitás (és a konvergencia) szükséges feltétele. Amennyiben , úgy a kezdeti feltételt leíró hullám minden időlépésben gerjesztődik és a diszkrét feladat nem lesz stabil. Ha , akkor a véges differencia séma fiktív csillapítást vezet be, azaz a kezdeti feltételhez tartozó hullám-megoldás az idővel folyamatosan csillapodik, szélsőséges esetben bizonyos hullámok teljesen el is tűnhetnek (ez bizonyos nagyfrekvenciás, nemkívánatos hullám-megoldásoknál előnyös lehet). Ideálisan pontosan 1-gyel egyenlő, ekkor a véges differencia séma nem változtatja meg a kezdeti feltételt leíró hullám amplitúdóját (neutrális).

Megjegyezzük, hogy a Neumann-módszernél is általánosabban alkalmazható az energia-módszer, például nem-lineáris problémák numerikus stabilitásának vizsgálatára. A módszert azonban a jegyzet keretében nem tárgyaljuk, az érdeklődőknek Durran (1999) írását ajánljuk.

A stabilitási feltétel különböző meteorológiai problémákra

Térjünk most vissza a CFL-kritérium általános alakjához, amihez tekintsük a lineáris gravitációs hullám (II.28.) egyenletét. A feladat átírható a következő alakban:

(II.45.)

ahol a jobb oldalon szereplő mátrix sajátértéke (a determináns alapján) , azaz a k hullámszámhoz két, ellentétes irányban fázissebességgel haladó hullám tartozik. A feladatban tehát a legnagyobb terjedési sebesség a gravitációs hullám alakú terjedési sebessége. A Neumann-módszer segítségével belátható (a Neumann-módszer alkalmazására az alfejezet végén a 2. feladat szolgáltat példát), hogy a feladat diszkretizációjára a fenti sémát alkalmazva a stabilitási kritérium:

(II.46.)

Becsüljük meg a napjaink numerikus modelljeiben alkalmazható időlépcsők nagyságát! A tropopauza átlagos magasságát (H-t) 10 km-nek, a gravitációs gyorsulás (g) átlagos értékét 10 ms-2-nak véve, ez a sebesség 300 ms-1 nagyságúnak adódik, s így a 10 km-es rácsfelbontásnál alkalmazható időlépcső nem nagyon haladhatja meg a 30 másodpercet. Ezzel szemben az advekciós feladatban szereplő advekciós sebesség még a nagy magasságokban sem nagyobb 100 ms-1-nál, ekkor az integrálási időlépés 10 km-es felbontás mellett 100 s is lehet. Látható tehát, hogy a most kapott (II.46.) stabilitási kritérium szigorúbb korlátot jelent az időlépcső hosszára, mint amit az advekciós egyenlet esetében kaptunk.

Megvizsgálva az előrejelzési feladatot reálisabban közelítő, advekciós és gravitációs hullám tagokat egyaránt tartalmazó (II.29.) feladatot is, hasonló módon belátható, hogy a fenti explicit sémát alkalmazva az időlépcsőnek az alábbi kritériumot kell kielégíteni:

(II.47.)

Tehát az adott felbontás esetén alkalmazható időlépcső hosszára az advekció és a gravitációs hullám terjedési sebességének összege szab korlátot. Mivel a korábbiakban láttuk, hogy az utóbbi lényegesen nagyobb az előbbinél, ezért kijelenthetjük, hogy a stabilitás szükséges feltételénél gyakorlatilag a gravitációs hullámok terjedési sebessége a meghatározó.

Hasonló eredményre vezet, ha az időbeli és térbeli deriváltak közelítésére leapfrog illetve centrált véges differencia módszereket alkalmazunk. A leapfrog-séma esetében egy további jelenséggel kell számolnunk, nevezetesen azzal, hogy a Neumann-módszerrel λk-ra két megoldás adódik. Az egyik a folytonos feladat fizikailag értelmes megoldásához tartozó ún. fizikai módusz, a másik viszont abból ered, hogy két kezdeti feltételt kell megadni a séma indításakor, s így jelenik meg a fizikai értelemmel nem bíró, ún. számítási módusz. Mivel a teljes numerikus megoldás a fizikai és a számítási módusz (lineáris) kombinációjaként áll elő, nagy jelentősége van annak, hogy a számítási módusz hogyan fejlődik az idővel. A fiktív móduszhoz általában negatív előjel társul, aminek következtében minden időlépcsőben előjelet vált, ezért célszerű valamilyen módszerrel kiszűrni. A szűrő megoldások a leapfrog-sémában egyébként gyengén csatolt páros és páratlan időlépcsők összekapcsolását szolgálják, a gyakorlatban a legelterjedtebb az Asselin-filter (1972).

Megállapíthatjuk tehát, hogy a tisztán explicit sémák segítségével az előrejelzési feladatot csak olyan időlépcső alkalmazásával tudjuk a stabilitás megőrzésével megoldani, amelyek kielégítik a fenti feltételt. Ez azonban, ahogy a nagyságrendi becslés során láttuk, nem kedvez a számítógépes megvalósítás hatékonyságának, mert a teljes primitív egyenletrendszer esetén csak korlátozott hosszúságú időlépést enged meg. Célunk tehát olyan nagypontosságú diszkretizációs módszerek alkalmazása, melyekkel a numerikus stabilitás gyakorlatilag tetszőlegesen hosszú integrálási időlépcső használata mellett sem sérül. Hogy ehhez milyen irányba kell elindulni, arra a fázishiba jelensége fog magyarázatot adni.

A hullám-megoldás amplitúdóján kívül a véges differencia séma a folytonos feladatban jellemző c fázissebességet is módosíthatja, ezt hívjuk fázishibának, amit a diszkrét feladatbeli c’ és a folytonos feladatbeli c fázissebességek arányaként definiálunk. Amennyiben c’ > c, a véges differencia séma gyorsítja a folytonos feladatbeli k hullámot; ha c’ < c, akkor a séma lassítja a k hullámot; s ha c’ = c, akkor a séma nem változtatja meg a hullám fázissebességét. Bár az utóbbi eset tűnik ideálisnak, mégis alkalmasan megválasztott véges differencia séma esetében bizonyos hullámok „lassúbbá torzítása” előnyös is lehet. Ennek illusztrálására tekintsük például a (II.28.) feladatot, amelynek a meteorológiai szempontból kevéssé releváns gyorsan terjedő gravitációs hullámok is megoldásai. Az előzőekben láttuk, hogyha erre explicit véges differencia sémát alkalmazunk, amely a hullámok fázissebességét lényegesen nem változtatja meg, akkor meglehetősen szigorú feltételt kell a numerikus stabilitás teljesítéséhez kielégíteni, amiben a vezető tag a feladat által leírt leggyorsabban terjedő mozgásforma (jelen esetben a gravitációs hullámok) fázissebessége. Implicit séma alkalmazásával azonban kedvezőbb feltétel nyerhető vagy éppen feltétel nélküli stabilitással tudunk dolgozni, mégpedig azáltal, hogy az implicit módszer ezeknek a nagyfrekvenciás hullámoknak a fázissebességét csökkenti. A megoldás fizikai értelmezésénél ez nem okoz problémát, mert ahogyan már említettük, ezek a hullámok a meteorológiai folyamatok szempontjából nem lényegesek, viszont a gyorsan terjedő hullámok lassításával stabilizálható a feladat ezekért felelős része.

Hatékony numerikus sémák

Az előző alfejezetben különböző véges differencia sémák numerikus stabilitását tekintettük át. Láttuk, hogy az explicit sémák alkalmazásával a stabilitás csak akkor teljesülhet, ha eleget teszünk a CFL-kritériumnak. Az implicit sémák esetében tetszőlegesen hosszú időlépcső választható, mégpedig azért, mert ezek a feladatban érvényes fázissebességet lassítják, s ezáltal a gyorsan terjedő hullám-megoldásokra is stabil megoldást biztosítanak. Az általuk okozott fázishiba ugyanakkor nem minden hullám esetében kívánatos, ezért, valamint bonyolult megvalósításuk miatt a meteorológiai gyakorlatban nem alkalmaznak tisztán implicit sémákat.

Szemi-implicit módszer

Alkalmazzunk most a (II.29.) feladatra ún. szemi-implicit diszkretizációt, ami az advekciós tagok esetében megtartja az explicit sémát, a gravitációs hullám tagokat viszont implicit módon kezeli (Robert, 1981). Például az előbbi tagok esetében az explicit középponti sémát használjuk, az utóbbi tagok esetében pedig centrált sémák átlagát vesszük két időszintre:

(II.48.)

A Neumann-módszer segítségével belátható, hogy a szemi-implicit módszer alkalmazása esetén a következő adódik:

(II.49.)

ez pedig minden esetben teljesül, amikor az advekciós egyenletnél kapott stabilitási feltétel érvényes, azaz amikor:

(II.50.)

Tehát a szemi-implicit módszer alkalmazásával elértük, hogy ne a feladat által leírt leggyorsabban terjedő mozgásforma, a gravitációs hullámok terjedési sebessége legyen a meghatározó az időlépcső megválasztásánál, hanem az annál jóval kisebb advekciós sebesség, lehetővé téve a nagyobb időlépésekben történő stabil modellintegrálást.

A szemi-implicit séma szépsége abban áll, hogy az implicit kezelést az egyenlet lineáris részére alkalmazza, s ezek a tagok felelősek egyben azokért a hullám-megoldásokért, melyek gyorsan terjednek, de jelentőségük meteorológiai szempontból kicsi. Az implicit módszer ezeket a hullámokat lelassítja, ezáltal stabilizálja a feladat lineáris részét, de érintetlenül hagyja a meteorológiailag releváns, lassúbb mozgásformákat. Általánosan a szemi-implicit módszer a következő alakban írható fel egy vektor-mennyiségre:

(II.51.)

ahol a teljes nem-lineáris modellt linearizáljuk egy referencia-állapot körül és az L operátor a modell linearizált része, az N operátor a nem-lineáris maradéktag, a 0, +, – indexek pedig az aktuális, a későbbi és a korábbi időlépcsőket reprezentálják, jelölve, hogy a szemi-implicit módszerben a linearizált tagokra implicit, a nem-lineáris tagokra explicit kezelést alkalmazunk. Hátránya az eljárásnak, hogy a linearizáláshoz használt referencia-állapot (izoterm, nyugvó légkör) gyakran távol esik a valós légköri állapottól. A nem-lineáris egyenletrendszer teljes implicit kezelése (ami biztosíthatná a feltétel nélküli stabilitást) nem lehetséges (és reális), mert ez egy nem-lineáris operátor invertálását igényelné. Ezért hogy az advekciós sebesség által meghatározott stabilitási feltételt tovább tudjuk enyhíteni, a szemi-implicit sémát ötvözni kell egy másik hatékony módszerrel.

Szemi-Lagrange módszer

A hidro-termodinamikai egyenletrendszerben a legdominánsabb nem-lineáris tag a nem-lineáris advekció, s ennek kezelésére hívjuk segítségül a Lagrange-módszert. A Lagrange-szemléletben nem egy térben rögzített koordináta-rendszer pontjaiban tekintjük az állapothatározók változását, hanem a részecskékhez (légelemekhez) rögzített lokális koordináta-rendszerekkel dolgozunk: a részecskék trajektóriáját követjük, ami mentén magukkal viszik kiindulási tulajdonságaikat. Adott φ tulajdonság trajektória-menti megmaradását fejezi ki az egydimenziós lineáris advekciós egyenlet alábbi Lagrange-alakja:

(II.52.)

Tehát az állapothatározók adott időpontban vett térbeli eloszlásának és pályájának ismeretében meghatározható jövőbeli eloszlásuk. A numerikus számítások során azonban szeretnénk az Euler-szemléletnek azt a kényelmes tulajdonságát megtartani, hogy a légköri változókat egy szabályos rács rácspontjaiban tekintjük (egyrészt mert egyéb műveleteket is rácson kell elvégezni, másrészt mert ezáltal egyenletes térbeli lefedettség biztosítható). Ez a tiszta Lagrange-módszer segítségével nem lehetséges, hiszen azok a részecskék, melyek kiinduláskor még szabályosan helyezkedtek el, egy időlépés után már szabálytalan és térben inhomogén elrendeződést ölthetnek. Ezért az ún. szemi-Lagrange módszerben minden időlépésben egy backward (visszafelé) trajektória számításával, majd térbeli interpoláció alkalmazásával állítjuk elő az állapothatározók értékét az általunk kívánt rácspontokban (II.14. ábra). Világos, hogy így az időbeli integrálás során nem ugyanazokat a részecskét követjük, hanem a rácspontok elhelyezkedése alapján minden időlépésben egy új részecskehalmazt definiálunk.

II.14. ábra. A szemi-Lagrange módszer sematikus rajza.

A trajektória kiindulási pontját (az ábrán * jelöli) az advekciós sebesség ismeretében tudjuk meghatározni, amit az egyszerűség kedvéért most tekintsünk állandónak: a részecske ekkor Δt idő alatt c sebességgel advektálódik, minek során távolságot tesz meg. A korábbi időszintről való indulás nem feltétlenül rácspontból történik, ezért a megtett utat két részre oszthatjuk: a Δt idő alatt bejárt rácspontok számát p-vel jelöljük, ahol p egész szám, a maradék út pedig a rácstávolság α-szorosa, ahol α értelemszerűen 0 és 1 közé eső törtszám, azaz:

(II.53.)

Az (+ 1)-edik időlépésben a φ állapothatározó értéke a konzervativitás miatt meg fog egyezni a φ  n-edik időlépésben és * pontban felvett értékével. Mivel az n-edik időszinten a φ eloszlását csak a rácspontokban ismerjük, ezért a *-gal jelölt pontban az értéket térbeli interpolációval tudjuk előállítani. Az egyszerűség kedvéért tekintsük most a legegyszerűbb lineáris interpolációt:

(II.54.)

A Neumann-módszer segítségével belátható, hogy

(II.55.)

ahol a korábbiak szerint k az adott hullám-módushoz tartozó hullámszám, λk a k hullám-számhoz tartozó és a numerikus sémától függő komplexértékű, amplitúdó-jellegű mennyiség. Ez alapján a séma stabilitásának kritériuma az, hogy α-ra fennálljon az feltétel. Ez azonban a séma konstrukciójából eredően teljesül, így a séma feltétel nélkül stabil. A szemi-Lagrange módszer alkalmazásakor tehát nem kell tekintettel lennünk a CFL-kritériumra, a kritérium által megengedettnél hosszabb időlépcsőt is választhatunk. Ez az időlépcső azonban továbbra sem lehet tetszőlegesen nagy, figyelembe kell venni a Lipschitz-feltételt, mely szerint a trajektóriák egy időlépcső alatt nem metszhetik egymást – ellenkező esetben nem lehetséges a részecskék pályáját egyértelműen meghatározni (Smolarkiewicz és Pudykiewicz, 1992). Mindazonáltal numerikus kísérletekkel igazolták, hogy mezoskálájú modellek esetén a szemi-Lagrange séma az időlépcső mintegy hatszoros növelését teszi lehetővé az euleri sémákkal szemben, továbbá a szemi-implicit sémával való kombinálása esetén az integrálás hatékonysága további hatszorosával növelhető (Staniforth és Côté, 1991). Hatékonyságuk miatt a szemi-Lagrange és a szemi-implicit módszereket, illetve kombinációjukat szívesen alkalmazzák az operatív időjárás-előrejelző modellekben, pl. az Országos Meteorológiai Szolgálatnál futtatott ALADIN és AROME modellekben is. Az ALADIN modellben a szemi-Lagrange séma 8 km-es horizontális felbontás mellett 300 másodperces időlépcső alkalmazását teszi lehetővé, ami gyors modellintegrálást eredményez, különösen ha összevetjük azzal, hogy a WRF modell euleri advekciókezelése hasonló, 7 km-es felbontáson csak 40 s hosszú időlépcsőt enged meg (André, 2012). Meg kell említeni, hogy vannak olyan kis időskálán zajló fizikai folyamatok, amelyek pontos leírásához nem elegendő egy hosszabb időlépcső. Például a csapadékképződés során a csapadékelemek néhány perc alatt kieshetnek az egységnyi rácstérfogatból, ezért ennek pontosabb leírásához a dinamikai időlépcsőt felosztják kisebb lépésekre.

A fentiekben a szemi-Lagrange séma alapelvét mutattuk be, a gyakorlati megvalósítás során általában a módszer bonyolultabb, két-időszintes változatát alkalmazzuk. A szemi-Lagrange sémával előállított megoldást összehasonlítva a folytonos megoldással, megállapítható, hogy a módszer olyan csillapítást vezet be, ami bizonyos esetekben a legrövidebb hullámok teljes kioltáshoz is vezethet (Dévényi et al., 1998). Ez a fiktív csillapítás a hosszabb szimulációknál (pl. éghajlati kísérlet esetében) tönkreteheti a megoldást, ezért elsősorban klímamodellezési célokra kifejlesztették a séma ún. nem-interpolációs változatát (Ritchie, 1986). Ekkor a trajektóriát lényegében egy a legközelebbi rácspontba mutató és egy „maradék” vektor összegére bontják fel, és csak a rácspontok közötti advekciót kezelik szemi-Lagrange módszerrel. Ezen túlmenően a szemi-Lagrange séma nem változtat a megoldás terjedési sebességén, ha a kiindulási pont éppen rácspontba esik, ha a részecske egy időlépés alatt nagy távolságot tesz meg, illetve a hosszúhullámok esetében. A fázishiba és a fiktív csillapítás mértéke a számítási költség növekedése árán javítható, ha lineáris interpoláció helyett magasabbrendű, pl. köbös spline interpolációt alkalmazunk (Dévényi et al., 1998).

Az elmúlt években a modellek felbontásának finomabbá válásával a fizikai parametrizációk számításigénye némileg megnövekedett: míg a 90-es években hozzávetőlegesen a modellintegrálás 30 %-át tették ki a parametrizációs számítások (Gustafsson és McDonald, 1996), addig ma ez már meghaladja a 40 %-ot (Miller et al., 2010). Ugyanakkor a modellintegrálásnak még mindig jelentős része fordítódik a dinamikai számítások végrehajtására, ezért továbbra is fontos a hatékony numerikus sémák használata és fejlesztése.

II.3.2. Spektrális módszerek

Horányi András

Szépszó Gabriella

A spektrális módszerek alapjai

A spektrális módszerek alkalmazása során a HTER-ben szereplő prognosztikai változókat (amelyik tipikusan a felszíni légnyomás, hőmérséklet, nedvesség és a szélsebesség két komponense) egy függvényrendszer elemei segítségével írjuk fel. A módszer leírásához tekintsük a következő egyenletet:

(II.56.)

ahol L egy differenciál-operátor, u a függő változó (a fentebb említett prognosztikai változók bármelyike) és f(x) adott kényszer függvény (tipikusan a HTER jobb oldalán szereplő kifejezések). Az u(x) ismeretlen változót a φ j(x) lineárisan független ortogonális bázis függvények lineáris kombinációjával közelítjük:

(II.57.)

ahol uj a j-edik bázis függvényhez tartozó együttható. A gyakorlatban természetesen csak véges összegzést tudunk elvégezni:

(II.58.)

A feladat megoldásához a konkrét feladathoz jól illeszkedő függvényrendszert célszerű választani. A spektrális módszer esetében, mint már említettük, a bázis függvények ortogonális rendszert alkotnak (ilyen például a trigonometrikus függvények rendszere). Az ilyen függvényeket szokás globális bázis függvényeknek is nevezni, megkülönböztetésül a Galjorkin-módszerek másik típusánál, a véges elem módszernél alkalmazott, kis területen értelmezett alacsony fokszámú polinomokból álló lokális (Kalnay, 2003) bázis függvényektől (például az ún. kalap függvények; Strang és Fix, 1973). A spektrális módszer végrehajtása során célunk az uj együtthatók meghatározása. Ezen együtthatók egyértelmű meghatározásához egy további feltételt is ki kell szabnunk (ezt nevezzük Galjorkin-feltételnek; Durran, 1999), amely szerint a közelítés hibájának ortogonálisnak kell lennie minden egyes bázis függvényre. Az ortogonalitás feltételét a következőképpen írhatjuk fel:

(II.59.)

ahol  vagy

(II.60.)

A hibát a (II.59.) ortogonalitási feltételben kifejtve az alábbi egyenletrendszerhez jutunk:

(II.61.)

Ez az összefüggés N darab ismeretlent (az uj együtthatókat) tartalmazó, N tagból álló algebrai egyenletrendszer. A spektrális módszer végrehajtása során ezt az egyenletrendszert oldjuk meg az ismeretlen uj együtthatókra. A spektrális módszer alkalmazása során elvben minden műveletet az így származtatott együtthatók segítségével az együtthatók terében, azaz az ún. spektrális térben végzünk el. A későbbiekben láthatjuk, hogy a spektrális módszer gyakorlati alkalmazása esetén ez nem feltétlenül van így.

A spektrális módszerben a bázis függvények optimális megválasztását az előrejelzési feladat geometriája, valamint a határfeltételek határozzák meg. Téglalap alakú tartomány és periodikus határfeltételek esetében kézenfekvő a komplex Fourier bázis függvények alkalmazása (ilyet alkalmaznak pl. a globális modellekben a szélességi körök mentén); nem-periodikus határfeltétellel bíró tartományok esetében használhatók a Csebisev-polinomok, míg a Legendre-függvények szférikus rendszerekben terjedtek el (pl. a globális modellekben a hosszúsági körök irányában alkalmazzák őket).

Egy egyszerű példa

A spektrális módszer gyakorlati részleteinek megvilágításához tekintsük az alábbi egyszerű példát:

(II.62.)

Ekkor .

Periodikus határfeltételeket adunk meg:

(II.63.)

Fourier bázis függvényeket alkalmazunk, a következő alakban:

(II.64.)

Ezek a bázis függvények a intervallumon ortogonálisak és kielégítik a határfeltételeket is. A (II.62.) egyenlet bal oldalának közelítése alkalmazva a bázis függvények második deriváltját:

(II.65.)

A közelítés hibájának és a bázis függvények ortogonalitásának feltételét ennek alapján a következő egyenlet fejezi ki:

(II.66.)

A (II.66.) egyenletben megjelennek a bázis függvények szorzatai, amelyek a bázis függvények konkrét alakját figyelembe véve az alábbi módon fejthetők ki:

(II.67.)

ahol a Kronecker delta (= 1, ha j és = 0, ha ≠ j). Az ortogonalitási feltétel alapján kapott egyenlet megoldása így az

(II.68.)

alakot veszi fel a keresett együtthatókra. Az együtthatók tehát a kényszer függvény véges Fourier-transzformációjával arányosak.

Transzformációs módszer

A spektrális módszerek egyik legfontosabb részlete, hogy nem minden művelet oldható meg hatékonyan az „ együtthatók terében” (azaz a spektrális térben), hanem a számolások egy részét a rácspontok terében kell végrehajtanunk (azaz abban a fizikai térben, ahol a modell változóit rácson reprezentáljuk). Ez a gyakorlatban annyit jelent, hogy a spektrálisnak nevezett meteorológiai modellek tulajdonképpen csak pszeudo-spektrálisak, hisz nem minden számolást valósítunk meg a spektrális együtthatók révén. Az e mögött rejlő gyakorlati részleteket az alábbiakban mutatjuk be röviden.

Tekintsük a barotróp örvényességi egyenletet, és alkalmazzuk a spektrális módszert a megoldására. A barotróp örvényességi egyenletet a következő alakban írhatjuk fel (Charney et al., 1950):

(II.69.)

ahol az áramfüggvény. Tételezzük fel, hogy a mezők a sík mindkét (x és y) irányában periodikusak, és teljesül, hogy

(II.70.)

ahol k és l hullámszámok. Ortogonális bázis függvényekként a mindkét irányban periodikus harmonikus függvényeket alkalmazzuk az alábbi formában:

(II.71.)

Az áramfüggvény a bázis függvények segítségével a következő módon közelíthető:

(II.72.)

ahol az összegzés természetesen véges. Megjegyezzük, hogy ahhoz, hogy a mennyiség valós legyen, az együtthatókra teljesülnie kell a feltételnek, ahol ( )* a komplex konjugálásra utal.

Vezessük be az mknlj hullámszám vektort, valamint az xyj hely vektort. A közelítése ekkor

(II.73.)

A barotróp örvényességi egyenletben megjelenő tagok a bázis függvények segítségével az alábbi alakban írhatóak fel:

(II.74.)

A H és L hullámszám vektorok azért kerültek bevezetésre, mert a későbbiekben az egyes mennyiségeket összeszorozzuk és átrendezzük, és olyankor mindig teljes összegeket kell összeszorozni. Ezek után a fenti kifejezéseket behelyettesítjük a barotróp örvényességi egyenletbe és felírjuk az egyenlet teljesülésének hibáját (eN):

(II.75.)

A spektrális módszer ortogonalitási feltétele a hiba bázis függvényekre való ortogonalitását írja fel. Mivel az ortogonalitás feltétele a skalárszorzat eltűnése, és komplex függvények esetén a függvények skaláris szorzata definíció szerint , az eN hibafüggvényt a bázisvektorokkal beszorozva nyerjük az ortogonalitási feltételre vonatkozó egyenleteket. Szorozzuk meg tehát az előző egyenletet -rel, majd integráljuk azt a periodikus tartományon és vegyük figyelembe az hiba minden egyes bázis függvényre vett ortogonalitását:

(II.76.)

minden egyes M-re az eredeti összegzésnek megfelelően. A periodikus határfeltétel miatt az exponenciális függvények minden határozott integrálja eltűnik, kivéve azt az esetet, amikor a kitevő zérussal egyenlő. Ennek alapján minden M-re a következő közönséges differenciálegyenlethez jutunk:

(II.77.)

A számolásnál az első két tagban az M esetben, míg az utolsó tagban az – H esetben léteztek nem nulla tagok.

Vegyük észre, hogy az egyenlet utolsó tagja az eredeti egyenletben megjelenő nem-lineáris advekciós tag miatt fellépő különböző hullámok közti kölcsönhatásokat írja le (kölcsönhatási tagnak is szokták nevezni). Ez annyit jelent, hogy az M. hullámra a H. és (– H). hullámok fejtenek ki hatást. Az is látható, hogy ez a kölcsönhatásokat leíró nem-lineáris tag igen bonyolult, és meghatározása a spektrális térben nehézségeket okoz. Ennek elkerülése érdekében vezették be a transzformációs módszert, ami annyit jelent, hogy az egyenletekben fellépő (ilyen és ehhez hasonló) nem-lineáris tagokat nem a spektrális együtthatók terében, a spektrális együtthatók révén határozzuk meg, hanem a rácsponti térben, pontosan ugyanúgy, ahogy azt a rácsponti modelleknél tesszük. Ily módon ötvözhetők a rácsponti és spektrális módszerek előnyös tulajdonságai: a lineáris műveletek elvégzése (tipikusan a deriváltak számolása) egyszerű és pontos a spektrális együtthatók segítségével (jól megválasztott bázis függvények segítségével), míg a nem-lineáris műveletekelvégzése (például az advekció származtatása) könnyebb rácsponti térben (hisz ott csak egyszerűen összeszorozzuk a nem-lineáris kifejezésben szereplő változókat). Mindazonáltal a transzformációs módszer csak abban az esetben alkalmazható hatékonyan, ha létezik olyan számítástechnikailag is hatékony transzformációs eljárás, ami gyorsan végrehajtható, mivel a spektrális modellek alkalmazása során minden integrálási időlépésben ingázunk a spektrális és a rácsponti terek között (a periodikus függvények esetében a Fourier-transzformáció megvalósítása az ún. gyors Fourier-transzformáció segítségével például ilyen algoritmus).

A transzformációs módszer gyakorlati szükségessége miatt a spektrális modellek végrehajtásának legfontosabb lépései az alábbiak:

  • Számítások a spektrális együtthatók terében (ez a kiindulás, mert a spektrális modellek esetében a spektrális együtthatók száma kisebb, mint a kapcsolódó rácspontok száma, s ezért érdemes az egyes meteorológiai változókat spektrális együtthatókként tárolni; ne felejtsük el, hogy a spektrális tér és a rácsponti tér között egy-egy értelmű megfeleltetés van): differenciál-operátorok számítása, valamint az egyenletekben szereplő lineáris műveletek (pl. a szemi-implicit séma) végrehajtása.

  • Inverz transzformáció a spektrális térből a rácsponti térbe, azaz a rendelkezésre álló spektrális együtthatók „átalakítása” rácspont térbeli fizikai változókká (példa: inverz gyors Fourier-transzformáció).

  • Az egyenlet nem-lineáris tagjainak származtatása a rácsponti térben (az egyes rácsponti értékek egyszerű összeszorzásával), és olyan további műveletek elvégzése, amelyek végrehajtása a spektrális együtthatók terében nehezebb, mint a rácsponti térben (pl. a fizikai parametrizációs számítások elvégzése), és végül az új időlépcsőre való áttérés (azaz a modell egy lépéssel való előreléptetése). Mivel az egyenletek legfontosabb nem-lineáris tagjaként az advekció számítása is itt történik, ezért amennyiben szemi-Lagrange sémát alkalmazunk, úgy arra is a rácsponti térben kerül sor.

  • Direkt transzformáció, azaz a rácsponti értékek spektrális reprezentációjának előállítása (példa: direkt gyors Fourier-transzformáció) és a fenti műveletek megismétlése most már a legfrissebb időlépcsőben.

A csonkítás fogalma

Már a bevezetésben említettük, hogy az ismeretlen változók bázis függvényekkel való felírásakor a bázis függvények lineáris kombinációja egy véges összegzést jelent (nem is jelenthet mást, hiszen nem vehetjük figyelembe az összes lehetséges spektrális együtthatót, azaz valahol határt kell szabnunk az összegzés hosszának). Ennek kapcsán jegyezzük meg azt is, hogy a csonkítás mértéke mindig kompromisszum eredménye az optimális pontosság és a rendelkezésre álló számítási kapacitás között. A csonkítás legfontosabb kérdése, hogy milyen logika szerint válasszuk ki azokat a spektrális együtthatókat, amelyeket az összegzés során figyelembe veszünk (A csonkítás elnevezése általában arra utal, hogy a kiválasztott együtthatók milyen alakot vesznek fel a horizontális hullámszámok terében).

A globális modellek esetében használt szférikus harmonikusok segítségével az ismeretlen változó az alábbi alakban írható fel:

(II.78.)

ahol a λ a földrajzi hosszúság, a μ a földrajzi szélesség szinusza, az Ym,n függvény λ irányában a Fourier-módusz, μ irányában pedig a Legendre-függvény, N és M a figyelembe vett spektrális együtthatók száma az egyes koordináta irányokban. Az összegzésben szereplő N(m) értéket többféle módon választhatjuk meg attól függően, hogy milyen feladatot szeretnénk megoldani, illetve milyen folyamatokat szeretnénk leírni. A globális modellekben leggyakrabban az M választással élünk, ami az ún trianguláris vagy háromszög alakú csonkítást eredményezi (II.15. ábra). Ez a csonkítás egységes térbeli felbontást biztosít a Föld teljes felszínén és invariáns bármely a Föld középpontja körüli elforgatásra. Ugyanakkor nem túl előnyös olyan folyamatok leírásánál, amelyek a földrajzi helytől függő változékonyságot mutatnak (erre példa a geopotenciál mező, amely jóval nagyobb változékonysággal bír a magasabb szélességeken, mint a trópusokon; Durran, 1999). Ezért a globális modellekben használatos még az ún. romboidális csonkítás, amelynél N-et N(m) = |m| + M-nek választjuk (II.15. ábra). Ennél a csonkításnál a közepes szélességeken a legnagyobb a térbeli felbontás, ami inkább a kis (alacsony) felbontású modellekben előnyös.

II.15. ábra. Háromszög alakú (bal) és romboidális (jobb) csonkítások (Riddaway, 2001).

A korlátos tartományú modelleknél alkalmazott, mindkét irányban periodikus Fourier-függvények esetében az ismeretlen függvény a bázis függvények segítségével az alábbi módon írható fel:

.

(II.79.)

Ez az összefüggés felfogható különböző frekvenciájú hullámok szuperpozíciójaként, ahol a kis n és m értékek nagy hullámhosszakat (hosszú hullámokat), míg a nagy n és m értékek kis hullámhosszakat (rövid hullámokat) jelentenek. Az ismeretlen változót tehát a különböző hullámhosszú hullámok összegeként írjuk fel. Természetesen annál pontosabban járunk el, ha minél több hullámot tudunk figyelembe venni, azaz N és M értékei minél nagyobbak.

A kétdimenziós Fourier-függvények esetében a téglalap alakú és az elliptikus csonkítás terjedt el, ekkor a zonális és a meridionális hullámszámok terében az összegzéshez felhasznált spektrális együtthatók egy téglalapot illetve ellipszist írnak le (és ezeket a téglalap illetve az ellipszis egyenlete alapján választjuk ki, II.16. ábra).

II.16. ábra. Az elliptikus és téglalap alakú csonkítások.

A nem-lineáris instabilitás

Itt érdemes kitérni a nem-lineáris instabilitás, idegen szóval az „aliasing” fogalmára is. Adott Δx felbontású ráccsal 2Δx hosszúságú hullámot tudunk leírni. Ugyanakkor az előrejelzési feladat megoldása során – tipikusan a nem-lineáris tagok számításakor – megjelennek olyan hullámok is, amelyek 2Δx-nél rövidebbek, és energiájuk hozzáadódik a rendszerhez. Ezeket a hullámokat a Δx felbontású ráccsal nem tudjuk leírni, ezért energiájuk a rácson még „felbontott” hullámok energiájához adódik hozzá (II.17. ábra; emiatt hívják ezt a jelenséget aliasingnek is: mert ezek a rövid hullámok „úgy tesznek”, mintha az adott felbontáson még leírt, hosszabb hullámok lennének). Mivel az energiatöbblet a spektrum jobb oldali végéhez, azaz a rövid hullámhosszú hullámokhoz adódik, ahol az energia egyébként kicsi, ez akár az energia megduplázódását is jelentheti, így instabilitást okozva. Ugyanez a jelenség a spektrális modellekben is fennáll, és jól illusztrálható a fentebb alkalmazott kétdimenziós Fourier-függvények segítségével. A fenti alakot tekintve látható, hogy két ismeretlen változó összeszorzása esetén olyan tagok is megjelennek, amelyek kitevője már nagyobb lesz akár az N vagy akár az M értékénél (az exponenciális függvények szorzata a kitevők összeadását jelenti, s így jelennek meg a „túlságosan” nagy hullámszámok). Más szavakkal a nem-lineáris tagok olyan információt hoznak létre, amelyet az adott csonkítás már nem tud kezelni, s megjelenésük instabilitást okoz. A jelenséget Philips (1956) vette észre, amikor az általános cirkulációt szimuláló numerikus kísérlete a gondosan megválasztott, stabil időlépcső ellenére 16 nap után „felrobbant”. Ő ismerte fel, hogy az instabilitás a nem-lineáris műveleteknek köszönhető és nevezte el a jelenséget nem-lineáris instabilitásnak (Philips, 1959).

II.17. ábra. A különböző hullámszámhoz (x-tengely) tartozó hullámok energiaspektruma (y tengely). Látható, hogy adott csonkításnál (π/Δx) nagyobb (a rácsfelbontás által le nem írt) (k1 + k2) hullámszámmal rendelkező hullámok energiája a rövid hullámhosszú, hullámszámú hullámok energiáját növeli.

A nem-lineáris instabilitást célszerű tehát elkerülni, ami a spektrális modellekben túlcsonkítással egyszerűen megtehető. Tekintsük egy dimenzióban a [0, 2π] tartományt, amin Δx rácsfelbontást alkalmazunk. A tartomány ekkor = 2π/Δ+ 1 rácspontot tartalmaz, és a Δx felbontással még leírható legrövidebb hullám hullámszáma = (– 1)/2 = π/Δx. Azaz ahhoz, hogy a spektrális és a rácsponti tér közötti transzformáció során ne veszítsünk információt, a csonkítási hullámszám (M) és az alkalmazott rácspontok (J) száma nem lehet független egymástól, közöttük a következő kapcsolatnak kell fennállnia:

(II.80.)

A nem-lineáris szorzatok miatt megjelenő, a (kk2) > K hullámszámú hullámok energiája „szimmetrikusan” a  = π/Δx – (kk– π/Δx) = 2π/Δx – (kk2) hullámszámú hullámok energiájához adódik hozzá (

II.17. ábra). A csonkítás mértékét tehát úgy kell megválasztani, hogy ez az energiatöbblet ne okozzon változást azokon a skálákon, amelyeket az adott felbontású rács leír, tehát legyen nagyobb K-nál. Ha a csonkítási hullámszám M, akkor (kk2) maximális értéke 2M, s ekkor a feltétel a következőképpen néz ki:

(II.81.)

Ha ezt tovább alakítjuk kihasználva, hogy 2π/Δx = – 1, akkor az alábbi összefüggést kapjuk:

(II.82.)

Tehát a rácspontok száma és a spektrális térben leírt hullámok száma közötti kapcsolatnak a fenti feltételt kell kielégítenie; ezzel a rendszert „túlcsonkítjuk”, azaz kevesebb hullámszámot veszünk figyelembe, mint amit a rácspontok száma amúgy megengedne. Orszag (1971) megmutatta, hogy ha ezt a feltételt betartjuk, akkor ezzel kiküszöbölhető a kvadratikus tagok által okozott instabilitás. A háromszoros nem-lineáris szorzatok esetében továbbra is fellép az aliasing, ez azonban már nem okoz számottevő további instabilitást.

A nem-lineáris instabilitás kiküszöbölésére a spektrális technikán kívül egyéb módszerek is használhatók. Például az advekciót mint legjellemzőbb nem-lineáris tagot szemi-Lagrange módszerrel kezelve szintén elkerülhetjük a problémát. Ezért a szemi-Lagrange módszer spektrál technikával való kombinálása esetén elegendő a feltételt kielégíteni.

Összefoglalás

A spektrális módszer előszeretettel alkalmazott eljárás a légkör hidro-termodinamikai egyenletrendszerének megoldására. Elsősorban globális modellek esetén alkalmazzák a Földre jól illeszkedő szférikus harmonikus bázis függvények használata miatt, ugyanakkor vannak korlátos tartományú alkalmazások is (Haugen és Machenhauer, 1993; Radnóti, 2003). A spektrális módszer alkalmazása során az ismeretlen változókat ortogonális bázis függvények lineáris kombinációjaként állítjuk elő, és azt tételezzük fel, hogy a közelítés hibája ortogonális a bázis függvényekre. Ezen feltétel alkalmazásával az eredeti parciális differenciálegyenlet a lineáris kombinációban szereplő bázis függvények együtthatóira (spektrális együtthatók) vonatkozó közönséges differenciálegyenlet-rendszerré alakítható át. Mindazonáltal, a spektrális módszert csak ritkán alkalmazzák tisztán, azaz csak a spektrális együtthatókkal dolgozva, hiszen a nem-lineáris műveletek továbbra is rácsponti térben kerülnek végrehajtásra. A spektrális és rácsponti tér közötti, időlépcsőnként megvalósuló átmenetet ún. transzformációk révén valósítjuk meg és hajtjuk végre (ezek számítástechnikailag hatékony eljárások annak érdekében, hogy a modell végrehajtása megfelelően gyors maradjon). A nem-lineáris instabilitás elkerülése érdekében túlcsonkítást alkalmazunk, azaz kevesebb hullámot veszünk figyelembe, mint amennyit a rácspontok lehetővé tennének. A spektrális modellek legfontosabb előnyei a következők:

  • Az egyenletekben megjelenő térbeli deriváltak pontosan számolhatóak, s nem lép fel lineáris fázishiba.

  • A nem-lineáris tagoknál a nem-lineáris instabilitás elkerülhető (túlcsonkítás révén).

  • A spektrális együtthatók száma kisebb, mint a rácspontok száma, azaz a változók tárolása gazdaságosabb, mint a rácsponti modellek esetében.

  • A globális modelleknél nem jelenik meg a pólus probléma (a meridiánok konvergenciája következtében a pólusoknál jóval rövidebb integrálási időlépcső engedhető csak meg a stabilitás megőrzése érdekében, mint az alacsony szélességeken). A rácsponti modelleknél ez komoly nehézséget okoz, ami a szférikus harmonikus bázis függvények alkalmazásával elegánsan kiküszöbölhető.

  • A szemi-implicit séma Helmholtz-egyenletének megoldása spektrális modellekben triviális (mert a bázis függvények sajátértékei a Helmholtz-egyenletnek), míg a rácsponti modellekben hatékony módszerek nélkül költséges lehet a számítása.

A fenti előnyök mellett a spektrális modelleknek vannak hátrányai is (jóllehet ezek többsége kiküszöbölhető):

  • Az alkalmazott sémák bonyolultakká és nagyon költségessé válhatnak (a transzformációs módszer révén ez a probléma áthidalható).

  • A végrehajtandó műveletek száma a felbontás növekedésével gyorsabban növekszik, mint a rácsponti modellek esetében (ez komoly problémát jelenthet a jövőben a rácsfelbontás további javulása esetén).

  • Alkalmazásuk a korlátos tartományú modellek esetében nehézkesebb (ugyanakkor a szférikus harmonikusok „átírhatóak” mindkét horizontális irányban periodikus függvényekké, például úgy, ahogyan az ALADIN modell kifejlesztése során történt).

  • A spektrális csonkítás fizikailag értelmetlen mennyiségekhez is vezethet, pl. negatív nedvességtartalom (ez orvosolható a szemi-Lagrange sémával, amennyiben az az interpoláció során megőrzi a nedvesség pozitív értékeit).

  • Az erős gradiensű mezők (pl. a domborzat) reprezentációja „zajos” lehet a spektrális modellekben (erre különböző szűrési technikákat fejlesztettek ki).

Ahogyan a fentiekből látható a spektrális modellek jelentős előnyei megtartásával a hátrányai kiküszöbölhetőek, ezért ezt a modellfajtát napjainkban előszeretettel alkalmazzák a numerikus prognosztikában. A Magyarországon használt két legfontosabb modell, a globális ECMWF IFS modell, valamint a korlátos tartományú ALADIN modell is spektrális technikát alkalmaz.

II.3.3. Feladatok

Feladat 1: A φ (x) mennyiségre vonatkozó derivált közelítésére a centrált sémát alkalmazzuk. Taylor-sorfejtés alkalmazásával határozzuk meg a közelítés konzisztenciájának rendjét!

Megoldás 1: A konzisztencia rendjét a csonkítási hiba adja meg, amihez a diszkrét és a folytonos feladat különbségét kell venni. Ehhez alkalmazzunk Taylor-sorfejtést a centrált séma egyes tagjaira a (jn) pont körül:

Ezt helyettesítsük be a centrált séma alakjába:

Az így kapott diszkrét differenciál-operátorból kivonva a folytonos deriváltat, a következőt kapjuk a csonkítási hibára: . Ez azt jelenti, hogy a centrált séma az x-szerinti derivált másodrendben pontos közelítése, azaz konzisztenciájának rendje 2.

Feladat 2: Lássuk be a Neumann stabilitásvizsgálati módszer alkalmazásával, hogy a φ (x, t) mennyiségre vonatkozó egydimenziós lineáris advekciós egyenlet időbeli és a térbeli deriváltjainak közelítésére használt forward illetve baloldali sémák alkalmazása esetén a stabilitási kritérium a következő: .

Megoldás 2: A (II.27.) diszkrét feladat megoldását Fourier-sor alakban keresve, a feladat a következő módon írható fel:

Az egyszerűsítések elvégzése után az alábbi összefüggést kapjuk λk-ra, illetve abszolút értékének négyzetére:

A Neumann-módszer alapján a stabilitáshoz |λk |2 nem lehet nagyobb 1-nél, azaz a fenti összefüggés utolsó tagja nem lehet pozitív. Ez pontosan akkor teljesül, ha a Courant-szám 0 és 1 közé esik, azaz > 0 mellett .