Preiswerte Energiegewinnung - Polynomiale … · Diskretisiere die Gleichung auf einem N-Punkt...

Preiswerte Energiegewinnung - PolynomialeEigenwertprobleme

Dominik Löchel

Betreuer: M. Hochbruck und M. Tokar

Graduiertenkolleg „Dynamik heißer Plasmen“

Mathematisches InstitutHeinrich-Heine-Universität Düsseldorf

Kompaktseminar Juli 2009

Gliederung

Physikalischer Hintergrund

Numerische Methoden

Energiegewinnung in der Zukunft

Energiegewinnung in der Zukunft

I Fossile Brennstoffe: gehen zu Neige, Treibhauseffekt

I Wasserkraftwerke: nicht weiter ausbaubar

I Solarkraftwerke: möglich, aber zu teuer

I Kernspaltungskraftwerk: radioaktiver Verseuchung, Uran begrenzt

I Fusionskraftwerk: preiswert, aber noch nicht realisierbar

Kernfusion – was ist das?

Kernfusion→

Kernspaltung←−

Fusion von Deuterium und Tritium zu Helium

E = mc2

Weshalb Kernfusion?Energiedichte

Energiegewinnung durch Kernfusion

I Fusion von Deuterium undTritium zu Helium

I Überwindung der CoulombBarriere→ starke Kernkraft

I hohe Temperatur undhohe Dichteausreichend lange Zeit

I PlasmaI magnetischer Einschluss





I Plasma

I magnetischer Einschluss





I PlasmaI magnetischer Einschluss

TEXTOR im Forschungszentrum Jülich

Tokamak Experiment for Technology Oriented Research

Bezeichnungen am Torus

ϕ ist der toroidale und θ der poloidale Winkel

Tokamak

I toroidale Magnetfeldspulen

I magnetische Feldlinien in der FlussoberflächeI PrimärspuleI Transformator-Eisenkern

Tokamak

I toroidale MagnetfeldspulenI magnetische Feldlinien in der Flussoberfläche

I PrimärspuleI Transformator-Eisenkern

Tokamak

~v~E×~B←−∇B

I toroidale MagnetfeldspulenI magnetische Feldlinien in der FlussoberflächeI PrimärspuleI Transformator-Eisenkern

Tokamak

I toroidale MagnetfeldspulenI magnetische Feldlinien in der FlussoberflächeI PrimärspuleI Transformator-Eisenkern

Drift-Instabilitäten


φ(θ, t) = φ(θ) · exp(ikr r + ik⊥y − iωt), =(ω) maximal


φ(θ, t) = φ(θ) ·exp(ikr r + ik⊥y− iωt)

0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen

I Störungswelle φ

I anomaler Transport Γ⊥


↓ φ(θ, t) = φ(θ) ·exp(ikr r + ik⊥y− iωt)

0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen↑I Störungswelle φ



→φ(θ, t) = φ(θ) ·exp(ikr r + ik⊥y− iωt)

0

0.5

1|φ|2

<(φ)

=(φ)




↑


0

0.5

1|φ|2

<(φ)

=(φ)




←φ(θ, t) = φ(θ) ·exp(ikr r + ik⊥y− iωt)

0

0.5

1|φ|2

<(φ)

=(φ)




↓ φ(θ, t) = φ(θ) ·exp(ikr r + ik⊥y− iωt)

0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen ↑I Störungswelle φ




0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen

I Störungswelle φ




0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen

I Störungswelle φI anomaler Transport Γ⊥


0

5

10

15

20

25

T /

eV

0

5

10

15

20n

/ 1019

m−

3

Tn

p = T (θ) n(θ)

innen θ außen


0

0.5

1|φ|2

<(φ)

=(φ)

innen θ außen

I Störungswelle φI anomaler Transport Γ⊥


0

5

10

15

20

25

T /

eV

0

5

10

15

20n

/ 1019

m−

3

Tn

p = T (θ) n(θ)

innen θ außen

Eigenwert Gleichung

∂2φ

∂θ2 =ω(β + zγ3µK 2

⊥)− (1 + λ)βK⊥ + izγ3CK 2⊥

γ31z

(K⊥ − ω

(1 + zγ3(1 + α)K 2

⊥))

·(

(1 + α)γ2

γ3LB(1− zγ3ωK⊥) + zω(ω + αK⊥)

)φ

generische Form:[ω3a3 + ω2a2 + ω

(a1 + b1

∂2

∂θ2

)+ a0 + b0

∂2

∂θ2

]φ = 0

aj ,bj , φ: [0,2π[ 7→C, 2π-periodische glatte Funktionen in θ

gesucht: Eigenpaar (ω, φ) mit maximaler Anwachsrate =(ω)

Diskretisierung[ω3a3 + ω2a2 + ω

(a1 + b1

∂2

∂θ2

)+ a0 + b0

∂2

∂θ2

]φ = 0

Diskretisiere die Gleichung auf einem N-Punkt Gitter:I [0,2π[ 7→ ~θ

I a(θ)→ Diag(a(~θ))

I ∂2

∂θ2 → D2 finite Differenzen oder Pseudo-Spektral-Methode

∂k f∂θk (θ0) ≈

m∑j=0

〈ψj , f 〉∂kψj

∂θk (θ0)

P(ω)~φ :=(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


(a1 + b1

∂2

∂θ2

)+ a0 + b0

∂2

∂θ2

]φ = 0

Diskretisiere die Gleichung auf einem N-Punkt Gitter:

I [0,2π[ 7→ ~θ


I ∂2


∂k f∂θk (θ0) ≈

m∑j=0


∂θk (θ0)

P(ω)~φ :=(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


(a1 + b1

∂2

∂θ2

)+ a0 + b0

∂2

∂θ2

]φ = 0

Diskretisiere die Gleichung auf einem N-Punkt Gitter:I [0,2π[ 7→ ~θ


I ∂2


∂k f∂θk (θ0) ≈

m∑j=0


∂θk (θ0)

P(ω)~φ :=(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0

Einschub: EigenwertproblemeStandard Eigenwertproblem

Ax = λx ⇔ (A− λI)x = 0

Standardlöser: QR-Algorithmus

Generalisiertes Eigenwertproblem

Ax = λBx ⇔ (A− λB)x = 0

Standardlöser: QZ -Algorithmus

Quadratisches Eigenwertproblem(λ2M + λC + K

)x = 0

Polynomiales Eigenwertproblemd∑

j=0

λjMjx = 0

Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0

Eine Linearisierung:

⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.

⇔ ωBx = −Ax (generalisiertes Eigenwertproblem)

I dreifache GrößeI 3N Eigenpaare (ω,~φ)

I Eigenvektoren ~φ linear abhängig

Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I dreifache Größe

I 3N Eigenpaare (ω,~φ)


Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I dreifache GrößeI 3N Eigenpaare (ω,~φ)


allgemeine Eigenschaften PEP

Polynomiales Eigenwertproblem:

P(λ) =d∑

j=0

λjMj~x = 0, M ∈ CN×N

Eigenschaften:

1. Matrizen A,B ∈ CdN×dN

2. P besitzt d · N Eigenpaare (λ,~x)

3. d ≥ 2: Eigenvektoren ~x linear abhängig.4. |λ| =∞ möglich, falls Kern(Md) 6= {~0}

Bsp.: (M2λ2 + M1λ+ M0)~x = 0, λ = α

β , α 6= 0

⇒ (M2

(αβ

)2+ M1

αβ + M0)~x = 0

⇒ (M2α2 + M1βα+ β2M0)~x = 0

⇒ (M2 + M1β + β2M0)~x = 0 (o.B.d.A. α = 1)



P(λ) =d∑

j=0


Eigenschaften:1. Matrizen A,B ∈ CdN×dN



Bsp.: (M2λ2 + M1λ+ M0)~x = 0, λ = α

β , α 6= 0

⇒ (M2

(αβ

)2+ M1

αβ + M0)~x = 0

⇒ (M2α2 + M1βα+ β2M0)~x = 0

⇒ (M2 + M1β + β2M0)~x = 0 (o.B.d.A. α = 1)



P(λ) =d∑

j=0




3. d ≥ 2: Eigenvektoren ~x linear abhängig.

4. |λ| =∞ möglich, falls Kern(Md) 6= {~0}

Bsp.: (M2λ2 + M1λ+ M0)~x = 0, λ = α

β , α 6= 0

⇒ (M2

(αβ

)2+ M1

αβ + M0)~x = 0

⇒ (M2α2 + M1βα+ β2M0)~x = 0

⇒ (M2 + M1β + β2M0)~x = 0 (o.B.d.A. α = 1)



P(λ) =d∑

j=0





Bsp.: (M2λ2 + M1λ+ M0)~x = 0, λ = α

β , α 6= 0

⇒ (M2

(αβ

)2+ M1

αβ + M0)~x = 0

⇒ (M2α2 + M1βα+ β2M0)~x = 0

⇒ (M2 + M1β + β2M0)~x = 0 (o.B.d.A. α = 1)

Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I N = 1024. QZ -Algorithmus: 1 Stunde.I Simulationen mit 100 000 Eigenwertgleichungen: 10 JahreI nur ein Eigenpaar gesucht→ iterativer Löser, speziell

Jacobi-Davidson-Verfahren

Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I N = 1024. QZ -Algorithmus: 1 Stunde.

I Simulationen mit 100 000 Eigenwertgleichungen: 10 JahreI nur ein Eigenpaar gesucht→ iterativer Löser, speziell


Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I N = 1024. QZ -Algorithmus: 1 Stunde.I Simulationen mit 100 000 Eigenwertgleichungen: 10 Jahre

I nur ein Eigenpaar gesucht→ iterativer Löser, speziellJacobi-Davidson-Verfahren

Lösung des PEP

(ω3M3 + ω2M2 + ωM1 + M0

)~φ = 0


⇔

ωM3

II

+

M2 M1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0.


I N = 1024. QZ -Algorithmus: 1 Stunde.I Simulationen mit 100 000 Eigenwertgleichungen: 10 JahreI nur ein Eigenpaar gesucht→ iterativer Löser, speziell


Jacobi-Davidson Verfahren

(0.) Wähle Suchraum V ∈ CN×k , k � N

loop(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.

(6.) Erweitere den Suchraum zu [V ,~t ].end loop


(0.) Wähle Suchraum V ∈ CN×k , k � Nloop

(1.) Orthonormalisiere V

(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.

(6.) Erweitere den Suchraum zu [V ,~t ].

end loop



(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0

(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.


end loop



(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).

(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.


end loop



(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .

if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.


end loop



(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if

(5.) Löse (näherungsweise) (z.B. mit GMRES)(I −

~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.


end loop



(1.) Orthonormalisiere V(2.) Berechne Eigenpaare (ν, ~y) von V HP(ν)V~y = 0(3.) Wähle Ritz Paar (ν, ~u := V~y).(4.) Berechne Residuum ~r := P(ν)~u. Es gilt ~r ⊥ V .if ‖~r‖ klein genug do Stopp end if(5.) Löse (näherungsweise) (z.B. mit GMRES)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u, ~w := P ′(ν)~u.

(6.) Erweitere den Suchraum zu [V ,~t ].end loop

Jacobi-Davidson Verfahren(5.) Berechne (näherungsweise)(

I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u

Ausgehend von (ν, ~u) einen Newton-Schritt

F (λ,~x) :=

(P(λ)~x~xH~x − 1

)= 0.

andere Auflösung: Ein-Schritt-Approximation [Sleijpen 1998]


I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u


F (λ,~x) :=

(P(λ)~x~xH~x − 1

)= 0.


~t =~uHQ(ν)~r

~uHQ(ν)P ′(ν)~uQ(ν)P ′(ν)~u −Q(ν)~r , Q(ν)P (ν) = I

Q(ν)~r = ~z ⇔ P (ν)~z = ~r , Q(ν)P ′(ν)~u = ~s ⇔ P (ν)~s = P ′(ν)~u

P (ν) = =


I −~w~uH

~uH ~w

)P(ν)(I − ~u~uH)~t = −~r , ~t ⊥ ~u


F (λ,~x) :=

(P(λ)~x~xH~x − 1

)= 0.


~t =~uHQ(ν)~r

~uHQ(ν)P ′(ν)~uQ(ν)P ′(ν)~u −Q(ν)~r , Q(ν)Pf (ν) = I

Q(ν)~r = ~z ⇔ Pf (ν)~z = ~r , Q(ν)P ′(ν)~u = ~s ⇔ Pf (ν)~s = P ′(ν)~u

Pf (ν) = =


(3.) Wähle Ritz Paar (ν, ~u := V~y).

Physik: Anwachsrate =(ω) maximal.

V = [~v1], ~v1 Zufallsvektor

I gesuchtes Eigenpaar wird i.A. nicht gefunden.



Möglichkeiten:I Zufallsvektor→Wahl des richtigen Ritzpaares fraglichI V =

(exp(ij~θ)

)j=−m,...,m→ solange ~φ entsprechend glatt

Idee: Approximation an ~φ verschaffen.Umsetzung: Rechnung auf gröberem GitterVorteil:

I gröbstes Gitter (z.B. N = 8) QZ , Eigenpaar auswählen.I hohe Genauigkeit durch Pseudo-Spektral-Methode oder optimierte

Gitterpositionen.I Eigenpaar auf nächst feinerem Gitter mit Jacobi-Davidson verbessern.I Wahl des Ritzpaares anhand der Ähnlichkeit zur Approximation.




(exp(ij~θ)


Idee: Approximation an ~φ verschaffen.

Umsetzung: Rechnung auf gröberem GitterVorteil:






(exp(ij~θ)


Idee: Approximation an ~φ verschaffen.Umsetzung: Rechnung auf gröberem Gitter

Vorteil:I gröbstes Gitter (z.B. N = 8) QZ , Eigenpaar auswählen.I hohe Genauigkeit durch Pseudo-Spektral-Methode oder optimierte





(exp(ij~θ)



I gröbstes Gitter (z.B. N = 8) QZ , Eigenpaar auswählen.

I hohe Genauigkeit durch Pseudo-Spektral-Methode oder optimierteGitterpositionen.

I Eigenpaar auf nächst feinerem Gitter mit Jacobi-Davidson verbessern.I Wahl des Ritzpaares anhand der Ähnlichkeit zur Approximation.




(exp(ij~θ)




Gitterpositionen.

I Eigenpaar auf nächst feinerem Gitter mit Jacobi-Davidson verbessern.I Wahl des Ritzpaares anhand der Ähnlichkeit zur Approximation.




(exp(ij~θ)




Gitterpositionen.I Eigenpaar auf nächst feinerem Gitter mit Jacobi-Davidson verbessern.

I Wahl des Ritzpaares anhand der Ähnlichkeit zur Approximation.




(exp(ij~θ)





Multilevel Jacobi-Davidson Verfahren

Beispiel:

0

1

|φ|2 N = 8

außen θ innen


Beispiel:

0

1

|φ|2 N = 16

außen θ innen


Beispiel:

0

1

|φ|2 N = 32

außen θ innen


Beispiel:

0

1

|φ|2 N = 64

außen θ innen


Beispiel:

0

1

|φ|2 N = 128

außen θ innen


Beispiel:

0

1

|φ|2 N = 256

außen θ innen


Beispiel:

0

1

|φ|2 N = 512

außen θ innen


Beispiel:

0

1

|φ|2 N = 1024

außen θ innen


Beispiel:

0

1

|φ|2 N = 2048

außen θ innen


Beispiel:

0

1

|φ|2 N = 4096

außen θ innen

Ergebnis

Beispiel 1

0 1 2 30

0.2

0.4

0.6

0.8

1

|φ|2

<(φ)

=(φ)

θ/π

Beispiel 2

0 1 2 30

0.2

0.4

0.6

0.8

1

|φ|2

<(φ)

=(φ)

θ/π

Ergebnis

Beispiel 1

0 1 2 30

0.2

0.4

0.6

0.8

1

|φ|2<(φ)=(φ)

θ/πω = −0.074 + 0.332i

0 1 2 30

0.2

0.4

0.6

0.8

1

|φ|2<(φ)=(φ)

θ/πω = −0.001 + 0.135i

Beispiel 2

8 16 32 128 256 512 1024 2048

1e−8

1e−6

1e−4

1e−2

1e0

|ωN − ω4096|

N

Aufwand der Rechnung

ω∗ exakter Eigenwert, ωN Näherung auf N-Punkt Gitter

Genauigkeitsforderung |ω∗ − ωN | ≤ 10−4|ω∗|

Beispiel 1

Beispiel 2

N dim t/sec

dim t/sec

4096 6 44.78

16 113.29

2048 6 11.17

16 28.35

1024 7 3.14

14 6.21

512 6 0.78

9 1.07

256 6 0.16

9 0.22

128 5 0.05

8 0.08

64 5 0.04

8 0.05

32 5 0.02

8 0.04

16 5 0.04

7 0.04

Summe 51 59

95 147




Beispiel 1 Beispiel 2N dim t/sec dim t/sec

4096 6 44.78 16 113.292048 6 11.17 16 28.351024 7 3.14 14 6.21

512 6 0.78 9 1.07256 6 0.16 9 0.22128 5 0.05 8 0.08

64 5 0.04 8 0.0532 5 0.02 8 0.0416 5 0.04 7 0.04

Summe 51 59 95 147

Skalierung

I Eigenwert skalieren

P(λ) = P(αω), α fest, ω skalierter Eigenwert

I löse SP(λ)~x = 0, S = Diag(~s), sj 6= 0 ∀jI optimale Skalierung zur Berechnung von (~φl , ω, ~φr ):

löse SlP(λ)Sr~y = 0, ~x = Sr~y , Sl = Diag(|~φl |), Sr = Diag(|~φr |)

I Sl , Sr durch Prolongation der Grobgitterlösung

I ~φl kann aus (ω,~φr ) bestimmt werden.

Skalierung



Auswirkung auf die LinearisierungλM3

II

+

M2 M1 M0−I

−I

λ2~φ

λ ~φ

~φ

= 0

ωα3M3

II

+

α2M2 αM1 M0−I

−I

ω2~φ

ω ~φ

~φ

= 0





Skalierung



I löse SP(λ)~x = 0, S = Diag(~s), sj 6= 0 ∀j

I optimale Skalierung zur Berechnung von (~φl , ω, ~φr ):




Skalierung



I löse SP(λ)~x = 0, S = Diag(~s), sj 6= 0 ∀jr := P(λ)x ⇒ Sr = SP(λ)x

I optimale Skalierung zur Berechnung von (~φl , ω, ~φr ):




Skalierung







Skalierung




löse SlP(λ)Sr~y = 0, ~x = Sr~y , Sl = Diag(|~φl |), Sr = Diag(|~φr |)I Sl , Sr durch Prolongation der Grobgitterlösung


Skalierung






∂2φ

∂θ2 =ω(β + zγ3µK 2

⊥)− (1 + λ)βK⊥ + izγ3CK 2⊥

γ31z

(K⊥ − ω

(1 + zγ3(1 + α)K 2

⊥))

·(

(1 + α)γ2


)φ

Skalierung






∂2φ

∂θ2 =Zähler(ω)

Nenner(ω)φ

Skalierung






Nenner(ω)∂2φ

∂θ2 = Zähler(ω)φ

Skalierung






Nenner(ω)∂2φ


~φl = Diag(Nenner(ω)

)−1~φr

Skalierung






Nenner(ω)∂2φ


~φl = Diag(Nenner(ω)

)−1~φr

V := S−1r~φr




Beispiel 1 Beispiel 2Skalierung ohne mit ohne mit

N dim t/sec dim t/sec dim t/sec dim t/sec4096 6 44.78

2 17.33

16 113.29

2 17.45

2048 6 11.17

5 9.15

16 28.35

5 9.79

1024 7 3.14

6 2.80

14 6.21

7 3.09

512 6 0.78

6 0.74

9 1.07

9 1.14

256 6 0.16

5 0.13

9 0.22

11 0.30

128 5 0.05

4 0.04

8 0.08

13 0.17

64 5 0.04

3 0.02

8 0.05

14 0.13

32 5 0.02

3 0.03

8 0.04

15 0.14

16 5 0.04

4 0.11

7 0.04

9 0.05

Summe 51 59

38 30

95 147

85 32





N dim t/sec dim t/sec dim t/sec dim t/sec4096 6 44.78 2 17.33 16 113.29

2 17.45

2048 6 11.17 5 9.15 16 28.35

5 9.79

1024 7 3.14 6 2.80 14 6.21

7 3.09

512 6 0.78 6 0.74 9 1.07

9 1.14

256 6 0.16 5 0.13 9 0.22

11 0.30

128 5 0.05 4 0.04 8 0.08

13 0.17

64 5 0.04 3 0.02 8 0.05

14 0.13

32 5 0.02 3 0.03 8 0.04

15 0.14

16 5 0.04 4 0.11 7 0.04

9 0.05

Summe 51 59 38 30 95 147

85 32





N dim t/sec dim t/sec dim t/sec dim t/sec4096 6 44.78 2 17.33 16 113.29 2 17.452048 6 11.17 5 9.15 16 28.35 5 9.791024 7 3.14 6 2.80 14 6.21 7 3.09

512 6 0.78 6 0.74 9 1.07 9 1.14256 6 0.16 5 0.13 9 0.22 11 0.30128 5 0.05 4 0.04 8 0.08 13 0.17

64 5 0.04 3 0.02 8 0.05 14 0.1332 5 0.02 3 0.03 8 0.04 15 0.1416 5 0.04 4 0.11 7 0.04 9 0.05

Summe 51 59 38 30 95 147 85 32

Wellenzahl

Wellenzahl K⊥

∂2φ

∂θ2 =ω(β + zγ3µK 2

⊥)− (1 + λ)βK⊥ + izγ3CK 2⊥

γ31z

(K⊥ − ω

(1 + zγ3(1 + α)K 2

⊥))

·(

(1 + α)γ2


)φ

Γ⊥ =

∫IC(T ,p,K⊥)

=(ω(K⊥))3

|ω(K⊥)|2|φ(K⊥)|2 dK⊥

I Diskrete Repräsentation des Wellenzahl-Intervalls.I Definiere Moden (ω(K⊥), φ(K⊥))

I Verfolgen der Moden und Auffinden der maximalen Anwachsrate

Wellenzahl

Wellenzahl K⊥

∂2φ

∂θ2 =ω(β + zγ3µK 2

⊥)− (1 + λ)βK⊥ + izγ3CK 2⊥

γ31z

(K⊥ − ω

(1 + zγ3(1 + α)K 2

⊥))

·(

(1 + α)γ2


)φ

Γ⊥ =

∫IC(T ,p,K⊥)

=(ω(K⊥))3

|ω(K⊥)|2|φ(K⊥)|2 dK⊥

I Diskrete Repräsentation des Wellenzahl-Intervalls.

I Definiere Moden (ω(K⊥), φ(K⊥))


Wellenzahl

Wellenzahl K⊥

∂2φ

∂θ2 =ω(β + zγ3µK 2

⊥)− (1 + λ)βK⊥ + izγ3CK 2⊥

γ31z

(K⊥ − ω

(1 + zγ3(1 + α)K 2

⊥))

·(

(1 + α)γ2


)φ

Γ⊥ =

∫IC(T ,p,K⊥)

=(ω(K⊥))3

|ω(K⊥)|2|φ(K⊥)|2 dK⊥

I Diskrete Repräsentation des Wellenzahl-Intervalls.I Definiere Moden (ω(K⊥), φ(K⊥))


WellenzahlK⊥ = 0.4, K⊥ = 0.5

−1

0

1

<(φ)

innen θ außen


−1

0

1

<(φ)

innen θ außensim(ω, ω) = exp(−|ω − ω|) sim(~φ, ~φ) = |~φH~φ|


−1

0

1

<(φ)


sim((ω,~φ), (ω,~φ)

):= sim(ω, ω) sim(~φ, ~φ)


−1

0

1

<(φ)


sim((ω,~φ), (ω,~φ)

):= sim(ω, ω) sim(~φ, ~φ)

sim(ω, ω)=

0.99 0.99 1.000.97 1.00 1.001.00 0.96 0.98

, sim(~φ, ~φ)=

0.47 0.87 0.990.28 0.99 0.871.00 0.34 0.45

Wellenzahl

−0.2 −0.1 0 0.1 0.2

0.02

0.04

0.06

0.08

Wellenzahl

−0.2 −0.1 0 0.1 0.2

0.02

0.04

0.06

0.08

0.2 0.4 0.6

0.02

0.04

0.06

0.08

K⊥

=(ω)

Wellenzahl

−0.2 −0.1 0 0.1 0.2

0.02

0.04

0.06

0.08

0.2 0.4 0.6

0.02

0.04

0.06

0.08

K⊥

=(ω)

0

0.2

0.4

0.6

0.8

1

|φ|2

innen θ außen

Neue MöglichkeitenEigenwertgleichung:

I früher: Mathieu-Gleichungmit 〈T 〉, 〈n〉

I jetzt: inhomogene Profile

0 0.5 10

0.1

0.2

=(ω)

K⊥

Γ⊥ = C(T ,p,K⊥)=(ω)3

|ω|2 |φ|2

0

5

10

15

20

25

T /

eV

innen θ außen

0

1

2

3

4

Γ ⊥ /

1023

m−

2 s−

1

innen θ außen

Selbstkonsistente Rechnung

I dynamische DämpfungI Trust region (T , p), exponentielle Terme (Ln, Ei)


Dämpfung in der FixpunktiterationI dynamische DämpfungI Trust region (T , p), exponentielle Terme (Ln, Ei)

Selbstkonsistente RechnungSimulation

〈n〉/1019 m−3 3 4 5 6 7

TemperaturprofilIterationen 220 282 36 97 549Iterationen mit Eigen-wertberechnung

19 20 8 17 38

Eigenwertgleichungen 1881 1980 792 1683 3762dimV N = 32 1.62 2.00 2.00 2.84 3.00

N = 64 1.25 1.81 2.00 2.00 3.00N = 128 1.00 1.04 1.44 2.00 3.00N = 256 1.00 1.00 1.00 1.00 3.00N = 512 1.00 1.00 1.00 1.00 3.00N = 1024 1.00 1.00 1.00 1.00 3.00

Dauer / Minuten 2:42 3:10 1:16 2:55 16:58


0

10

20

30

40

50

T /

eV

HFS θ LFS

45

810

20

304050

80

n / 1

019 m

−3

Neutralteilchenquelle bei— HFS and - - LFS— 〈n〉 = 4 · 1019 m−3

0

2

4

6

8

Γ ⊥ /

1021

m−

2 s−

1

HFS θ LFS


0

10

20

30

40

50

T /

eV

HFS θ LFS

45

810

20

304050

80

n / 1

019 m

−3


— 〈n〉 = 5 · 1019 m−3

0

2

4

6

8

Γ ⊥ /

1021

m−

2 s−

1

HFS θ LFS


0

10

20

30

40

50

T /

eV

HFS θ LFS

45

810

20

304050

80

n / 1

019 m

−3


— 〈n〉 = 5 · 1019 m−3

— 〈n〉 = 6 · 1019 m−3

0

2

4

6

8

Γ ⊥ /

1021

m−

2 s−

1

HFS θ LFS


0

10

20

30

40

50

T /

eV

HFS θ LFS

45

810

20

304050

80

n / 1

019 m

−3


— 〈n〉 = 5 · 1019 m−3

— 〈n〉 = 6 · 1019 m−3

— 〈n〉 = 7 · 1019 m−3

0

2

4

6

8

Γ ⊥ /

1021

m−

2 s−

1

HFS θ LFS

Simulation: magnetische GeometrieEinfluss der magnetischen Geometrie auf den anomalen Transport

0 50 100150

−200

−100

0

100

200

D = 0, E = 1,2,3

0 50 100150

−200

−100

0

100

200

D = 0.2, E = 1,2,30 50 100150

−200

−100

0

100

200

D = 0.4, E = 1,2,3


0 50 100150

−200

−100

0

100

200

D = 0, E = 1,2,30 50 100150

−200

−100

0

100

200

D = 0.2, E = 1,2,3

0 50 100150

−200

−100

0

100

200

D = 0.4, E = 1,2,3


0 50 100150

−200

−100

0

100

200

D = 0, E = 1,2,30 50 100150

−200

−100

0

100

200

D = 0.2, E = 1,2,30 50 100150

−200

−100

0

100

200

D = 0.4, E = 1,2,3


0 50 100150

−200

−100

0

100

200

HFS LFS0

5

10

15

E=3E=2E=1

D = 0|φ|2

0 50 100150

−200

−100

0

100

200

HFS LFS0

5

10

E=3E=2E=1

D = 0.2|φ|2

0 50 100150

−200

−100

0

100

200

HFS LFS0

2

4

6

E=3E=2E=1

D = 0.4|φ|2

MAST parameters

Simulation: magnetische Geometrie

1 2 3

0.6

0.7

0.8

ℑ(ω

)

E

1 2 3

−0.23

−0.22

−0.21

−0.2

ℜ(ω

)

E

1 2 30.42

0.43

0.44

0.45

0.46

E

K⊥

1 2 30

5

10

x 1019

E

Γm

ax

D = 0

D = 0.2

D = 0.4

X-Punkt Geometrie

Bisher: Limiter Maschine

Besser: Divertor Maschine

X-Punkt Geometrie

Bisher: Limiter Maschine Besser: Divertor Maschine

Limiter und Divertor Maschine

Randbedingung

φ(θ) = φ(θ + 2πn), n ∈ Z

φ(θ) = cnφ(θ + 2πn), c = exp(−2πikq), n ∈ Z, k ∈ N

φ(θ) = φ(θ + 2πn), n ∈ Z

φ(θ) = cnφ(θ + 2πn), c = exp(−2πikq), n ∈ Z, k ∈ N

Irrationale q:

Torus Koordinaten

0 0.5 1 1.5 2−1

−0.5

0

0.5

1

|φ|2

<(φ)

=(φ)

θ/π

Magnetfeld Koordinaten

0 0.5 1 1.5 2−1

−0.5

0

0.5

1|φ|2

<(φ)

=(φ)

θ/π

Zusammenfassung

I multilevel Jacobi-Davidson Eigenwertlöser• Approximation des Suchraumes auf grobem Gitter• Korrekturgleichung über preiswerte LR-Zerlegung• Verbesserung der Kondition durch Skalierung SlP(λ)Sr

I Verfolgen der Eigenmoden (Wellenzahl)

I selbstkonsistente Rechnung: Fixpunktiteration mitproblemangepasster Dämpfung

I physikalisches Modell bedarf der Erweiterung

I Auswirkung der magnetischen Geometrie auf den Verlust

I Ausblick: X-Punkt-Geometrie in der divertor Maschine

Preiswerte Energiegewinnung - PolynomialeEigenwertprobleme

Energiekosten

2009 Zeit

P(ω)~φ = 0

Preiswerte Energiegewinnung - Polynomiale … · Diskretisiere die Gleichung auf einem N-Punkt...

Documents

Transcript of Preiswerte Energiegewinnung - Polynomiale … · Diskretisiere die Gleichung auf einem N-Punkt...