305
pages

Voir plus
Voir moins

Vous aimerez aussi

Fakult¨at fur¨ Mathematik und Physik

Solving the System of Radiation

Magnetohydrodynamics

for solar physical simulations in 3d

Andreas Dedner

Dissertation zur Erlangung des Doktorgrades der Fakult¨at fur¨ Mathematik und

Physik der Albert-Ludwigs-Universitat¨ Freiburg im Breisgau

Betreuer: Prof. Dr. Dietmar Kroner¨

Abteilung fur¨ Angewandte Mathematik

Freiburg im Breisgau, April 2003Dekan : Prof. Dr. Rolf Schneider

Referenten : Prof. Dr. Dietmar Kr¨oner

: Prof. Dr. Gerald Warnecke, Universit¨at Magdeburg

Datum der Promotion : 22. September 2003

Picture on title page:

simulation of a circular slip stream. In the purely hydrodynamic setting the interface is unstable with

respect to perturbations (Kelvin–Helmholtz instability). By means of a suﬃciently strong magnetic

ﬁeld that is tangential to the interface, this instability can be suppressed. For the simulation shown

here, the magnetic ﬁeld is not yet strong enough so that the development of the Kelvin–Helmholtz

instabilities can be clearly seen. The large picture shows the density of the ﬂuid. The two smaller

pictures show the locally adapted grid (top) and the third component of the magnetic ﬁeld (bottom)

in a small section of the domain in the vicinity of the interface.Abstract

In this study we present a ﬁnite–volume scheme for solving the equa-

tions of radiation magnetohydrodynamics in two and three space di-

mensions. Among other applications this system is used to model the

plasmainthesolarconvectionzoneandinthesolarphotosphere. Itis

a non–linear system of balance laws derived from the Euler equations

of gas dynamics and the Maxwell equations; the energy transport

through radiation is also included in the model. The starting point of

our presentation is a standard explicit ﬁrst and second order ﬁnite–

volume scheme on both structured and unstructured grids. We ﬁrst

study the convergence of a ﬁnite–volume scheme applied to a scalar

model problem for the full system of radiation magnetohydrodynam-

ics. We then present modiﬁcations of the base scheme. These make

itpossibletoapproximatethesystemofmagnetohydrodynamicswith

anarbitraryequationofstate; theyreduceerrorsduetoaviolationof

the divergence constraint on the magnetic ﬁeld, and they lead to an

improved accuracy in the approximation of solution near an equilib-

rium state. These modiﬁcations signiﬁcantly increase the robustness

oftheschemeandareessentialforanaccuratesimulationofprocesses

in the solar atmosphere. For simulations in the solar photosphere, we

have to take the radiation intensity into account. A scheme for solv-

ing the radiation transport equation is a further focus of this study.

We present both analytical results and numerical tests, comparing

our scheme with some standard schemes found in the literature. We

conclude our presentation with a study of the parallelization strategy

for distributed memory computers that we use in our 3d code.iiIntroduction

Numerical simulations have become an important tool for studying many diﬀerent

physical and technical problems. Ranging from the formation of galaxies to weather

forecasts to the design for parts of complex machinery, the applications are numerous.

On the one hand, numerical simulations serve as a tool for the veriﬁcation of physical

theories deduced from observation; on the other hand, they play an important role

in reducing development cost in manufacturing. Although the range of applications

is extremely broad, the methods used for solving problems numerically have many

features in common. This is due to the fact that the physical models used have similar

properties. For example, ﬂuid ﬂow in the atmosphere of stars or in car engines can be

modeled by very similar systems of equations and can be simulated using very similar

numerical methods.

In this study we investigate numerical schemes that can be used to simulate the evo-

lution of a compressible ﬂuid. The governing system of partial diﬀerential equations

is based on the Euler equations of gas dynamics. Over the last centuries, this system

has been the focus of both analytical and numerical studies. A large number of dif-

ferent schemes have been developed and tested using this system. One very successful

approach turned out to be the ﬁnite–volume framework, and many diﬀerent schemes

basedonthisapproachhavebeenpresented. Thesamemethodshavealsobeenapplied

to diﬀerent extensions of the basic system of Euler equations, including, for example,

reactive ﬂow and magnetohydrodynamics. The latter will be the main focus of our

study.

Solar physical applications

The material presented in the following is part of a project ﬁnanced by the Deutsche

Forschungsgemeinschaft(DFG)aimedatderivingandanalyticallyjustifyingnumerical

methods for studying ﬂuid ﬂow in the solar atmosphere. The development of many of

the methods is a direct consequence of the interaction between members of our group

here in Freiburg (Dietmar Kroner,¨ Christian Rohde, Matthias Wesenberg, and myself)

and solar physicists (Manfred Schuss¨ ler and Peter Vollm¨oller from the Max–Planck

Institute for Aeronomie in Kattlenburg–Lindau), whose ideas greatly inﬂuenced our

work. Many of the problems discussed here occur only if the methods are applied

not to academic test cases, but to realistic settings. Therefore, the discussions with

the solar physicists and their help in developing and testing the numerical methods

inﬂuenced the direction in which our work progressed.

Although a variation in solar activity has a strong impact on life here on earth, a

thorough understanding of the physical processes behind these phenomena is still the

subject of research all over the world:

iiiiv INTRODUCTION

http://image.gsfc.nasa.gov/poetry/storm0/black1.html http://sec.noaa.gov/SWN/index.html

Storms are usually responsible for the losses of Early records of sunspots indicate that the Sun

electricity we endure, but did you know that went through a period of inactivity in the late

”storms” as far away as the sun are capable 17thcentury. Veryfewsunspotswereseenonthe

of knocking out large areas of electric service? Sun from about 1645 to 1715. [...] This period

Amazingly, the sun is capable of not only dis- of solar inactivity also corresponds to a climatic

rupting electrical power, but also short wave period called the ”Little Ice Age” when rivers

radio, television and telegraph signals, naviga- that are normally ice-free froze and snow ﬁelds

tional equipment (GPS and LORAN), defense remained year-round at lower altitudes. There

(military) early warning radar systems, the cli- is evidence that the Sun has had similar peri-

mate,andcanevenknockoutourcommunication ods of inactivity in the more distant past. The

satellites in space. connection between solar activity and terrestrial

http://image.gsfc.nasa.gov/poetry/storm0/black1.html climate is an area of on-going research.

http://sec.noaa.gov/SWN/index.html

Since the possibilities for direct observation of physical processes below the solar sur-

face are limited, numerical simulations play an important role in obtaining a clearer

understanding of solar phenomena. A further example of a solar phenomena not yet

fully understood is the eleven year cycle in which the number of sun spots on the solar

surfaceincreaseanddecrease. Onediﬃcultyisthattheﬁlamentsattheboundaryofthe

5sun spots are made up of magnetic ﬂuxtubes that are formed about 2·10 kilometers

below the solar surface in the lower convection zone of the sun. In this region direct

observation is hardly possible so that the formation and evolution of the ﬂuxtubes has

to be studied by means of numerical simulations. Although the presentation here is far

more general and such solar phenomena are not the immediate focus, the application

of our method to problems in solar physics has been a constant motivation.

Mathematical model

Themathematicalmodelconsistsofasystemofbalancelawscombiningtheequationsof

magnetohydrodynamics (MHD) and the radiation transport (RT) equation. The MHD

equations are a non–linear system of eight conservation laws; the energy transport

through radiation leads to an additional source term that is non–local in space.

The MHD equations describe the evolution of an electrically conductive plasma in the

presence of magnetic ﬁelds and combine the Euler equations of gas dynamics and theINTRODUCTION v

Maxwellequations. Thelatteralsointroduceaconstraintequationonthedivergenceof

the magnetic ﬁeld. In the solar atmosphere the force of gravity plays an important role

and is included in our model via source terms. To perform the simulations, we have to

prescribe suitable initial conditions for the ﬂuid. These often consist of a perturbation

of a stratiﬁed and static background atmosphere. One intrinsic problem of simulations

with this type of initial data is the size of the computational domain. The setting

allows for no physical boundaries, and the construction of suitable artiﬁcial boundary

conditions that can be used in numerical simulations is no easy undertaking.

The main diﬃculty in approximating the radiation ﬁeld is the high dimensionality

of the problem and the propagation speed of the radiation, which is several orders

of magnitude above the speed of the ﬂuid. In our model we deal with the second

problem by assuming an instantaneous radiation equilibrium. We have thus removed

the diﬀerent time scales, but we introduce a non–local dependency into our problem,

which we have to cope with in our numerical scheme. The high dimensionality of the

radiation intensity — it depends on space, time, propagation angle, and frequency —

forces us to construct a very eﬃcient solver to compute the radiation ﬁeld.

Numerical scheme

We use a ﬁrst and second order ﬁnite–volume scheme on locally adapted structured

and unstructured grids. We have implemented this method in one, two, and three

space dimensions, using both Cartesian and triangular grids in 2d and hexahedral

and tetrahedral grids in 3d. To increase the eﬃciency of the scheme, we make use

of parallelization strategies including distributed memory parallelization with dynamic

loadbalancing. Mostofthemethodspresentedinthisstudyare,however,notrestricted

tousewithaﬁnite–volumeschemeandhavebeenconstructedtobeeasilyaddedtoany

existing method for solving the system of magnetohydrodynamics. We have already

pointed out that the application to solar physical problems serves as a motivation for

the development and test of the scheme; the presentation, however, is kept at a far

moregenerallevel. Forexample, thecorrectionmethodusedtocomputesolutionsnear

anequilibriumstatecanbeusedformanytypesofatmosphericﬂow,orevenfortotally

diﬀerent applications where the problem of balancing source terms and ﬂux gradients

plays a crucial role.

An important consideration for the development of our methods is their simple im-

plementation within the framework of an existing numerical scheme, which we modify

as little as possible. The complexity of our applications also requires an eﬃcient so-

lution algorithm. Consequently, none of the modiﬁcations should lead to an increase

in the computational cost. Furthermore, we try to reduce the number of free param-

eters as much as possible; if available, we use an analytically motivated choice for the

parameters, otherwise we try to ﬁnd suitable values by means of numerical tests.

Analytical justiﬁcation

There are very few analytical results for complex non–linear coupled systems of the

type studied here. Even for the MHD system without radiation very little is knownvi INTRODUCTION

concerning the existence and uniqueness of solutions for general initial data. Conse-

quently, convergence analysis for numerical schemes is not yet available. One approach

often used in the analysis of complex systems is to reduce the complexity (often down

to a scalar balance law), taking care to retain the important characteristic features of

the original system. A multitude of analytical results are available for scalar balance

laws, ranging from existence and uniqueness results to the convergence of numerical

schemes in higher space dimensions. We employ this approach to justify the use of a

ﬁnite–volume scheme to solve the coupled system of radiation magnetohydrodynam-

ics. We carefully derive a scalar balance law that includes a non–local operator. This

source term, which models the radiation transport, is the novel feature of our model

problem. We ﬁrst study the inﬂuence of this non–local term on the solution of the

model problem; then we prove the convergence of a ﬁnite–volume scheme including an

explicit approximation of the non–local operator.

Numerical tests

Themathematicalmodelconsistsoftwoparts, onedescribingtheevolutionoftheﬂuid

andtheothertheradiationﬁeld. Theconstructionofournumericalschemesisbasedon

this splitting, which leadstoa MHDandaRT“module”. Thesemodulesare discussed

and tested separately since the coupling of the radiation and the ﬂuid ﬂow occurs only

on a source term level. We thereby assume that the performance of the full scheme

can be measured by the performance of both contributing modules. This indirect test

of our algorithm is necessary since we are not aware of any simple test cases for the

full coupled system. A rigorous test of the full algorithm is very diﬃcult and requires

a detailed understanding of the underlying physical processes in the solar atmosphere;

this is beyond the scope of this presentation.

The main focus of our study is a comparison of the eﬃciency of diﬀerent numerical

schemes. We compare often used approaches from the literature with newly developed

schemes. We measure the eﬃciency of a scheme by studying the error to runtime ratio.

Since the complexity of our problems (especially of our simulations in 3d) leads to a

high demand on computational cost, the runtime eﬃciency of the numerical scheme

has to be the essential aspect of our study.

Hardware and software used

The numerical scheme is implemented in C++. We used many diﬀerent computer

systems for our numerical tests, including single processor Linux PC, a shared memory

SGI computer system (Origin with 46 processors), the IBM RS/6000 SP computer

at the Rechenzentrum in Karlsruhe, and the IBM Regatta at the Rechenzentrum in

Freiburg. BothGnuPlotandthegraphicslibraryGraPEwereusedforthevisualization

of the data. Detailed references to all software packages used are given later.INTRODUCTION vii

Outline of the thesis

In the ﬁrst chapter we derive the relevant system of equations for our solar physical

applications. The physical derivation of the full system is not discussed in detail,

only the relevant notation is introduced. The main part of our study is divided into

three parts, each of which is preceded by an overview chapter and concluded with a

summary. In the ﬁrst part we outline a very general numerical scheme for solving the

system of radiation magnetohydrodynamics. We justify the numerical scheme through

the analysis of a simpliﬁed setting. In the second and third parts, we extend this basic

numericalscheme, treatingthepartsfortheﬂuidandtheradiationseparately. Wenow

describe the three parts in more detail.

In the ﬁrst part (Chapters 2–5) we present a standard ﬁnite–volume scheme for

solving the system of radiation magnetohydrodynamics. At this stage we describe only

the standard building block as can be found in the literature. The scheme described

in Chapter 3 does not yet include all aspects of our mathematical model and in its

basic form it is not suitable for use in challenging applications. It serves rather as the

skeleton for the extensions described in the second and third parts. Before we study

the necessary modiﬁcations of the scheme, we justify the general approach with an

analyticalstudyofasimpliﬁedmodelproblem. Thefactthatthecentralnon–standard

aspect of the system is the non–local eﬀect of the radiation source term dictates the

choice of material presented in Chapter 4. The model problem consists of a scalar

balance law with a right hand side including a non–local integral operator. We

ﬁrst study the properties of special solutions to the model scalar balance law and then

present a general convergence proof for ﬁnite–volume schemes in 2d.

In the second part (Chapters 6–10) we present modiﬁcations of that part of the numer-

ical scheme in which the evolution of the ﬂuid variables is computed. In the overview

Chapter6wepresentanumberofchallengesthatournumericalschememustmeetand

also describe approaches found in the literature, approaches we then use as comparison

schemes for our own solution technique. The comparison methods are chosen in accor-

dancewiththeguidelineswesetupforourownmodiﬁcationasdiscussedabove(simple

extension of existing scheme and no additional computational cost). Chapters 7–9 are

devoted to the description of the methods and numerical tests for three central chal-

lenges: we ﬁrst study a relaxation approach that allows us to extend a solver for a

perfect gas to approximate the MHD equations with a general equation of state;

thenwepresentageneralframeworkinwhichthedivergenceconstraintonthemag-

netic ﬁeld is coupled with the evolution equation for the magnetic ﬁeld. Based on this

approach we derive a number of diﬀerent correction mechanisms for reducing errors

in the divergence of the magnetic ﬁeld. Finally we study a modiﬁcation of the base

scheme that facilitates the accurate approximation of solutions near an equilibrium

state.

Thethirdpartofourinvestigation(Chapters11–14)isdevotedtothepresentationand

study of numerical schemes for solving the radiation transport equation. Again

we start with an overview in which we discuss the central aspects of this part of the

numerical scheme and present a standard solver found in the literature that we use as

a reference method. In Chapter 12 we derive a numerical scheme for approximating

the radiation intensity for a ﬁxed propagation direction. This is the central buildingviii INTRODUCTION

block used for the approximation of the radiation source term that enters into the

balance law for the total energy. We present a convergence proof for our method

and, after presenting numerical tests for ﬁxed propagation directions, we conclude

our investigation of the radiation transport module in Chapter 13 by studying the

approximation of the radiation source term itself.

In the last Chapter we then present some results using our3d MHD code, including

a simulation for a problem from solar physics.

The enclosed CD ROM contains a pdf version of this thesis and the sources of our

2d and 3d MHD code. Furthermore, we have included the web pages of our project.

TheCDROMalsocontainsadditionalmaterialincludingmoviesandpostersthatwere

produced during the project. The ﬁle (MHD.html) in the root directory of the CD

ROM gives speciﬁc details of the layout.

Acknowledgments

I would like to express my thanks to all my colleagues here at the IAM in Freiburg

— especially to Matthias Wesenberg and Christian Rohde for the discussions of the

methods,oftheresults,andoflife,theuniverse,andeverything. Thereareinfactmany

reasons for thanking Matthias Wesenberg with whom I have been working together

since the beginning of our studies. All my colleagues helped me by supplying such

important raw material as coﬀee and tea; special thanks, however, should go to the

diﬀerent people who, over the last few years, invested a huge amount of their time to

keeping our computer system in working order.

Let me continue by saying thank you to my supervisor Dietmar Kr¨oner for suggesting

andplanningthisprojectandforthemanyfruitfuldiscussionsandsuggestions. Iwould

also like to thank him for introducing me to several other scientists who also inﬂuenced

theworkpresentedhere, butofwhomIcanmentiononlyafewhere. Ontheonehand,

these include our project partners from the Max–Planck Institute of Aeronomie, Man-

fred Schu¨ssler and Peter Vollm¨oller. Many of the problems studied here were brought

toourattentiononlythroughtheircontinuingdesiretoapplyourmethodstoproblems

from solar physics; the development of the new radiation transport solver was one of

the results of this cooperation. On the other hand, I want to mention Ivan Sofronov

from the Keldysh Institute of Applied Mathematics RAS in Moscow, with whom we

developed transparent boundary conditions for our atmospheric ﬂow problems, and

Claus–Dieter Munz from the Institut fur¨ Aerodynamik und Gasdynamik, Universit¨at

Stuttgart, who suggested extending his divergence cleaning technique derived for the

Maxwell equations to the MHD system. The analytical results for the ﬁnite–volume

scheme were proven in cooperation with Christian Rohde, large parts of the numerical

scheme were implemented together with Matthias Wesenberg, and the grid concept

goes back to Bernhard Schupp, whom I would also like to thank for saving me from

the tedious task of having to design the hierarchical grid concept.

My thanks for proofreading this thesis go to Christian Rohde and to my mother. To

my wife Sabine Voigt go my very special thanks for reasons of which she is well aware.

This study was supported by the Deutsche Forschungsgemeinschaft (DFG) as part of

thepriorityresearchprogramAnalysisandNumericsforConservationLaws (ANumE).