med level set metoden

INOM

EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP

STOCKHOLM SVERIGE 2020 ,

Hantering av rörliga geometrier

med level set metoden

ISABEL WEDHOLM

MARCUS RICKENLUND

Abstract

This report describes what the level set method is and why it was developed, with a brief description of some suitable use cases.

The method is then described through an explanation of the overall concept, review of the mathematical basics and an account of how to implement the method numerically.

There are many numerical methods one can use to implement the level set method. In this project, we tested three such methods. The first one is named Lax–Friedrichs method, the second is a first-order Upwind scheme and the last is a second-order Upwind ENO scheme.

These methods are tested and the result shows that Lax–Friedrichs method is highly diffusive, and that the first order Upwind scheme is a more practical method.

The report also investigates the concept of reinitialization, which is first done by explaining what reinitialization means and why it can be useful. The concept of reinitialization is illustrated using some examples to gain a deeper understanding.

It is shown that applying reinitialization can reduce the error of the numerical solution. However, performing reinitialization also introduces some error by moving the zero contour. Therefore, an accurate reinitialization-procedure is very important when applying reinitialization on a problem.

The methods are implemented in Matlab, and the report contains some pseudo-code to explain the implementation steps.

Finally, it turns out that the method has flaws such as a high computational cost. What

it gives in return is that events such as separation and merging of different geometries are

handled automatically by the implicit nature of the level set method, a desirable trait in many

situations.

Sammanfattning

I denna rapport beskrivs vad level set metoden är och varför den kommit till. Även några lämpliga användningsområden tas upp och diskuteras.

Metoden beskrivs genom att förklara övergripande koncept, genomgång av de matematiska grunderna och redogörelse för hur man går till väga vid själva implementeringen.

Det finns flertalet numeriska metoder som kan tillämpas i level set metoden och vi redogör för tre stycken av dessa i denna rapport, vilka är Lax–Friedrichs metod, Upwind order 1 och Upwind med ENO order 2. Dessa metoder testas och resultatet visar på att Lax–Friedrichs metod är mycket diffus och att Upwind schemat är en mera praktisk metod att använda.

Rapporten redogör även för begreppet återinitialisering. Detta görs genom att till att börja med förklara vad en återinitialisering innebär och varför den kan komma att behövas. Begrepp illustreras även med hjälp av några exempel för att få fördjupad förståelse. Det visas att återinitialisering kan användas för att minska felet i en numerisk lösning. Dock inför

återinitialiseringen alltid ett litet fel då den flyttar på nollkonturen. Därför är en noggrann återinitialiserings-procedur viktig när man behöver tillämpa återinitialisering på ett problem.

De numeriska metoderna implementeras i Matlab och rapporten innehåller pseudokod som förklarar implementeringen.

Avslutningsvis visar det sig att level set metoden har vissa svagheter så som hög beräkningskostnad.

Det som dock är bra med den är att hantering av objektseparation och sammanfogning sker

utan svårigheter vilket är en önskvärd egenskap i många situationer.

Innehåll

1 Inledning 1

2 Teoribakgrund level set metoden 2

2.1 Övergripande koncept . . . . 2

2.2 Härledning av akvektionsekvationen . . . . 4

3 Implementering 5 3.1 Rumsdiskretisering . . . . 5

3.2 Tidsdiskretisering . . . . 7

3.3 Lax–Friedrichs metod . . . . 8

3.4 Upwind . . . . 9

3.5 ENO . . . . 11

3.6 Återinitiering . . . . 12

4 Numeriska resultat 18 4.1 Problem 1, expanderande cirkel . . . . 19

4.2 Problem 2, snurrande cirkel . . . . 22

4.3 Problem 3, applicering av återinitilitisering . . . . 25

4.4 Problem 4, kombination av olika typer av hastighetsfält . . . . 28

5 Slutsats och diskussion 32

6 Litteraturförteckning 33

1 Inledning

Det finns många numeriska metoder för att hantera rörliga geometrier. En rörlig geometri kan vara gränsskiktet mellan två olika ytor. En tillämpning där rörliga geometrier används är simulering av deformerbara objekt. Ett exempel på detta är att försöka återskapa hur vatten rör sig i en behållare. En annan tillämpning är bildbehandling [1].

Det som dessa metoder har gemensamt är att de gör det möjligt att simulera komplexa fysiska fenomen. Ett exempel på detta är fysikalisk modelleringen av en eldsflamma [2].

Ett sätt att tackla problem som innefattar rörliga geometrier är att använda explicita metoder [3].

Då definieras gränsskiktet av en mängd punkter och en parametrisering av dessa. Gränsskiktet kan då flyttas direkt genom att flytta punkterna. Ett annat sätt är att representera gränsskiktet implicit som en nivåkurva till en funktionsyta, detta innebär att en funktion i högre dimension än gränsskiktet måste definieras. De implicita metoderna medför således större beräkningskostnad, men ger i gengäld att kurvan inte behöver parametriseras. Detta gör att det blir mycket enklare att hantera topologiska förändringar dvs. när två olika geometrier går ihop och blir till en sammanfogad geometri. Ett exempel på när detta händer är när två vattendroppar slås samman till en enda. Level set är en av de implicita metoderna och den metod som kommer att granskas i denna rapport.

Level set metoden togs fram på 1980 talet av de amerikanska matematikerna Stanley Osher och James Sethian [2]. Det är en väl beprövad metod som blivit populär vid bildbehandling, datorgrafik och fysikalisk modellering. Level set metoden kan delas in i olika varianter, och den som hanteras i denna rapport kallas för standard level set metoden. Med denna variant definieras geometrin, som fortsättningsvis kommer kallas för Γ, implicit som nollkonturen till en avståndsfunktion φ där denna har värdet noll. En annan variant på metoden är en vid namn Conservative level set [4]. Där representeras funktionsytan istället av en generaliserad stegfunktion, dvs. en funktion som är noll i området utanför geometrin och ett innanför. Den ändras sedan ”snällt” mellan dessa värden.

I denna rapport kommer vi beskriva vad level set metoden är, varför den kommit till och berätta om några användningsområden. Detta kommer göras dels genom att redogöra för dess teoribakgrund med hjälp av en övergripande konceptuell beskrivning och genomgång av bakomliggande matematiska härledningar. Vi kommer även redogöra för hur level set metoden kan implementeras. Detta kommer göras med hjälp av Matlab [5]. Den största svårigheten har varit att en partiell differentialekvation i flera dimensioner behövts lösas för att kunna flytta på Γ. Denna ekvation kallas för advektionsekvationen och är en hyperbolisk PDE. Om inte detta moment genomförs noggrant kan egenskaper som att φ är en avståndsfunktion gå förlorad. För att få ytterligare förståelse för hur man tillämpar metoden kommer även fyra olika exempel problem att beskrivas, lösas och diskuteras. De olika problemen varierar brett för att ge en inblick i olika sätt att implementera metoden numeriskt.

Syftet med projektet är att få en grundläggande förståelse för hur metoden fungerar, vilka

fördelar och nackdelar som finns, utvecklingsmöjligheter och viktiga tillämpningsområden.

2 Teoribakgrund level set metoden

2.1 Övergripande koncept

Level set metoden är ett verktyg för att representera olika ytor och former i simuleringar.

Fördelen med denna metod är att den gör det möjligt att utföra numerisk beräkning för kurvor och ytor på en fixerad kartesisk grid, utan att behöva parametrisera dessa objekt.

Anta att det finns en punkt i R

. Med level set metoden modelleras punkten som en linje och punktens position ges implicit av den punkt där linjen skär y axeln. I två dimensioner, vilket hanteras i denna rapport, blir funktionen en yta och randen blir en kurva i planet. Randen representeras som en nivåkurva till funktionsytan. Ytan, planet och dess skärningspunkter dvs.

randen illustreras i Figur 1 nedan.

Figur 1: De tre övre bilderna visar ett plan som skär funktionsytan i olika nivåer och de tre undre bilderna visar samma sak fast sett från sidan så att hela funktionsytan blir synlig.

Level set metoden medför både fördelar och nackdelar. Jämfört med en explicit variant är metoden mer kostsam. Det blir mer kostsamma numeriska beräkningar jämfört med den explicita metoden på grund av att man representerar geometrin via en funktion som finns i en dimension högre. Fördelen med en implicit metod som level set är att hantering av separationen och kombinationen av objektet blir trivialt. Med en parametriserad rand uppstår ett stort problem då ytor kommer i kontakt med varandra. Level set metoden lämpar sig således bra för

simuleringsområden där detta är en vanligt förekommande händelse så som modelleringen av en eldflamma eller en droppe av olja som flyter på vattnet.

För att illustrera konceptet tittar vi på ett enkelt exempel där geometrin Γ är en cirkel med

radie R centrerad i origo, och ställer oss frågan hur ser funktionen φ verkligen ut? Svaret är att

den i teorin kan se ut hur som helst så länge nollkonturen inte ändras. Detta betyder således

att φ inte är unik. Däremot är det vanligt att φ definieras som en så kallad avståndsfunktion,

vilket innebär att värdet i en punkt (x, y) är minsta avståndet till Γ. Det är även praktiskt att definiera avståndsfunktionen så att tecknet för φ definierar huruvida en punkt befinner sig innanför eller utanför konturen.

Levelsetfunktionen φ(x, y, t) definieras som minsta avståndet d till geometrin Γ. Låt Ω

definieras som området innanför Γ. Detta kan uttryckas enligt:

φ(x, y, t) =



 

 

−d, (x, y) ∈ Ω

0, (x, y) ∈ Γ d, (x, y) / ∈ Ω

.

(1)

Detta innebär att beloppet av gradienten har värdet ett överallt dvs. |∇φ|= 1. I vårt exempel, det stationära fallet då Γ är en cirkel, fås att funktionen som uppfyller dessa krav ges av:

φ(x, y) = p

x

+ y

− R. (2)

Geometrin Γ definieras nu implicit som nollkonturen enligt:

Γ = n

(x, y) : φ(x, y) = 0 o

. (3)

I Figur 2 nedan visas φ(x, y) och nollkonturen Γ som är utmarkerad i blått.

Figur 2: Levelsetfunktionen φ(x, y) från ekv.(2) med R = 1 och dess nollkontur, en cirkel med radie R centrerad i origo, markerat med blått. Innanför den blå cirkeln har φ negativt värde och utanför ett positivt.

Genom att representera geometrin på detta vis kan geometriska kvantiteter som normalvektorer samt krökning räknas ut från φ. Ekvationen för normalvektorn ges av:

n = ∇φ

|∇φ| ∀x ∈ Γ, (4)

där n är normalvektorn som pekar utåt från kurvan med längden ett och ∇φ är gradienten till φ som ges enligt:

∇φ = ( ∂φ

∂x , ∂φ

∂y ). (5)

Implicita krökning dvs. måttet på nollkonturen ges av:

κ = ∇ · n = φ

φ

− 2φ

φ

φ

+ φ

φ

(φ

+ φ

)

∀(x, y) ∈ Γ, (6) där κ är krökningen, φ

, φ

, φ

, φ

, φ

är första, andra och blandade partiella derivator av φ. Ekv.(4) och ekv.(6) definierar en extension av normalen och krökningen till hela begränsnings området.

Det stationära fallet är dock inte speciellt intressant. Vi vill flytta på geometrin, och för att göra det måste φ ändras med tiden.

2.2 Härledning av akvektionsekvationen

Anta att φ är definierad enligt ekv.(1). Hela φ kan beskrivas som en mängd av nivåkurvor φ(x, y) = C, av vilka nollkonturen Γ är en delmängd. Här är C någon godtycklig konstant.

Frågan är hur punkterna som definierar φ ska ändras i tiden då Γ flyttas med någon hastighet U . Inför tidsvariabeln t som en parameter i φ. Alla nivåkurvor representeras då av:

φ(x(t), y(t), t) = C. (7)

För att lösa hur funktionen utvecklats med tiden differentieras funktionen med avseende på t.

Kedjeregeln ger då:

∂φ(x, (t), y(t), t)

∂t = ∂φ

∂x(t)

∂x(t)

∂t + ∂φ

∂y(t)

∂y(t)

∂t + ∂φ

∂t = 0. (8)

Hastighet och gradient kan faktoriseras till en skalärprodukt, vilket ger akvetionsekvationen:

∂φ

∂t = −V · ∇φ. (9)

Där hastigheten V = U på Γ. Om hastighetsfältet är givet som en normalhastighet V

till nivå kurvan är ändring till lokala koordinater enkel då skalärprodukten kan skrivas om som hastighetsvektorns längd i gradientvektorns, dvs. normalens riktning enligt:

∂φ

∂t = −|V |cos(θ)|∇φ|= −V

|∇φ|. (10)

Ett exempel på ett hastighetsfält i lokala koordinater kan vara att krökningen från ekv.(6) används som ett normalhastighetsfält. Krökningen ger då ett värde vid varje punkt. Detta värde beror på hur krökt konturen är, vilket kan tolkas som normalkomponenten av hastighetsfältet.

Resultatet har vissa konsekvenser. Problemet som level set metoden försöker lösa är evolutionen

av en kurva, eller en yta i tre dimensioner, men hastigheten V är nu ett fält som är definierat

och verkar över hela domänen. För att kurvan ska följa den korrekta lösningen är det enda kravet

att hastighetsfältet är definierat enligt problemuppställningen i ett område nära kurvan. Längre

ifrån kurvan är det möjligt att V kan vara nånting annat. I de flesta fall appliceras samma

hastighetsfält överallt, men att hitta ett annat fält som kan appliceras utanför kurvan är ett

aktivt forskningsområde [2].

3 Implementering

Teorin bakom advektionsekvationen är till sedes ganska enkel men numerisk implementering tillför en viss problematik. Vid en naiv implementering där gradienten av φ räknas ut med centrala differanser blir lösningen instabil. Den numeriska instabiliteten i advektionen av ytan sker på grund av att PDEn, (advektionsekvationen), är hyperbolisk. För att lösa detta används olika numeriska stenciler.

3.1 Rumsdiskretisering

Representation av funktionsytan

I en dator måste φ och dess gradient diskritiseras. Antag att ett område [a,b]x[c,d] i R

ska diskritiseras uniformt i N − 1 intervall i x-led och M − 1 interval i y-led. Detta ger M xN diskreta punkter i en grid. Upplösningen i x-led och y-led är ∆x =

respektive ∆y =

. Låt φ

beteckna φ(x

, y

) där i och j är heltals-index på intervallet (1 ≤ i ≤ N ) och (1 ≤ j ≤ M ).

Figur 3: Här visas hur diskretiseringen ska tolkas geometriskt.

Diskretisering av hastighetsfältet

Utförs på samma sätt som för funktionsytan φ. Detta görs så att matrisindexen i, j representerar samma koordinat som i funktionsytan.

Hur representeras då de partiella rumsderivatorna?

Här används finita differenser, ett sätt att numeriskt räkna ut derivator. I det enklaste exemplet räknas derivatan av en funktion f (x) enligt:

df

dx = f (x + ∆x) − f (x)

∆x + O(∆x). (11)

Denna approximation av derivatan har noggrannhetsordning 1. En annan viktig egenskap hos

denna approximation är att den är frammåtriktad, dvs. att information om funktionen i och till

höger om punkten x används. På samma vis är en funktion som använder information i och till vänster om x bakåtriktad.

I det här fallet är funktionen en yta som beror på två diskreta rumsparametrar x

och y

, där φ

är värdet för φ vid den diskreta punkten (x

, y

). Finita differenser appliceras i dessa punkter.

Vi får då följande stenciler för de partiella derivatorna av φ.

Derivata Stenciler

Följande stenciler approximerar

upp till noggranhetordning 2.

Framåtdifferans och bakåtdifferens med noggrannhetsoring 1 ges av:

φ

(x

, y

) = φ

− φ

∆x + O(∆x), (12)

φ

(x

, y

) = φ

− φ

∆x + O(∆x). (13)

Skeva differensformler följt av centraldifferens med noggrannhetsoring 2 ges av:

φ

(x

, y

) = −3φ

+ 4φ

− φ

2∆x + O(∆x

), (14)

φ

(x

, y

) = 3φ

− 4φ

+ φ

2∆x + O(∆x

), (15)

φ

(x

, y

) = φ

− φ

2∆x + O(∆x

). (16)

För differentiering i y-led görs samma sak, men i är konstant och indexet j ändras. Steglängden

∆y används istället för ∆x.

Finita differenser används även för att approximera högre ordningens partiella derivator. I denna rapport används enbart centrerade stenciler för den andra partiella derivatan i x- och y-led samt för blandderivatan. Dessa ges av:

φ

(x

, y

) = φ

− 2φ

+ φ

∆x

+ O(∆x

), (17)

φ

(x

, y

) = φ

− 2φ

+ φ

∆y

+ O(∆x

), (18)

φ

(x

, y

) = φ

− φ

− φ

+ φ

4∆x∆y + O(∆x

). (19)

Då dessa approximationer kommer användas mycket i den numeriska implementeringen av

advektionsekvationen införs D

φ

, D

φ

och D

φ

för att indikera olika finita differenser

av φ i punkten x

och y

, där index p visar på vilken partiell derivata som approximeras och

potensen redogör om det är framåt(+), bakåt(-) eller centrerad(ingen potens) differens.

3.2 Tidsdiskretisering

Funktionsytan måste även diskretiseras med avseende på tiden. Liksom rumsdiskretiseringen indikeras varje steg i tiden med ett heltal, n och längden mellan tidstegen med ∆t. För stegning i tiden används explicita stegmetoder. Den första metoden är Euler framåt. Låt φ

och F (φ

, t

) beteckna funktionsytan samt dess tidsderivata vid det diskreta tidsteget t

. Med Euler fram definieras diskretiseringen av φ

= F (φ, t) som:

φ

= φ

+ ∆tF (φ

, t

). (20)

Då rumsderivator av andra noggranhetsordningen appliceras så används också en RK metod av ordning 2, TVD-RK2 (Total Variation Diminishing Runge Kutta), en vanlig metod i dessa problem [6]. Funktionsytan vid steget n + 1 ges av medelvärdet mellan två virtuella Euler steg, φ ˆ

och ˆ φ

och startvärdet är φ

. Tidsteget φ

ges av:



 

 

φ ˆ

= φ

+ ∆tF (φ

, t

) φ ˆ

= φ

+ ∆tF ( ˆ φ

, t

) φ

= (φ

+ ˆ φ

)/2.

(21)

Algorithm 1 Lösningsprocess

1:

procedure Initiering

2:

[~ x, ~ y] ← N xM punkter i ett nät

3:

for alla punkter (x

, y

) do

4:

φ(i,j) ← φ

(x

, y

)

5:

∆x ← x

− x

6:

T ← Tid då lösningen ska sluta stega

7:

∆t ← ∆x · (En konstant)

8:

n ← round(T /∆t)

9:

t ← 0

10:

∆t ← T /n

11:

F (φ, t) ← definiera reducerad funktion från: F (φ, f

, f

, f

, t, ∆t, ∆x, ∆y)

12:

procedure Tidstegning

13:

for n steg y

do

14:

φ ← output från tidsstegfunktion(φ, F, t, ∆t)

15:

t ← t + ∆t

1:

procedure tidstegning Euler

2:

input ← φ, F, t, ∆t

3:

output ← φ + ∆t · F (φ, t, ∆t)

1:

procedure tidstegning RK2

2:

input ← φ, V

, V

, V

, K, ∆x, ∆t

3:

φ ← output från Eulersteg(φ, F, t, ∆t) ˆ

4:

φ ← output från Eulersteg( ˆ ˆ φ, F, t, ∆t)

5:

output ← (φ + ˆ φ)/2

f

,f

, f

är godtyckliga funktioner för hastighetsfälten V

, V

och V

som i varje punkt beror

på ytan φ, tid t, eller om fälten är konstanta enbart på rumskordinaterna x

och y

. ’Reducerad

funktion’ på rad 11 innebär här enbart att antalet inputs till F reduceras till enbart de som

ändras i varje steg.

3.3 Lax–Friedrichs metod

Lax–Friedrichs metod är en av de numeriska metoderna som används för att lösa en hyperbolisk partiell differential ekvation och är baserad på finita differenser. Metodens tidsdiskretisering liknar Euler framåt, med en liten men viktig ändring. För att hantera den numeriska instabiliteten av advektionsekvationen, φ

= −V

· ∇φ − V

|∇φ|, där V

är ett externt hastighetsfält och V

är ett normalhastighetsfält, så medelvärdesbildas varje punkte φ

från omkringliggande punkter enligt följande:

φ

= φ

+ φ

+ φ

+ φ

4 − ∆t

 V

i,j

· D

φ

D

φ

 + V

i,j

D

φ

D

φ

 . (22) Här används centrerade finita differenser av ordning 2. Se kap 3.1 för djupare beskrivning.

En uppenbar konsekvens av Lax–Friedrichs metod är att om rumsupplösningen inte ökas när tidsteget minskar så kommer felet i lösningen öka på grund av att fler medelvärdesbildningar utförs under samma tidsintervall. Lax–Friedrichs metod har förväntad noggrannhetsordning 1 i både rum och tid.

Pseudokod: För att slippa extrapolera φ två gånger i kodimplementeringen har vi valt att stega Lax–Friedrichs med samma procedur som Euler fram i algoritm 1. Medelvärdesbildningen i stencilen utförs fortfarande, men gör i derivatafunktionen F , vilket är förklaringen till φ

på rad 16 och 17 i algoritmen för ax–Friedrichs metod som implementeras på följande vis:

Algorithm 2 Lax-Fredrich metod

1:

procedure F funktion

2:

input ← φ, f V

(φ, t), f V

(φ, t), f V

(φ, t), t, ∆t, ∆x, ∆y

3:

φ

← definiera matris 2 högre och bredare än φ

4:

φ

← Inre punkter = φ. Interpolera yttre punkter

5:

V

← f V

(φ, t)

6:

V

← f V

(φ, t)

7:

V

← f V

(φ, t)

8:

D

, D

, |D|, φ

← definiera matriser av samma storlek som φ

9:

N ← antalet kolumner i φ

10:

M ← antalet rader i φ

11:

for i = 1 to N do

12:

for j = 1 to M do

13:

D

(i, j) ← (φ

(i + 1, j) − φ

(i − 1, j))/(2∆x)

14:

D

(i, j) ← (φ

(i, j + 1) − φ

(i, j − 1))/(2∆y)

15:

|D|(i, j) ← pD

(i, j)

+ D

(i, j)

16:

φ

(i, j) ← (φ

(i − 1, j) + φ

(i + 1, j) + φ

(i, j − 1) + φ

(i, j + 1))/∆t

17:

φ

(i, j) ← φ

(i, j) − φ

(i, j)/∆t

18:

output ← V

 D

+ V

 D

+ V

 |D| + φ

Här är  Hadamard produkten, dvs. elementvis multiplikation. Funktionerna f V

, f V

och

f V

kan bero på φ som i fallet då ett krökningsdrivet hastighetfällt appliceras. Notera att

indexeringen φ

egentligen bör vara en större för både i och j, men här ignoreras detta för att

hålla pseudokoden lättläst. För beskrivning termerna D se kap 3.1.

3.4 Upwind

Upwind metoden löser stabilitetsproblemet genom att använda framåt och bakåt finita differenser beroende på vilken riktning hastighetsfältet har. Om hastighetsfältet är positivt i en punkt används framåtriktad differens och tvärt om då hastighetsfältet är negativt. Detta kallas Upwind [2].

Betrakta det endimensionella fallet av advektionsekvationen:

φ

= −V

φ

, (23)

där V

är hastigheten. Låt φ

och φ

notera rumsderivatan i uppåt, respektive nedåt riktning.

Applicering av regeln om Upwind differentiering ger:

φ

=

( −V

φ

, V

< 0

−V

φ

, V

> 0. (24)

För att använda upwind schemat i level set metoden måste principen i ekv.(24) appliceras i ett större antal rumsdimentioner. Eftersom hastighetsfältet är givet som ett externt fält i x och y riktning blir lösningen enkel då skalärprodukten resulterar i att hastigheten i en dimension enbart multipliceras med rumsderivatan i samma riktning. Då fås:

V · ∇φ = V

φ

+ V

φ

. (25)

Detta gör att upwind schemat kan hanteras separat i varje dimension och sedan summeras. Med max och min funktioner kan ekv.(25) skrivas kompakt enligt:

V · ∇φ = max(V

, 0)D

+ min(V

, 0)D

+ max(V

, 0)D

+ min(V

, 0)D

, (26) här är max och min funktionerna definierade så att operationen görs elementvis. Argumentet kommer även vara två matriser av samma storlek, på samma sätt som Hadamard produkten . Se kap 3.1 för definition av finita differens operatorn D. Hantering av hastighet i normalriktningen blir dock mer komplicerad. Problemet är att lista ut hur höger och vänster derivator ska kombineras på rätt sätt. Se återigen advektionsekvationen i ekv.(9). Om hastighetsvektorn V är normalhastigheten ges den av:

V = V

∇φ

|∇φ| . (27)

Det som är intressant är tecknet på hastigheten V i x och y riktning. Anta först att V

från det givna normalhastighetsfältet är positivt i en punkt. Applicering av upwind regeln ger fyra olika fall enligt:

V

|∇φ|=



 

 

 



 

 

 

 V

q

(φ

)

+ (φ

)

, φ

< 0 och φ

< 0 V

q

(φ

)

+ (φ

)

, φ

< 0 och φ

>= 0 V

q

(φ

)

+ (φ

)

, φ

>= 0 och φ

< 0 V

q

(φ

)

+ (φ

)

, φ

>= 0 och φ

>= 0

V

>= 0. (28)

Notera att framåt och bakåt differenser oftast kommer ge samma tecken och då kan alla fall slås ihop till en enda rad. Undantagsfallen där framåt och bakåt differenser ger olika tecken hanteras genom att ta det största värdet enligt:

V

|∇φ|= V

q

max(min(D

, 0)

, max(D

, 0)

) + max(min(D

, 0)

, max(D

, 0)

), V

>= 0.

(29)

I det fallet då V

är negativt byter samtliga olikheter i ekv.(28) riktning, vilket leder till att min blir max och max blir min för varje differens i ekv.(30). Notera dock att det största bidraget används i undantagsfallet:

V

|∇φ|= V

q

max(max(D

, 0)

, min(D

, 0)

) + max(max(D

, 0)

, min(D

, 0)

), V

< 0.

(30) När både ett externt hastighetsfält samt ett normalhastighetsfält är givet summeras resultatet från ekv.(26), ekv.(29) och ekv.(30).

Pseudokod för implementeringen av Upwind. Här är G

, G

och G

approximationer av V

|∇φ|, V

φ

respektive V

φ

.

Algorithm 3 Upwind

1:

procedure F funktion

2:

input ← φ, f V

(φ, t), f V

(φ, t), f V

(φ, t), t, ∆t, ∆x, ∆y

3:

φ

← definiera matris 2 högre och bredare än φ

4:

φ

← Inre punkter = φ. Interpolera yttre punkter

5:

V

← max(f V

(φ, t), 0)

6:

V

← max(f V

(φ, t), 0)

7:

V

← max(f V

(φ, t), 0)

8:

V

← min(f V

(φ, t), 0)

9:

V

← min(f V

(φ, t), 0)

10:

V

← min(f V

(φ, t), 0)

11:

G

, G

, H, ← definiera matriser av samma storlek som φ

12:

D

, D

, D

, D

← definiera matriser av samma storlek som φ

13:

N ← antalet kolumner i φ

14:

M ← antalet rader i φ

15:

for i = 1 to N do

16:

for j = 1 to M do

17:

D

(i, j) ← (φ

(i + 1, j) − φ

(i, j))/∆x

18:

D

(i, j) ← (φ

(i, j + 1) − φ

(i, j))/∆y

19:

D

(i, j) ← (φ

(i, j) − φ

(i − 1, j))/∆x

20:

D

(i, j) ← (φ

(i, j) − φ

(i, j − 1))/∆y

21:

G

← V

 D

+ V

 D

22:

G

← V

 D

+ V

 D

23:

H

← q

max(max(D

, 0)

, min(D

, 0)

) + max(max(D

, 0)

, min(D

, 0)

)

24:

H

← q

max(min(D

, 0)

, max(D

, 0)

) + max(min(D

, 0)

, max(D

, 0)

)

25:

G

← V

 H

+ V

 H

26:

output ← G

+ G

+ G

3.5 ENO

Det är möjligt att nå högre noggranhetsordning genom att applicera högre ordningens finita differenser i upwind schemat. Detta resulterar dock i minskad stabilitet. ENO, (essentially non oscillatory), är en metod som försöker lösa problemet genom att i varje punkt välja en av flera finita differenser för att approximera rumsderivatorna.

Vi har valt att implementera en simpel variant från [6]. Schemat blir mycket likt Upwind schemat förutom en liten ändring. Framåt och bakåt finita differenser för första partiella rumsderivatorna definieras om på följande vis:

D

φ

= φ

− φ

∆x − minmod(D

φ

, D

φ

), (31) D

φ

= φ

− φ

∆x + minmod(D

φ

, D

φ

), (32) D

φ

= φ

− φ

∆y − minmod(D

φ

, D

φ

), (33) D

φ

= φ

− φ

∆y + minmod(D

φ

, D

φ

), (34) där D

φ

betecknar andra ordningens centrala finita differensapproximation av φ

i punkten x

, y

och tidsteget t

.

Funktionen minmod definieras på följande vis:

minmod(a, b) =



 

 

0, ab < 0

min(a, b), a > 0 & b > 0 max(a, b), a ≤ 0 & b ≤ 0,

(35)

där a och b är reela tal.

3.6 Återinitiering

När funktionsytan advekteras kommer fel i beräkningarna resultera i att avståndsfunktionen inte bevaras, dvs. att |∇φ|= 1 inte längre gäller. Väldigt stora eller små värden på ∇φ leder till stora fel i representationen av geometrin samt beräkningen av normal och krökning.

Ett sätt att bibehålla avståndsfunktionen är att utföra en återinitialiserings procedur där syftet är att återställa avståndsfunktionen utan att flytta nollkonturen. Detta kan göras t.ex. genom att lösa följande PDE tills φ inte ändras längre:

φ

= sgn(φ

)(|∇φ|−1), (36)

här är φ

levelsetfunktionen, som inte längre är någon avståndsfunktion, den förvrängda ytan och sgn(φ

) är den vanliga teckenfunktionen som ger antingen -1, 0 eller 1 då φ

är negativt, 0 eller positivt. För att implementera detta fält med ENO stencilerna kan sgn(φ

) tolkas som en normalhastighet. Ett problem som uppstår vid implementering är numeriska fel från den abrupta ändringen i hastighetsfältet. För att lösa detta modifieras teckenfunktionen till en snällare variant. Den vanliga teckenfunktionen och två snälla varianter ges av:

sgn(φ

) =



 

 

1, φ(i, j) > 0 0, φ(i, j) = 0

−1, φ(i, j) < 0,

(37)

sgn(φ

) = φ

q

φ

+ ∆x

, (38)

sgn(φ

) = φ

q

φ

+ ∆x

|∇φ|

. (39)

Ett problem är att punkter nära nollkonturen är fria att flytta sig. Ett sätt att förbättra detta är att införa gridpunkter vid nollkonturen mellan dom diskreta punkterna där levelsetfunktionen φ byter tecken. När funktionsytan advekteras används dessa virtuella punkter vid uträkningen av gradienten, och på så sätt låses funktionsytan fast [6].

För att bibehålla stabilitet ändrades advektionshastigheten vid de kritiska punkterna proportionellt mot kvoten av ˆ ∆x/∆x, där ˆ ∆x ger det minsta avståndet till randen vid dessa punkter.

Pseudokod: När återinitialisering utförs utan att låsa nollkonturen appliceras bara någon av sgn(φ

) funktionerna som normalhastighetsfält. Lösningen stegas med ENO metoden i rum och TVD-RK2 metoden i tid. För att låsa nollkonturen måste först positionen av dessa nollpunkter hittas. Eftersom lösningen av återinitialiseringen beror på finita differenser i x och y riktning undersöks enbart nollpunkterna längsmed nätlinjerna i rumsdiskritiseringen. Detta görs genom att ta fram de diskreta punkter där φ byter tecken.

En nollpunkt ˆ x finns mellan φ

och φ

då φ

φ

< 0.

Så samma vis finns en nollpunkt ˆ y mellan φ

och φ

då φ

φ

< 0.

När en sådan punkt hittas approximeras avståndet ˆ ∆x och ˆ ∆y från φ

till nollpunkten med linjär approximation enligt:

∆x = − ˆ φ

∆x φ

− φ

, (40)

∆y = − ˆ φ

∆y

φ

− φ

. (41)

Indexen för dessa punkter i och j, samt tillhörande avstånd till nollkonturen läggs in i två listor, en för alla ˆ x och en för alla ˆ y. När dessa punkter används vid uträkningen av gradienten kommer mindre tidsteg krävas för att få en stabil lösning. Detta är på grund av stabilitetsvilkoret i ekv.(45). Eftersom hastighetsfältet, dvs. sgn(φ

), har max belopp 1 beror detta enbart på hur

∆x och ˆ ˆ ∆y förhåller sig till ∆x respektive ∆y.

Istället för att minska tidsteget ändrar vi istället hastighetsfältets storlek i de relevant punkterna.

Vi definierar en matris R där alla punkter har värdet 1, förutom punkterna ˆ x och ˆ y. I dessa punkter har istället värdet av kvoten mellan avståndet till nollpunkten och ∆x eller ∆y. Notera att en punk φ

kan samtidigt ha nollpunkter till höger, vänster samt uppåt och ner. I dessa punkter där flera fall kolliderar används den minsta kvoten av samtliga fall. När detta är utfört används matrisen R för att reducera värdet på hastighetsfältet så att:

V

= sign(φ

)R

. (42)

Algoritmen för ovannämnda process implementeras på följande vis.

Algorithm 4 Reinit hitta kritiska punkter

1:

procedure criticalPoints funktion

2:

input ← φ, ∆x, ∆y

3:

xlist ← tom lista för alla kritiska punkter i x led ˆ

4:

ylist ← tom lista för alla kritiska punkter i y led ˆ

5:

for alla punkter (x

, y

) förutom sista do

6:

if φ

φ

< 0 then

7:

xlist ← append(i, j, ˆ ˆ ∆x)

8:

if φ

φ

< 0 then

9:

ylist ← append(i, j, ˆ ˆ ∆y)

10:

R ← matris med ettor av samma storlek som φ

11:

for alla i, j, ˆ ∆x i ˆ xlist do

12:

if R(i, j) > ˆ ∆x/∆x then

13:

R(i, j) ← ˆ ∆x/∆x

14:

if R(i + 1, j) > (∆x − ˆ ∆x)/∆x then

15:

R(i + 1, j) ← (∆x − ˆ ∆x)/∆x

16:

for alla i, j, ˆ ∆y i ˆ ylist do

17:

if R(i, j) > ˆ ∆y/∆y then

18:

R(i, j) ← ˆ ∆y/∆y

19:

if R(i, j + 1) > (∆y − ˆ ∆y)/∆y then

20:

R(i, j + 1) ← (∆y − ˆ ∆y)/∆y

21:

output ← ˆ xlist, ˆ ylist, R

Detta används sedan i återinitialiseringsfunktionen, som gör n antal steg med TVD-RK2, (Runga Kutta 2). Stegningen görs med funktionen ˆ F , en modifierad implementering av derivatafunktionen F från ENO implementeringen.

Algorithm 5 Återinitialisering

1:

procedure reinit funktion

2:

input ← φ, ∆x, ∆y, n

3:

xlist, ˆ ˆ ylist, R ← output från criticalP oints(φ, ∆x, ∆y)

4:

V

← R  sign(φ)

5:

∆τ ← 0.6min(∆x, ∆y)

6:

for n antal steg do

7:

φ ← φ − ∆τ ˆ ˆ F (φ, V

, ∆x, ∆y, ˆ xlist, ˆ ylist)

8:

φ ← ˆ ˆ φ − ∆τ ˆ F (φ, V

, ∆x, ∆y, ˆ xlist, ˆ ylist)

9:

φ ← (φ + ˆ φ)/2

10:

output ← φ

Funktionen ˆ F är modifierad från ENO implementeringen så att den först räknar ut alla derivatamatriser D

, D

, D

samt D

som vanligt. Därefter går algoritmen igenom listorna ˆ x och ˆ y. För varje i,j och ˆ ∆x i ˆ x listan görs följande

D

φ

= − (∆x

− ˆ ∆x

)φ

+ ∆x ˆ

φ

∆x ˆ ∆x(∆x − ˆ ∆x) , (43)

D

φ

= (∆x

− ˜ ∆x

)φ

+ ∆x ˜

φ

∆x ˜ ∆x(∆x − ˜ ∆x) , (44)

där ˜ ∆x = ∆x − ˆ ∆x. Dessa är uppåt och nedåt riktade finita differenser av nogranhetsordning två då φ i punkten (x

+ ˆ ∆x, y

) antas vara 0. En liknande uträkning av D

φ

och D

φ

görs med avsende på alla ˆ y listan. Som tidigare nämnt tenderar enbart högre ordningens finita

differenser i upwind schema försämra stabiliteten, men det verkade inte ha någon större effekt i

detta fall.

För att testa hur en återinitilitisering fungerar kan avsiktligt avståndsfunktion försämras genom att den multipliceras med en annan ekvation så att nollkonturen fortfarande är den samma. Då behöver en återinitiering införas för att återställa avståndsfunktion från nollkonturen. För att testa ett sådant problem replikeras ett problem från Chohong Mins artikel om ämnet [6] där två cirklar skär varandra. Observera den korrekta avståndsfunktionen jämfört med den förvrängda som visas i Figur 4 nedan.

Korrekt avståndsfunktion

-2 -1 0 1 2

x -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

y

Förstörd avståndsfunktion

-2 -1 0 1 2

x -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

y

Figur 4: Vänstra bilden visar korrekt avståndsfunktion med 20 nivåkurvor utmarkerade och

bilden till höger visar den förstörda avståndsfunktionen, dvs. då den multiplicerats med en

annan funktion med 80 nivåkurvor. Problemet i återinitiering är att gå från den högra figuren

till den vänstra utan att ändra på den blåa nollkonturen.

Vi tillämpar de fyra ovannämnda varianterna på återinitiering på den förstörda avståndsfunktionen från Figur 4, dvs. den vanliga sign funktionen enligt ekv.(37), de två snällare varianterna som visas i ekv.(38) och ekv.(39) samt metoden då nivåkurvan låses fast. I Figur 5 här nedan kan det ses att de olika varianterna ger varierande resultat på hur nollkonturen påverkas.

Standard sign function

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

sgn( ₀) = ₀/(dx² + ₀²)

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

sgn( ₀) = ₀/(dx²| ₀|² + ₀²)

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5

2 Standard sign med konturlåsning

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5 2 Startkontur Återinitialiserad

Figur 5: Figuren visar hur nollkonturen ser ut efter återinitiering med de olika metoderna. Här

kördes 512 iterationer så att φ inte längre märkbart ändras vid vidare iterationer. Resultatet

visar utan djupare analys att de första tre varianterna märkbart ändrar på nollkonturen i det

här fallet.

Återinitialiserings-metoden med låsning av kritiska punkter hade betydligt mycket bättre resultat än de utan. Denna hade mycket litet fel och återskapar nästan helt den förstörda avståndsfunktionen på inte allt för många iterationer. Processen illustreras i Figur 6 nedan.

0 steg

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5

2 25 steg

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

50 steg

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5

2 75 steg

-2 -1 0 1 2

y -2

-1.5 -1 -0.5 0 0.5 1 1.5 2

Figur 6: Här visas utveckling efter tillämpning av återinitieringen. Längst upp till vänster visas

ursprungsläget, figurerna efter visar hur konturen ändras när man löser återinitialiseringen,

ekv.(36), med 25, 50 och 75 steg

4 Numeriska resultat

Stencilerna för Upwind, Upwind med ENO samt Lax–Friedrichs metod, som hädanefter kommer att kallas LAX, appliceras först på två enkla typ exempel, problem 1 och problem 2. Första problemet gäller en expanderade cirkel och det andra en cirkel som roterar kring en mittpunkt.

Då dessa problem har relativt enkla analytiska lösningar för φ(x, y, t) som kan härledas intuitivt kommer resultatet från den numeriska akvektionen jämföras med de analytiska lösningarna för att ta fram nogrannhetsorningen. Stabilitetsvillkoren för metoderna testates genom att pröva olika ∆t och se vart det blir instabilt. I teorin bör samtliga metoder ha stabilitetsvilkoret:

∆t < cf ∆x

max(V ) (45)

där cf är en en liten positiv konstant som är något mindre än ett, och max(V ) är största värdet i hastighetsfältet. Tidsteget ∆t kommer anges som produkten av en konstant och ∆x.

Därefter testas återinitilitisering i ett tidsvariant problem, problem 3, där meningen är att se om återinitilitisering kan förbättra resultatet på en lösning. Det problem som testas ett vanlig problem i från litteraturen [7] där ett normalhastighetsfält som varierar i tid och rum får stora fel i resultatet på grund av uppskattningen av normalen.

I problem 4 testas avslutningsvis några olika kombinationer av hastighetsfält. Undersökningen görs genom att titta på vad som händer när följande fall läggs till på akvektionsekvationen:

• Fall 1, Krökningsdrivet fält på en stjärnformat nollkontur

• Fall 2, Separation av en cirkel.

• Fall 3, Normalexpansion i ett geometriskt begränsat utrymme.

Samtliga deformationer redovisas och diskuteras. I fall 1 och 2 appliceras även ett krökningsdrivet hastighetsfält. Ett problem som uppstår är mycket höga värden nära singulariteter. Om φ är t.ex avståndsfunktionen till en cirkel är krökningen 1/r, vilket ökar obegränsat nära cirkelns mitt. För att undvika detta används en liten positiv konstant  så att krökningen istället blir 1/(r + ). Krökningen ˆ κ som appliceras som hastighetsfält räknas då ut som:

ˆ

κ = κ

|κ|+1 (46)

4.1 Problem 1, expanderande cirkel

Levelsetfunktionen φ representerar initialt en cirkel med centrum i origo och radien r

. Cirkeln expanderar längs normalen n med konstant hastighet V

. För detta problem löses således följande ekvation:

φ

+ V

|∇φ|= 0, (47)

där V

är ett konstant hastighetsfält i normalriktingen. Här är både r

och V

givna till 0.3.

Den exakta lösningen till problemet alltså det som de numeriska lösningarna kommer jämföras emot ges enligt:

φ(x, y, t) = p

x

+ y

− r(t) = p

x

+ y

− (r

+ V

t). (48) I Figur 7 nedan visas hur cirkeln expanderar under tiden T när de olika stencilerna LAX, Upwind och Upwind med ENO tillämpas för att lösa problemet numeriskt. Även den korrekta lösningen finns representerad. Antalet gridpunkter ges av N

och här är N satt till 17. Tidssteget ∆t har valt till

.

T = 0.25

-0.5 0 0.5

x -0.6

-0.4 -0.2 0 0.2 0.4 0.6

y

T = 0.5

-0.5 0 0.5

x -0.6

-0.4 -0.2 0 0.2 0.4 0.6

y

T = 0.75

-0.5 0 0.5

x -0.6

-0.4 -0.2 0 0.2 0.4 0.6

y

T = 1

-0.5 0 0.5

x -0.6

-0.4 -0.2 0 0.2 0.4 0.6

y

Lax Upwind ENO exakt

Figur 7: De fyra olika figurerna visar hur LAX, Upwind och Upwind med ENO förhåller sig mot

den exakta lösningen vid olika tidpunker under cirkelns expansion. Här är N = 17 och ∆t har

valts till

I Figur 8, precis här under, jämförs återigen de olika metoderna mot den korrekta lösningen.

Resultatet visas för tre olika val av N där N

är totala antalet gridpunkter.

N = 17

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 x

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

y

N = 17

0.4 0.42 0.44 0.46

x 0.39

0.4 0.41 0.42 0.43 0.44 0.45 0.46

y

N = 33

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 x

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

y

N = 33

0.4 0.42 0.44 0.46

x 0.39

0.4 0.41 0.42 0.43 0.44 0.45 0.46

y

N = 65

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 x

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

y

N = 65

0.4 0.42 0.44 0.46

x 0.39

0.4 0.41 0.42 0.43 0.44 0.45 0.46

y

Lax Upwind ENO exakt

Figur 8: Här visas samtliga tre numeriska stenciler vid t = 1 samt den exakta lösningen för olika val av N och ∆t =

N

. De högra bilderna är inzoomningar av de vänstra.

I Figur 9, som återfinns här nedan, visas felet i position och area vid tiden t = 1 som funktion av N .

10² 10³

N 10^-6

10^-5 10^-4 10^-3 10^-2 10^-1

fel

Avståndsfel loglog plot

10² 10³

N 10^-6

10^-5 10^-4 10^-3 10^-2 10^-1

areafel

Areafel loglog plot Lax

Upwind ENO ref: 1/N ref: 1/N²

Figur 9: Visar en studie av noggrannhetsordingen för de tre olika metoderna. I figuren finns även guidelinjer representerade för att underlätta avläsningen, och noggrannhetsordningen är som förväntat 1 för Lax–Friedrichs och Upwind samt 2 för Upwind med ENO.

När N ökas ses att felet avtar och noggrannhetsordningen för metoderna hamnar på den förväntade. För stabilitet visade det sig att ∆t =

n

var litet nog för samtliga metoder.

4.2 Problem 2, snurrande cirkel

En cirkel med konstant radie roterar kring origo med vineklhastigheten ω.

När t = 0 har cirkeln radien 0.3 och är centrerad i x = −0.4, y = 0.

För detta problem läggs ett externt hastighetsfält på och advektionsekvationen ges då av:

φ

+ V · ∇φ = 0, (49)

där V är det pålagda externa hastighetsfältet som i varje position ges av:

V

(x, y) V

(x, y)



= −2πy 2πx



. (50)

Den exakta levelsetfunktionen ges av:

φ(x, y, t) = p

(x + L cos(t2πω))

+ (y + L sin(t2πω))

− r. (51) Resultatet från level set metoden med Upwind, Upwind med ENO och den exakta lösningen vid tidpunkterna 0.25, 0.5, 0.75 och 1 redovisas i Figur 10 nedan. Antalet gridpunkter är N

och här är N satt till 51.

-0.5 0 0.5

y -0.8

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

T = 0.25

Exakt Lax Upwind ENO

-0.5 0 0.5

y -0.8

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

T = 0.5

-0.5 0 0.5

y -0.8

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

T = 0.75

-0.5 0 0.5

y -0.8

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

T = 1

Figur 10: Advektion av nollkonturen vid olika tidpunkter med hastighetsfältet i ekv.(50). Här är antalet gridpunkter N = 51 och ∆t = ∆x/(3 √

2π)

.

Då LAX stencilen är väldigt diffus medförde det att cirkeln nästan försvann direkt. Detta syns tydligt i Figur 8. Således utesluts den stencilen vid kommande beräkningar.

I Figur 11, precis här under, jämförs återigen de olika metoderna mot den korrekta lösningen.

Resultatet visas för tre olika val av N där N

är totala antalet gridpunkter.

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 y

-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

N = 33

-0.3 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 x

-0.43 -0.42 -0.41 -0.4 -0.39 -0.38 -0.37

y

N = 33

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 y

-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

N = 65

-0.3 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 x

-0.43 -0.42 -0.41 -0.4 -0.39 -0.38 -0.37

y

N = 65

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 y

-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

N = 129

-0.3 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 x

-0.43 -0.42 -0.41 -0.4 -0.39 -0.38 -0.37

y

N = 129 Exakt

Upwind ENO

Figur 11: Visar Upwind, Upwind med ENO samt den korrekta lösningen. Resultatet visas för tre olika val av N där N

är totala antalet gridpunkter och ∆t = ∆x/(3 √

2π). De högra bilderna

är inzoomningar av de vänstra.

I Figur 12 här nedan visas en loglog plot av felet i position och area vid tiden t = 0.25 som funktion av N.

10² 10³

fel 10^-5

10^-4 10^-3 10^-2 10^-1 10⁰

Avståndsfel loglog plot

10² 10³

areafel 10^-5

10^-4 10^-3 10^-2 10^-1 10⁰

Areafel loglog plot

Lax Upwind ENO ref: 2/N ref: 10/N²

Figur 12: Visar en studie av noggrannhetsordingen för de tre olika metoderna. I figuren finns även guidelinjer representerade för att underlätta avläsningen, och noggrannhetsordningen är som förväntat 1 för Lax–Friedrichs och Upwind samt 2 för Upwind med ENO.

Som i första testet visar resultatet attnoggrannhetsordningen för metoderna hamnar på den förväntade. För stabilitet visade det sig att ∆t =

=

3√

2π

var litet nog för samtliga

metoder.

4.3 Problem 3, applicering av återinitilitisering

I det här problemet undersöker vi hur återinitilitisering kan förbättra resultatet av en lösning i de fall då stor förvrängning av avståndsfunktionen orsakar stora fel. Ett vanligt test problem är att applicera ett tidsberoende normalhastighetsfält som varierar med riktningen från origo [7].

Detta fält ges av:

V

= cos(8 arctan(y/x)) sin(2πt/T ). (52) Normalhastighetsfältet appliceras på en cirkel, med radien 3, i beräkningsdomänet [-6,6]x[-6,6].

Perioden T väljs till 6 och tidssteget sätts till ∆x/5. Den korrekta lösningen vid tiden t = T är samma cirkel som i begynnelsetillståndet.

Den numeriska metoden som används är Upwind då ENO implementeringen i detta fall resulterade i numerisk instabilitet. Det är oklart varför. Dock används TVD-RK2 metoden för att diskritisera i tiden. Detta ökar inte noggrannhetsordningen men återinitialisering kan enbart reducera felet från rumsdiskritiseringen och då felet från tidsdiskritisering minimeras blir detta tydligare.

Lösningen utan återinitialisering visas vid olika tider nedan i Figur 13.

T = 0

-6 -4 -2 0 2 4 6

x -6

-4 -2 0 2 4 6

y

T = 2

-6 -4 -2 0 2 4 6

x -6

-4 -2 0 2 4 6

y

T = 4

-6 -4 -2 0 2 4 6

x -6

-4 -2 0 2 4 6

y

T = 6

-6 -4 -2 0 2 4 6

x -6

-4 -2 0 2 4 6

y

Figur 13: Lösningen efter advektion med hastighetsfältet i ekv.(52) utan återinitialisering.

Lösningen har stort fel på grund av fel i beräkning av gradienten. Den korrekta lösningen visas

i blått. Här är N = 256 och ∆t = ∆x/5

Samma lösning körs, fast nu med 30 återinitialiseringssteg med jämna mellanrum och N/8 steg i varje återinitialisering. Resultatet visar en tydlig reduktion av felet i lösningen. Detta illustreras i Figur 14 här nedan.

-3 -2 -1 0 1 2 3

x -3

-2 -1 0 1 2 3

y

T = 6

Figur 14: Nollkonturen av level set funktionen vid tiden T = 6. Advektion med hastighetsfält i ekv.(52) och regelbunden återinitialisering. N = 256 och ∆t = ∆x/5.

Även om det är uppenbart att återinitialiseringen reducerar felet är det intressant att se hur

snabbt felet avtar när N ökar. Felet i position och area som funktion av N visa i Figur 15, som

återfinns på nästa sida.

10² 10³ N

10^-2 10^-1 10⁰

fel

Avståndsfel loglog plot

10² 10³

N 10^-2

10^-1

areafel

Areafel loglog plot

Upwind ref: 1/N

Figur 15: Loglog plot för areafelet och maximalt fel i radie då återinitialisering används. Arean

visar sämre konvergens än avståndsfelet. En möjlig förklaring kan vara att arean har både

positivt och negativt fel och att areafelet t.ex. råkade bli litet vid N = 128. Lutningen verkar

öka igen efteråt så det är möjligt att noggranhetsordningen skulle konvergera mot 1. På grund

av den extremt snabbt ökande beräknignskostnaden var det dock inte möjligt att undersöka

större N .

4.4 Problem 4, kombination av olika typer av hastighetsfält

Här används en levelsetfunktion som ger en stjärnformad nollkontur:

φ(x, y) = p

x

+ y

− (0.5 + 0.2 cos(5(y/x))). (53)

Funktionsytan är ej en nivåkurva, så återinitalisering körs innan advektionen påbörjas. Advektionsekvationen som löses är:

φ

= 0.02 κ

|κ|+1 |∇φ|, (54)

där κ är krökningen och konstanten  = 0.1. Lösningen körs till tiden 2.5 med ∆t = ∆x/4, och nätet har upplösningen 121x121. Resultatet från det krökningsdrivna fältet är att stjärnan kollapsar till en cirkel, som visas i Figur 16 här nedan.

T = 0

-1 -0.5 0 0.5 1

x

-0.5 0 0.5 1

y

T = 0.84

-1 -0.5 0 0.5 1

x

-0.5 0 0.5 1

y

T = 1.7

-1 -0.5 0 0.5 1

x

-0.5 0 0.5 1

y

T = 2.5

-1 -0.5 0 0.5 1

x

-0.5 0 0.5 1

y

Figur 16: Visar hur stjärnan förändras under tiden T när ekv.(54) löses genom tillämpning av

Upwind med ENO. Här är N = 121 och ∆t =

.