Outflow

Roland Haas <roland.haas@physcis.gatech.edu>

August 15 2009

Abstract

Outflow calculates the flow of rest mass density across a SphericalSurface, eg. and apparent horizon or a sphere at “infinity”.

1 Introduction

Hydrodynamic simulations should conserve eg. the total rest mass in the system (outside he event horizon, that is). This thorn allow to measure the flux of rest mass across a given SphericalSurface.

2 Physical System

The Valencia formalism defines its \(D = \rho W\) variable such that one can define a total rest mass density as: \begin {displaymath} M = \int \sqrt {\gamma } D d^3x \end {displaymath} and from the EOM of \(D\) (see [1]) \begin {displaymath} \frac {\partial }{\partial x^0} (\sqrt {\gamma }D )+ \frac {\partial }{\partial x^i} (\sqrt {\gamma } \alpha D (v^i- \beta ^i/\alpha )) = 0 \end {displaymath} one obtains:

\begin {eqnarray} \dot M &=& - \int _V \frac {\partial }{\partial x^i} (\sqrt \gamma \alpha D(v^i - \beta ^i/\alpha )) d^3x \\ &=& -\int _{\partial V} \sqrt \gamma \alpha D (v^i-\beta ^i/\alpha ) d\sigma _i \end {eqnarray}

where \(\sigma _i\) is the ordinary flat space directed surface element of the enclosing surface, eg. \begin {displaymath} d\sigma _i = \hat r_i r^2 \sin \theta d\theta d\phi \end {displaymath} with \(\hat r_i = [\cos \phi \sin \theta , \sin \phi \sin \theta , \cos (\theta )]\) for a sphere of radius r.

For a generic SphericalSurface parametrized by \(\theta \) and \(\phi \) one has:

\begin {eqnarray} x &=& \bar r(\theta , \phi ) \cos \phi \sin \theta \\ y &=& \bar r(\theta , \phi ) \sin \phi \sin \theta \\ z &=& \bar r(\theta , \phi ) \cos \theta \end {eqnarray}

where \(\bar r\) is the isotropic radius. Consequently the surface element is

\begin {eqnarray} d \sigma _i &=& \left ( \frac {\partial \vec {\bar r}}{\partial \theta } \times \frac {\partial \vec {\bar r}}{\partial \phi } \right )_i \\ &=& \bar r^2 \sin \theta \hat {\bar r}_i - \frac {\partial \bar r}{\partial \theta } \bar r \sin \theta \hat \theta _i - \frac {\partial \bar r}{\partial \phi } \bar r \hat \phi _i \end {eqnarray}

where \(\hat {\bar r}\), \(\hat \theta \) and \(\hat \phi \) are the flat space standard unit vectors on the sphere [3].

3 Numerical Implementation

We implement the surface integral by interpolating the required quantities (\(g_{ij}\), \(\rho \), \(v^i\), \(\beta ^i\), \(\alpha \)) onto the spherical surface and then integrate using a fourth order convergent Newton-Cotes formula.

For the \(\theta \) direction SphericalSurfaces defines grid points such that \begin {displaymath} \theta _i = -(n_\theta - 1/2) \Delta _\theta + i \Delta _\theta \qquad 0 \le i < N_\theta -1 \end {displaymath} where \(N_\theta \) is the total number of intervals (sf_ntheta), \(n_\theta \) is the number of ghost zones in the \(\theta \) direction (nghosttheta) and \(\Delta _\theta = \frac {\pi }{N_\theta -2n\theta }\). Since with this \begin {displaymath} \theta _i \in [\Delta _\theta /2, \pi - \Delta _\theta /2] \qquad \mbox {for $i \in [n_\theta , N_\theta -n_\theta -1]$} \end {displaymath} we do not have grid points at the end of the interval \(\{0,\pi \}\) we derive an open extended Newton-Cotes formula from Eq. 4.1.14 of [2] and a third order accurate extrapolative rule (see Maple worksheet). We find

\begin {eqnarray} \int _{x_0}^{x_{N-1}} f(x) dx &\approx & h \Bigl \{ \frac {13}{12} f_{1/2} + \frac {7}{8} f_{3/2} + \frac {25}{24} f_{5/2} + f_{7/2} + f_{9/2} + \cdots + f_{N-1-7/2} \nonumber \\ && + f_{N-1-9/2} + \frac {25}{24} f_{N-1-5/2} + \frac {7}{8} f_{N-1-3/2} + \frac {13}{12} f_{N-1-1/2} \Bigr \} + O(1/N^4) \end {eqnarray}

.

For the \(\phi \) direction SphericalSurfaces defines grid points such that \begin {displaymath} \phi _i = -n_\phi \Delta _\phi + i \Delta _\phi \qquad 0 \le i < N_\phi -1 \end {displaymath} where \(N_\phi \) is the total number of intervals (sf_nphi), \(n_\phi \) is the number of ghost zones in the \(\phi \) direction (nghostphi) and \(\Delta _\phi = \frac {\pi }{N_\phi -2n\phi }\). With this \begin {displaymath} \phi _i \in [0, 2 \pi - \Delta _\phi ] \qquad \mbox {for $i \in [n_\phi , N_\phi -n_\phi -1]$} \end {displaymath} we use a simple extended trapezoid rule to achieve spectral convergence due to the periodic nature of \(\phi \) (note: \(x_N = x_0\))

\begin {eqnarray} \int _{x_0}^{x_N} f(x) dx \approx & h \sum _{i=0}^{N-1} f_i \end {eqnarray}

.

The derivatives of \(\vec {\bar r}\) along \(\theta \) and \(\phi \) are obtained numerically and require at least two ghost zones in \(\theta \) and \(\phi \).

4 Using This Thorn

Right now surface can only be prescribed by SphericalSurfaces, the flux through each surface is output in in a file outflow_det_%d.asc

4.1 Interaction With Other Thorns

Takes care to schedule itself after SphericalSurfaces_HasBeenSet.

4.2 Examples

See the parameter file in the test directory. For spherical symmetric infall the flux through all detectors should be equal (since rest mass must be conserved).

References

[1]   José A. Font, “Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity”, Living Rev. Relativity 11, (2008), 7. URL (cited on August 15. 2009): http://www.livingreviews.org/lrr-2008-7

[2]   Press, W. H. (2002). Numerical recipes in C++: the art of scientific computing. Array Cambridge, UK: Cambridge University Press.

[3]   Charles W. Misner, Kip S. Thorne, and John Archibald Wheeler (1973). “Gravitation”. New York: W.H. Freeman and Company.

5 Parameters

compute_every
Scope: private  INT
Description: How frequently to compute the outflow
Range   Default: 1
1:*
number of iterations

compute_every_det
Scope: private  INT
Description: How frequently to compute the outflow
Range   Default: -1
-1
take from compute_every
do not compute
1:*
number of iterations

coord_system
Scope: private  STRING
Description: What is the coord system?
Range   Default: cart3d
.*
Any smart string will do

extra_variables
Scope: private  STRING
Description: extra (scalar) variables to project onto the spherical surface
Range   Default: (none)
.*
Any Cactus variables

interpolator_name
Scope: private  STRING
Description: Which interpolator should I use
Range   Default: Hermite polynomial interpolation
.+
Any nonempty string

interpolator_pars
Scope: private  STRING
Description: Parameters for the interpolator
Range   Default: order=3
.*
”Any string that Util_TableSetFromStr ing() will take”

num_detectors
Scope: private  INT
Description: how many detectors do we have
Range   Default: (none)
0:*
number of detectors

num_thresholds
Scope: private  INT
Description: how many thresholds to use
Range   Default: (none)
ignore thresholds, output only one flux
1:20
how many to use

out_format
Scope: private  STRING
Description: Which format for Scalar floating-point number output
Range   Default: .19g
see [1] below
output with given precision in exponential / floating point notation

[1] 

\^({\textbackslash}.[1]?[0-9])?[EGefg]\$

output_2d_data
Scope: private  BOOLEAN
Description: output 2d data on spherical surface
  Default: no

output_relative_coordinates
Scope: private  BOOLEAN
Description: output coordinates on spheroid relative to sf_centroid, absolute coordinates otherwise
  Default: no

override_radius
Scope: private  BOOLEAN
Description: do not take radius values from spherical surface
  Default: no

rad_rescale
Scope: private  REAL
Description: Factor for scaling radius for this detector surface (e.g. for buffer around AH surface)
Range   Default: 1
(0:*
Anything positive

radius
Scope: private  REAL
Description: spherical radius for this detector surface
Range   Default: -1
-1
Invalid
(0:*
any positive value

surface_index
Scope: private  INT
Description: the indices of the sphericalsurfaces from which we take the surfaces to compute the outflow on
Range   Default: -1
-1
Invalid
0:*
valid surface

surface_name
Scope: private  STRING
Description: the indices of the sphericalsurfaces from which we take the surfaces to compute the outflow on
Range   Default: (none)
use surface_index
.*
any string

threshold
Scope: private  REAL
Description: compute fluxes with Lorentz factors / specific energy at infinity above these levels only
Range   Default: (none)
*:*
Any real number

threshold_on_var
Scope: private  KEYWORD
Description: Which variable to use as threshold
Range   Default: w_lorentz
eninf
(Approximate) specific energy at infinity, - u_t - 1
w_lorentz
Lorentz factor

verbose
Scope: private  INT
Description: How verbose do you want it?
Range   Default: 1
0:*
0 is nothing, 1 is interesting, >1 is crazy

maxnphi
Scope: shared from SPHERICALSURFACE INT

maxntheta
Scope: shared from SPHERICALSURFACE INT

nghostsphi
Scope: shared from SPHERICALSURFACE INT

nghoststheta
Scope: shared from SPHERICALSURFACE INT

sphericalsurfaces_nsurfaces
Scope: shared from SPHERICALSURFACE INT

6 Interfaces

General

Implements:

outflow

Inherits:

admbase

hydrobase

sphericalsurface

Grid Variables

6.0.1 PRIVATE GROUPS
  Group Names    Variable Names    Details   
outflow_flux outflow_flux   compact0
  descriptionflux of mass through the detectors
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  tagscheckpoint=”no”
  timelevels1
 vararray_sizenum_detectors
 variable typeREAL
fluxdens_projected fluxdens_projected   compact0
  description2D (theta
    descriptionphi) grid arrays for flux density
  dimensions2
  distributionCONSTANT
  group typeARRAY
  sizeSPHERICALSURFACE::MAXNTHETA
    sizeSPHERICALSURFACE::MAXNPHI
  tagscheckpoint=”no”
  timelevels1
 vararray_sizenum_detectors
 variable typeREAL
w_lorentz_projected w_lorentz_projected   compact0
  description2D (theta
    descriptionphi) grid arrays for Lorentz factor
  dimensions2
  distributionCONSTANT
  group typeARRAY
  sizeSPHERICALSURFACE::MAXNTHETA
    sizeSPHERICALSURFACE::MAXNPHI
  tagscheckpoint=”no”
  timelevels1
 vararray_sizenum_detectors
 variable typeREAL
eninf_projected eninf_projected   compact0
  description2D (theta
    descriptionphi) grid arrays for the specific energy at infinity
  dimensions2
  distributionCONSTANT
  group typeARRAY
  sizeSPHERICALSURFACE::MAXNTHETA
    sizeSPHERICALSURFACE::MAXNPHI
  tagscheckpoint=”no”
  timelevels1
 vararray_sizenum_detectors
 variable typeREAL
surface_projections   compact0
surface_projection_0   description2D (theta
    descriptionphi) grid arrays for points on the spherical surfaces
surface_projection_1   dimensions2
surface_projection_2   distributionCONSTANT
surface_projection_3   group typeARRAY
surface_projection_4   sizeSPHERICALSURFACE::MAXNTHETA
    sizeSPHERICALSURFACE::MAXNPHI
surface_projection_5   tagscheckpoint=”no”
surface_projection_6   timelevels1
surface_projection_7  vararray_sizenum_detectors
surface_projection_8  variable typeREAL

7 Schedule

This section lists all the variables which are assigned storage by thorn EinsteinAnalysis/Outflow. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

 

Always:  
outflow_flux  
fluxdens_projected  
w_lorentz_projected 
eninf_projected  
surface_projections  
   

Scheduled Functions

CCTK_ANALYSIS

  outflow

  compute outflow

 

 Language:c
 Options: global
 Type: function

CCTK_STARTUP

  outflow_setup

  set up global ompute data structures

 

 Language:c
 Options: global
 Type: function