RAMSES_poisson_data.hh

/*
		This file is part of libRAMSES++
			a C++ library to access snapshot files
			generated by the simulation code RAMSES by R. Teyssier

    Copyright (C) 2008-09  Oliver Hahn, ojha@gmx.de

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef __RAMSES_POISSON_DATA_HH
#define __RAMSES_POISSON_DATA_HH

#include <fstream>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>

#include "FortranUnformatted_IO.hh"
#include "RAMSES_info.hh"
#include "RAMSES_amr_data.hh"

namespace RAMSES{
namespace POISSON{


/*!
 * @class data
 * @brief encapsulates hydro data from a RAMSES simulation snapshot
 *
 * This class provides low-level read access to RAMSES hydro_XXXXX.out files.
 * Data from a given list of computational domains can be read and is
 * stored in internal datastructures.
 * Access to cell position and threaded tree structure of the cell is provided
 * through the member functions of class RAMSES_amr_level.
 * @sa RAMSES_amr_level
 */
template< typename TreeType_, typename Real_=double >
class data : public RAMSES::HYDRO::proto_data<TreeType_,Real_>{

public:

	//! the poisson file header structure
	struct header{
		unsigned ncpu;		//!< number of CPUs in simulation
		unsigned ndim;		//!< number of spatial dimensions
		unsigned nlevelmax;	//!< maximum allowed refinement level
		unsigned nboundary;	//!< number of boundary regions
	};


	std::string 	m_fname;		//!< the file name	
	struct header	m_header;	 	//!< header meta data

protected:
	//! generates a poisson_XXXX filename for specified cpu
	std::string gen_fname( int icpu );

	//! generate poisson_XXXXX filename from info filename
	std::string rename_info2poisson( const std::string& info );

	//! generate poisson_XXXXX filename from amr filename
	std::string rename_amr2poisson( const std::string& info );

	//! read header data containing meta information
	void read_header( void );

	//! perform read operation of poisson variables (internal use)
	/*!
	 *
	 */
	void read( void );

public:

	//! constructor for poisson solver gravity data
	/*!
	 * @param AMRtree underlying AMR hierarchical tree data structure
	 */
	explicit data( TreeType_& AMRtree )
	: proto_data<TreeType_,Real_>( AMRtree ),
	  m_fname( rename_amr2poisson(AMRtree.m_fname) )
	{
		read_header();

		if( this->m_cpu > m_header.ncpu || this->m_cpu < 1 )
			throw std::runtime_error("RAMSES::POISSON::data : expect to read from out of range CPU.");

		if( this->m_minlevel < 0 || this->m_maxlevel >= m_header.nlevelmax )
			throw std::runtime_error("RAMSES::POISSON::data : requested level is invalid.");
	}
	
	//! access the value of the cells associated with the oct designated by the iterator
	/*!
	 * @param it the grid iterator pointing to the current oct
	 * @param ind index of the child cell of the current oct (0..7)
	 * @param idir the cartesian direction of the force vector (0..2)
	 */
	inline ValueType_& cell_value( const typename TreeType_::iterator& it, int ind, int idir )
	{
		unsigned ipos   = it.get_absolute_position();
		unsigned ilevel = it.get_level();//-m_minlevel;
		return (m_var_array[ilevel])[3*(m_twotondim*ipos+ind)+idir];
	}
	
	//! access the value of the cells associated with the oct designated by the iterator
	/*!
	 * a call to this function will always fail since the poisson files contain only forces so far
	 * @param it the grid iterator pointing to the current oct
	 * @param ind index of the child cell of the current oct (0..7)
	 */
	inline ValueType_& cell_value( const typename TreeType_::iterator& it, int ind )
	{
		throw std::runtime_error("You should not call this two variable cell_value(.,.) function for forces!");
		return (m_var_array[ilevel])[3*(m_twotondim*ipos+ind)];
	}

	//! access the value of the cells associated with the oct designated by the iterator
	/*!
	 * @param it the grid iterator pointing to the current oct
	 * @param ind index of the child cell of the current oct (0..7)
	 * @param idir the cartesian direction of the force vector (0..2)
	 */
	inline ValueType_& operator()( const typename TreeType_::iterator& it, int ind, int idir )
	{	return cell_value(it,ind,idir); }

};


template< typename TreeType_, typename Real_ >
void data<TreeType_,Real_>::read_header( void )
{
	FortranUnformatted ff( gen_fname( this->m_cpu ) );
	
	//-- read header data --//
	ff.read( m_header.ncpu );
	ff.read( m_header.ndim );
	ff.read( m_header.nlevelmax );
	ff.read( m_header.nboundary );
}


template< typename TreeType_, typename Real_ >
void data<TreeType_,Real_>::read( void )
{
	this->m_var_array.clear();
	//int twotondim = (int)(pow(2,m_header.ndim)+0.5);

	FortranUnformatted ff( gen_fname( this->m_cpu ) );

	//.. skip header entries ..//
	ff.skip_n_from_start( 4 ); //.. skip header

	this->m_var_array.clear();

	for( unsigned ilvl = 0; ilvl<=this->m_maxlevel; ++ilvl ){

		this->m_var_array.push_back( std::vector<Real_>() );

		for( unsigned icpu = 0; icpu<m_header.ncpu+m_header.nboundary; ++icpu ){

				unsigned file_ilevel, file_ncache;
				ff.read(file_ilevel);
				ff.read(file_ncache);

				if( file_ncache == 0 )
						continue;

			if( ilvl >= this->m_minlevel ){
				if( file_ilevel != ilvl+1 )
					throw std::runtime_error("RAMSES::POISSON::data::read : corrupted file " \
						 "or file seek failure in file \'"+m_fname+"\'.");


				std::vector<float> tmpx,tmpy,tmpz;
				for( unsigned i=0; i<this->m_twotondim; ++i )
				{
					ff.read<double>( std::back_inserter(tmpx) );	//.. x-force
					ff.read<double>( std::back_inserter(tmpy) );	//.. y-force
					ff.read<double>( std::back_inserter(tmpz) );	//.. z-force
				}
				//.. reorder array to increase data locality..//
				this->m_var_array.reserve( 3*tmpx.size() );
				for( unsigned i=0; i<file_ncache; ++i ){
					for( unsigned j=0; j<this->m_twotondim; ++j ){
						//... interleave directional data ...//
						(this->m_var_array.back()).push_back( tmpx[i+j*file_ncache] );
						(this->m_var_array.back()).push_back( tmpy[i+j*file_ncache] );
						(this->m_var_array.back()).push_back( tmpz[i+j*file_ncache] );
					}
				}
			}else{
				ff.skip_n( this->m_twotondim*3 );
			}
		}
	}
}


template< typename TreeType_, typename Real_ >
std::string data<TreeType_,Real_>::gen_fname( int icpu )
{
	std::string fname;
	char ext[32];
	fname = m_fname;
	fname.erase(fname.rfind('.')+1);
	sprintf(ext,"out%05d",icpu);
	fname.append(std::string(ext));
	return fname;
}


template< typename TreeType_, typename Real_ >
std::string data<TreeType_,Real_>::rename_info2hydro( const std::string& info )
{
	std::string amr;
	unsigned ii = info.rfind("info");
	amr = info.substr(0,ii)+"poisson" + info.substr(ii+4, 6) + ".out00001";
	return amr;
}


template< typename TreeType_, typename Real_ >
std::string data<TreeType_,Real_>::rename_amr2hydro( const std::string& info )
{
	std::string amr;
	unsigned ii = info.rfind("amr");
	amr = info.substr(0,ii)+"poisson" + info.substr(ii+3, 6) + ".out00001";
	return amr;
}


} // namespace POISSON
} // namespace RAMSES


#endif // .. __RAMSES_POISSON_DATA_HH