IceModelVec2T Class Reference

#include <iceModelVec2T.hh>

Inheritance diagram for IceModelVec2T:

Public Member Functions
	IceModelVec2T ()
	IceModelVec2T (const IceModelVec2T &other)
virtual void	set_n_records (unsigned int N)
	Sets the number of records to store in memory. Call it before calling create().
virtual PetscErrorCode	create (IceGrid &mygrid, const char my_short_name[], bool local, int width=1)
virtual PetscErrorCode	init (string filename)
virtual PetscErrorCode	update (double t_years, double dt_years)
	Read some data to make sure that the interval (t_years, t_years + dt_years) is covered.
virtual PetscErrorCode	set_record (int n)
	Sets the record number n to the contents of the (internal) Vec v.
virtual PetscErrorCode	get_record (int n)
	Sets the (internal) Vec v to the contents of the nth record.
virtual double	max_timestep (double t_years)
	Given the time t_years and the current selected time-step dt_years, determines the maximum possible time-step this IceModelVec2T allows.
virtual PetscErrorCode	interp (double t_years)
	Extract data corresponding to t_years using linear interpolation.
virtual PetscErrorCode	interp (int i, int j, int N, PetscScalar times, PetscScalar results)
	Gets an interpolated time-series out. Has to be surrounded with begin_access() and end_access().
virtual PetscErrorCode	average (int i, int j, double t_years, double dt_years, double &result)
	Finds the average value at i,j over the interval (t_years, t_years + dt_years) using trapezoidal rule.
virtual PetscErrorCode	begin_access ()
	Checks if an IceModelVec is allocated and calls DAVecGetArray.
virtual PetscErrorCode	end_access ()
	Checks if an IceModelVec is allocated and calls DAVecRestoreArray.
Protected Member Functions
virtual PetscErrorCode	destroy ()
virtual PetscErrorCode	get_array3 (PetscScalar ***&a3)
virtual PetscErrorCode	update (int start)
	Update by reading at most n_records records from the file.
virtual PetscErrorCode	discard (int N)
	Discard the first N records, shifting the rest of them towards the "beginning".
Protected Attributes
vector< double >	times
	all the times available in filename
vector< double >	T
	times stored in memory
string	filename
	file to read (regrid) from
LocalInterpCtx *	lic
DA	da3
Vec	v3
	a 3D Vec used to store records
void ***	array3
int	n_records
	maximum number of records to store in memory
int	first
	in-file index of the first record stored in memory

Detailed Description

A class for storing and accessing 2D time-series (for climate forcing)

This class was created to read time-dependent and spatially-varying climate forcing data, in particular snow temperatures and precipitation.

It allocates a number (given as an argument to the create() method) of records and reads them from a file if necessary.

If requests (calls to update()) go in sequence, every records should be read only once.

Note that this class is optimized for use with a PDD scheme -- it stores records so that data corresponding to a grid point are stored in adjacent memory locations.

IceModelVec2T is always global (i.e. has no ghosts).

Both versions of interp() use linear interpolation and extrapolate (by a constant) outside the available range.

Usage example:

  // initialization:
  char filename[] = "climate_inputs.nc";
  IceModelVec2T v;
  v.set_n_records(config.get("climate_forcing_buffer_size"));
  ierr = v.create(grid, "snowtemp", false); CHKERRQ(ierr);
  ierr = v.set_attrs("climate_forcing", "snow surface temperature", "K", ""); CHKERRQ(ierr);
  ierr = v.init(filename); CHKERRQ(ierr);

  // actual use:

  // ask it how long a time-step is possible:
  double t = 0, max_dt = 1e16;
  ierr = v.max_timestep(t, max_dt); CHKERRQ(ierr);

  // update (this might read more data from the file)
  ierr = v.update(t, max_dt); CHKERRQ(ierr);

  // get a 2D field (such as precipitation):
  ierr = v.interp(t + max_dt / 2.0); CHKERRQ(ierr);
  // at this point v "looks" almost like an IceModelVec2S with data we need (and
  // can be type-cast to IceModelVec2S)

  // get a time-series for every grid location:
  int N = 21;
  vector<double> ts(N), values(N);
  double dt = max_dt / (N - 1);

  for (int j = 0; j < N; ++j) {
    ts[j] = t + dt * j;
  }

  // is is OK to call update() again, it will not re-read data if at all possible
  ierr = v.update(t, max_dt); CHKERRQ(ierr);

  ierr = v.begin_access(); CHKERRQ(ierr);
  for (PetscInt i=grid->xs; i<grid->xs+grid->xm; ++i)
    for (PetscInt j=grid->ys; j<grid->ys+grid->ym; ++j) {
      ierr = v.interp(i, j, N, &ts[0], &values[0]); CHKERRQ(ierr);
      // do more
    }
  ierr = v.end_access(); CHKERRQ(ierr);

  // compute an average value over a time interval at a certain grid location:

  double average;
  ierr = v.begin_access(); CHKERRQ(ierr);
  ierr = v.average(grid.xs, grid.ys, t, max_dt, average); CHKERRQ(ierr);
  ierr = v.end_access(); CHKERRQ(ierr);