Obtaining anisotropic second order derivatives

Discretization of the diffusion term in the RT transfer equation in SPH poses some diffi-culties. We are basically confronted with the task to find an efficient and accurate approx-imation to terms of the form

∂²Qαβ

∂xs∂xk

, (2.41)

where Qαβ is the product of the local Eddington tensor h and the photon density nγ . Simply differentiating a kernel interpolant twice is not a good solution, as this becomes very noisy because the kernel-interpolant of SPH is only second-order accurate. On the other hand, the discretization of the Laplacian discussed by Jubelgas et al. (2004) does not work either, as it only works for the isotropic case.

We now describe the solution we have found for this problem, which basically consists of the task to approximate the second order partial derivatives of an element Qαβ(x) of the matrix Q(x) with a kernel-interpolant. We consider a Taylor-series for Q_αβ(x_j) in the proximity of Qαβ(xi), i.e.

Q_αβ(x_j)−Q_αβ(x_i) =∇Q_αβ

_x

i

·(x_j−x_i) +1 2

X

s,k

∂²Q_αβ

∂xs∂xk

_x

i

(x_j−x_i)_s(x_j−x_i)_k+O(x_j−x_i)³. (2.42) Let us use the short-hand notationx_ij =x_j−x_i and Wij =W(|x_j−x_i|), whereW(r) is the SPH smoothing kernel. Neglecting higher order terms, we multiply the above expansion with

(xij)l Wij,m

|x_ij|² , (2.43)

and integrate over all x_j. Here (xij)l is the l-th component of the vector x_ij, and Wij,m

is the partial derivative of W_ij with respect to the m-component of x_i. In particular, this means we have

Wij,m= (∇ⁱWij)m =−(∇^jWij)m = ∂W(|x_ij|)

∂(x_i)_m =−W^′(|x_ij|)(xij)m

|x_ij| . (2.44)

2.4 Numerical Implementation 27

We now find that

Z (xij)k(xij)lWij,m

|x_ij|² dxj = 0, (2.45)

for all combinations of k, l and m. This is because there is always at least one single component ofx_ij left so that the integral vanishes by symmetry. As a result, the first order term of our integrated Taylor expansion drops out.

We now consider the second order term, where we encounter the expression Tsklm =

Z (xij)s (xij)k (xij)l (xij)m W^′(|x_ij|)

|x_ij|³ dxj. (2.46)

There are a number of different cases. If l and m are equal, then s and k must also be equal, otherwise the integral vanishes. So here we would have three possible contributions to a s, k-sum, corresponding to the three coordinates thats and k can assume. If l and m are unequal, then we must either have s = l and k = m, or have s = m and k = s. So here there are two contributions to a s, k-sum in this case. Evaluating the integral T_sklm for these cases gives:

Tsklm=































−³₅ if l=m and s=k =l =m,

−¹5 if l=m and s=k, but s6=l,

−¹₅ if l6=m, and s=l, k =m or s=m, k=l,

0 in all other cases.

(2.47)

Note that we can pick l and m freely when we multiply the Taylor expansion with the term (2.43) and integrate over it. In particular, we can also use several different choices one after the other and then form a linear combination of the results. This can in fact be used to isolate any of the second derivatives of the Hessian matrix of Qαβ. Let us assume for example that we want to calculate the second derivative of Q_αβ with respect to x₀. Choosing l = m = 0, then the three choices k =s = 0, k =s = 1 and k =s = 2 all give terms that contribute to the integral over the expansion. These are:

2

Z Q(xj)−Q(xi)

|x_ij|² (xij)0Wij,0dxj = 3 5

∂²Q

∂x²₀ +1 5

∂²Q

∂x²₁ +1 5

∂²Q

∂x²₂ . (2.48) HereQ is to be understood asQ=Qαβ for brevity. Based on this, we can now isolate the desired partial derivative by forming a linear combination:

∂²Q

∂x²₀ = 2

Z Q(xj)−Q(xi)

|x_ij|² ×

2(xij)0Wij,0− 1

2(xij)1Wij,1− 1

2(xij)2Wij,2

dxj. (2.49) In a similar fashion, we can obtain a mixed partial derivative in the following way:

∂²Q

∂x₀∂x₁ = 2

Z Q(xj)−Q(xi)

|x_ij|² ×

5

4(xij)0Wij,1+5

4(xij)1Wij,0

dxj. (2.50)

Formulae for all other second-order partial derivatives can be obtained from these ex-pressions by cyclic permutation. Also, they are valid for each of the matrix elements Qαβ. Using these results, we can now turn to obtaining an expression for the sum of the second derivatives, as needed in the anisotropic diffusion equation. Based on the above, we can write the desired expression in the compact form:

∂²Qαβ

∂xα∂xβ

_x

i

= 2

Z x^T_ij[Q(x˜ j)−Q(x˜ i)]∇ⁱWij

|x_ij|² dx_j. (2.51)

Here we defined a new matrix Q˜ through the matrix elements of the original matrix Q= (Qαβ), in the following way:

Q˜ = 5 2Q−1

2Tr(Q)I. (2.52)

Inspection of this result highlights one interesting issue that could potentially become a numerical stability problem in certain situations. The matrix Q˜ is not guaranteed to correspond to a positive definite quadratic form when it is used in the SPH discretization form of equation (2.51). If the radiation transfer is very anisotropic, the matrix Q˜ can contain negative diagonal elements and thus gives rise to an ‘anti-diffusive’ behavior in the discretized radiation transfer equation, where radiation is transported from a particle of lower radiation intensity to one with higher radiation intensity. It is not clear right away whether this will lead to numerical stability problems of the radiative diffusion treatment, but it could.

In case this is a problem, one way to avoid it would be to somehow suppress transport of radiation opposite to the direction of the gradient of the radiation intensity between a particle pair. Another way is to add in an isotropic component toQ˜ such that

Q˜^∗ =αQ˜ + (1−α)I

3. (2.53)

Here the idea is to makeQ˜^∗ slightly more isotropic, such thatQ˜^∗ becomes positive definite again. In order to guarantee this, we need to assign α = ²₅. Form equations (2.52) and (2.53) we can see that this ‘anisotropy-limited’ matrix is then actually Q˜^∗ = Q. One interpretation of this result is that the unmodified matrixQ mediates diffusion which is a mix of 2/5 of the ‘correct’ anisotropic diffusion and 3/5 of isotropic diffusion. In some of our tests we will compare results from both formulations of the matrix. We will refer to Q˜^∗ as the ‘anisotropy-limited’ tensor and toQ˜ as the ‘fully-anisotropic’ tensor.

Note that in the case where Q is diagonal and proportional to the identity matrix, equation (2.51) reduces to the isotropic result for a scalar function derived by Jubelgas et al.

(2004) for the thermal conduction problem. The important point about equation (2.51) is that it involves only a first order derivative of the kernel function. As a result, it can be discretized straightforwardly in the usual SPH way, where the integration is replaced by a sum over all neighboring SPH particles within the kernel volume.

2.4 Numerical Implementation 29

Im Dokument Numerical Radiative Transfer and the Hydrogen Reionization of the Universe (Seite 44-47)