Roland Haas <>

August 17 2020


A FileReader like thorn that uses InterpolateLocalUniform to interpolate the data read in onto the new grid.

1 Introduction

IOUtil [1] defines a file reading interface via the IO parameters filereader_ID_dir, filereader_ID_files, filereader_ID_vars parameters. IOUtil uses the checkpoint recovery routines to read grid variable data from disk. This requires that the grid resolution used to write the files agrees with the resolution using when reading in data and that, for each refinement level, the to be read in region is fully contained in the corresponding refinement level in the source files.

ReadInterpolate relaxes this restriction by allowing arbitrary changes in resolution and grid structure in between source files and simulation. It uses Cactus local interpolation operators CCTK_InterpLocalUniform to interpolate data using the highest resolution source data available for each target grid point.

It provides parameters files and only_these_datasets that mirror IOUtil’s parameters to select files and datasets in IO::filereader_ID_dir.

Auxiliary parameters are provided to fine tune what is read and to control the interpolator.

2 Physical System

ReadInterpolate takes input data ψ on a grid patch described by its origin O, grid spacing Δx and shape N. Data is interpolated to the coordinates of the simulation grid x using Cactus’ interpolation operators: ψ(x) = L(ψ,x;O,Δx,N).

3 Numerical Implementation

For a target point x multiple datasets in the source files may contain this physical coordinate. When this happens the following selection algorithm applies:

Once a source data set has been selected its data and the target coordinate triplet are passed to CCTK_InterpLocalUniform for interpolation.

4 Using This Thorn

This thorn is intended as an initial data thorn, scheduling itself in the CCTK_INITIAL time bin of Cactus.

4.1 Obtaining This Thorn

This thorn is part of the Einstein Toolkit in the EinsteinInitialData repository.

4.2 Basic Usage

To use this thorn, set IO::filereader_ID_dir to the directory containing the data files, set ReadInterpolate’s files to the base names of the file to read and only_these_datasets to a comma separated list of regular expressions to match the datasets to read. To obtain the list of data sets present in the file the h5ls tool, which is part of HDF5, can be used. The parameter max_number_of_read_variables must be set to the expected number of grid variables to be read.

If desired read in data can be moved to a different position on the grid using shift_read_datasets_by[3] which relocates the origin of the coordinates in the read in data to the given coordinate triplet.

4.3 Special Behavior

ReadInterpolate attempts to avoid one sided interpolation as much as possible, yet tried to allow reading in data at the outer boundaries of the grid. A heuristic is used to allow one sided derivatives there. This requires changing interpolator_pars to allow off centered interpolation stencils.

A number of parameters are provided to select which datasets to read:

ignore any datasets with refinement level less than this,
ignore any datasets with refinement level higher than this,

Some Carpet version contain a bug and do not set origin correctly for cell centered grids (it is off by half a grid point, i.e. it does not take the staggering into account). This option fixes this issue by offsetting all datasets that should be offset. Setting fix_cell_centered_origins adds an offset to correct for this. This is incorrect for files where Carpet already used the correct origin. This parameter only affects cell centered input data.

Interpolation parameters ReadInterpolate uses interpolator_name as the local interpolator to interpolate read in data to simulation grid points. It passes interpolator_vars to the interpolator. The default values select AEILocalInterp’s Hermite interpolator and forbid off centered interpolation stencils or extrapolation.

ReadInterpolate uses epsilon as a fuzziness parameter when determining if a point is contained in a read in dataset. This may require changes for very large grids.

When deciding if a read in dataset is large enough to allow for interpolation without off centered stencils ReadInterpolate requires knowledge about the stencil width of the chosen interpolation operator. This can be explicitly specified as interpolator_half_width or alternatively assumed to be cctk_nghostzones if this parameter is set to -1, which is the default.

This makes sure that the interpolator does not off center the interpolation stencil if there are insufficient points to interpolate, which can happen if there are insufficient ghost-zones for the interpolation method used, and can lead to processor-number dependent results.

ReadInterpolate is no aware of any symmetry conditions used in either the simulation used to produce the data files the current simulation. Thus it will interpolated into symmetry points in the same manner as for outer boundary points. In these situations the parameter enforce_symmetries_after_reading to apply symmetry boundary conditions after having read in all data. This does not enable ReadInterpolate to read in data produced on, for example, an octant grid and provide data on a full grid.

4.4 Interaction With Other Thorns

ReadInterpolate assumes that data files were written by CarpetIOHDF5 and may not work with files written by PUGH’s IOHDF5.

This thorn is intended as an initial data thorn, however since it is generic it does not extend ADMBase’s or HydroBase’s initial_XXX parameters.

ReadInterpolate can be used to read data into curvilinear grids set up by Llama but cannot read data from curvilinear grids.

4.5 Examples

The test suite file admhydro.par shows how to use ReadInterpolate to read in metric and fluid variables:

ActiveThorns = "ReadInterpolate"  
# ReadInterpolate ignore these but they are required to allocate storage  
ADMBase::initial_data    = "Cartesian Minkowski"  
ADMBase::initial_lapse   = "one"  
ADMBase::initial_shift   = "zero"  
ADMBase::initial_dtlapse = "zero"  
ADMBase::initial_dtshift = "zero"  
HydroBase::initial_hydro = "zero"  
IO::filereader_ID_dir = "tov_write"  
ReadInterpolate::files = "alp rho"  
ReadInterpolate::max_number_of_read_variables = 2  
ReadInterpolate::shift_read_datasets_by[0] = +3.14  
ReadInterpolate::shift_read_datasets_by[1] = -42  
ReadInterpolate::shift_read_datasets_by[2] = -6

4.6 Support and Feedback

Please use the Einstein Toolkit ticket system to report issues.

5 History

This thorn was originally inspired by SpEC’s [2] ReadTensorsFromDiskWithMap item and developed to read in single neutron star data produced by THC [3] into GRHydro [4]. It has been found useful to in combination with CT_MultiLevel since then.

5.1 Acknowledgments

David Radice and Philipp Mösta provided bug reports and comments during the time this thorn was developed.



[2]   Spectral Einstein Code,

[3]   WhiskyTHC: The General-Relativistic Templated Hydrodynamics Code,

[4]   GRHydro is an evolution code for a general-purpose 3D relativistic hydrodynamics,

6 Parameters

Scope: private  BOOLEAN

Description: apply symmetry boundary conditions after reading data

  Default: no

Scope: private  REAL

Description: Tolerance when deciding if a grid point lies within an overlap region

Range   Default: 1e-12
any positive number, should be similar to AEILocalInterp’s tolerance

Scope: private  STRING

Description: List of basenames of files to read in as initial data (e.g. omit the .file_XXX.h5 part)

Range   Default: (none)
Space-separated list of initial data filenames
An empty string for not recovering initial data

Scope: private  BOOLEAN

Description: offset rl>1 cell centered input data by half a grid step to fix broken HDF5 files

  Default: no

Scope: private  INT

Description: maximum number of points to one side of the interpolated location the interpolator uses

Range   Default: -1
this is similar to cctk_nghostzones, 1 for linear and constant interpolation
2nd for second order and 3rd interpolation
use cctk_nghostzones

Scope: private  STRING

Description: Which interpolator should I use

Range   Default: Hermite polynomial interpolation
Any nonempty string

Scope: private  STRING

Description: Parameters for the interpolator

Range  Default: order=3 boundary_off_centering_tolerance={1e12 1e12 1e12 1e12 1e12 1e12} boundary_extrapolation_tolerance={1e-12 1e-12 1e-12 1e-12 1e-12 1e-12}
”Any string that Util_TableSetFromStr ing() will take”

Scope: private  INT

Description: how many variables will be read in?

Range   Default: 1
any positive number

Scope: private  INT

Description: ignore any datasets with refinement level higher than this

Range   Default: 1000
maximum refinement level to keep

Scope: private  INT

Description: ignore any datasets with refinement level less than this

Range   Default: (none)
minimum refinement level to keep

Scope: private  STRING

Description: comma separated list of dataset regular expression patterns to consider for reading

Range   Default: (none)
any valid Cactus regular expression
see [1] below
refinement level 6, all metric
metric and rho everywhere



Scope: private  REAL

Description: add to all read in coordinates

Range   Default: 0.0
what used to be origin in datasets appear here

Scope: private  KEYWORD

Description: control how test data is used

Range   Default: no
do not use test data
generate test data
compare read data to test data

Scope: private  INT

Description: how much diagnostic output to produce

Range   Default: 1
only errors
warn about suspicous states
more and more detailed debug output (values up to 10 are defined. 10 is very verbose.)

Scope: shared from IO STRING

7 Interfaces






Grid Variables


  Group Names    Variable Names    Details   

reflevelseen reflevelseen   compact0
  descriptionhighest refinement level in source data that overlaps this point so far
  group typeGF
  tagsCheckpoint=”no” Prolongation=”none”
 variable typeINT

interpthispoint interpthispoint   compact0
  descriptionfill this point this time around
  group typeGF
  tagsCheckpoint=”no” Prolongation=”none”
 variable typeINT

interp_coords   compact0
interp_x   descriptioncoordinate arrays for the interpolator
interp_y   dimensions3
interp_z   distributionDEFAULT
interp_data   group typeGF
  tagsCheckpoint=”no” Prolongation=”none”
 variable typeREAL

test_values test_values   compact0
  descriptionstore values for test
  group typeGF
  tagsCheckpoint=”yes” ProlongationOperator=”none”
 variable typeREAL

test_results test_results   compact0
  descriptiondifference values for test
  group typeGF
  tagsCheckpoint=”no” ProlongationOperator=”none”
 variable typeREAL

Uses header:


8 Schedule

This section lists all the variables which are assigned storage by thorn EinsteinInitialData/ReadInterpolate. 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.



  test_values[3] test_results

Scheduled Functions

CCTK_PARAMCHECK (conditional)


  sanity check given parameters


 Options: global
 Type: function

CCTK_INITIAL (conditional)


  read in datasets from disk


 After: admbase_initialdata
 Type: group

ReadInterpolate_ReadData (conditional)


  read in datasets


 Options: level
 Storage: reflevelseen
 Type: function



  free memory used for dataset caches


 Options: global
 Type: function

ReadInterpolate_ReadData (conditional)


  enforce symmeries if desired


 After: readinterpolate_read
 Options: level
 Type: function

ReadInterpolate_ReadData (conditional)


  apply symmetry conditions to read in variables



CCTK_INITIAL (conditional)


  generate polynomial test data


 Type: function



  compare to polynomial test data


 Type: function

Aliased Functions


Alias Name:        Function Name:
ApplyBCs ReadInterpolate_ApplyBCs