Abstract
Liouville’s equation on phase space in geometrical optics describes the evolution of an energy distribution through an optical system, which is discontinuous across optical interfaces. The discontinuous Galerkin spectral element method is conservative and can achieve higher order of convergence locally, making it a suitable method for this equation. When dealing with optical interfaces in phase space, non-local boundary conditions arise. Besides being a difficulty in itself, these non-local boundary conditions must also satisfy energy conservation constraints. To this end, we introduce an energy conservative treatment of optical interfaces. Numerical experiments are performed to prove that the method obeys energy conservation. Furthermore, the method is compared to the industry standard ray tracing. The numerical experiments show that the discontinuous Galerkin spectral element method outperforms ray tracing by reducing the computation time by up to three orders of magnitude for an error of \(10^{-6}\).
Similar content being viewed by others
1 Introduction
Illumination optics deals with the design of optical components for various applications, like LED lighting [28] and automotive headlamps [10, 36]. The design of these optical components requires a different approach than used in imaging optics, as imaging effects are highly undesirable [9]. Light propagation through an optical system is usually computed using ray tracing [14]. This means computing the evolution of many rays through an optical system, where upon hitting an optical interface Snell’s law of refraction or the law of specular reflection have to be applied. Typical optical interfaces are lenses or mirrors. Ray tracing is commonly employed to directly obtain the illuminance or intensity on a target. Although, forward quasi-Monte Carlo ray tracing converges rather slowly with rates close to \(\mathcal {O} \left( N_{RT} ^{-1} \right) \) where \(N_{RT}\) denotes the number of rays, it is the industry standard. For details on quasi-Monte Carlo integration, see e.g. [24].
The analysis of light propagation using a phase space description provides a new approach to understanding optical systems [17, 29, 35]. Phase space, being defined as the collection of all positions and direction coordinates of rays, provides a complete description of the spatial and angular distribution of light. A point in phase space evolves according to a Hamiltonian, describing the evolution of one single light ray, whenever the refractive index is smooth. When a ray hits an optical interface the laws for refraction or reflection have to be applied. In [11,12,13] new ray tracing methods are presented based on the phase space description. These methods allow for tracing of fewer rays to achieve the same accuracy as classical ray tracing.
An alternative approach to ray tracing is based on directly obtaining an energy distribution on phase space, rather than its integrated quantities such as the illuminance or luminous intensity. The propagation of an energy distribution, related to the luminance, through an optical system is governed by Liouville’s equation for geometrical optics on phase space. Recently, numerical schemes for Liouville’s equation were developed that incorporate the optical interfaces. In [33] van Lith et al. derived a first-order upwind finite difference scheme and in [32] van Lith et al. introduced a third-order active flux finite volume scheme on moving meshes. The third-order active flux scheme was proved to be faster and more accurate compared to classical ray tracing for obtaining an energy distribution on phase space. Additionally in [31] van Lith made use of the discontinuous Galerkin spectral element method to solve Liouville’s equation.
The discontinuous Galerkin spectral element method (DGSEM) discussed by Kopriva in [21] is a collocation scheme for the semi-discretisation of the spatial domain for conservation laws, leaving the time-like variable continuous. In terms of Liouville’s equation this entails discretising phase space. The phase space domain is partitioned into elements with each element having interior nodes placed at collocation points. The solution is approximated using a polynomial, where the polynomial degree determines the number of interior nodes used for each element. Consequently the method has an extraordinary flexibility as it is an hp-method, where h refers to the mesh size and p to the polynomial degree, i.e. the accuracy can be increased by decreasing the mesh size or by increasing the polynomial degree. Additionally, the method does not enforce continuity across the boundary of each element, making it particularly suitable for the discontinuous solutions across optical interfaces.
At an optical interface Liouville’s equation is not valid. Instead, Snell’s law of refraction or the law of specular reflection describe the discontinuous change in the direction coordinate, i.e., a jump in phase space. In what follows, we will see that this results in non-local boundary conditions for the energy distribution in phase space. Our contribution consists of describing the treatment of these optical interfaces so that they obey energy conservation. In the DGSEM the elements communicate using numerical fluxes. Snell’s law and the law of specular reflection are incorporated in these numerical fluxes at an optical interface. In addition to the discontinuous change in the direction coordinate described by these laws, a single element before the optical interface might contribute to multiple elements after the optical interface. This connection to multiple elements is similar to fully non-conforming geometries when using subdomain refinement [3, 4]. Kopriva et al. outlined such a strategy for the DGSEM in [23]. In [5] an analysis of this method is presented by Bui-Thanh and Ghattas. Across an optical interface the numerical fluxes are discontinuous and therefore we have to take a different approach. Inspired by [23], we present a method that directly incorporates the laws of optics and obeys energy conservation.
The article is outlined as follows: in Sect. 2 we discuss the conserved quantities in an optical system and Liouville’s equation, and in Sect. 3 we discuss the DGSEM. In Sect. 4 we discuss the energy conservative treatment of the optical interfaces, and in Sect. 5 we present numerical experiments proving energy conservation for two examples. The first example features a smooth refractive index field, while the second example is a test case featuring a discontinuity in the refractive index described by van Lith et al. in [33]. In the latter example, we compare the DGSEM for solving Liouville’s equation to quasi-Monte Carlo ray tracing for obtaining the illuminance. Finally we present our conclusions in Sect. 6.
2 Conserved Quantities and Liouville’s Equation
In optics we consider the transfer of luminous flux between surfaces. A source emits a beam of radiation or light, carrying a finite amount of luminous flux denoted by \(\varPhi \). In the absence of losses by absorption or scattering in an optical system, the total flux \(\varPhi \) is conserved, i.e., energy throughout the optical system is conserved. A related quantity is the luminance denoted by \(\rho ^*\), which is defined as [9, 25]
where \(\mathrm {d}\varPhi \) is an infinitesimal amount of flux carried by an infinitesimal beam, \(\mathrm {d}\omega \) is an element of solid angle subtended at the center of the source by the area at the detector \(\mathrm {d}A\), and \(\mathrm {d}A \cos \theta \) the projected area perpendicular to the beam, i.e., \(\mathrm {d}A\) is an element of the surface area of the detector and \(\theta \) describes the angle between the normal of the detector and the beam. The solid angle describes a cone on the unit sphere with the center of the source as its vertex and \(\mathrm {d}\omega \) the area on the unit sphere subtended by the cone.
Another important quantity that is also conserved in an optical system is étendue, which is defined by [9]
where n is the refractive index in which the beam is immersed. This allows us to write the luminance as
When a beam is propagating through a homogeneous medium the luminance \(\rho ^*\) is conserved, as is implied by conservation of energy and conservation of étendue. When a beam of light strikes an optical interface, e.g., a lens or a mirror, the beam is subject to Snell’s law of refraction or the law of specular reflection. In the case where the beam is refracted, e.g., a transition from a medium with refractive index \(n_1\) to a medium with refractive index \(n_2\), the luminance is not conserved. Assuming no Fresnel reflections, and applying conservation of energy and étendue, we obtain [25, 26]
where \(\rho _\text {i}^*\) and \(\rho _\text {t}^*\) describe the incident and transmitted luminance, respectively. The quantity \(\rho ^* / n^2\), known as basic luminance is conserved for refractions, cf. (4). A similar result can be derived for reflections, where the refractive indices are equal. Relation (4) will be referred to as basic luminance invariance. For a complete derivation including Fresnel reflections see [9, 25, 26].
The definitions of luminance, étendue and basic luminance invariance described above hold for three-dimensional optics, whereas in two-dimensional optics the definitions are slightly altered. For more details, see [9]. In summary, we denote the basic luminance for both two- and three-dimensional optics by \(\rho \), which is defined by
where the étendue \(\mathrm {d}\mathcal {U}\) for two- and three-dimensional systems reads [9]
where for 2D optics \(\theta \) denotes the angle between the normal of the detector and the beam, and \(\mathrm {d}\theta \) an element of angle subtended at the center of the source by the infinitesimal line segment at the detector \(\mathrm {d}l\). The basic luminance is related to the luminance by \(\rho = \rho ^* / n^2\) for three-dimensional optics, whereas for two-dimensional optics \(\rho = \rho ^* / n\).
2.1 Liouville’s Equation
In geometrical optics the evolution of light rays in a beam of light can be cast in a Hamiltonian system, where we denote with \({\varvec{q}} \in \mathbb {R}^d\) the position and \({\varvec{p}} \in \mathbb {R}^d\) the momentum coordinates [35]. For two-dimensional optics \(d = 1\), while for three-dimensional optics \(d = 2\). Both terms together describe a point \(({\varvec{q}}, {\varvec{p}})\) in phase space, where the phase space \(\mathcal {P}\) is defined as the collection of all positions \({\varvec{q}}\) and momenta \({\varvec{p}}\) at a certain position along the optical axis denoted by the z-coordinate. A point in phase space evolves when we move along the optical axis.
The momentum \(\mathbf {{p}} = (\varvec{p}, p_z) \in \mathbb {R}^3\) is restricted to Descartes’ sphere \(\left| \mathbf {{p}} \right| = n(z, {\varvec{q}})\) where n is the refractive index field as a function of the three-dimensional position coordinates \(\mathbf {{q}} = (\varvec{q}, z)\) [35]. This restriction invites us to use spherical coordinates to represent the momentum vector \(\mathbf \mathbf{{p}}\) as
where \(\theta \) represents the polar angle, describing the angle between the z-axis and \(\mathbf {{p}}\) measured from the z-axis, and \(\varphi \) the azimuthal angle for describing the direction in the \(\varvec{q}\)-plane. Therefore, at a given position \(z_0\) along the optical axis, one can visualise the phase space coordinates on the screen that is perpendicular to the z-axis and intersects the z-axis at \(z_0\), where \({\varvec{q}}\) is the position on the screen and \({\varvec{p}}\) describes the projection of \(\mathbf {{p}}\) on the screen [32]. The restriction of the momentum for physical rays \(\mathbf {{p}}\) to Descartes’ sphere also implies that the two-dimensional momentum vector is restricted by \(\left| {\varvec{p}} \right| \le n\), describing a region known as Descartes’ disc [35].
The phase space coordinates of a light ray evolve as a function of the distance along the optical axis according to Hamilton’s equations, which read
with \(h = h(z, \varvec{q}, \varvec{p})\) the optical Hamiltonian given by
Here, \(\sigma \in \{-1, 0, +1\}\) denotes the direction of the light ray travelling along the optical axis, with \(\sigma = 0\) being marginal rays that travel perpendicular to the optical axis [35]. For simplicity, we assume that all rays travel in the positive z-direction given by \(\sigma = +1\).
Hamilton’s Eqs. (8) hold for a single light ray, however, this may be generalised to a beam of light carrying a finite amount of energy in terms of luminous flux. The flow generated by Hamilton’s equations describes canonical transformations, otherwise known as symplectic transformations, on phase space. These transformations preserve the symplectic structure of phase space [1]. In other words the phase space volume element \(\mathrm {d}q_1 \mathrm {d}q_2 \mathrm {d}p_1 \mathrm {d}p_2\) is constant. In the context of optics this has the equivalent meaning of étendue conservation. In fact \(\mathrm {d} \mathcal {U} = \mathrm {d}q_1 \mathrm {d}q_2 \mathrm {d}p_1 \mathrm {d}p_2\), which can be obtained from the first two components of the three-dimensional momentum vector described by expression (7). The Jacobian determinant of \(\varvec{p}\) with respect to the polar and azimuthal angles \(\theta \) and \(\varphi \), can be computed as
with \(\mathrm {d}\omega = \sin \theta \mathrm {d}\theta \mathrm {d}\varphi \) an element of solid angle. Noting that the differential area on a screen can be written as \(\mathrm {d}A = \mathrm {d}q_1 \mathrm {d}q_2\) and substituting (10) into relation (6) for 3D optics, we obtain
The basic luminance invariance (4) implies that \(\rho \) remains constant if we move an arbitrary distance \(\varDelta z\) along the optical axis, i.e.,
Note that this relation also holds when a beam of light is reflected or refracted. If the solution is sufficiently smooth, one can derive Liouville’s equation by subtracting the right-hand side of (12) from its left-hand side and dividing by \(\varDelta z\) and subsequently taking the limit \(\varDelta z \rightarrow 0\), resulting in
Here, we have already applied Hamilton’s equations (8). The advective form of Liouville’s equation (13) may be written in conservative form by assuming that h is twice differentiable, upon which we obtain
with \(\nabla = ( \frac{\partial }{\partial \varvec{q}}, \frac{\partial }{\partial \varvec{p}} )^\mathrm{T}\) and the flux vector \(\varvec{f} = \varvec{f}(\varvec{q}, \varvec{p})\) defined as
where we have used that the velocity field \(\varvec{u}\) is divergence-free, and the superscript T denotes transpose. Note that an optical interface causes the Hamiltonian to be discontinuous. Therefore, at an optical interface, Liouville’s equation does not hold and we must apply (12) together with Snell’s law and/or the law of specular reflection, in the limit \(\varDelta z \rightarrow 0\).
Solving Liouville’s equation on phase space at any point z along the optical axis tells us how the basic luminance changes when light propagates through an optical system, allowing us to compute at each z-coordinate the related integral quantities such as luminous flux on the screen. The total luminous flux \(\varPhi \) in the optical system at \(z = \mathrm {const}\) reads
Here, the phase space dependence on the z-coordinate is denoted explicitly, since the momentum domain is restricted according to Descartes’ disc. Assuming the optical system is lossless the total luminous flux should be constant, i.e., \(\varPhi (z) = \varPhi (0)\).
An infinitesimal element of illuminance E is defined by
Applying definition (5) for the basic luminance and relation (11), \(\mathrm {d}E\) (16) can be rewritten as
Next, integrating over momentum space results in the illuminance \(E(z, \varvec{q})\), i.e.,
where \(\varvec{p} = (p_1, p_2) \in P(z)\) in which P(z) denotes the two-dimensional momentum space at a certain position z along the optical axis. Alternatively, an infinitesimal element of luminous intensity I is defined by
Applying again definition (5) for the basic luminance and definition (6) for the étendue, \(\mathrm {d}I\) (18) can be written as
Subsequently using the relation for \(p_z\) defined in (7) and \(\mathrm {d}A = \mathrm {d}q_1 \mathrm {d}q_2\), followed by integration over the position coordinates on the screen, denoted by \(\varvec{q} = (q_1, q_2) \in Q(z)\), we obtain
With these definitions, the main quantities of interest in optics can thus be easily computed from the basic luminance, satisfying Liouville’s equation. In the next section, we explore a method for solving Liouville’s equation.
3 Derivation of DGSEM
In what follows, we restrict ourselves to two-dimensional optical systems, hence reducing the complexity from a four-dimensional phase space to a two-dimensional phase space with position coordinate q and momentum coordinate p, and the distance along the optical axis denoted by z. Note that the position and momentum are now scalar quantities, therefore we omit the bold-face notation for these quantities. Next, we outline a spatial semi-discretisation of Liouville’s equation, leaving only z continuous. For the semi-discretisation we apply the discontinuous Galerkin spectral element method (DGSEM) described by Kopriva in [21], to the two-dimensional Liouville equation for \(\rho = \rho (z, q, p)\) in conservative form
where \(\nabla = \left( \frac{\partial }{\partial q}, \frac{\partial }{\partial p} \right) ^\text {T}\) and the flux vector \(\varvec{f}\) now reads
The Hamiltonian h for two-dimensional optics reduces to
and consequently the velocity \(\varvec{u}\) reads
For phase space discretisation, the two-dimensional phase space domain \(\mathcal {P}\) is covered with straight-sided quadrilaterals \(\varOmega ^k \subset \mathcal {P}\) with k the index of the element. In a more general discretisation, the boundaries of quadrilaterals are allowed to be curved, such that curved boundaries from physical constraints can be modelled appropriately. In fact, when the refractive index field changes continuously as a function of \({\varvec{q}}\), then the maximum allowed momentum varies as a function of q due to the restriction of \({\varvec{p}}\) to Descartes’ sphere. This restriction can be accommodated by curved boundaries when solving Liouville’s equation, see [31]. For a discussion on DGSEM with curved quadrilateral elements, see for example [8, 18, 21, 22]. In this paper we only consider straight-sided quadrilaterals.
Each quadrilateral \(\varOmega ^k\) has four vertices \(\{ \varvec{x}_1, \varvec{x}_2, \varvec{x}_3, \varvec{x}_4 \}\) labelled in counter-clockwise direction where \(\varvec{x} = \left( q, p \right) ^\text {T}\) and we have omitted the element index (superscript k), see Fig. 1. For ease of computation, the reference square \(\chi = [-1, 1]^2\) is mapped to each quadrilateral \(\varOmega ^k\), transforming a point in the reference domain \((\xi , \eta ) \in \chi \) to a point in physical space \(\varvec{x}(\xi , \eta ) \in \mathcal {P}\) using the following bilinear transformation
The Jacobian of the transformation is given by \(\frac{\partial (q, p)}{\partial (\xi , \eta )} = \left( \frac{\partial \varvec{x}}{\partial \xi }, \frac{\partial \varvec{x}}{\partial \eta } \right) \), where the columns read
The divergence term in (20a) can be rewritten by applying the chain rule resulting in
where \(\mathcal {J} = \frac{\partial q}{\partial \xi } \frac{\partial p}{\partial \eta } - \frac{\partial q}{\partial \eta } \frac{\partial p}{\partial \xi }\) denotes the Jacobian determinant, \(\nabla _{\varvec{\xi }} = \left( \frac{\partial }{\partial \xi }, \frac{\partial }{\partial \eta } \right) ^\text {T}\) and \(\widetilde{\varvec{f}}\) is an auxiliary flux defined by the product of the adjoint Jacobian matrix and the flux \(\varvec{f}\), i.e.,
Mapping from reference square \(\chi \) to a quadrilateral \(\varOmega ^k\)
Applying the transformation (25) to Liouville’s Eq. (20a), we obtain
where \(\rho = \rho (z, \xi , \eta )\).
The weak formulation of Liouville’s equation is obtained by first multiplying the PDE (27) by the Jacobian determinant \(\mathcal {J}\) and by a smooth test function \(\phi \), and subsequently integrating over the reference domain \(\chi \). This results in
The second term is rewritten by applying the product rule and Gauss’s theorem, so that
where \(\hat{\varvec{n}}\) is the outward unit normal on \(\partial \chi \) and the orientation of the closed curve \(\partial \chi \) is counter-clockwise. Using this, we obtain the weak formulation of Liouville’s equation on the reference domain
Note that for strong solutions we require the flux to be differentiable, hence, h(z, q, p) should be twice differentiable. However, the DGSEM uses the weak form of the solution and only requires the flux to be continuous, therefore, h(z, q, p) being once continuously differentiable is sufficient. For typical optical interfaces this is not sufficient since the refractive index field is discontinuous and, therefore, h(z, q, p) and also the flux are discontinuous. In particular, for these interfaces we require a special treatment of the fluxes which we will discuss in Sect. 4.
3.1 Tools for Approximating the Solution
The solution \(\rho \) in Eq. (29) is approximated by an expansion in basis functions [21]. We choose one-dimensional basis-functions \(\varphi _i\), \(i = 0, ..., N\), for which \(\varphi _i(\xi _j) = \delta _{ij}\) holds for chosen points \(\xi _j\). Moreover, we require that the basis-functions form an orthogonal basis with respect to the standard \(L^2\)-inner product. A suitable choice are the Lagrange polynomials defined on Gauss-Legendre nodes. In the following, we will replace \(\rho \) by the approximation
where \(\rho _{ij}(z)\) are the expansion coefficients for the chosen basis. In this paper we will restrict ourselves to using the same basis functions in both directions with an equal number of expansion coefficients, although in general this restriction is not necessary.
The quadrature rule defined by Gauss-Legendre nodes \(\{ \xi _i \}_{i=0}^N\) and corresponding weights \(\{ w_i \}_{i=0}^N\) allows us to approximate the integral of any function g, i.e.,
with \(-1< \xi _i < 1\) and \(w_i > 0\). The quadrature rule on the reference domain has a tensor product structure, hence,
Thus, we place nodes inside the reference domain at the Gauss-Legendre nodes \((\xi _i, \eta _j)\). Furthermore, for one-dimensional integrals the Gauss-Legendre quadrature gives exact integration for at least all polynomials of degree \(2N + 1\).
Focusing on one dimension, the Lagrange polynomials on the Gauss-Legendre nodes read
which satisfy the Kronecker property, i.e.,
It can be readily verified, that the Lagrange polynomials defined on the Gauss-Legendre nodes are orthogonal with respect to the standard \(L^2\)-inner product, i.e.,
where the integration is exact since \(\ell _i(\xi ) \ell _j(\xi )\) is a polynomial of degree at most 2N.
A very useful property of Lagrange polynomials is that the polynomial interpolation of any function g is rather easy, i.e.,
where \(\mathcal {I}_N\) is the polynomial interpolation operator using \(N + 1\) nodes. The approximation of the derivative of g is defined as the derivative of the interpolant, i.e.,
Note that if g is a polynomial of degree N or less, both the interpolant and the derivative are exact. In what follows, we require the derivative at the nodes, i.e.,
where \(D_{ij} = \frac{\mathrm {d} \ell _j}{\mathrm {d} \xi } (\xi _i)\). The elements of the differentiation matrix \(\varvec{D} = ( D_{ij} )\) can be found by differentiating the Lagrange polynomial followed by evaluation at the node, and thus read
and the diagonal elements read
due to the fact that the derivative of a constant function vanishes.
3.2 Approximating the Solution
To derive an approximation of the solution, we expand both the solution and the flux in Lagrange polynomials. The expansions read
The coefficients, indicated by the index-subscript ij, in each expansion are related to the position of an element’s interior node \((q_{ij}, p_{ij})\), by \(\rho _{ij}(z) = \rho (z, q_{ij}, p_{ij})\) and \(\widetilde{\varvec{f}}_{ij}(z) = \widetilde{\varvec{f}}(z, q_{ij}, p_{ij})\). The auxiliary flux coefficients \(\widetilde{\varvec{f}}_{ij}(z)\) are related to \(\rho \) by
with \(\widetilde{\varvec{u}}_{ij}\) the transformed velocity, similarly defined to (26). Here the velocity \(\widetilde{\varvec{u}}_{ij}(z) = \widetilde{\varvec{u}}(z, q_{ij}, p_{ij})\) depends on z if the refractive index n depends on z. In the following, we omit \(\widetilde{\varvec{f}}_{ij}\)’s dependence on z for ease of notation.
Next, we have to approximate the integrals in Eq. (29). The test function \(\phi \) is chosen to be in the same basis as the solution \(\rho \), resulting in a Galerkin method. Therefore, taking
allows us to derive \((N+1)^2\) equations for the \((N+1)^2\) coefficients \(\rho _{ij}\). Combining this together with the approximations (40a) and (40b) for \(\rho \) and \(\widetilde{\varvec{f}}\) we can approximate the integrals using the Gauss-Legendre quadrature rules. Therefore, substituting the approximation (40a) in the first term of (29), we obtain
Applying the Kronecker property (34) of the Lagrange polynomials, the sums reduce to
where \(\mathcal {J}_{ij} := \mathcal {J}(\xi _i, \eta _j)\). Note that the integral is exact for the given combination of a bilinear mapping \(\varvec{x}(\xi , \eta )\) and Lagrangian polynomials, since then the integrand is a polynomial of degree \(2N+1\) in \(\xi \) and in \(\eta \). The Gauss-Legendre quadrature rule is exact for this bivariate polynomial.
For the third term in (29), we substitute the approximation (40b) and denote \(\widetilde{\varvec{f}} = (\widetilde{f}, \widetilde{g})\), resulting in
where we have used the definition of the differentiation matrix (39). Furthermore, we introduce the following auxiliary matrix
for ease of computation. The third term then reads
In what follows, we will replace the flux appearing in the boundary integral from Eq. (29) with a numerical flux \(\widetilde{\varvec{F}} = \left( \widetilde{F}, \widetilde{G} \right) \). The boundary integral can be split into four parts and evaluated for each boundary segment, see Fig. 2. Along each segment the numerical flux \(\widetilde{\varvec{F}}\) is described by an Nth degree polynomial at the boundary nodes shown in the figure. For the bottom part, with \(\eta = -1\), the integral can be exactly evaluated using Gauss-Legendre quadrature, such that we obtain
Similarly, we can compute the other components and the result for the full boundary integral reads
In the discontinuous Galerkin spectral element method the elements communicate by fluxes through the faces of each element. The solution on the boundary between two elements is allowed to be discontinuous, thus the limit towards the boundary of an element can have two values, one for each element it touches. The flux on the boundary must be replaced by a numerical flux so that the neighbouring elements can communicate. The numerical flux depends on the values of \(\rho \) just left and right of the boundary, i.e., \(\rho _\text {L}\) and \(\rho _\text {R}\), which are computed by evaluating the interior solution (40a) at the boundary. For the numerical flux we take the upwind flux. Due to the transformation (26) of the flux to the reference domain, the physical upwind flux F is scaled at an edge by \(\varDelta l / 2\) with \(\varDelta l\) the length of the edge, such that the upwind flux \(\widetilde{F} = \tfrac{1}{2} \varDelta l F\) over an edge reads
where \(\hat{\varvec{n}}\) denotes the outward normal vector w.r.t. the element left of the boundary.
Reference square with polynomial degree \(N = 4\). The normals are denoted by arrows, interior nodes are denoted by filled circles and boundary points by open squares
Next, we substitute expressions (43), (45) and (47) in equation (29), so that we obtain the semi-discrete ODE system for the expansion coefficients \(\rho _{ij}(z)\):
with the numerical fluxes \(\widetilde{\varvec{F}} = \big ( \widetilde{F}, \widetilde{G} \big )\) given by (48). This ODE system can be solved using any numerical time integrator, e.g., the classical fourth order Runge–Kutta method. Other popular choices in the literature are explicit low-storage Runge-Kutta methods, see [6, 19, 34].
The discontinuous Galerkin spectral element method approximates the exact solution by an Nth degree polynomial, so the global spatial error e for a typical mesh size \(\varDelta x\) behaves as
Furthermore, the scheme is restricted by stability in terms of a CFL condition. For discontinuous Galerkin methods on quadrilaterals there is no direct known bound for the CFL condition. For triangular grids the relation between the Courant number and the shape of the triangles is studied in [7, 30].
4 Optical Interfaces
In the phase space representation, the flow of \(\rho \) describes a beam of light propagating through an optical system. When the beam hits an optical interface, the momentum p changes discontinuously according to the law of specular reflection or Snell’s law of refraction. Furthermore, from the discussion in Sect. 2 we know that the total luminous flux should remain constant throughout the optical system. The numerical solution should respect the actual physics, therefore, the discontinuity at optical interfaces coupled with conservation of energy should be incorporated into the DGSEM when we solve Liouville’s equation.
In the DGSEM the solution is allowed to be discontinuous across the boundary connecting two or multiple elements. Therefore, the mesh in phase space is aligned such that elements adjacent to the interface have edges that coincide with the optical interface, across which the solution is discontinuous [31]. The elements in the DGSEM communicate through numerical fluxes, hence, we have to incorporate both Snell’s law and the energy conservation constraint in the numerical flux when integrating (49). In particular, for a beam of light moving towards an optical interface, i.e., the velocity is directed towards the interface, we have to leave \(\rho \) free, whilst for a beam moving away from an optical interface, i.e., the velocity is directed away from the interface, we have a Dirichlet boundary condition for \(\rho \) due to Snell’s law or the law of specular reflection [31].
Refraction or reflection causes the elements to be connected in a non-trivial manner at the optical interface. For example, one single element, on the side where light is moving towards the interface, can contribute to multiple elements on the other side. This occurs because both Snell’s law and the law of reflection are non-linear in the momentum p. Therefore, this requires special treatment of the numerical fluxes to ensure that the scheme obeys energy conservation.
Van Lith et al. present Snell’s function in [33], which is an explicit version of Snell’s law and the law of specular reflection combined on phase space. Snell’s function \(\mathcal {S}\) relates the momentum p of an incident ray to the outgoing momentum \(\bar{p}\) of the ray. Let \(n_1\) be the refractive index of the incident medium and \(n_2\) the index of the transmitted medium. Then for a generic two-dimensional optical interface in the (q, z)-plane with surface unit normal \({\mathbf {\nu }} = ({{\nu }_q}, {{\nu }_z})\) directed towards the incident medium, Snell’s function reads
with the auxiliary variables \(\delta \) and \(\psi \) defined by
In the expression for \(\psi \) the plus sign should be taken for rays that propagate in the positive z-direction, while the minus sign should be taken for rays that propagate in the negative z-direction. Furthermore, the sign of \(n_2\), i.e., \({{\,\mathrm{sgn}\,}}(n_2)\), in the first case of (51a) can be used to accommodate embedded mirrors in a medium of refractive index \(n_1 \ge 1\) by taking \(n_2 = -n_1\), see [32]. Note that there is a so-called critical momentum \(p_\text {c}\) when \(n_1 \ge n_2\) so that \(\delta = 0\). For \(\delta < 0\) all light will be reflected, referred to as total internal reflection (TIR), while for \(\delta \ge 0\) light will be refracted. The outgoing momentum is computed as \(\bar{p} = \mathcal {S}(p; n_1, n_2, \mathbf {\nu })\), for which we will frequently use the shorthand notation \(\bar{p} = \mathcal {S}(p)\) and take the other parameters as given. Furthermore, the inverse of Snell’s function will also be frequently used, for which we will use the shorthand notation \(p = \mathcal {S}^{-1}(\bar{p})\). This means, find the momentum p such that \(\bar{p} = \mathcal {S}(p)\). For example, for refraction the inverse reads [33]
Snell’s function (51) combined with (12) results in [33]
where \(p^+ = \mathcal {S}(p^-; n_1, n_2, \mathbf {\nu })\) and the ± denote one-sided limits towards the optical interface. This relation allows us to relate the basic luminance \(\rho \) on both sides of the interface.
To elaborate the energy conservation constraint, we consider the following flat optical interface parallel to the z-axis
Note that the optical interface in phase space is represented by a line parallel to the p-axis. The optical interface has two sides where on one side the normal in phase space is directed towards \(q < q_0\) and describes the part with refractive index \(n_1\), whereas on the other side the normal is directed towards \(q > q_0\) and describes the part with refractive index \(n_2\). The normal in phase space is given by \(\hat{\varvec{n}} = ( \pm 1, 0 )\) with the plus sign for the direction towards \(q > q_0\). Since the optical interface is represented by a line parallel to the p-axis at some constant q-value, only the q-component of the flux (20b) needs to be considered, i.e.,
cf. (22).
In what follows, we assume that light is initially in the medium with refractive index \(n_1\). We partition the optical interface, represented in phase space, into line segments for both sides of the optical interface. The partitioning is based on whether light is moving towards or away from the optical interface.
The line segment on the side of the optical interface with velocity directed towards the optical interface is denoted L. The line segment L describes the incoming momentum of light from the medium with refractive index \(n_1\), due to the assumption of light being initially in this medium.
The line segments just on the optical interface with velocity directed away from the optical interface is denoted R. The line segments R describe the outgoing momentum of light, and is further split into two parts, i.e., \(R = R_\text {R} \cup R_\text {T}\). The line segment denoted \(R_\text {R}\) represents the momentum of light after total internal reflection, and is part of the optical interface where the refractive index is \(n_1\), while the line segment denoted \(R_\text {T}\) is in the medium with refractive index \(n_2\), and represents the momentum after transmission.
To distinguish the momentum taken from either side of the optical interface, we write \(p \in L\) and \(\bar{p} \in R\). First, consider the integral of the flux entering an arbitrary momentum interval \([\bar{p}_1, \bar{p}_2] \subseteq R_\text {T}\). The integral reads
where \(q_0^+\) denotes the limit towards the optical interface from the line segment \(R_\text {T}\). Relation (53) implies
where \(q_0^-\) denotes the limit towards the optical interface from the line segment L. Subsequently, we transform the integral using \(\bar{p} = \mathcal {S}(p) = \mathcal {S}(p; n_1, n_2, \mathbf {\nu })\) resulting in
where \(\bar{p}_i = \mathcal {S}(p_i)\) for \(i = 1, 2\), and \([p_1, p_2] \subseteq L\).
The relation for reflection can be derived similarly by considering the integral of the flux entering an arbitrary momentum interval \([\bar{p}_3, \bar{p}_4] \subseteq R_\text {R}\). We obtain the relation
with \(\bar{p}_i = \mathcal {S}(p_i)\) for \(i = 3, 4\), and \([p_3, p_4] \subseteq L\). The relations (58) and (59) describe how the fluxes leaving L are related to the fluxes entering \(R_\text {T}\) or \(R_\text {R}\), respectively. Henceforth, they are known as energy conservation constraints.
For the particular optical interface (54) the constraints can be simplified. Since, we assumed that light was initially in the medium with refractive index \(n_1\), the optical interface normal on the full position space is equal to \(\mathbf {\nu } = ({{\nu }_q}, {{\nu }_z}) = (-1, 0)\). Snell’s function (51) reduces for this flat interface to
with \(p_\text {c} = \sqrt{n_1^2 - n_2^2}\). Then, the energy conservation constraint for \(R_\text {T}\), given by (58), can be simplified by noting that for refraction the following relations hold
Hence, we obtain
and similarly for reflection the constraint for \(R_\text {R}\), given by (59), reduces to
The balances (61) have to be combined with relation (53) to ensure the scheme conserves energy. However, the coupling between the line segments L and R is not straightforward, which will be discussed in the next section.
Conservative handling of fluxes. The incident and transmitted momenta are related by \(\bar{p} = \mathcal {S}(p)\)
4.1 Conservative Handling of Fluxes
First, consider only the refractive part of the optical interface. Elements adjacent to the optical interface in phase space are shown in Fig. 3a. These elements have edges on the optical interface and these edges are denoted by \(L_i\) (\(i = 1, 2\)) and \(R_j\) (\(j = 1, 2, 3, 4\)). Due to refraction, the value of \(\rho \) in the elements that contain \(R_j\) as an edge is determined by the flow through the elements that contain \(L_1\) and \(L_2\). In fact, taking a closer look at how Snell’s function connects the line segments from L to R in momentum space at the optical interface, we obtain for example Fig. 3b. In Fig. 3b \(L_i\) and \(R_j\) denote line segments along the optical interface. The basic luminance \(\rho \) along the line segments \(L_i\) and \(R_j\) are represented by their inner-element solution evaluated at the optical interface. To simplify notation, we denote these polynomials along the optical interface by \(\rho ^{L_i}(p)\) with \(i = 1, 2\) and \(\rho ^{R_j}(p)\) with \(j = 1, 2, 3, 4\). For example:
where \(\zeta = \zeta (p)\) denotes the line segment’s local reference coordinate along the interface.
In Fig. 3b also virtual line segments \(\bar{L}_i\) are shown. The virtual line segment \(\bar{L}_i\) is the image of \(L_i\) under \(\mathcal {S}\), i.e.,
Hence, the endpoints of these line segments are found applying Snell’s function to the endpoints of \(L_i\), i.e., \(\bar{p}_i^L = \mathcal {S}(p_i^L)\). Note that due to Snell’s function, the line segments \(L_i\) are stretched or compressed in the momentum direction. Computing the endpoints \(\bar{p}_i^L\) allows us to determine which line segments before the optical interface contribute to a single line segment after the optical interface. From the figure we see that part of \(L_1\) contributes to \(R_2\) (the blue coloured region). Therefore, a relation connecting \(\rho ^{L_1}(p)\) and \(\rho ^{R_2}(p)\) on opposite sides of the optical interface must be found. Hence, as a first step applying relation (53) to a polynomial on \(L_i\), allows us to find the corresponding \(\rho \) on \(\bar{L}_i\), i.e.,
with \(\bar{p} = \mathcal {S}(p)\).
The coupling between line segments that do not exactly match, as shown in Fig. 3b, is similar to what is known as a geometrically non-conforming mesh [20, 23]. In [23] the authors describe a discontinuous Galerkin method for non-conforming meshes, applied to Maxwell’s equations that form a hyperbolic system of PDEs. In their approach to treating non-conforming interfaces the solutions are first transferred to an intermediate construct called a ‘mortar’, and on this mortar the numerical fluxes are computed and transferred back to the corresponding elements. The transfer of the solutions and numerical fluxes is done using a least-squares matching, with integrals evaluated using Gauss-Legendre quadrature [5].
We will take a slightly different approach since in Liouville’s equation for optics the flux \(\varvec{f}\) is discontinuous across an optical interface. Instead, relation (53) and Snell’s function together describe how \(\rho \) transforms across an optical interface, cf. (64), therefore, a least-squares matching of the polynomials describing \(\rho \) along either side of the interface is used with Snell’s function directly incorporated and an additional constraint is used to satisfy energy conservation.
For the reflective part of a flat interface that is parallel to the z-axis, Snell’s function reduces to \(\bar{p} = -p\), see (60). The conservative treatment of these types of optical interfaces is easily accommodated by choosing a mesh such that the elements and nodes are symmetric with respect to the line \(p = 0\) and, therefore, the constraint (61b) is easily satisfied. Due to this choice of mesh each node \(\bar{p}_j \in R_\text {R}\) will exactly correspond to \(-\bar{p}_j = p_j \in L\) and a point-by-point transfer of \(\rho \) can be made. Henceforth, the following exposition of the method describes the method considering only refraction.
4.2 Contribution from One Element
From Fig. 3b we see that the line segment \(R_2\) only depends on the solution in \(L_1\). The polynomial \(\rho ^{R_2}\) must thus be computed from the polynomial \(\rho ^{L_1}\) with the additional constraint of energy conservation. That is, the integral of the flux within the blue interval on either side of the optical interface should be equal analogous to equation (61a). Therefore, the constrained least-squares approximation reads
Here, \([p_2^R, p_3^R] \subseteq L_1 = [p_1^L, p_2^L]\) and the momenta on both sides are related by \(p_i^R := \mathcal {S}^{-1}(\bar{p}_i^R)\), see Fig. 3b. Furthermore, the numerical fluxes are defined as expansions in the Lagrange polynomial basis on Gauss-Legendre nodes, similar to (62), with flux coefficients \(F_j := u_j \rho _j\). The minimisation of the integral in (65a) requires finding a polynomial that matches in the least-squares sense, while the constraint (65b) ensures that the scheme conserves energy.
The integrals in the constrained minimisation problem (65) are transformed to reference line segments, so that we can compute the integrals using Gauss-Legendre quadrature. To be more specific, the integral on the LHS of (65b) and the integral in (65a) are transformed to the reference line segment along \(R_2\), while the integral on the RHS of (65b) is transformed to the reference line segment along \(L_1\). Omitting the element’s subscripts, applying relation (64) and introducing an auxiliary function \(\varXi \) for ease of notation, we obtain
where \(\varDelta \bar{p}^R := \bar{p}_3^R - \bar{p}_2^R\) and \(\varDelta p^L := p_2^L - p_1^L\). Furthermore, the coefficients \(\sigma ^L \in [-1, 1]\) and \(\lambda ^L \in [0, 2]\) denote the offset and scaling in \(L_1\)’s reference frame, such that \(p(\sigma ^L) = p_2^R\) and \(p(\sigma ^L + \lambda ^L) = p_3^R\) in \(L_1\). Finally, the auxiliary function \(\varXi \) reads
This function relates the reference frame coordinates for a momentum interval \([\bar{p}^R, \bar{p}^R + \varDelta \bar{p}^R]\) past the optical interface to the reference frame coordinates on a momentum interval \([p^L, p^L + \varDelta p^L]\) before the optical interface.
Next, we write the constrained minimisation problem (66) in terms of a Lagrange function \(\mathcal {L}\) with a Lagrange multiplier \(\mu \) for the energy conservation constraint, i.e.,
The coefficients \(\rho ^R_j\) for the polynomial \(\rho ^R \in \mathbb {P}_N\) can then be computed by solving
Recalling that both \(F^L\) and \(F^R\) are written as expansions in the Lagrange polynomial basis on Gauss-Legendre nodes, we obtain for the energy conservation constraint (69b)
where we have used \(F_j = u_j \rho _j\).
To evaluate the second integral, we transform it to the reference interval \([-1, 1]\) using \(\zeta (\xi ) = \sigma ^L + \frac{1}{2} \lambda ^L \left( 1 + \xi \right) \). Next, we replace both integrals with Gauss-Legendre quadrature to find the exact values, since both integrands are at most Nth degree polynomials, therefore we obtain
with \(\xi _k\) and \(w_k\) the Gauss-Legendre nodes and weights, respectively.
Recalling again that the polynomials \(\rho ^L\) and \(\rho ^R\) are written as an expansion in a Lagrange polynomial basis, cf. (62), we can rewrite the equations (69a) to
For the evaluation of the first integral, we introduce a generic auxiliary variable \(S_{ij}\), given by
with \(\varXi \) defined in (67). The integral is evaluated by transforming to the reference interval and subsequently applying Gauss-Legendre quadrature. The integration interval \([\sigma ^R, \sigma ^R + \lambda ^R]\) of \(S_{ij}\) depends on how large \(R_j\) is compared to \(\bar{L}_i\). In this case, the entire line segment \(R_2\) fits in \(\bar{L}_1\), therefore, the integral over \(R_2\) is transformed to a reference line segment \([-1, +1]\), corresponding to \(\sigma ^R = -1\) and \(\lambda ^R = 2\).
The integrals in (72) are evaluated using Gauss-Legendre quadrature resulting in
with
and \(S_{ij}\) given by (73). The integral \(M_{ij}\) describes the orthogonality of the Lagrange polynomials on Gauss-Legendre nodes, and therefore, is easily evaluated, cf. (35), to be
The coefficients \(M_{ij}\) and \(S_{ij}\) are elements of the matrices \(\varvec{M}, \varvec{S} \in \mathbb {R}^{(N+1) \times (N+1)}\). The matrix \(\varvec{M}\) is simply a diagonal matrix containing the Gauss-Legendre quadrature weights, i.e., \(\varvec{M} = {{\,\mathrm{diag}\,}}( \varvec{w} )\) with \(\varvec{w} = (w_0, w_1, ..., w_N)^\text {T}\).
By defining
as the components of the vectors \(\varvec{\alpha }^R\) and \(\varvec{\beta }^L\), we can write the linear system given by (74) and (71) for \(\varvec{\rho }^R = (\rho _0^R, \rho _1^R, ..., \rho _N^R)^\text {T}\) and \(\mu \) compactly in matrix-vector form:
where we take the arguments for \(\varvec{S}\) as understood. Furthermore, \(\circ \) denotes the Hadamard product for vectors \(\varvec{a} = \left( a_0, a_1, ..., a_N \right) ^\text {T}\) and \(\varvec{b} = \left( b_0, b_1, ..., b_N \right) ^\text {T}\), i.e., \(\varvec{a} \circ \varvec{b} = \left( a_0 b_0, a_1 b_1, ..., a_N b_N \right) ^\text {T}\). Let us denote the matrix on the LHS by \(\varvec{A}\):
The determinant of this matrix reads
see Appendix A for details. Therefore, the matrix is only singular if all coefficients \(\alpha _i^R = 0\), or equivalent if all velocities \(u_i = 0\), which would mean no flux can enter the element from that side. Hence, we can safely assume that the matrix \(\varvec{A}\) is regular.
An analytical inverse for the matrix \(\varvec{A}\) is derived in Appendix A, and reads
where the coefficients of the matrix \(\varvec{B}\) read
Now, we can directly obtain an expression for the Dirichlet boundary condition values \(\varvec{\rho }^R\) in terms of \(\varvec{\rho }^L\), i.e.,
Note that for problems where the refractive index n does not depend on z, the coefficient matrix \(\varvec{C}\) relating \(\varvec{\rho }^R\) and \(\varvec{\rho }^L\) can be pre-computed and re-used during integration along the z-axis.
4.3 Contributions from Multiple Elements
From Fig. 3b we see that the element \(R_3\) depends on both \(\bar{L}_1\) and \(\bar{L}_2\). The idea remains the same, i.e., to use a least-squares matching with a constraint to ensure that the scheme is energy conservative. The constrained least-squares problem for \(R_3\) reads
where \(p_{34}^R := \mathcal {S}^{-1}(\bar{p}_{34}^R)\) and \(\bar{p}_{34}^R\) is the momentum value where the intervals \(\bar{L}_1\) and \(\bar{L}_2\) meet, see Fig. 3b. Furthermore, \(\rho ^{\bar{L}}\) contains the contributions from \(\rho ^{L_1}\) and \(\rho ^{L_2}\), and is defined by
The integrals in (83) are transformed to their respective line segments, e.g., the integral on the LHS of (83b) and the integral in (83a) are transformed to the reference interval \([-1, 1]\) along \(R_3\), such that we obtain
with
where we write R instead of \(R_3\) for brevity, and \(\kappa \) is defined such that \(p(\kappa ) = \bar{p}_{34}^R\) in \(R_3\) and \(\varDelta \bar{p}^R := \bar{p}_4^R - \bar{p}_3^R\), \(\varDelta p^{L_1} := p_2^L - p_1^L\) and \(\varDelta p^{L_2} := p_3^L - p_2^L\). Note that \(\sigma ^{L_1} + \lambda ^{L_1} = 1\) and \(\sigma ^{L_2} = -1\), however, for illustration purposes we will keep using the variables rather than these values. The Lagrange function \(\mathcal {L}\) for this constrained minimisation problem reads
The coefficients \(\rho _j^R\) for the polynomial \(\rho ^R \in \mathbb {P}_N\) can be found by solving
Following the same steps as in Sect. 4.2, we obtain the following system of equations
with
The linear system described by (87) can once again be assembled into a matrix-vector form:
where we have used the shorthand notation \(\varvec{S}^{L_1} = \big ( S_{ij}( -1, 1 + \kappa ) \big )\) and \(\varvec{S}^{L_2} = \big ( S_{ij}( \kappa , 1 - \kappa ) \big )\). Note that the matrix on the LHS is exactly the same as the matrix obtained in the previous section, except for possibly different values for \(\alpha _j^R\). Therefore, we can again solve the linear system explicitly for the Dirichlet boundary condition values \(\varvec{\rho ^R}\), resulting in
cf. (82). This result can of course be generalised to K elements contributing to \(\varvec{\rho }^R\), resulting in
4.4 Overview
To summarise, during a z-step the numerical fluxes over the optical interface are evaluated as follows. First, the elements are identified that have an edge on the optical interface. Those elements are separated into elements with velocities directed towards the optical interface, denoted L, and elements with velocities directed away from the optical interface, denoted R. For the elements from L the solution is evaluated at edges on the optical interface. The numerical flux over the edges for the elements L can be directly computed as there is no constraint on \(\rho \). For each element from R there is a Dirichlet boundary condition on the edge at the optical interface given by (53), that is incorporated into the numerical flux.
The value for the Dirichlet boundary condition is determined from the elements L, as follows. To determine which elements from L contribute to a single R element, Snell’s function is applied to the momentum boundaries of the elements L. Subsequently, the geometric quantities relating the element sizes are computed. Next, the momenta p at the quadrature nodes, for evaluation of the integral \(S_{ij}\), are determined. Subsequently, we apply \(\mathcal {S}^{-1}\) to these nodes, and compute \(\varXi \) using (67). Hereafter, the integrals \(S_{ij}\) are evaluated and the coefficients \(\alpha _j^L\), \(\alpha _j^R\) are computed. Finally, the values for the Dirichlet boundary condition can be found from their contributing L-elements by applying (90).
5 Results
Numerical experiments were performed for two examples. The first example features light propagating through a gradient-index medium. The smooth refractive index field of the medium fits naturally into the DGSEM for solving Liouville’s equation. For such optical systems ray-tracers usually have to resort to difficult to obtain closed-form expressions for the trajectory of the rays [2], or use symplectic integrators to solve Hamilton’s equations for every ray [27]. Solving Liouville’s equation with the DGSEM provides directly the energy distribution, i.e., the basic luminance \(\rho \) for the optical system. Furthermore, the method conserves energy by design.
The second example features a single optical interface. The problem exhibits both total internal reflection and refraction. At the optical interface we apply the strategy outlined in Sect. 4. Furthermore, a comparison is made between solving Liouville’s equation using the DGSEM and applying quasi-Monte Carlo ray tracing [14]. The illuminance is solved using both methods and the performance of both methods is tested.
5.1 Elliptic Waveguide
As a first example, we consider the elliptic waveguide [35] which features a smooth refractive index field given by
The parameters \(n_0\) and \(\kappa \) are taken to be \(n_0 = 1.4\) and \(\kappa = \sqrt{n_0^2 - 1}\). The refractive index field and several rays are shown in Fig. 4. We observe that the elliptic waveguide contains light much like an optical fibre. Hamilton’s equations (8) for rays inside the elliptic waveguide read
Since the refractive index field does not depend on z, the Hamiltonian h remains constant for each ray. The solution of (92) reads
where the initial conditions are given by \(( q(0), p(0) ) = (q_0, p_0)\). Note that from the refractive index field n and the Hamiltonian h we obtain [35]
where the right-hand side is constant when we move along the z-axis. We can readily see that the trajectories follow an elliptical path in phase space, hence, the name elliptic waveguide.
Elliptic waveguide: background colour indicates the refractive index value n(q), and the solid lines represent ray trajectories. Arrows indicate the direction of the ray, i.e., the momenta \((p_z, p)\)
Let the function \(\varphi _m\) be defined as
which is a \(C_0^m\)-function, meaning its first m derivatives are continuous and has compact support. The function \(\varphi _m\) is plotted in Fig. 5 for \(m = 7, 28\). We solve Liouville’s equation (20a) with the following initial condition
at \(z = 0\) and on the boundary of the domain we leave \(\rho \) free whenever the velocity field is pointing out of the domain, otherwise we prescribe \(\rho = 0\). In (96) we take \(m = 7\), \(\sigma _q = 0.25\) and \(\sigma _p = 0.1\).
Function \(\varphi _m\) for \(m = 7\) and \(m = 28\)
The ODE system (49) is integrated using the low-storage 4th order Runge-Kutta method by Zingg and Chisholm [37]. The numerical solution is integrated from \(z = 0\) to \(z = Z = 3\). The result using a 6th degree polynomial (\(N = 6\)) and \(K = 16 \times 16 = 256\) elements is shown in Fig. 6, together with the initial condition. The numerical solution at \(z = Z\) has roughly the same phase space area and is approximately a rotation of the initial condition. The scheme is also energy conservative up to machine precision, as can be seen in the plot of Fig. 7. Here, the relative error in the total luminous flux is plotted as a function of z. The total luminous flux is computed according to its definition (15) by applying a quadrature rule in agreement with the chosen polynomial degree.
Elliptic waveguide: basic luminance distributions \(\rho (z, q, p)\). Parameters are \(N = 6\), \(K = 256\), \(Z = 3\)
Elliptic waveguide: relative error in luminous flux
Furthermore, a convergence test is performed for this example by changing the number of elements K and varying the polynomial degree from \(N = 1, 2, ..., 6\). The numerical solution is compared to the exact solution, which can be found from the trajectory of the rays given by expressions (93). The expressions describe the evolution of a ray, given the initial conditions of the ray. For the analytical solution to Liouville’s equation at an arbitrary z we want to know where the ray originated from, since the phase space coordinates on the mesh are known. This amounts to tracing the ray backwards starting from an arbitrary z to \(z = 0\), resulting in
with q(z) and p(z) given in (93).
Using the exact solution (97) we can evaluate the discretisation error for which we take the \(L^1\)-norm, i.e.,
where \(\rho _\text {DG}\) denotes the numerical solution and \(\rho \) denotes the exact solution (97). The integrals in (98) are evaluated using Gauss-Legendre quadrature with \(N + 3\) nodes. The convergence order \(\gamma _\text {DG}\) is estimated from the empirical relation
with \(C_\text {DG} > 0\) an arbitrary constant.
The convergence data is shown in Table 1. The spatial discretisation is done using an Nth degree polynomial, and therefore the spatial order of accuracy is \(N + 1\). The temporal discretisation is done using a 4th order explicit Runge-Kutta method, where we choose \(\varDelta z\) to be the maximum allowable step such that the temporal integration is stable. Furthermore, a uniform rectangular mesh is used, where upon mesh refinement the mesh size in each direction is halved and similarly \(\varDelta z\) is halved to ensure stability.
The global error depends on whether the spatial or temporal discretisation errors dominate. From Table 1, we observe that the spatial discretisation error dominates for the polynomial degrees \(N = 1\) to \(N = 6\). Choosing a smaller \(\varDelta z\)-step in the numerical experiments did not influence the discretisation error. The results show that we obtain the expected \(N + 1\) order of convergence.
5.2 Bucket of Water
To illustrate that the strategy outlined in Sect 4 for handling optical interfaces is energy conservative, we apply it to a test case. The test case ‘bucket of water’ introduced by van Lith et al. [32, 33] is a suitable choice. The refractive index for this problem is given by
where we take \(n_1 = 1.4\) and \(n_2 = 1\). Using an initial basic luminance \(\rho _0\) that has support \(\mathcal {D}\) where \(q < 0\) and \(p > 0\) for all \((q, p) \in \mathcal {D}\), the solution features both refraction and total internal reflection in two separate quadrants of phase space. The exact solution reads [33]
where \(p_\text {c} = \sqrt{n_1^2 - n_2^2}\), \(\bar{p} = -\mathcal {S}(-p; n_2, n_1, -\mathbf {\nu })\) with \(\mathbf {\nu } = (-1, 0)\), and
The region described by \(\{q < 0, p \ge 0\}\) features propagation through the medium with refractive index \(n_1\). The region \(\{q< 0, -p_\text {c}< p < 0\}\) describes light that was reflected at the optical interface, and the region \(\{q > 0, p \ge 0\}\) describes light that was refracted.
As an initial condition we use
with \(\varphi _m\) defined in (95) and on the part of the boundary of the domain that is not on the optical interface, we prescribe \(\rho = 0\) whenever the velocity field is pointing into the domain, otherwise we leave \(\rho \) free. Since the q-position is restricted to \(q \in [-1, 1]\), this means that at \(q = \pm 1\) we place virtual detectors that capture any luminous flux leaving the system. For the parameters in (102), we take \(q_0 = -0.35\), \(\sigma _q = 0.25\), \(p_0 = 0.45\), \(\sigma _{p,0} = 0.45\), \(p_1 = \tfrac{1}{2} (1.3 + p_\text {c} )\) and \(\sigma _{p,1} = 1.3 - p_1\). Furthermore, we take \(m = 7\) unless specified otherwise.
Again, the explicit 4th order RK method from the previous example is used with a constant \(\varDelta z\)-step as determined by the stability of the temporal integration. The numerical solution is integrated from \(z = 0\) to \(z = Z = 0.7\) and \(z = 2Z\), and is shown in Fig. 8, together with the used initial condition. The result was obtained using a 6th degree polynomial (\(N = 6\)) and \(K = 480\) elements. The mesh uses only rectangular elements and is almost uniform. To easily treat the optical interface, we have shifted the elements below and above the critical momentum \(p_\text {c} \approx 0.98\) in the p-direction such that the critical momentum is aligned with the edges of these elements. The mesh spacings for \(K = 480\) are \(\varDelta q = 0.1\) and \(\varDelta p \approx 0.1\).
In Fig. 8 the quadrants featuring reflection and refraction can be clearly distinguished, while the solution is, as expected, perfectly discontinuous along the optical interface. Furthermore, at \(z = 2Z\) some of the solution has passed \(q = 1\), meaning some energy has hit the detectors. We observe that a total 7.5 % of the initial luminous flux has hit the detectors at \(z = 2Z\). Taking into account the luminous flux on the detectors, we compute the relative error in the total luminous flux as a function of z which is plotted in Fig. 9a. The plot shows that the method obeys energy conservation up to machine precision.
Bucket of water: basic luminance distributions \(\rho (z, q, p)\). Parameters are \(N = 6\), \(K = 480\), \(Z = 0.7\)
Furthermore, to show that the optical interface treatment does not incur any penalty on the convergence order, we compute the discretisation error for this example defined in (98). The convergence data for \(N = 1, 2, ..., 6\) is shown in Table 2. Also for this example, we observe that the spatial discretisation error is dominant and choosing smaller \(\varDelta z\)-steps did not result in different discretisation errors. Moreover, the expected spatial order of convergence \(N + 1\) is obtained.
Next, we verify the exponential convergence of DGSEM by increasing the polynomial degree, whilst keeping the number of elements fixed to \(K = 1920\) and choosing \(m = 28\) in (102). For temporal integration a fixed number of \(2 \cdot 10^4\) z-steps are performed, chosen such that the temporal integration error does not interfere with the convergence test. The result is shown in Fig. 9b and exponential convergence is observed.
Bucket of water
5.2.1 Comparison with Ray Tracing
We compare forward ray tracing with bin-counting to solving Liouville’s equation using the DGSEM. Solving Liouville’s equation already has two advantages, i.e., it conserves energy and provides a more complete picture due to solving the luminance instead of its integrated quantities, the illuminance or luminous intensity. The latter advantage also comes at a price of having to solve a two-dimensional problem in phase space followed by integration to compute these quantities. Ray tracing on the other hand can directly use bins on a one-dimensional grid to compute either the illuminance or luminous intensity.
For a fair comparison, we compute the illuminance E, defined by (17), for this test case using both quasi-Monte Carlo ray tracing and the DGSEM. For quasi-Monte Carlo ray tracing we fix the number of bins to \(B = 1000\) and employ a uniform grid on \(q \in [-1, 1]\), i.e.,
with \(\varDelta q = \frac{2}{B}\). The jth bin is defined by \([Q_j, Q_{j+1}]\) with midpoint \(q_j = \tfrac{1}{2}(Q_j + Q_{j+1})\). The global error for quasi-Monte Carlo integration using a 2D Sobol sequence behaves as \(\mathcal {O}( \log (M)^2 / M )\) with M the number of 2D points [11]. The 2D points are in our case the initial phase space coordinates \((q_i, p_i) \in \mathcal {P}\) of each ray. For more details on quasi-Monte Carlo integration, see [24]. In the bucket of water example \(M = N_\text {RT}\) denotes the number of rays and we use a fixed number of bins.
For the DGSEM we compute the luminance followed by integration such that we obtain the illuminance. Ray tracing defines an average illuminance on each bin, hence, for a fair comparison we also average the illuminance for the DGSEM when computing the discretisation error. For the discretisation error we take the \(L^1\)-norm and compare the numerical solution to the exact illuminance, which is computed by integrating the exact luminance (101) numerically up to machine precision.
Once again we take the initial condition (102) and (95) with \(m = 7\). The illuminance computed using ray tracing with \(N_\text {RT} = 0.64 \cdot 10^6\) rays and the illuminance obtained with DGSEM on a mesh with \(K = 480\) elements and \(N = 4\) are shown in Fig. 10a, together with the exact solution. The ray tracing (RT) solution is noisy, which is inherent in the method due to the quasi-random Monte Carlo process, while the DGSEM solution is almost indistinguishable from the exact solution.
The discretisation error for ray tracing for increasing number of rays is shown in Table 3, while the results for the DGSEM with increasing number of elements K is shown in Table 4. In the tables \(e_\text {RT}\) and \(e_\text {DG}\) denote the errors for ray tracing and solving Liouville’s equation using the DGSEM, respectively, while \(t_\text {RT}\) and \(t_\text {DG}\) denote their respective computation times using only a single core. Furthermore, \(\gamma _\text {RT}\) is estimated from the empirical relation
while \(\gamma _\text {DG}\) is estimated from the empirical relation (99).
From the tables we observe that ray tracing uses \(2.62 \cdot 10^9\) rays taking almost an hour and a half, while the DGSEM achieves roughly the same accuracy in only 8.0 seconds when using 1920 elements. Varying the polynomial degrees results in the performance graph shown in Fig. 10b. It can be observed that the DGSEM always achieves a better accuracy for \(N \ge 1\), compared to quasi-Monte Carlo ray tracing in the same amount of time. The DGSEM significantly outperforms ray tracing and, moreover, can achieve high accuracies in reasonable time.
Left: the illuminance computed using ray tracing (RT) with \(N_\text {RT} = 0.64 \cdot 10^6\) rays and the DGSEM with \(K = 480\) and \(N = 4\) (DG). Right: the error as a function of the computation time for both methods. Both results were computed at \(z = Z = 0.7\)
6 Conclusions
We have solved Liouville’s equation for geometrical optics, using the discontinuous Galerkin spectral element method. For smooth refractive index fields the scheme obeys energy conservation by design. At optical interfaces Snell’s law of refraction or the law of specular reflection needs to be applied. Together with the basic luminance invariance, this corresponds to non-local boundary conditions in phase space. A method was presented to treat the non-local boundary conditions along the optical interface such that the scheme remains energy conservative in the presence of optical interfaces.
Energy conservation is verified in an example. Moreover, in the same example the scheme was compared with forward ray tracing when computing the illuminance. Ray tracing uses bins on a one-dimensional grid to compute the illuminance, while the DGSEM has to solve a two-dimensional problem followed by integration. This still resulted in a better performance compared to ray tracing. In particular, for a fourth degree polynomial, the DGSEM has a computation time of 8.0 seconds, while ray tracing took 1 hour and 26 minutes to achieve almost the same accuracy.
At the moment Fresnel reflections [15, 16] were not included. Therefore, an obvious next step will be to include this in the method by modifying the basic luminance invariance over an optical interface (53), see [26]. Moreover, only two-dimensional optics was considered in this paper. Hence, for future research we intend to extend the method to a three-dimensional optics settings. This requires a four-dimensional phase space together with the propagation along the z-coordinate, making it a five-dimensional problem. Despite the increased computational complexity due to the high dimensionality of the problem, the high convergence rates of DGSEM will still be an advantage over ray tracing, where in theory the DGSEM might achieve a better performance [31].
Availability of Data
Data will be made available on reasonable request.
Change history
24 January 2022
A Correction to this paper has been published: https://doi.org/10.1007/s10915-021-01754-y
References
Arnold, V.I.: Mathematical Methods of Classical Mechanics, vol. 60. Springer Science & Business Media, Berlin (2013)
Bahrami, M., Goncharov, A.V.: Geometry-invariant GRIN lens: finite ray tracing. Opt. Express 22(23), 27797–27810 (2014)
Bernardi, C., Maday, Y., Landriani, G.S.: Noncorming matching conditions for coupling spectral and finite element methods. Appl. Numer. Math. 6(1–2), 65–84 (1989)
Bernardi, C., Maday, Y., Patera, A. T.: Domain Decomposition by the Mortar Element Method. In: Asymptotic and Numerical Methods for Partial Differential Equations with Critical Parameters, pages 269–286. Springer (1993)
Bui-Thanh, T., Ghattas, O.: Analysis of an hp-nonconforming discontinuous Galerkin spectral element method for wave propagation. SIAM J. Numer. Anal. 50(3), 1801–1826 (2012)
Carpenter, M. H., Kennedy, C. A.: Fourth-order 2N-storage Runge-Kutta schemes. NASA TM 109112 (1994)
Chalmers, N., Krivodonova, L.: A robust CFL condition for the discontinuous Galerkin method on triangular meshes. J. Comput. Phys. 403, 10905 (2020)
Chan, J., Hewett, R.J., Warburton, T.: Weight-adjusted discontinuous Galerkin methods: curvilinear meshes. SIAM J. Sci. Comput. 39(6), A2395–A2421 (2017)
Chaves, J.: Introduction to Nonimaging Optics. CRC Press, Boca Raton (2017)
Cvetkovic, A., Dross, O., Chaves, J., Benitez, P., Miñano, J.C., Mohedano, R.: Etendue-preserving mixing and projection optics for high-luminance LEDs, applied to automotive headlamps. Opt. Express 14(26), 13014 (2006)
Filosa, C.: Phase Space Ray Tracing for Illumination Optics. PhD thesis, Eindhoven University of Technology (2018)
Filosa, C., ten Thije Boonkkamp, J.H.M., Ijzerman, W.L.: Ray tracing method in phase space for two-dimensional optical systems. Appl. Op. 55(13), 3599–3606 (2016)
Filosa, C., ten Thije Boonkkamp, J.H.M., Ijzerman, W.L.: Phase space ray tracing for a two-dimensional parabolic reflector. Math. Stat. 5(4), 135–142 (2017)
Glassner, A.S.: An Introduction to Ray Tracing. Elsevier, Amsterdam (1989)
Griffiths, D. J.: Introduction to Electrodynamics (2005)
Hecht, E., et al.: Optics, vol. 4. Addison Wesley, San Francisco (2002)
Herkommer, A.M.: Phase space optics: an alternate approach to freeform optical systems. Opt. Eng. 53(3), 031304 (2013)
Hesthaven, J.S., Warburton, T.: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer Science & Business Media, Berlin (2007)
Ketcheson, D.I.: Runge-Kutta methods with minimum storage implementations. J. Comput. Phys. 229(5), 1763–1773 (2010)
Kopriva, D.A.: A conservative staggered-grid Chebyshev multidomain method for compressible flows. II. A semi-structured method. J. Comput. Phys. 128(2), 475–488 (1996)
Kopriva, D.A.: Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers. Springer Science & Business Media, Berlin (2009)
Kopriva, D.A., Gassner, G.J.: Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable. Appl. Math. Comput. 272, 274–290 (2016)
Kopriva, D.A., Woodruff, S.L., Hussaini, M.Y.: Computation of electromagnetic scattering with a non-conforming discontinuous spectral element method. Int. J. Numer. Meth. Eng. 53(1), 105–122 (2002)
Leobacher, G., Pillichshammer, F.: Introduction to Quasi-Monte Carlo Integration and Applications. Springer, Berlin (2014)
McCluney, W.R.: Introduction to Radiometry and Photometry. Artech House, New York (2014)
Nicodemus, F.E.: Radiance. Am. J. Phys. 31(5), 368–377 (1963)
Ohno, H.: Symplectic ray tracing based on Hamiltonian optics in gradient-index media. JOSA A 37(3), 411–416 (2020)
Pelka, D.G., Patel, K.: An overview of LED applications for general illumination. Des. Eff. Illum. Syst. 5186, 15–26 (2003)
Rausch, D., Rommel, M., Herkommer, A.M., Talpur, T.: Illumination design for extended sources based on phase space mapping. Op. Eng. 56(6), 065103 (2017)
Toulorge, T., Desmet, W.: CFL conditions for Runge-Kutta discontinuous Galerkin methods on triangular grids. J. Comput. Phys. 230(12), 4657–4678 (2011)
van Lith, B. S.: Principles of Computational Illumination Optics. PhD thesis, Eindhoven University of Technology (2017)
van Lith, B.S., ten Thije Boonkkamp, J.H.M., IJzerman, W.L.: Active flux schemes on moving meshes with applications to geometric optics. J. Comput. Phys. X 3, 100030 (2019)
van Lith, B.S., ten Thije Boonkkamp, J.H.M., IJzerman, W.L., Tukker, T.W.: A novel scheme for Liouvilles equation with a discontinuous Hamiltonian and applications to geometrical optics. J. Sci. Comput. 68(2), 739–771 (2016)
Williamson, J.: Low-storage Runge-Kutta schemes. J. Comput. Phys. 35(1), 48–56 (1980)
Wolf, K.B.: Geometric Optics on Phase Space. Springer Science & Business Media, Berlin (2004)
Zhu, X., Zhu, Q., Wu, H., Chen, C.: Optical design of LED-based automotive headlamps. Op. Laser Technol. 45, 262–266 (2013)
Zingg, D.W., Chisholm, T.T.: Runge-Kutta methods for linear ordinary differential equations. Appl. Numer. Math. 31(2), 227–238 (1999)
Acknowledgements
This work is part of the research programme NWO-TTW Perspectief with project number P15-36, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO).
Funding
This work is funded as a part of the research programme NWO-TTW Perspectief with project number P15-36, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO).
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Code Availability
Custom code was written in MATLAB and will not be made available online.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The original article was revised: Modifications have been made to some equations and table 2.
Appendix
Appendix
The matrix \(\varvec{A} \in \mathbb {R}^{(N+2) \times (N+2)}\) defined in equation (78) has a certain structure that allows for analytical computation of its determinant and inverse. Omitting the superscripts from (78) the matrix \(\varvec{A}\) reads
which can be rewritten as
with \(\varvec{Q}\) defined by
The determinant of \(\varvec{A}\) is equal to product of the determinant of the diagonal matrix and the determinant of \(\varvec{Q}\). The determinant of \(\varvec{Q}\) can be found by Laplace (cofactor) expansion along the first row, i.e.,
The second term on the RHS can be easily evaluated using a cofactor expansion along the first column, since it has all zeros except for \(\alpha _0 w_0\) and the remaining minor is the determinant of an identity matrix. Therefore, we obtain
The first term on the RHS of (107) can again be expanded along the first row, resulting in
Repeating these steps we obtain the expression for r, defined in (80), i.e.,
so that
The inverse of \(\varvec{A}\) can readily be found if the inverse of \(\varvec{Q}\) is known. The derivation of the inverse of \(\varvec{Q}\) is briefly outlined for a \(3 \times 3\) matrix, as it can easily be extended to the generic case. We start with the augmented matrix
First, we subtract \(\alpha _i w_i\) times the \((i+1)\)th row from the last row, resulting in
Next, we subtract \(\alpha _i / r\) times the \((i+1)\)th row from the first two rows, and divide the last row by r, such that we obtain
The inverse of \(\varvec{Q}\) now allows for easy evaluation of the inverse of \(\varvec{A} \in \mathbb {R}^{3 \times 3}\), by taking inverse of \(\varvec{A}\) in (106a) resulting in
Similarly, for the general case we obtain
with r defined in (109) and the coefficients of the matrix \(\varvec{B} = (B_{ij})\) read
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
van Gestel, R.A.M., Anthonissen, M.J.H., Boonkkamp, J.H.M.t.T. et al. An Energy Conservative hp-method for Liouville’s Equation of Geometrical Optics. J Sci Comput 89, 27 (2021). https://doi.org/10.1007/s10915-021-01612-x
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10915-021-01612-x