Radiosity equation

$
\def\x{{\bf x}}
\def\y{{\bf y}}
\def\out{{\rm out}}
\newcommand{\Spro}[2]{\langle {#1},{#2} \rangle}
$

Die Radiosity-Gleichung beschreibt den Strahlungsaustausch zwischen diffusen Oberflächen und ist die Grundlage für Surface-to-Surface-Strahlungsmodelle. Bei großen Modellen wird die numerische Behandlung allerdings schnell anspruchsvoll, da die entstehenden Gleichungssysteme einen hohen Speicher- und Rechenaufwand verursachen. Ein zentraler Grund dafür ist die dicht besetzte View-Factor-Matrix. Um solche Probleme effizient behandeln zu können, sind daher Verfahren erforderlich, die nicht auf ihrer vollständigen Berechnung und Speicherung beruhen. Auf dieser Seite werden die zugrunde liegende Gleichung und ein numerischer Ansatz vorgestellt, der auf hierarchischen Datenstrukturen und iterativen Lösungsverfahren basiert.

Mathematisches Modell

Die radiosity equation ist eine Integralgleichung, die den Strahlungs- bzw. Energieaustausch zwischen diffusen grauen Oberflächen modelliert. Sie lautet
\begin{equation} \label{IGLqout2}
q_\out (\x) = \epsilon \, \sigma \, T^4(\x) + \rho \, \int_\Gamma k(\x,\y) \, q_\out (\y) \, d\y,
\end{equation}
für einen Punkt $\x$ auf der Oberfläche $\Gamma$. Die radiosity equation ist eine vereinfachte Version der rendering equation für diffuse Oberflächen.

Die Gleichung \eqref{IGLqout2} repräsentiert keinen Zusammenhang zwischen diskreten Größen, sondern wird im Rahmen eines Funktionenraums betrachtet. Das bedeutet, dass $q_\out$, $T$ und der Kern $k$ Funktionen sind, die bestimmte (Regularitäts-) Eigenschaften haben. Dies führt zur Fredholm-Theorie, in der die Gleichung üblicherweise mit einem Integraloperator $\mathcal{K}$ dargestellt wird:
$$
( I – {\mathcal K} ) \, q_{\rm out} (\x) = \epsilon \, \sigma \, T^4(\x) , \quad \forall \x \in \Gamma .
$$

Ausgehend von den Temperaturen auf den beteiligten Oberflächen wird die entsprechende Helligkeit, d.h. die Ausstrahlung oder radiance, durch $\eqref{IGLqout2}$ festgelegt. Eine besondere Rolle spielt hierbei die Sichtbarkeitsfunktion. Obwohl sie auf den ersten Blick harmlos erscheint, entwickelt sich ihre Berechnung zu einem ernsthaften Problem hinsichtlich des erforderlichen Rechenaufwands, siehe Visibility.

Bezeichnungen

Bezeichnungen (klicke hier)

$\Gamma$ eine Oberfläche in $\mathbb{R}^3$
$T(\x)$ die Temperatur im Oberflächenpunkt $\x$,
$q_{\rm out}$ die Helligkeit eines Paneels
$\epsilon$ der Emissionsgrad, $\rho=1-\epsilon$ der Reflexionsgrad und $\sigma$ die Stefan-Boltzmann Konstante,
$ $
\begin{equation} \label{kern3D}
k(\x,\y) := \frac{\cos \phi_\x \, \cos \phi_\y }{\pi \; \|\x-\y\|^2} \, \chi(\x,\y), \qquad \x,\y \in \Gamma,
\end{equation}
der Kern des Integraloperators $\mathcal{K}$
\begin{equation} \label{Kdef}
({\mathcal K} \, q )(\x) := \rho \, \int_\Gamma k(\x,\y) \, q (\y) \; d\y, \quad \x \in \Gamma,
\end{equation}
$ \chi(\x,\y)$ die visibility-Funktion, siehe hier.

Numerische Lösung

Für die numerische Lösung auf einem Computer muss die kontinuierliche Gleichung \eqref{IGLqout2} in ein diskretes System überführt werden. Dies kann durch verschiedene Methoden erfolgen, darunter das Kollokations- und das Galerkin-Verfahren. Die Diskretisierung der radiosity equation mit dem Kollokationsverfahren führt auf ein lineares Gleichungssystem der Form
\begin{equation} \label{DCR}
(I_n – K_n) \, q_n = E_n
\end{equation}
mit

$n$ der Anzahl der Paneele (geometrische Diskretisierung der Oberfläche),
der vollbesetzten(!) Matrix $K_n \in \mathbb{R}^{n \times n}$, die Viewfactor-Matrix (Cohen^[1]),
dem $n$-dimensionalen Vektor $q_n \in \mathbb{R}^n$, die gesuchten Wärmestromdichten (Helligkeit) bzgl. der Paneele
und der bekannten rechte Seite (die Quellen) $E_n \in \mathbb{R}^n$ .

Für feine Diskretisierungen, also für großes $n$, steigt der Aufwand stark an. Sowohl die Rechenzeit zur Bestimmung als auch der Speicherbedarf für die $n^2$ Koeffizienten der Matrix $K_n$ sind dann nicht mehr akzeptabel. So sind bei einfacher Floating-Point-Genauigkeit $4 \, n^2$ Byte zur Speicherung erforderlich. Die folgende Tabelle enthält einige Werte zum Speicherplatzbedarf der Viewfactor-Matrix:

$n$	1000	100000	250000	1000000	5000000
$4\, n^2$ byte	3.81 [MB]	38.1 [GB]	238.4 [GB]	3.8 [TB]	95.3 [TB]

Ist ein Modell mit 5 Millionen Oberflächenelementen aufgelöst, so benötigt die Berechnung der vollen Viewfactor-Matrix eine Speicherplatzbelegung von 95,3 TB. Wenn man davon ausgeht, dass die Systemmatrix auf einem Cluster verteilt ist, entspricht dies etwa 186 Workstations mit jeweils 512 GB RAM-Speicher.

Iteratives Verfahren

Ein iteratives Lösungsverfahren für das Gleichungssystem (4) erfordert nicht zwingend die explizite Kenntnis der Matrix $K_n$. Stattdessen genügt es, wenn ein Algorithmus für die Matrix-Vektor-Multiplikation
\begin{equation} \label{MatVek}
(K_n \, q_n)_i = \sum_{j=1}^n q_j \int_{t_j} k(\xi_i, \y) \, d\y, \quad i=1,\ldots,n ,
\end{equation}
vorhanden ist. Die Idee besteht nun darin, nicht mehr die Matrix $K_n$ zu berechnen, sondern eine sinnvolle Approximation der Matrix-Vektor-Multiplikation $K_n q_n$ zu finden. Dies wird durch eine spezielle Nahfeld-Fernfeld-Zerlegung, das sogenannte panel clustering, ermöglicht.

panel clustering

Die Entwicklung eines funktionierenden Verfahrens zur numerischen Lösung der Gleichung \eqref{IGLqout2} mittels panel clustering^[2] ist anspruchsvoll und erfordert einen hohen Arbeitsaufwand. Insbesondere sind die folgenden Punkte zu beachten:

Als notwendige Voraussetzung zum Panel-Clustering ist aus der geometrischen Diskretisierung ein sogenannter Clusterbaum zu bestimmen.
Es bedarf einer Datenstruktur für den “Panel-Clustering Operator”, um die (approximierte) Matrix-Vektor-Multiplikation effizient zu berechnen.
Das Auftreten von schwach singulären Integralen im Nahfeld muss berücksichtigt werden.
Es ist erforderlich, die Sichtbarkeiten zu bestimmen und entsprechende Datenstrukturen zu implementieren,
und schließlich braucht man noch ein passendes iteratives Verfahren, um die diskrete Gleichung \eqref{DCR} zu lösen.

Das Panel-Clustering-Verfahren bietet deutliche Vorteile gegenüber klassischen Methoden, insbesondere im Vergleich zur Verwendung der vollen Viewfactor-Matrix:

Durch das Panel-Clustering wird der Speicherplatzbedarf erheblich reduziert, was die Durchführung der numerischen Lösung auf einer Workstation oder wenigen Clusterknoten erst ermöglicht.
Schnelle Matrix-Vektor-Multiplikation und damit auch eine schnellere Lösung des diskreten GLeichungssystems.

Bei der Wahl des iterativen Lösers kann man eine positive Eigenschaft der radiosity equation nutzen, deren Darstellung hier den Rahmen sprengen würde. Eine Variante des BiCG-Verfahrens^[3] hat sich als äußerst erfolgreich erwiesen. Dabei zeigt sich ein sehr interessantes Phänomen: Die Konvergenzrate verschlechtert sich nicht mit einer geometrischen Verfeinerung, also mit einer höheren Dimension der Diskretisierungsmatrix. Sie ist unabhängig von der Dimension des Gleichungssystems.

Das Verfahren der GWR GmbH ist an realen Problemstellungen erprobt. Es ermöglicht daher eine effiziente Berechnung sowohl auf einer Workstation als auch auf einem HPC-System.

Beispiel

Der nachfolgende Plot zeigt den Verlauf der Iterationen für ein Testproblem mit 155640 Freiheitsgraden. Bereits nach 7 Iterationen wurde eine hinreichend genaue Lösung ermittelt, für die die eingesetzte Workstation (Ryzen 9 5950X-Prozessor mit 16 Kernen) lediglich 1.4 Sekunden benötigte.
Im Vergleich zum klassischen Verfahren, das etwa 90 GB Speicher für die Berechnung der vollen Viewfactor-Matrix erfordern würde, kam das gesamte Programm nur etwa 4,5 GB (einschließlich der Wärmeleitungsberechnung) aus. Davon entfielen nur ca. 1 GB auf das Panel-Clustering-Verfahren.

Konvergenz des Strahlungslösers — Konvergenz: Reduktion des Residuum; man beachte den logarithmischen Maßstab für die y-Achse

Bestrahlungsstärke (Irradiance)

Die im Punkt $\x$ eingestrahlte Energiestromdichte $q_{\rm in}$ ist die Summe aller $q_{\rm out} (\y)$, die auf $\x$ “fallen”. Mit einer geometrischen Betrachtung über den Raumwinkel findet sich ein analytischer Ausdruck mit dem die eingestrahlte Energie (-stromdichte) in Abhängigkeit von $q_{\rm out}$ angegeben werden kann (Siegel^[4]):
\begin{equation} \label{IGL1}
q_{\rm in} (\x) = \int_\Gamma k(\x,\y) \, q_{\rm out}(\y) \; d\y = \mathcal{K} q_{\rm out}(\x), \quad \forall \x \in \Gamma .
\end{equation}
Die Einstrahlung in einem Punkt $\x$ der Oberfläche ergibt sich durch Anwendung des Integraloperators $\mathcal{K}$ auf die Lösung der Gleichung $\eqref{IGLqout2}$.

Arbeiten Sie an großen Strahlungsmodellen?

Wenn Sie einen Anwendungsfall mit vielen Strahlungsflächen oder hohem Rechenaufwand haben, freue ich mich über den fachlichen Austausch. Gern prüfe ich, ob das Verfahren für Ihren Fall geeignet ist. Kontakt aufnehmen.

Quellen

Michael F. Cohen and John R. Wallace (1993): Radiosity and Realistic Image Synthesis. 1993, ISBN: 0121782700.
Wolfgang Hackbusch and Zenon Paul Nowak (1989): On the fast matrix multiplication in the boundary element method by panel clustering. In: Numerische Mathematik, Bd. 54, Nr. 4, S. 463–491, 1989, ISSN: 0945-3245.
H. A. Van der Vorst (1992): BiCGSTAB: A fast and smoothly converging variant of Bi–CG for the solution of nonsymmetric linear systems.. In: SIAM J. Sci. Statist. Comput., Bd. 13, Nr. 2, S. 631–644, 1992.
R. Siegel and Howell (2002): Thermal Radiation Heat Transfer. 4th Edition, Taylor & Francis, New York, London, 2002.