Artificial viscosity in comoving curvilinear coordinates: towards a differential geometrically consistent implicit advection scheme

Abstract

We propose a modification for the tensor of artificial viscosity employable for generally comoving, curvilinear grids. We present a strong conservation form for the equations of radiation hydrodynamics for studying nonlinear pulsations of stars. However, the modification we propose is of general mathematical nature. We study a differential geometrically consistent artificial viscosity analytically and visualize a comparison of our approach to previous implementations by applying it to a simple self-similar velocity field which has a direct application in stars as the fundamental mode of pulsation is radial. We first give a general introduction to artificial viscosity and motivate its application in numerical computations. We then show how a tensor of artificial viscosity has to be designed when going beyond common static Eulerian or Lagrangian comoving rectangular grids. We derive and state the modified equations which include metrical terms that adjust the isotropic (pressure) part of the tensor of artificial viscosity.

1 Background

In astrophysics a multitude of systems and configurations are described with concepts from hydrodynamics, often combined with gravitation, radiation and/or magnetism. Mathematically radiation hydrodynamics (RHD) and magnetohydrodynamics (MHD) are described by systems of coupled nonlinear partial differential equations. The Euler equations of hydrodynamics, the Maxwell equations as well as radiative transport equations are hyperbolic PDEs that connect certain densities and fluxes via conservation laws. The numerical solutions of these equations essentially need to comprise this quality. Today there exists a wide range of numerical schemes for conservation laws that ensure the conservation of mass, momentum, energy etc. if applied properly. Multiple fields in physics and astrophysics have adopted these sophisticated numerical methods for studying various applications.

Standard numerical methods for partial differential equations are established under the assumption of classical differentiability. Routine finite difference schemes of first order usually smear or smoothen the solution in the vicinity of discontinuities as they come with intrinsic numerical viscosity. Standard second-order methods often suffer from the Gibbs phenomenon, where oscillations around shocks emerge. In the past decades so-called high-resolution methods have been developed in order to achieve proper accuracy and resolution for nonlinear, discontinuous problems as they appear also in RHD or MHD. High order explicit Gudonov schemes have been dominating numerical applications and literature for the past decades. This trend was amplified by the enormous advancements made in parallel computing over the past decades.

In this work we consider physical configurations and problems that demand implicit advection schemes due to stability requirements. The parallelization of implicit nonlinear advection schemes, however, is still merely partially possible (solvers such as GMRES (Griewank and Walther [2008]) solve the linear sub-problems in parallel but iteratively in a sequential fashion for each time step). Hence it is desired to minimize the number of grid cells particularly for implicit schemes, where usually the inversion of a nonlinear matrix consumes a major part of the computational power needed. We suggest the adoption of problem oriented grids in combination with a modified artificial viscosity (which we will motivate in Section 1.1) for curvilinear coordinates. In higher-dimensional problems this artificial viscosity emerges as a tensorial quantity, which we demonstrate in Section 1.3. The result present in this paper can be seen as a tensor analytical consequence of the artificial viscosity in general curvilinear coordinates when using consistent metric tensors. In Section 2 we propose a correction for the commonly used tensor of artificial viscosity for curvilinear grids.

This correction is motivated by astrophysical applications where we consider time-dependent comoving nonlinear coordinates represented by non-conformal (non-angle preserving) maps from spherical coordinates. The authors are currently investigating the generation of grids that are asymptotically spherical but which allow certain asymmetries that can be found in rotating configurations as well as nonlinear pulsation processes in stars. This new approach to grid-based astrophysical simulation techniques will be addressed extensively with numerical applications in a future paper.

As an example of non-conformal two-dimensional coordinates, Figure 1 shows a grid that corresponds to the map $(x,y)\to (\xi ,\eta )$,

$$\begin{array}{rl}x=& \xi cos\eta ,\\ y=& (a_1\xi +a_2{\xi }^2)(1+\frac{a_3{\pi }^3-16a_2\xi +a_2{a}_3{\pi }^3\xi }{4\pi (1+a_2\xi )}\eta \\ +\frac{4a_2\xi -a_3{\pi }^3-a_2{a}_3{\pi }^3\xi }{{\pi }^2(1+b_2\xi )}{\eta }^2+a_3{\eta }^3)sin\eta \end{array}$$

(1)

which yields standard orthogonal polar coordinates for the choice of parameters $(a_1,a_2,a_3)=(1,0,0)$. In such a nonorthogonal grid the metric tensor is no longer diagonal and one has to consider a consistent differential geometric approach to the formulation of the governing equations of RHD and MHD, and also to the mathematical formulation of the artificial viscosity, which will be stressed in Section 2.

Figure 1

Example of a non-conformal non-steady 2D grid with oblateness governed by three choices of parameters $\mathbf{(}{\mathit{a}}_{\mathbf{1}}\mathbf{,}{\mathit{a}}_{\mathbf{2}}\mathbf{,}{\mathit{a}}_{\mathbf{3}}\mathbf{)}$ as defined in equation ( 1 ). (a) corresponds to $(a_1,a_2,a_3):(1,0,0)$, (b): $(0.8,0,0)$, (c): $(0.8,0,0.5)$ and (c): $(0.8,0.2,0.5)$.

The benefit of the consistent formulation is especially evident when considering time-dependent grids, e.g. when using time-dependent parameters $(a_1,a_2,a_3)$ in (1). We refer to the Appendix for a depiction of the system of equations of RHD for generally comoving curvilinear coordinates with time-dependent metrics.

1.1 Brief introduction to conservation laws

For the sake of stringency we recapitulate some important results from the theory and numerics of conservation laws and thereby introduce a few mathematical terms needed to motivate artificial viscosity. We refer to LeVeque ([1991]) and Richtmyer and Morton ([1994]) for the complete picture.

The equations of RHD and MHD form a system of hyperbolic conservation laws that describe the interaction of a density function $\mathbf{d}(\mathbf{x},t):{\mathbb{R}}^n\times [0,\mathrm{\infty })\to {\mathbb{R}}^m$ and its flux $\mathbf{f}(\mathbf{d}):{\mathbb{R}}^m\to {\mathbb{R}}^{m\times n}$. Equation (6) shows how a concrete choice for the density and the flux field can look like in a given coordinate system.

The temporal change of the integrated density in a connected set $\mathrm{\Omega }\subset {\mathbb{R}}^n$ then equals the flux over the boundary ∂ Ω, i.e.,

$${\partial }_t{\int }_{\mathrm{\Omega }}\mathbf{d}\mathrm{d}V+{\int }_{\partial \mathrm{\Omega }}\mathbf{f}\cdot \mathbf{n}\mathrm{d}S=\mathbf{0}\text{for all }t>0,$$

(2)

where n is the outward oriented normal of the surface.

The system is called hyperbolic if the Jacobian matrix ${\mathrm{\nabla }}_{\mathbf{d}}\mathbf{f}$ associated with the fluxes has real eigenvalues and if there exists a complete set of eigenvectors. In case of MHD and RHD this property has a direct physical relevance (Pons et al. [2000]).

Assuming f to be a continuously differentiable function, equation (2) can be rewritten via the divergence theorem as

$$\begin{array}{r}{\int }_t{\int }_{\mathrm{\Omega }}({\partial }_t\mathbf{d}+{div}_{\mathbf{x}}\mathbf{f}(\mathbf{d}))\mathrm{d}V\mathrm{d}t\\ =\mathbf{0}\text{for all }t>0,\mathrm{\Omega }\subset {\mathbb{R}}^n,\end{array}$$

(3)

which gives a system of partial differential equations for the density function d:

$${\partial }_t\mathbf{d}+{div}_{\mathbf{x}}\mathbf{f}(\mathbf{d})=\mathbf{0}\text{for all }t>0,\mathbf{x}\in {\mathbb{R}}^n.$$

(4)

With an initial condition $\mathbf{d}(\mathbf{x},0)={\mathbf{d}}_0(\mathbf{x})$, $\mathbf{x}\in {\mathbb{R}}^n$, this is called the Cauchy problem.

In order to illustrate the connection of hydrodynamical applications to this formalism, we express the Euler equations in the form (4). The appearing variables are the gaseous density $\rho (\mathbf{x},t)$, the gas velocity $\mathbf{u}(\mathbf{x},t)$, the inner energy $ϵ(\mathbf{x},t)$ and the gaseous pressure tensor $\mathbf{P}(\mathbf{x},t)$. Considering the differential form (4), we recognize the continuity equation, the equation of motion and the energy equation as the components of the hyperbolic problem. In case of the most relevant problem, that of 3D hydrodynamics, the density and its flux are given as

$$\mathbf{d}=\left(\begin{array}{c}\rho \\ \rho \mathbf{u}\\ \rho ϵ\end{array}\right)\in {\mathbb{R}}^5,\mathbf{f}(\mathbf{d})=\left(\begin{array}{c}\rho {\mathbf{u}}^{\mathrm{T}}\\ \rho \mathbf{u}{\mathbf{u}}^{\mathrm{T}}+\mathbf{P}\\ \rho ϵ\mathbf{u}+{(\mathbf{P}\mathbf{u})}^{\mathrm{T}}\end{array}\right)\in {\mathbb{R}}^{5\times 3}.$$

(5)

For a given coordinate system with base vectors ${\mathbf{e}}_i$, the tensorial fields are given explicitly as (using the Einstein notation)

$$\mathbf{d}=\left(\begin{array}{c}\rho \\ \rho u^i{\mathbf{e}}_i\\ \rho ϵ\end{array}\right),\mathbf{f}(\mathbf{d})=\left(\begin{array}{c}\rho u^i{\mathbf{e}}_i^{\mathrm{T}}\\ (\rho u^i{u}^j+P^{ij}){\mathbf{e}}_i{\mathbf{e}}_j^{\mathrm{T}}\\ (\rho ϵu^i+P_j^i{u}^j){\mathbf{e}}_i^{\mathrm{T}}\end{array}\right).$$

(6)

The gaseous pressure tensor can be assumed to be isotropic in most applications, which means that $\mathbf{P}=g^{ij}p$ where $p(\mathbf{x},t)$ is the scalar gas pressure and $g^{ij}={{\mathbf{e}}^i}^{\mathrm{T}}{\mathbf{e}}^j$ the contravariant metric tensor. In case of adaptive grids the base vectors are time-dependent as well, i.e., ${\mathbf{e}}_i={\mathbf{e}}_i(\mathbf{x},t)$.

Since even the simplest one-dimensional scalar conservation laws like the Burgers’ equation have classical solutions only in some special cases, one has to broaden the considered function space of possible solutions. For the so-called weak solutions, we appeal to generalized functions where the discontinuities are defined properly. The generalized concept of differentiation of distributions shifts the operations to test functions $\gamma :{\mathbb{R}}^n\times {\mathbb{R}}^{+}\supset G\to \mathbb{R}$ (G open) which are infinitely differentiable and have a compact support (meaning that for each γ there exists a closed and bounded subset K such that $\gamma (\mathbf{x},t)=0$ for all $x\in G\setminus K$). We denote this space of test functions by $D(G)$. In this generalized space of solutions the Cauchy problem (3) is written as

$${\int }_{t\ge 0}{\int }_{{\mathbb{R}}^n}({\partial }_t\mathbf{d}+{div}_{\mathbf{x}}\mathbf{f}(\mathbf{u}))\gamma \mathrm{d}V\mathrm{d}t=0\text{for all }\gamma \in D(G).$$

The weak formulation of the conservation law (3) is obtained by shifting the derivatives to the test functions by partial integration, and by using the compactness of the support. We get that the following has to hold for each $\gamma \in D(G)$:

$$\begin{array}{r}{\int }_{t\ge 0}{\int }_{{\mathbb{R}}^n}(\mathbf{d}{\partial }_t\gamma +\mathbf{f}(\mathbf{d}){\mathrm{\nabla }}_{\mathbf{x}}\gamma )\mathrm{d}V\mathrm{d}t\\ =-{\int }_{{\mathbb{R}}^n}\gamma (\mathbf{x},0){\mathbf{d}}_0(\mathbf{x})\mathrm{d}V.\end{array}$$

(7)

The function $\mathbf{d}\in L^{\mathrm{\infty }}$ is called a weak solution of the PDE (4), if it satisfies (7) and $\mathbf{d}\in U$ with ${\mathbf{d}}_0\in L^{\mathrm{\infty }}$. However, there is a small drawback. This weak solution is not necessarily unique and usually further constraints have to be imposed in order to guarantee its uniqueness. This leads us to the actual topic of this paper.

1.2 Introduction to artificial viscosity

For most physical problems it is sufficient to look for weak solutions from the function space of piecewise continuously differentiable functions. Constraining the space of solutions in this way, we call the physical variables d weak solutions of the Cauchy problem (4), if they are classical solutions wherever they are continuously differentiable, and if at discontinuities (shocks) they satisfy additional conditions in order to be physically reasonable (we elaborate on these conditions below).

The mathematical theory provides several techniques to distinguish physically valuable solutions out of a manifold of mathematically possible. One method is to add an artificial viscosity term to the right-hand side of (4), to get the equation:

$${\partial }_t\mathbf{d}+{div}_{\mathbf{x}}\mathbf{f}(\mathbf{d})=\epsilon \nu \mathrm{\Delta }\mathbf{d},\epsilon >0$$

(8)

and then consider the limiting case $\epsilon \to 0$. This idea is motivated by physical diffusion which broadens sincere discontinuities to differentiable steep gradients at the (microscopic) length scale of the mean free path of the particles. The physical solution of the weakly formulated problem is thus the zero diffusion limit of the diffusive problem. However, in practice this limit is difficult to calculate analytically and hence simpler conditions have to be found. A common technique to do this is motivated by continuum physics as well. Here an additional conservation law is set to hold for another quantity - the entropy of the fluid flow - as long as the solution remains smooth. Moreover, it is known that along admissible shocks this physical variable never decreases, and the conservation law for the entropy can be formulated as an inequality.

We denote the (scalar valued) entropy function by $\sigma (\mathbf{d})$ and the entropy flux function by $\varphi (\mathbf{d})$, and they satisfy

$${\partial }_t\sigma (\mathbf{d})+{div}_{\mathbf{x}}\varphi (\mathbf{d})=0.$$

(9)

Assuming the functions to be differentiable, we may rewrite this conservation law via the chain rule and the equation (4) as

$${\mathrm{\nabla }}_{\mathbf{d}}\sigma (\mathbf{d}){div}_{\mathbf{x}}\mathbf{f}(\mathbf{d})={div}_{\mathbf{x}}\varphi (\mathbf{d}),$$

(10)

where in higher-dimensional case the appearing matrices of gradients have to fulfill further constraints, see e.g. Godlewski and Raviart ([1992]). For scalar equations, it is always possible to find an entropy function of this kind. Furthermore it is assumed that the entropy function is convex, i.e.

$${\mathrm{\nabla }}_{\mathbf{d}}^2\sigma >0,\text{for all }\mathbf{d}\in U.$$

(11)

To get our actual entropy condition, we first rewrite our entropical conservation law (9) in the viscous form

$${\partial }_t\sigma (\mathbf{d})+{div}_{\mathbf{x}}\varphi (\mathbf{d})=\epsilon {\mathrm{\nabla }}_{\mathbf{d}}\sigma (\mathbf{d})\mathrm{\Delta }\mathbf{d}.$$

(12)

Integrating over an arbitrary time interval $[t_0,t_1]$ and a connected set $\mathrm{\Omega }\subset {\mathbb{R}}^n$, and using partial integration, we find that

$$\begin{array}{r}{\int }_{t_0}^{t_1}{\int }_{\mathrm{\Omega }}({\partial }_t\sigma (\mathbf{d})+{div}_{\mathbf{x}}\varphi (\mathbf{d}))\mathrm{d}V\mathrm{d}t\\ =\epsilon {\int }_{t_0}^{t_1}{\int }_{\partial \mathrm{\Omega }}({\mathrm{\nabla }}_{\mathbf{x}}\mathbf{d}{\mathrm{\nabla }}_{\mathbf{d}}\sigma (\mathbf{d}))\cdot \mathbf{n}\mathrm{d}S\mathrm{d}t\\ -\epsilon {\int }_{t_0}^{t_1}{\int }_{\mathrm{\Omega }}{\mathrm{\nabla }}_{\mathbf{x},i}\mathbf{d}{\mathrm{\nabla }}^2\sigma (\mathbf{d}){\mathrm{\nabla }}_{\mathbf{x},i}\mathbf{d}\mathrm{d}V\mathrm{d}t,{\mathrm{\nabla }}^2\sigma (\mathbf{d})>0.\end{array}$$

(13)

When we now consider the non-diffusive limit $\epsilon \to 0$, the first term on the right-hand side vanishes without further restriction whereas the second term has to remain non-positive. Then, using partial integration and the divergence theorem we obtain the entropy condition

$$\begin{array}{rl}{\int }_{\mathrm{\Omega }}\sigma (\mathbf{d}(\mathbf{x},t_1))\mathrm{d}V\le & {\int }_{\mathrm{\Omega }}\sigma (\mathbf{d}(\mathbf{x},t_0))\mathrm{d}V\\ -{\int }_{t_0}^{t_1}{\int }_{\partial \mathrm{\Omega }}\varphi (\mathbf{d})\cdot \mathbf{n}\mathrm{d}S\mathrm{d}t.\end{array}$$

(14)

For bounded, continuous pointwise solutions ${\mathbf{d}}^{\ast }$ of (12) such that ${\mathbf{d}}^{\ast }\to \mathbf{d}$ for $\epsilon \to 0$, the vanishing viscosity solution d is a weak solution of the initial value problem (3) and fulfills entropy condition (14). Generally spoken, applying the entropy condition to systems with shock solutions unveils those propagation velocities that ensure that no characteristics rise from discontinuities which would be non-physical. For detailed motivation, stringent argumentation and proofs to mathematical techniques presented in this section we refer to Harten et al. ([1976]) and LeVeque ([1991]).

1.3 Numerical artificial viscosity

As mentioned we are looking for high-resolution methods for nonlinear PDEs derived from hyperbolic conservation laws. In the past decades major efforts have been made in developing numerical methods for these problems that are at least of second order. One patent attempt to finding such a high-resolution method is to adapt a well-known high-order method for linear problems for nonlinear problems (such as the Lax-Wendroff scheme (Lax and Wendroff [1960])).

As illustrated above we can add an artificial viscosity term to the conservation law in a way that the entropy condition is satisfied and non-physical solutions are excluded. The viscosity term has to be designed in such a manner that it affects sincere discontinuities but vanishes sufficiently elsewhere so that the order of accuracy can be maintained in those regimes where the solution is smooth. The idea of numerical artificial viscosity was inspired by physical dissipation mechanisms and dates back more than half a century to von Neumann and Richtmyer ([1950]).

We denote as customary the approximate solution of the exact density $\mathbf{d}(x,t)$ at discrete grid points $\mathbf{d}(x_j,t_n)$ by ${\mathbf{D}}_j^n$, and set $\mathbf{D}={[{\mathbf{D}}_1\cdots {\mathbf{D}}_k]}^{\mathrm{T}}$, where k is the total number of grid points. The numerical representation of the flux function $\mathbf{f}(\mathbf{d})$ is denoted respectively by $\mathbf{F}(\mathbf{D})$, where ${[\mathbf{F}(\mathbf{D})]}_j=\mathbf{f}({\mathbf{D}}_j)$. The numerical flux function gets modified by a an artificial viscosity $\mathbf{Q}{[\mathbf{D}]}_j$ for instance in the following way:

$${[{\mathbf{F}}_{\mathrm{visc}}(\mathbf{D})]}_j={[\mathbf{F}(\mathbf{D})]}_j-h{[\mathbf{Q}(\mathbf{D})]}_j({\mathbf{D}}_{j+1}-{\mathbf{D}}_j),$$

(15)

where h denotes the resolution of the spatial discretization, $h=x_{j+1}-x_j$. Since the original design of the additional viscous pressure in the scalar form, $Q=c\rho {(\mathrm{\Delta }\mathbf{u})}^2$, $c\in \mathbb{R}$, as suggested in von Neumann and Richtmyer ([1950]) for one-dimensional advection ${\partial }_t\mathbf{d}+a{\partial }_x\mathbf{d}=Q{\partial }_{xx}\mathbf{d}$, it has undergone a number of modifications and generalizations. It has turned out to be numerically preferable to add a linear term (see Landshoff ([1955])) in order to control oscillations. Generalizations to multi-dimensional flows mostly retain the original analogy to physical dissipation and reformulate the velocity term accordingly, see e.g. Wilkins ([1980]).

The artificial viscosity broadens shocks to steep gradients at some characteristic length scale, but should not cause too large smearing. The concrete composition and implementation of this artificial viscosity coefficient Q depends on the application. As an example we discuss the following form of the tensor of the artificial viscosity in higher-dimensional RHD numerics. Similar forms of artificial viscosity can be found also in pure hydrodynamics and MHD calculations in 2D and 3D.

Tscharnuter and Winkler ([1979]) pointed out that the viscous pressure in 3D radiation hydrodynamics has to unravel the normal stress, quantified by the divergence of the velocity field and the shear stress, which is expressed by the symmetrized gradient of the velocity field according to the general theory of viscosity. It is designed to switch on only in case of compression (${div}_{\mathbf{x}}\mathbf{u}<0$), and this is all ensured by the form

$$\mathbf{Q}=-q_2^2{l}_{\mathrm{visc}}^2\rho max(-{div}_{\mathbf{x}}\mathbf{u},0)({[\mathrm{\nabla }\mathbf{u}]}_s-\frac{1}{3}\mathbf{e}{div}_{\mathbf{x}}\mathbf{u}),$$

(16)

where the symmetrization rule is defined componentwise for the lower indices as

$${\left({[\mathrm{\nabla }\mathbf{u}]}_s\right)}_{ij}=\frac{1}{2}({\mathrm{\nabla }}_i{u}_j+{\mathrm{\nabla }}_j{u}_i).$$

2 Numerical artificial viscosity in curvilinear coordinates

We introduced artificial viscosity in form of a three dimensional viscous pressure tensor (16) in Section 1.3 along the lines of Tscharnuter and Winkler ([1979]). In this section we want to point out, how such a definition must be adapted for curvilinear coordinates in order to ensure tensor analytical consistency.

When formulating PDEs derived from hyperbolic conservation laws on a curvilinear grid, the tensorial equations (3) have to be transformed to the according coordinate system. Not only the vectorial and tensorial quantities have to be transformed but also the differentiation operators, in particular the divergence operator in our case. The appropriate framework to do this is provided by differential geometry. Like the gradient of a scalar is natively a covector, there are rules for co- and contravariant indices of tensors such as the one we are interested in. The crucial term in (16) is the symmetrized velocity gradient ${[\mathrm{\nabla }\mathbf{u}]}_s$ that accounts for shear stresses, and one sees that the form (16) comes into conflict with the demand of vanishing trace ($Tr\mathbf{Q}=0$) when the divergence term is simply of the form $\mathbf{e}{div}_{\mathbf{x}}\mathbf{u}$, as we find it commonly in several MHD and RHD grid codes.

Proposition

The correct form of the viscous pressure tensor (16) in general coordinates is

$$\mathbf{Q}=-q_2^2{l}_{\mathrm{visc}}^2\rho max(-{div}_{\mathbf{x}}\mathbf{u},0)({[\mathrm{\nabla }\mathbf{u}]}_s-\frac{1}{3}\mathbf{g}{div}_{\mathbf{x}}\mathbf{u}),$$

(17)

where g is the fundamental metric tensor. We show that this tensor has the desired properties. The viscous pressure tensor must be symmetric by definition, i.e. $Q_{ij}=Q_{ji}$, which can be easily verified from (17). Also, the trace of the tensor has to vanish ($Tr\mathbf{Q}=Q_i^i=0$), as pointed out in von Neumann and Richtmyer ([1950]). To show that this holds for (17), we first consider its native covariant components

$$\begin{array}{rl}{Q}_{ij}=& -{\mu }_{2}max(-{div}_{\mathbf{x}}\mathbf{u},0)\\ \times (\frac{1}{2}({\mathrm{\nabla }}_i{u}_j+{\mathrm{\nabla }}_j{u}_i)-\frac{1}{3}{g}_{ij}{div}_{\mathbf{x}}\mathbf{u}),\end{array}$$

(18)

where we have renamed $q_2^2{l}_{\mathrm{visc}}^2\rho ={\mu }_2$. Next, we need to raise an index with the metric,

$$\begin{array}{rl}{Q}_j^i=& Q_{lj}{g}^{li}\\ =& -{\mu }_{2}max(-{div}_{\mathbf{x}}\mathbf{u},0)\\ \times (\frac{1}{2}{g}^{li}({\mathrm{\nabla }}_l{u}_j+{\mathrm{\nabla }}_j{u}_l)-\frac{1}{3}{\delta }_j^i{div}_{\mathbf{x}}\mathbf{u}),\end{array}$$

and use the essential identity $g^{li}{g}_{lj}=g_j^i={\delta }_j^i$. The Ricci lemma ${\mathrm{\nabla }}_i{g}_{jk}={\partial }_i{g}_{jk}-{\mathrm{\Gamma }}_{ij}^l{g}_{lk}-{\mathrm{\Gamma }}_{ik}^l{g}_{jl}=0$ for the fundamental tensor naturally also holds for the contravariant components and we can permute g into the derivatives ${\mathrm{\nabla }}_l{g}^{li}{u}_j$ and ${\mathrm{\nabla }}_j{g}^{li}{u}_l$ which yields two times the divergence ${\mathrm{\nabla }}_i{u}^i={div}_{\mathbf{x}}\mathbf{u}$ when we conduct the contraction $j\to i$. We use the summation ${\delta }_i^i=3$ and obtain the desired result

$$Q_i^i=\cdots =(\frac{1}{2}(2{div}_{\mathbf{x}}\mathbf{u})-\frac{1}{3}3{div}_{\mathbf{x}}\mathbf{u})=0.$$

The commonly used (see e.g. Dorfi ([1999]) in RHD, Iwakami et al. ([2008]) in MHD, Fryxell et al. ([2000]) in MHD) form of Q (16) is not compatible with these requirements since the symmetrization is only defined for lower indices, whereas the unit tensor e of a metric space is only defined for mixed indices, meaning there is no such thing as ${\delta }_{ij}$. However, the above mentioned and other authors have neglected this inconsistency since they have considered mixed indices from the start when concentrating on Cartesian or affine coordinates. Nonlinear corrections have been suggested in Benson and Schoenfeld ([1993]) albeit this approach does not explicitly concern curvilinear coordinates and is based on a TVD approach.

In several hydro- and MHD-codes that include non-Cartesian grids such as Pluto (Mignone et al. [2007]), the geometric source terms are coded explicitly for several geometries (polar, cylindrical, spherical), and not only for the artificial viscosity flux. The suggestions made e.g. in Vinokur ([1974]) lead to geometrical source terms that correct curvilinear grid effects. However the strong conservation form as elaborated in Warsi ([1981]) would need to appeal to our differential geometrically consistent approach in order to deal with the viscosity in an intrinsically consistent way. Especially when the metric tensor itself is not only a function of space but also time-dependent (as discussed in Section 1), the latter approach reaches its limits. Our correction affects curvilinear coordinates in multiple dimensions, whereas it is not necessary that the coordinates are orthogonal, i.e., the metric tensor does not need to be diagonal. Our initial motivation to study more general coordinates comes from the idea to generate problem-oriented coordinate systems for astrophysical numerical calculations. In a following paper we want to present some feasible approaches to grid generation under certain physical restrictions. Such nonlinear grids that are adaptive in multiple dimensions have time-dependent metric tensors and thus benefit directly from our consistent definition. On the contrast, when using adaptive mesh refinement, the metric tensor remains geometrically constant in time.

In order to support the theoretical results in this work, in the upcoming section we study as an example a simple velocity field with a non-vanishing divergence and visualize the according artificial viscosities for the two presented cases.

3 Application and visualization

The most common application of curvilinear coordinates in 3D is the map $(x,y,z)\to (r\in {\mathbb{R}}^{+},\vartheta \in [0,\pi ],\phi \in [0,2\pi ])$ with $x=rsin\vartheta cos\phi$, $y=rsin\vartheta sin\phi$ and $z=rcos\vartheta$ as spherical coordinates. The corresponding diagonal covariant metric tensor in this simple orthogonal case is given by $\mathrm{diag}(1,r,rsin\vartheta )$.

3.1 Toy model velocity field

As the presented considerations for artificial viscosity on non-steady curvilinear coordinates stem from astrophysical applications, we want to exemplarily address a velocity field with a certain practice in RHD. We study a toy model of a self-similar fluid flow solution, namely the velocity field given by

$${\mathbf{u}}_{\text{Ex}}=\frac{\mathbf{x}}{\sqrt{x^2+y^2+z^2}}=\frac{\mathbf{x}}{r}$$

(19)

here in Cartesian coordinates. Such self-similar solutions appear in idealized spherical models of stars, for example as shocks driven by radial stellar pulsations.

This vector field is obviously symmetric with respect to the origin and has a non-vanishing divergence ${div}_{\mathbf{x}}{\mathbf{u}}_{\text{Ex}}=2/r$. The covariant components of this vector field are given in any other coordinates by scalar product with the base vectors, i.e. $u_{\text{Ex},i}={\mathbf{u}}_{\text{Ex}}\cdot {\mathbf{e}}_i$. This leads to the covariant components $(1,0,0)$ in spherical coordinates. The nonzero covariant components of the tensorial part $({[\mathrm{\nabla }\mathbf{u}]}_s-\frac{1}{3}\mathbf{e}{div}_{\mathbf{x}}\mathbf{u})$ of the artificial viscosity (17) are given for this field by

$$Q_{rr}=-\frac{2}{3r},Q_{\theta \theta }=\frac{r}{3},Q_{\varphi \varphi }=\frac{r{sin}^2\theta }{3}.$$

(20)

In the following section we visualize the tensor of artificial viscosity (TAV) for the velocity field (19). A uniform distribution of the leading eigenvalues over the whole domain is expected due to the symmetry of the vector field.

One easily verifies the identity $Q_i^i=0$ summing over the mixed components $Q_r^r=-2/3r$, $Q_{\theta }^{\theta }=1/3r$, $Q_{\varphi }^{\varphi }=1/3r$.

With the previous version of the TAV from Tscharnuter and Winkler ([1979]) we would get the following covariant components of the artificial viscosity tensor:

$$\begin{array}{r}{Q}_{rr}=-\frac{2}{3r},Q_{\theta \theta }=-\frac{2}{3r}+r,\\ Q_{\varphi \varphi }=-\frac{2}{3r}+r{sin}^2\theta .\end{array}$$

(21)

The visualization of this non-metric version of the artificial viscosity for the symmetric velocity field (19) shows obviously a field unequal in strength and direction over the whole domain. In numerical computations this would clearly lead to artificial anisotropies in the flux of the density field and destroy all efforts in constructing a higher-order conservative numerical scheme with artificial viscosity.

3.2 Visualization of scalar and tensor fields

We can see the incorrect vs. correct behavior immediately by even just displaying the major eigenvalue or the trace of the viscosity tensor. However, since the major eigenvalue represents only one degree of freedom out of the six available in the tensor field, a measure involving all six components is a more objective metric for validation.

We used the Vish Visualization Shell (Benger et al. [2007]) to numerically sample (20) and (21) on a uniform grid for analyzing scalar fields (Figures 2 and 3) and a radial sampling distribution (Figure 4) for the full tensor field.

Figure 2

Major eigenvalue the incorrect (a) and correct (b) viscosity, shown along the XZ plane.

Figure 3

Trace the incorrect (a) and correct (b) viscosity, shown along the XZ plane. Note the range of values in Figure 3(b) - what we see is just numerical noise, the trace is perfectly zero within numerical precision.

Figure 4

Reynold glyphs (Moore et al.[1995]) of the incorrect and correct viscosity tensor. The incorrect tensor shows a strongly negative component, whereas the correct tensor is balanced, indicating zero trace.

Figure 2 displays a structure of the eigenvalue corresponding to the major eigenvector of the TAV on the XZ plane (i.e., in the plane $y=0$), evidently showing some asymmetric coordinate-dependent behavior of the incorrect TAV, whereas the correct TAV is radially symmetric, as desirable. The tensor field is symmetric and of rank two, however it is not positive definite and may exhibit vanishing trace. The correct TAV is trace-free in the entire domain, whereas the trace of the incorrect TAV ranges through positive values (for large radial distances) to large negative values close to the coordinate origin, as depicted in Figure 3. A direct visualization method depicting the full six degrees of freedom of the tensor field is thus favorable, or even rather essential. However, many direct visualization methods for tensor fields require positive definiteness and are thus not applicable to this data. Only very few methods are suitable for general tensors. Figure 4 shows so-called Reynold glyphs (Moore et al. [1995]) for the TAV field. A Reynold glyph is the surface generated by mapping a tangential vector $v(P)$ at each data sample point P as

$$P\to Q(v(P),v(P)).$$

Such glyphs shown at each sampling point provide a direct visualization of the full six degrees of freedom of the tensor properties. Reynold glyphs are also able to depict negative definite tensors, whereas a quadric surface (ellipsoids representing $P\to 1/\sqrt{Q(v(P),v(P))}$ becomes hyperbolic for negative eigenvalues and problematic for visualization purposes. The Reynold glyph directly shows the ‘directional value’ $Q(v(P),v(P))$ of the tensor field Q in the direction v around a the sampling point P - the resulting surface is intersecting the sampling point P whenever $Q(v(P),v(P))=0$, which is the case for points where the tensor is degenerating and not positive definite. In such areas the glyph will visually appear like two intersecting surfaces corresponding to the isopotential surfaces of second order spherical harmonic functions. These both surface components represent positive and negative eigenvalues of the tensor field - if both positive and negative component ‘counter-balance’ themselves they therefore indicate vanishing trace, which is the sum of the eigenvalues. If only one surface component is visible, then the tensor is either positive or negative definite on that certain point.

We used a radial sampling for the direct visualization of the tensor field in order to minimize coordinate artifacts. As depicted in Figure 4(b) and Figure 5(b) for the correct TAV all ‘modes’ are equivalently represented, indicating vanishing trace of the tensor field while being radially aligned with the underlying coordinate system. In contrast, the incorrect TAV, Figure 4(a) and Figure 5(a), exhibits a dominantly negative trace.

Figure 5

Detailed Reynold glyphs (Moore et al.[1995]) of the incorrect and correct viscosity tensor. Glyphs of the correct tensor show spherical harmonics that are balanced in their positive and negative half-widths.

4 Conclusions

We studied a generalization of the tensor of numerical artificial viscosity for curvilinear coordinates and compared our result to previous definitions found in literature. We analyzed a symmetric toy velocity field and visualized its viscosities. Clearly, the non-metric version of the TAV as used by many authors of hydro- (HD) respectively magneto-hydro- (MHD) or radiation hydrodynamic-codes (RHD) leads to incorrect results in curvilinear coordinates, whereas our suggestion for the numerical artificial viscosity gives geometrically consistent results.

Appendix: Strong conservation form

As mentioned in Section 1, the benefit of the strong formulation of Cauchy problems arising from physical conservation laws will be discussed in more detail here.

Following the ideas of Warsi ([1981]), in non-steady coordinates the geometrically consistent strong form of tensorial conservation laws (2) is given as

$${\partial }_t\left(\sqrt{|\mathbf{g}|}\mathbf{d}\right)+{\partial }_i(\sqrt{|\mathbf{g}|}\mathbf{f}(\mathbf{d})\cdot {\mathbf{e}}^i)=0,$$

(22)

where d and f are decomposed according to their native tensorial components. The scalar multiplication with the i th contravariant base vectors ${\mathbf{e}}^i$ yields a projection on the contravariant coordinate lines which in case of generalized grids can differ in direction and length of their covariant counterparts. Equation (22) gives also the integral form of the conservation law, which should be treated numerically in a correct way for non-steady coordinates in any finite-volume discretization. The important difference to the componentwise structure $\mathrm{\nabla }(\cdot )=\partial (\cdot )+\mathrm{\Gamma }(\cdot )$, where Christoffel symbols account for the geometry is that in this case undifferentiated terms arise (see Vinokur ([1974])) which act like geometric sources in the equations and destroy conservativeness. A comprehensive proof of Vinokurs theorem using differential forms can be found in Bridges ([2008]).

With non-steady curvilinear grids not only the nonlinearity of the metric tensor but also its time-dependence has to be taken into account numerically. The motion of the grid itself and its implications on the formulation of the set of equations respectively the calculation of the occurring fluxes is discussed in the following section.

A.1 Adaptive grids

In fluid dynamics we distinguish two main reference systems that suit unequally for various applications. The Eulerian frame is the fixed reference system of an external observer in which the fluid moves with velocity u whereas the Lagrangian approach describes the physics in the rest frame of the fluid. Between these two systems, the transformation of an advection term for a density d (that moves with a relative velocity u) is given via the material derivative $D_t\mathbf{d}={\partial }_t\mathbf{d}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{d}$.

Hence, when we work with comoving frames, the coordinate system respectively the computational grid is time-dependent. There is a number of purposes where strict Eulerian or Lagrangian grids are suboptimal and thus we need to consider the generalized the concept of the comoving frames.

We next state the strong conservation form for time-dependent general coordinate systems. The time derivative of a density d in the coordinate system ${\mathrm{\Sigma }}_{(\beta )}$ relative to a (e.g. static) coordinate system ${\mathrm{\Sigma }}_{(\alpha )}$ is given by ${\partial }_t{\mathbf{d}}_{(\beta )}={\partial }_t{\mathbf{d}}_{(\alpha )}+{\mathrm{\nabla }}_{(\alpha )}\mathbf{d}{\partial }_t{\mathbf{x}}_{(\beta )}$ and from the view point of system ${\mathrm{\Sigma }}_{(\alpha )}$, the time derivative is given as

$${\partial }_t{\mathbf{d}}_{(\alpha )}={\partial }_t{\mathbf{d}}_{(\beta )}-\dot{\mathbf{x}}\cdot {\mathrm{\nabla }}_{(\alpha )}{\mathbf{d}}^{\mathrm{T}},$$

(23)

where $\dot{\mathbf{x}}$ denotes the grid velocity. The second term on the right side we call grid advection, see e.g. Warsi ([1981]) and Thompson et al. ([1985]). An inhomogeneous advective term of a conservation law $\mathbf{K}={\mathbf{d}}_t+div({\mathbf{u}}^{\mathrm{T}}\mathbf{d})$ (in a fixed coordinate system) is given in the case of moving grid as

$$\mathbf{K}={\mathbf{d}}_t-\dot{\mathbf{x}}\cdot \mathrm{\nabla }{\mathbf{d}}^{\mathrm{T}}+div\left({\mathbf{u}}^{\mathrm{T}}\mathbf{d}\right).$$

(24)

When we apply the appropriate transformation prescriptions to the spatial derivatives, we obtain the following form

$$\mathbf{K}={\mathbf{d}}_t-\dot{\mathbf{x}}\cdot \frac{1}{\sqrt{|\mathbf{g}|}}{\partial }_i\left(\sqrt{|\mathbf{g}|}{\mathbf{e}}^i\mathbf{d}\right)+\frac{1}{\sqrt{|\mathbf{g}|}}{\partial }_i(\sqrt{|\mathbf{g}|}\mathbf{u}\cdot {\mathbf{e}}^i\mathbf{d}).$$

(25)

This is not yet geometrically conservative, since it is not of an integral structure. Following the idea of the Reynolds transport theorem we consider the temporal derivative of the volume, i.e., the determinant of the time-dependent metric tensor $\sqrt{|\mathbf{g}(\mathbf{x},\mathbf{t})|}$ in order to study the conservation of a density function in variable volumes and obtain the strong conservation form for time-dependent coordinate systems:

$$\sqrt{|\mathbf{g}|}\mathbf{K}={\partial }_t\left(\sqrt{|\mathbf{g}|}\mathbf{d}\right)+{\partial }_i\left(\sqrt{|\mathbf{g}|}{U}^i\mathbf{d}\right),$$

(26)

where $U^i$ denotes the contravariant velocity components relative to the moving grid, $U^i={\mathbf{e}}^i\cdot (\mathbf{u}-\dot{\mathbf{x}})$. For the full analytic derivation we refer to Thompson et al. ([1985]) again.

A.2 Set of RHD equations in strong conservation form

We exhibit the system of equations of radiation hydrodynamics in a somewhat simplified formulation. The following system has been basis for a number of implicit RHD computations (see e.g. Dorfi ([1999])). All the astrophysical assumptions, implications and simplifications can be found in Mihalas and Mihalas ([1984]). In this paper we only want to emphasize the structural form of such a set of equations in a strong conservation form for comoving curvilinear coordinates. Note that for scalar equations the only effectively remaining geometric term inside the derivatives is the volume element $\sqrt{|\mathbf{g}|}$. The vectorial equations contain however also the time-dependent base vectors.

The equations presented in the upcoming sections as well as the coordinate transformation 1 were partially generated using the computer algebra system Mathematica (Wolfram Research [2007]). The source code can be downloaded at https://bitbucket.org/.

A.2.1 Continuity equation

The strong conservation form of the continuity equation

$${\partial }_t\rho +div(\rho \mathbf{u})=0$$

is given for time-dependent coordinates by

$${\partial }_t\left(\sqrt{|\mathbf{g}|}\rho \right)+{\partial }_i(\sqrt{|\mathbf{g}|}{\mathbf{e}}^i\cdot (\mathbf{u}-\dot{\mathbf{x}})\rho )=0.$$

A.2.2 Equation of motion

The equation of motion that we consider in radiation hydrodynamics contains the conservation of moment of plain fluid dynamics (5), the radiative flux as a coupling term (H, with the specific Rosseland opacity of the fluid ${\kappa }_R$), gravitational force (G) and the artificial viscosity term (Q from (17)):

$${\partial }_t(\rho \mathbf{u})+div(\rho {\mathbf{u}}^{\mathrm{T}}\mathbf{u}+\mathbf{P}+\mathbf{Q})+\rho \mathbf{G}-\frac{4\pi }{c}{\kappa }_R\rho \mathbf{H}=0.$$

We inspect the elements of the energy stress tensor a little closer before we give the consistent strong conservation form. We define an effective tensor of gaseous momentum R that accounts for the motion of the coordinates as

$$\mathbf{R}=r^{ij}{\mathbf{e}}_i{\mathbf{e}}_j=\rho (\mathbf{u}-\dot{\mathbf{x}}){\mathbf{u}}^{\mathrm{T}}.$$

The isotropic gas pressure tensor P is defined by the scalar pressure and the metric tensor as $\mathbf{P}=p\mathbf{g}$. The viscous pressure tensor Q has to be modified in the way we suggested in (17). Since in most applications of RHD with self-gravity involved, the gravitational force G is determined by solving the Poisson equation for the potential Φ, we substitute $\mathbf{G}=-\mathrm{\nabla }\mathrm{\Phi }$. The equation of motion in strong conservation form is then written as

$$\begin{array}{r}{\partial }_t\left(\sqrt{|\mathbf{g}|}\rho \mathbf{u}\right)+{\partial }_i(\sqrt{|\mathbf{g}|}(\mathbf{R}+\mathbf{P}+\mathbf{Q})\cdot {\mathbf{e}}^i)\\ +\rho {\partial }_i\left(\sqrt{|\mathbf{g}|}\mathrm{\Phi }{\mathbf{e}}^i\right)-\frac{4\pi }{c}{\kappa }_R\sqrt{|\mathbf{g}|}\rho \mathbf{H}=0.\end{array}$$

The k th component of the strong conservation equation of motion is given by the k th Cartesian component of the unit vector, e.g. in spherical coordinates ${\mathbf{e}}^r=cos\phi sin\vartheta {\mathbf{e}}^x+sin\phi sin\vartheta {\mathbf{e}}^y+cos\vartheta {\mathbf{e}}^z$ and its derivatives. The projection of each physical tensor on the contravariant coordinate lines can be simplified by its contravariant components with respect to its covariant basis without losing strong conservation form, i.e. $\mathbf{f}\cdot {\mathbf{e}}^i=f^{i,j}{\mathbf{e}}_j$. We prefer this form with contravariant components since it meets the native design of the stress tensor R and the pressure tensor P. The tensor of artificial viscosity as given in (18) is brought to contravariant form by the summation

$$\begin{array}{rl}{Q}^{ij}& =Q_{lm}{g}^{li}{g}^{mj}=\cdots \\ =(g^{li}{g}^{mj}({\mathrm{\nabla }}_l{u}_m+{\mathrm{\nabla }}_j{u}_m)-g^{ij}\frac{1}{3}div\mathbf{u})\end{array}$$

and then the equation of motion is given by

$$\begin{array}{r}{\partial }_t\left(\sqrt{|\mathbf{g}|}\rho u^i{\mathbf{e}}_i\right)+{\partial }_i\left(\sqrt{|\mathbf{g}|}(r^{ij}+p^{ij}+Q^{ij}){\mathbf{e}}_j\right)\\ +\rho {\partial }_i\left(\sqrt{|\mathbf{g}|}\mathrm{\Phi }{\mathbf{e}}^i\right)-\frac{4\pi }{c}{\kappa }_R\sqrt{|\mathbf{g}|}\rho H^i{\mathbf{e}}_i=0.\end{array}$$

(27)

A.2.3 Equation of internal energy

The energy equation

$${\partial }_t(\rho ϵ)+div(\mathbf{u}\rho ϵ)+\mathbf{P}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}-4\pi {\kappa }_P\rho (J-S)+\mathbf{Q}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}=0$$

accounts for the thermodynamics of the fluid, namely the energy balance including kinetic and pressure parts as well as inner energy. Latter is a thermodynamic quantity which is associated with the equation of state. The specific inner energy (ϵ) is in case of an ideal fluid its thermic energy. Another term comes from the energy exchange with the radiation field ($(J-S)$-term) containing the specific Planck opacity ${\kappa }_P$ and viscous energy dissipation, expressed by the contraction of the viscosity with the velocity gradient $\mathbf{Q}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}$.

Since we assume isotropic gas pressure $\mathbf{P}=p\mathbf{g}$ we can reformulate its contribution via the Ricci lemma and obtain a very simple scalar expression

$$\mathbf{P}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}=g^{ij}p{\mathrm{\nabla }}_i{u}_j=p{\mathrm{\nabla }}_i{u}^i=pdiv\mathbf{u}.$$

The viscous energy dissipation is given by the contraction

$$\mathbf{Q}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}=Q^{ij}{\mathrm{\nabla }}_i{u}_j=Q^{ij}({\partial }_i{u}_j-{\mathrm{\Gamma }}_{ij}^k{u}_k)=:{\mathsf{E}}_{\mathrm{diss}},$$

and the strong conservative form of the energy equation is then given by

$$\begin{array}{r}{\partial }_t\left(\sqrt{|\mathbf{g}|}\rho ϵ\right)+{\partial }_i(\sqrt{|\mathbf{g}|}\rho ϵ{\mathbf{e}}^i\cdot (\mathbf{u}-\dot{\mathbf{x}}))+\sqrt{|\mathbf{g}|}pdiv\mathbf{u}\\ -4\pi \sqrt{|\mathbf{g}|}{\kappa }_P\rho (J-S)+\sqrt{|\mathbf{g}|}{\mathsf{E}}_{\mathrm{diss}}=0.\end{array}$$

A.2.4 Equation of radiation energy

We write a simplified frequency integrated radiation energy equation in the comoving frame as follows

$${\partial }_{t}J+div(\mathbf{u}J)+cdiv\mathbf{H}+\mathbf{K}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}+c{\chi }_P(J-S)=0.$$

For the scalar energy input of radiative pressure into the material we define a new coupling variable

$$\mathbf{K}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}=K^{ij}{\mathrm{\nabla }}_i{u}_j=:{\mathsf{P}}_{\mathrm{coup}}$$

and in strong conservation form the equation of radiation energy is given by

$$\begin{array}{r}{\partial }_t\left(\sqrt{|\mathbf{g}|}J\right)+{\partial }_i(\sqrt{|\mathbf{g}|}{\mathbf{e}}^i\cdot (J(\mathbf{u}-\dot{\mathbf{x}})+c\mathbf{H}))\\ +\sqrt{|\mathbf{g}|}{\mathsf{P}}_{\mathrm{coup}}+\sqrt{|\mathbf{g}|}c{\chi }_P(J-S)=0.\end{array}$$

A.2.5 Radiative flux equation

Another vectorial conservation law, the radiative flux equation, is given by

$${\partial }_t\mathbf{H}+div(\mathbf{u}\mathbf{H})+cdiv\mathbf{K}+\mathbf{H}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}+c{\chi }_R\mathbf{H}=0.$$

We define an effective radiative flux tensor L analogously to the effective tensor of gaseous momentum:

$$(\mathbf{u}-\dot{\mathbf{x}})\mathbf{H}=:\mathbf{L}=l^{ij}{\mathbf{e}}_i{\mathbf{e}}_j$$

and for the contribution of radiative momentum to the material $\mathbf{H}\cdot \mathrm{\nabla }{\mathbf{u}}^{\mathrm{T}}$ we define another coupling variable F with components

$${\mathsf{F}}_{\mathrm{coup}}^j=H^i{\mathrm{\nabla }}_i{u}^j.$$

The geometrically conservative form of the radiative flux equation in non-steady coordinates is then written as

$$\begin{array}{r}{\partial }_t\left(\sqrt{|\mathbf{g}|}\mathbf{H}\right)+{\partial }_i(\sqrt{|\mathbf{g}|}{\mathbf{e}}^i\cdot (\mathbf{L}+c\mathbf{K}))\\ +\sqrt{|\mathbf{g}|}{\mathsf{F}}_{\mathrm{coup}}+\sqrt{|\mathbf{g}|}{\kappa }_R\rho \mathbf{H}=0.\end{array}$$