Modify output to add read functionality that will be used to power absorption. (edd49ce7) · Commits · Cianciosa, Mark / graph_framework

graph_driver/xrays.cpp

+47 −0

Original line number	Diff line number	Diff line
		@@ -190,6 +190,48 @@ void trace_ray(const size_t num_times,
		}
		}

		//------------------------------------------------------------------------------
		/// @brief Calculate absorption.
		///
		/// @tparam T Base type of the calculation.
		/// @tparam SAFE_MATH Use safe math operations.
		//------------------------------------------------------------------------------
		template<typename T, bool SAFE_MATH=false>
		void calculate_power(const size_t num_times,
		const size_t sub_steps,
		const size_t num_rays) {
		std::vector<std::thread> threads(std::max(std::min(static_cast<unsigned int> (jit::context<T, SAFE_MATH>::max_concurrency()),
		static_cast<unsigned int> (num_rays)),
		static_cast<unsigned int> (1)));

		const size_t batch = num_rays/threads.size();
		const size_t extra = num_rays%threads.size();

		for (size_t i = 0, ie = threads.size(); i < ie; i++) {
		threads[i] = std::thread([num_times, sub_steps, num_rays, batch, extra] (const size_t thread_number,
		const size_t num_threads) -> void {

		const size_t num_steps = num_times/sub_steps;
		const size_t local_num_rays = batch
		+ (extra > thread_number ? 1 : 0);

		auto omega = graph::variable<T, SAFE_MATH> (local_num_rays, "\\omega");
		auto kx = graph::variable<T, SAFE_MATH> (local_num_rays, "k_{x}");
		auto ky = graph::variable<T, SAFE_MATH> (local_num_rays, "k_{y}");
		auto kz = graph::variable<T, SAFE_MATH> (local_num_rays, "k_{z}");
		auto x = graph::variable<T, SAFE_MATH> (local_num_rays, "x");
		auto y = graph::variable<T, SAFE_MATH> (local_num_rays, "y");
		auto z = graph::variable<T, SAFE_MATH> (local_num_rays, "z");
		auto t = graph::variable<T, SAFE_MATH> (local_num_rays, "t");

		}, i, threads.size());
		}

		for (std::thread &t : threads) {
		t.join();
		}
		}

		//------------------------------------------------------------------------------
		/// @brief Main program of the driver.
		///
		@@ -206,7 +248,12 @@ int main(int argc, const char * argv[]) {
		const size_t sub_steps = 100;
		const size_t num_rays = 100000;

		const bool use_safe_math = true;

		trace_ray<double> (num_times, sub_steps, num_rays);
		calculate_power<std::complex<double>, use_safe_math> (num_times,
		sub_steps,
		num_rays);

		std::cout << std::endl << "Timing:" << std::endl;
		total.print();

graph_framework.xcodeproj/project.pbxproj

+2 −0

Original line number	Diff line number	Diff line
		@@ -272,6 +272,7 @@
		/* End PBXCopyFilesBuildPhase section */

		/* Begin PBXFileReference section */
		C70AA9792AB48FE2002DAF1C /* absorption.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = absorption.hpp; sourceTree = "<group>"; };
		C70B705629F4F86A00098AA0 /* piecewise.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = piecewise.hpp; sourceTree = "<group>"; };
		C70D93132A30FF4E006A4227 /* special_functions.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = special_functions.hpp; sourceTree = "<group>"; };
		C713425B29413F0500672AD4 /* jit.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = jit.hpp; sourceTree = "<group>"; };
		@@ -518,6 +519,7 @@
		C723210222DC0D0A006BBF13 /* arithmetic.hpp */,
		C77E6DF522DD64E700469621 /* trigonometry.hpp */,
		C71C1FF527F5379D006997C2 /* dispersion.hpp */,
		C70AA9792AB48FE2002DAF1C /* absorption.hpp */,
		C71C1FF627F5F5A8006997C2 /* solver.hpp */,
		C7B677D829E45C9500D3ADC6 /* backend.hpp */,
		C70D93132A30FF4E006A4227 /* special_functions.hpp */,

graph_framework/absorption.hpp

0 → 100644

+25 −0

Original line number	Diff line number	Diff line
		//------------------------------------------------------------------------------
		/// @file absoprtion.hpp
		/// @brief Base class for a dispersion relation.
		///
		/// Defines functions for computing power absorbtion.
		//------------------------------------------------------------------------------

		#ifndef absorption_h
		#define absorption_h

		namespace absorption {
		//******************************************************************************
		// Root finder.
		//******************************************************************************
		//------------------------------------------------------------------------------
		/// @brief Class interface for the root finder.
		///
		/// @tparam DISPERSION_FUNCTION Class of dispersion function to use.
		//------------------------------------------------------------------------------
		template<class DISPERSION_FUNCTION>
		class root_finder {
		};
		}

		#endif /* absorption_h */

graph_framework/dispersion.hpp

+1 −2

Original line number	Diff line number	Diff line
		@@ -52,8 +52,7 @@ namespace dispersion {
		//------------------------------------------------------------------------------
		/// @brief Class interface to build dispersion relation functions.
		///
		/// @tparam T Base type of the calculation.
		/// @tparam SAFE_MATH Use safe math operations.
		/// @tparam DISPERSION_FUNCTION Class of dispersion function to use.
		//------------------------------------------------------------------------------
		template<class DISPERSION_FUNCTION>
		class dispersion_interface {

graph_framework/output.hpp

+270 −56

Original line number	Diff line number	Diff line
		@@ -18,10 +18,7 @@ namespace output {

		//------------------------------------------------------------------------------
		/// @brief Class representing a netcdf based output file.
		///
		/// @tparam T Base type of the calculation.
		//------------------------------------------------------------------------------
		template<typename T>
		class result_file {
		private:
		/// Netcdf file id.
		@@ -30,26 +27,12 @@ namespace output {
		int unlimited_dim;
		/// Number of rays dimension.
		int num_rays_dim;
		/// Dimension of ray. 1 for real, 2 for complex.
		int ray_dim;
		/// Number of rays.
		const size_t num_rays;

		//------------------------------------------------------------------------------
		/// @brief Struct map variables to a gpu buffer.
		//------------------------------------------------------------------------------
		struct variable {
		/// Variable id.
		int id;
		/// Pointer to the gpu buffer.
		T *buffer;
		};
		/// Variable list.
		std::vector<variable> variables;
		size_t num_rays;

		public:
		//------------------------------------------------------------------------------
		/// @brief Construct a result file.
		/// @brief Construct a new result file.
		///
		/// @params[in] filename Name of the result file.
		/// @params[in] num_rays Number of rays.
		@@ -64,13 +47,23 @@ namespace output {

		nc_def_dim(ncid, "time", NC_UNLIMITED, &unlimited_dim);
		nc_def_dim(ncid, "num_rays", num_rays, &num_rays_dim);
		if constexpr (jit::is_complex<T> ()) {
		nc_def_dim(ncid, "ray_dim", 2, &ray_dim);
		nc_def_dim(ncid, "num_rays", num_rays*2, &num_rays_dim);
		} else {
		nc_def_dim(ncid, "ray_dim", 1, &ray_dim);
		nc_def_dim(ncid, "num_rays", num_rays*1, &num_rays_dim);
		sync.unlock();
		}

		//------------------------------------------------------------------------------
		/// @brief Open a new result file.
		///
		/// @params[in] filename Name of the result file.
		//------------------------------------------------------------------------------
		result_file(const std::string &filename="") {

		sync.lock();
		nc_open(filename.c_str(), NC_WRITE, &ncid);

		nc_inq_dimid(ncid, "time", &unlimited_dim);
		nc_inq_dimid(ncid, "num_rays", &num_rays_dim);
		nc_inq_dimlen(ncid, num_rays_dim, &num_rays);
		nc_redef(ncid);
		sync.unlock();
		}

		@@ -81,84 +74,305 @@ namespace output {
		nc_close(ncid);
		}

		//------------------------------------------------------------------------------
		/// @brief End define mode.
		//------------------------------------------------------------------------------
		void end_define_mode() const {
		sync.lock();
		nc_enddef(ncid);
		sync.unlock();
		}

		//------------------------------------------------------------------------------
		/// @brief Get ncid.
		///
		/// @returns The netcdf file id.
		//------------------------------------------------------------------------------
		int get_ncid() const {
		return ncid;
		}

		//------------------------------------------------------------------------------
		/// @brief Get the number of rays.
		///
		/// @returns The number of rays.
		//------------------------------------------------------------------------------
		size_t get_num_rays() const {
		return num_rays;
		}

		//------------------------------------------------------------------------------
		/// @brief Get the number of rays dimension.
		///
		/// @returns The number of rays dimension.
		//------------------------------------------------------------------------------
		int get_num_rays_dim() const {
		return num_rays_dim;
		}

		//------------------------------------------------------------------------------
		/// @brief Get unlimited dimension.
		///
		/// @returns The unlimited dimension.
		//------------------------------------------------------------------------------
		int get_unlimited_dim() const {
		return unlimited_dim;
		}

		//------------------------------------------------------------------------------
		/// @brief Get unlimited size.
		///
		/// @returns The size of the unlimited dimension.
		//------------------------------------------------------------------------------
		size_t get_unlimited_size() const {
		size_t size;
		sync.lock();
		nc_inq_dimlen(ncid, unlimited_dim, &size);
		sync.unlock();

		return size;
		}

		//------------------------------------------------------------------------------
		/// @brief Sync the file.
		//------------------------------------------------------------------------------
		void sync_file() const {
		sync.lock();
		nc_sync(ncid);
		sync.unlock();
		}
		};

		//------------------------------------------------------------------------------
		/// @brief Class representing a netcdf dataset.
		///
		/// @tparam T Base type of the calculation.
		//------------------------------------------------------------------------------
		template<typename T>
		class data_set {
		private:
		/// Dimension of ray. 1 for real, 2 for complex.
		int ray_dim;
		/// Dimension of a data item.
		std::array<int, 3> dims;
		/// Data sizes.
		std::array<size_t, 3> count;
		/// Get the ray dimension size.
		const size_t ray_dim_size = 1 + jit::is_complex<T> ();
		/// The NetCDF type.
		const nc_type type = jit::is_float<T> () ? NC_FLOAT : NC_DOUBLE;

		//------------------------------------------------------------------------------
		/// @brief Struct to map variables to a gpu buffer.
		//------------------------------------------------------------------------------
		struct variable {
		/// Variable id.
		int id;
		/// Pointer to the gpu buffer.
		T *buffer;
		};
		/// Variable list.
		std::vector<variable> variables;

		//------------------------------------------------------------------------------
		/// @brief Struct to map references to a gpu buffer.
		//------------------------------------------------------------------------------
		struct reference {
		/// Variable id.
		int id;
		/// Pointer to the gpu buffer.
		T *buffer;
		/// Count stride.
		size_t ray_dim_size;
		/// Stride length.
		std::ptrdiff_t stride;
		};
		/// References list.
		std::vector<reference> references;

		public:
		//------------------------------------------------------------------------------
		/// @brief Construct a dataset.
		///
		/// @params[in] result A result file reference.
		//------------------------------------------------------------------------------
		data_set(const result_file &result) {
		sync.lock();
		if constexpr (jit::is_complex<T> ()) {
		if (NC_NOERR != nc_inq_dimid(result.get_ncid(), "ray_dim_cplx", &ray_dim)) {
		nc_def_dim(result.get_ncid(), "ray_dim_cplx", ray_dim_size, &ray_dim);
		}
		} else {
		if (NC_NOERR != nc_inq_dimid(result.get_ncid(), "ray_dim", &ray_dim)) {
		nc_def_dim(result.get_ncid(), "ray_dim", ray_dim_size, &ray_dim);
		}
		}
		sync.unlock();

		dims = {
		result.get_unlimited_dim(),
		result.get_num_rays_dim(),
		ray_dim
		};

		count = {
		1,
		result.get_num_rays(),
		ray_dim_size
		};
		}

		//------------------------------------------------------------------------------
		/// @brief Create a variable.
		///
		/// @tparam SAFE_MATH Use safe math operations.
		///
		/// @params[in] result A result file reference.
		/// @params[in] name Name of the variable.
		/// @params[in] node Node to create variable for.
		/// @params[in] context Context for the gpu.
		//------------------------------------------------------------------------------
		template<bool SAFE_MATH=false>
		void create_variable(const std::string &name,
		void create_variable(const result_file &result,
		const std::string &name,
		graph::shared_leaf<T, SAFE_MATH> &node,
		jit::context<T, SAFE_MATH> &context) {
		variable var;
		const std::array<int, 3> dims = {unlimited_dim, num_rays_dim, ray_dim};
		sync.lock();
		if constexpr (jit::is_float<T> ()) {
		nc_def_var(ncid, name.c_str(), NC_FLOAT, dims.size(),
		dims.data(), &var.id);
		} else {
		nc_def_var(ncid, name.c_str(), NC_DOUBLE, dims.size(),
		dims.data(), &var.id);
		}
		nc_def_var(result.get_ncid(), name.c_str(), type,
		static_cast<int> (dims.size()), dims.data(),
		&var.id);
		sync.unlock();

		var.buffer = context.get_buffer(node);

		variables.push_back(var);
		}

		//------------------------------------------------------------------------------
		/// @brief End define mode.
		/// @brief Load reference.
		///
		/// @tparam SAFE_MATH Use safe math operations.
		///
		/// @params[in] result A result file reference.
		/// @params[in] name Name of the variable.
		/// @params[in] node Node to create variable for.
		/// @params[in] context Context for the gpu.
		//------------------------------------------------------------------------------
		void end_define_mode() const {
		template<bool SAFE_MATH=false>
		void reference_variable(const result_file &result,
		const std::string &name,
		graph::shared_leaf<T, SAFE_MATH> &node,
		jit::context<T, SAFE_MATH> &context) {
		reference ref;
		nc_type type;
		std::array<int, 3> ref_dims;
		size_t ref_dim_size;

		sync.lock();
		nc_enddef(ncid);
		nc_inq_varid(result.get_ncid(), name, &ref.id);
		nc_inq_var(result.get_ncid(), ref.id, NULL, &type,
		NULL, ref_dims.data());
		nc_inq_dimlen(result.get_ncid(), ref.id, &ref.ray_dim_size);
		sync.unlock();

		assert(ref_dim_size <= ray_dim_size &&
		"Context variable too small to read reference.");

		ref.stride = ref.ray_dim_size < ray_dim_size ? 2 : 1;
		ref.buffer = context.get_buffer(node);
		references.push_back(ref);
		}

		//------------------------------------------------------------------------------
		/// @brief Write step.
		///
		/// @params[in] result A result file reference.
		//------------------------------------------------------------------------------
		void write() {
		size_t size;
		sync.lock();
		nc_inq_dimlen(ncid, unlimited_dim, &size);
		sync.unlock();
		const std::array<size_t, 3> start = {size, 0, 0};
		void write(const result_file &result) {
		const std::array<size_t, 3> start = {
		result.get_unlimited_size(), 0, 0
		};

		for (variable &var : variables) {
		sync.lock();
		if constexpr (jit::is_float<T> ()) {
		if constexpr (jit::is_complex<T> ()) {
		const std::array<size_t, 3> count = {1, num_rays, 2};
		nc_put_vara_float(ncid, var.id, start.data(), count.data(),
		nc_put_vara_float(result.get_ncid(), var.id,
		start.data(), count.data(),
		reinterpret_cast<float *> (var.buffer));
		} else {
		const std::array<size_t, 3> count = {1, num_rays, 1};
		nc_put_vara_float(ncid, var.id, start.data(), count.data(),
		nc_put_vara_float(result.get_ncid(), var.id,
		start.data(), count.data(),
		var.buffer);
		}
		} else {
		if constexpr (jit::is_complex<T> ()) {
		const std::array<size_t, 3> count = {1, num_rays, 2};
		nc_put_vara_double(ncid, var.id, start.data(), count.data(),
		nc_put_vara_double(result.get_ncid(), var.id,
		start.data(), count.data(),
		reinterpret_cast<double *> (var.buffer));
		} else {
		const std::array<size_t, 3> count = {1, num_rays, 1};
		nc_put_vara_double(ncid, var.id, start.data(), count.data(),
		nc_put_vara_double(result.get_ncid(), var.id,
		start.data(), count.data(),
		var.buffer);
		}
		}
		sync.unlock();
		}

		result.sync_file();
		}

		//------------------------------------------------------------------------------
		/// @brief Read step.
		///
		/// @params[in] result A result file reference.
		//------------------------------------------------------------------------------
		void read(const result_file &result,
		const size_t index) {
		const std::array<size_t, 3> ref_start = {
		index, 0, 0
		};

		for (reference &ref : references) {
		const std::array<size_t, 3> ref_count = {
		1,
		result.get_num_rays(),
		ref.ray_dim_size
		};
		const std::array<std::ptrdiff_t, 3> stride = {
		1, 1, ref.stride
		};
		const std::array<std::ptrdiff_t, 3> map = {
		0, 0, 0
		};

		sync.lock();
		nc_sync(ncid);
		if constexpr (jit::is_float<T> ()) {
		if constexpr (jit::is_complex<T> ()) {
		nc_get_varm_float(result.get_ncid(), ref.id,
		ref_start.data(), ref_count.data(),
		stride.data(), map.data(),
		reinterpret_cast<float *> (ref.buffer));
		} else {
		nc_get_varm_float(result.get_ncid(), ref.id,
		ref_start.data(), ref_count.data(),
		stride.data(), map.data(), ref.buffer);
		}
		} else {
		if constexpr (jit::is_complex<T> ()) {
		nc_get_varm_double(result.get_ncid(), ref.id,
		ref_start.data(), ref_count.data(),
		stride.data(), map.data(),
		reinterpret_cast<double *> (ref.buffer));
		} else {
		nc_get_varm_double(result.get_ncid(), ref.id,
		ref_start.data(), ref_count.data(),
		stride.data(), map.data(), ref.buffer);
		}
		}
		sync.unlock();
		}
		}
		};
		}