Merge branch 'powfix' into 'main' (77d02d39) · Commits · Cianciosa, Mark / graph_framework

CMakeLists.txt

+1 −0

Original line number	Diff line number	Diff line
		@@ -328,6 +328,7 @@ endmacro ()
		add_subdirectory (graph_driver)
		add_subdirectory (graph_benchmark)
		add_subdirectory (graph_playground)
		add_subdirectory (graph_korc)

		#-------------------------------------------------------------------------------
		# Define macro function to register tests.

graph_framework.xcodeproj/project.pbxproj

+175 −0

Original line number	Diff line number	Diff line
		@@ -1374,6 +1374,93 @@
		);
		LOCALIZATION_PREFERS_STRING_CATALOGS = YES;
		MACOSX_DEPLOYMENT_TARGET = 14.6;
		OTHER_LDFLAGS = (
		"-lnetcdf",
		"-ld_classic",
		"-L/Users/m4c/Projects/graph_framework/build/_deps/llvm-build/lib",
		"-lz",
		"-lLLVMCoverage",
		"-lLLVMSupport",
		"-lLLVMDebugInfoCodeView",
		"-lLLVMRemarks",
		"-lLLVMJITLink",
		"-lLLVMLinker",
		"-lLLVMTextAPI",
		"-lLLVMRuntimeDyld",
		"-lLLVMOrcShared",
		"-lLLVMOrcDebugging",
		"-lLLVMOrcTargetProcess",
		"-lLLVMOrcJIT",
		"-lLLVMHipStdPar",
		"-lLLVMAggressiveInstCombine",
		"-lLLVMVectorize",
		"-lLLVMAsmParser",
		"-lLLVMOption",
		"-lLLVMLTO",
		"-lLLVMObject",
		"-lLLVMWindowsDriver",
		"-lLLVMDemangle",
		"-lLLVMIRReader",
		"-lLLVMIRPrinter",
		"-lLLVMInstCombine",
		"-lLLVMBinaryFormat",
		"-lLLVMCoroutines",
		"-lLLVMBitstreamReader",
		"-lLLVMBitReader",
		"-lLLVMBitWriter",
		"-lLLVMDebugInfoDWARF",
		"-lLLVMInstrumentation",
		"-lLLVMCFGuard",
		"-lLLVMObjCARCOpts",
		"-lLLVMipo",
		"-lLLVMGlobalISel",
		"-lLLVMExecutionEngine",
		"-lLLVMFrontendDriver",
		"-lLLVMFrontendHLSL",
		"-lLLVMFrontendOpenMP",
		"-lLLVMFrontendOffloading",
		"-lLLVMSelectionDAG",
		"-lLLVMProfileData",
		"-lLLVMAnalysis",
		"-lLLVMScalarOpts",
		"-lLLVMCodeGenTypes",
		"-lLLVMCodeGenData",
		"-lLLVMCodeGen",
		"-lLLVMTargetParser",
		"-lLLVMScalarOpts",
		"-lLLVMTarget",
		"-lLLVMTransformUtils",
		"-lLLVMPasses",
		"-lLLVMSupport",
		"-lLLVMMCParser",
		"-lLLVMMC",
		"-lLLVMCore",
		"-lLLVMAsmPrinter",
		"-lLLVMAArch64Utils",
		"-lLLVMAArch64Info",
		"-lLLVMAArch64Desc",
		"-lLLVMAArch64AsmParser",
		"-lLLVMAArch64CodeGen",
		"-lLLVMSandboxIR",
		"-lLLVMFrontendAtomic",
		"-lLLVMCGData",
		"-lclangFrontend",
		"-lclangBasic",
		"-lclangEdit",
		"-lclangLex",
		"-lclangDriver",
		"-lclangSerialization",
		"-lclangAST",
		"-lclangSema",
		"-lclangAnalysis",
		"-lclangASTMatchers",
		"-lclangSupport",
		"-lclangParse",
		"-lclangAPINotes",
		"-lclangCodeGen",
		"-rpath",
		/usr/local/lib,
		);
		PRODUCT_NAME = "$(TARGET_NAME)";
		};
		name = Debug;
		@@ -1388,6 +1475,93 @@
		GCC_C_LANGUAGE_STANDARD = gnu17;
		LOCALIZATION_PREFERS_STRING_CATALOGS = YES;
		MACOSX_DEPLOYMENT_TARGET = 14.6;
		OTHER_LDFLAGS = (
		"-lnetcdf",
		"-ld_classic",
		"-L/Users/m4c/Projects/graph_framework/build/_deps/llvm-build/lib",
		"-lz",
		"-lLLVMCoverage",
		"-lLLVMSupport",
		"-lLLVMDebugInfoCodeView",
		"-lLLVMRemarks",
		"-lLLVMJITLink",
		"-lLLVMLinker",
		"-lLLVMTextAPI",
		"-lLLVMRuntimeDyld",
		"-lLLVMOrcShared",
		"-lLLVMOrcDebugging",
		"-lLLVMOrcTargetProcess",
		"-lLLVMOrcJIT",
		"-lLLVMHipStdPar",
		"-lLLVMAggressiveInstCombine",
		"-lLLVMVectorize",
		"-lLLVMAsmParser",
		"-lLLVMOption",
		"-lLLVMLTO",
		"-lLLVMObject",
		"-lLLVMWindowsDriver",
		"-lLLVMDemangle",
		"-lLLVMIRReader",
		"-lLLVMIRPrinter",
		"-lLLVMInstCombine",
		"-lLLVMBinaryFormat",
		"-lLLVMCoroutines",
		"-lLLVMBitstreamReader",
		"-lLLVMBitReader",
		"-lLLVMBitWriter",
		"-lLLVMDebugInfoDWARF",
		"-lLLVMInstrumentation",
		"-lLLVMCFGuard",
		"-lLLVMObjCARCOpts",
		"-lLLVMipo",
		"-lLLVMGlobalISel",
		"-lLLVMExecutionEngine",
		"-lLLVMFrontendDriver",
		"-lLLVMFrontendHLSL",
		"-lLLVMFrontendOpenMP",
		"-lLLVMFrontendOffloading",
		"-lLLVMSelectionDAG",
		"-lLLVMProfileData",
		"-lLLVMAnalysis",
		"-lLLVMScalarOpts",
		"-lLLVMCodeGenTypes",
		"-lLLVMCodeGenData",
		"-lLLVMCodeGen",
		"-lLLVMTargetParser",
		"-lLLVMScalarOpts",
		"-lLLVMTarget",
		"-lLLVMTransformUtils",
		"-lLLVMPasses",
		"-lLLVMSupport",
		"-lLLVMMCParser",
		"-lLLVMMC",
		"-lLLVMCore",
		"-lLLVMAsmPrinter",
		"-lLLVMAArch64Utils",
		"-lLLVMAArch64Info",
		"-lLLVMAArch64Desc",
		"-lLLVMAArch64AsmParser",
		"-lLLVMAArch64CodeGen",
		"-lLLVMSandboxIR",
		"-lLLVMFrontendAtomic",
		"-lLLVMCGData",
		"-lclangFrontend",
		"-lclangBasic",
		"-lclangEdit",
		"-lclangLex",
		"-lclangDriver",
		"-lclangSerialization",
		"-lclangAST",
		"-lclangSema",
		"-lclangAnalysis",
		"-lclangASTMatchers",
		"-lclangSupport",
		"-lclangParse",
		"-lclangAPINotes",
		"-lclangCodeGen",
		"-rpath",
		/usr/local/lib,
		);
		PRODUCT_NAME = "$(TARGET_NAME)";
		};
		name = Release;
		@@ -1427,6 +1601,7 @@
		GCC_PREPROCESSOR_DEFINITIONS = (
		"DEBUG=1",
		"$(inherited)",
		USE_INPUT_CACHE,
		);
		MACOSX_DEPLOYMENT_TARGET = 13.3;
		OTHER_LDFLAGS = (

graph_framework/arithmetic.hpp

+32 −3

Original line number	Diff line number	Diff line
		@@ -490,12 +490,41 @@ namespace graph {
		}
		}

		// Handle cases like:
		auto pl = pow_cast(this->left);
		auto pr = pow_cast(this->right);

		// (ab)^c + (ad)^c -> a^c*(b^c + d^c)
		// (ba)^c + (ad)^c -> a^c*(b^c + d^c)
		// (ab)^c + (da)^c -> a^c*(b^c + d^c)
		// (ba)^c + (da)^c -> a^c*(b^c + d^c)
		if (pl.get() && pr.get() &&
		pl->get_right()->is_match(pr->get_right())) {
		auto plm = multiply_cast(pl->get_left());
		auto prm = multiply_cast(pr->get_left());
		if (plm.get() && prm.get()) {
		if (plm->get_left()->is_match(prm->get_left())) {
		return pow(plm->get_left(), pl->get_right())*
		(pow(plm->get_right(), pl->get_right()) +
		pow(prm->get_right(), pl->get_right()));
		} else if (plm->get_left()->is_match(prm->get_right())) {
		return pow(plm->get_left(), pl->get_right())*
		(pow(plm->get_right(), pl->get_right()) +
		pow(prm->get_left(), pl->get_right()));
		} else if (plm->get_right()->is_match(prm->get_left())) {
		return pow(plm->get_right(), pl->get_right())*
		(pow(plm->get_left(), pl->get_right()) +
		pow(prm->get_right(), pl->get_right()));
		} else if (plm->get_right()->is_match(prm->get_right())) {
		return pow(plm->get_right(), pl->get_right())*
		(pow(plm->get_left(), pl->get_right()) +
		pow(prm->get_left(), pl->get_right()));
		}
		}
		}

		// (a/y)^e + b/y^e -> (a^2 + b)/(y^e)
		// b/y^e + (a/y)^e -> (b + a^2)/(y^e)
		// (a/y)^e + (b/y)^e -> (a^2 + b^2)/(y^e)
		auto pl = pow_cast(this->left);
		auto pr = pow_cast(this->right);
		if (pl.get() && rd.get()) {
		auto rdp = pow_cast(rd->get_right());
		if (rdp.get() && pl->get_right()->is_match(rdp->get_right())) {

graph_framework/equilibrium.hpp

+208 −61

Original line number	Diff line number	Diff line
		@@ -17,6 +17,7 @@
		#include "piecewise.hpp"
		#include "math.hpp"
		#include "arithmetic.hpp"
		#include "newton.hpp"

		namespace equilibrium {
		/// Lock to syncronize netcdf accross threads.
		@@ -156,6 +157,17 @@ namespace equilibrium {
		graph::shared_leaf<T, SAFE_MATH> y,
		graph::shared_leaf<T, SAFE_MATH> z) = 0;

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// The characteristic field is equilibrium dependent.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) = 0;

		//------------------------------------------------------------------------------
		/// @brief Get the contravariant basis vector in the x1 direction.
		///
		@@ -358,6 +370,19 @@ namespace equilibrium {
		auto zero = graph::zero<T, SAFE_MATH> ();
		return graph::vector(zero, zero, zero);
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// To avoid divide by zeros use the value of 1.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		return graph::one<T, SAFE_MATH> ();
		}
		};

		//------------------------------------------------------------------------------
		@@ -469,6 +494,19 @@ namespace equilibrium {
		graph::shared_leaf<T, SAFE_MATH> z) final {
		return graph::vector(0.0, 0.0, 0.1*x + 1.0);
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		return graph::one<T, SAFE_MATH> ();
		}
		};

		//------------------------------------------------------------------------------
		@@ -585,6 +623,19 @@ namespace equilibrium {
		auto zero = graph::zero<T, SAFE_MATH> ();
		return graph::vector(zero, zero, graph::one<T, SAFE_MATH> ());
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		return graph::one<T, SAFE_MATH> ();
		}
		};

		//------------------------------------------------------------------------------
		@@ -700,6 +751,19 @@ namespace equilibrium {
		graph::shared_leaf<T, SAFE_MATH> z) final {
		return graph::vector(0.0, 0.0, 0.01*x + 1.0);
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		return graph::one<T, SAFE_MATH> ();
		}
		};

		//------------------------------------------------------------------------------
		@@ -813,6 +877,19 @@ namespace equilibrium {
		auto zero = graph::zero<T, SAFE_MATH> ();
		return graph::vector(graph::one<T, SAFE_MATH> (), zero, zero);
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		return graph::one<T, SAFE_MATH> ();
		}
		};

		//------------------------------------------------------------------------------
		@@ -959,30 +1036,19 @@ namespace equilibrium {
		/// Cached magnetic field vector.
		graph::shared_vector<T, SAFE_MATH> b_cache;

		/// Cached magnetic field vector.
		graph::shared_leaf<T, SAFE_MATH> psi_norm_cache;

		//------------------------------------------------------------------------------
		/// @brief Set cache values.
		/// @brief Build psi.
		///
		/// Sets the cached values if x and y do not match.
		///
		/// @param[in] x X position.
		/// @param[in] y Y position.
		/// @param[in] z Z position.
		/// @param[in] r_norm The normalized radial position.
		/// @param[in] z_norm The normalized z position.
		/// @returns The psi value.
		//------------------------------------------------------------------------------
		void set_cache(graph::shared_leaf<T, SAFE_MATH> x,
		graph::shared_leaf<T, SAFE_MATH> y,
		graph::shared_leaf<T, SAFE_MATH> z) {
		if (!x->is_match(x_cache) \|\|
		!y->is_match(y_cache) \|\|
		!z->is_match(z_cache)) {
		x_cache = x;
		y_cache = y;
		z_cache = z;

		auto r = graph::sqrt(xx + yy);

		auto r_norm = (r - rmin)/dr;
		auto z_norm = (z - zmin)/dz;

		graph::shared_leaf<T, SAFE_MATH>
		build_psi(graph::shared_leaf<T, SAFE_MATH> r_norm,
		graph::shared_leaf<T, SAFE_MATH> z_norm) {
		auto c00_temp = graph::piecewise_2D(c00, num_cols, r_norm, z_norm);
		auto c01_temp = graph::piecewise_2D(c01, num_cols, r_norm, z_norm);
		auto c02_temp = graph::piecewise_2D(c02, num_cols, r_norm, z_norm);
		@@ -1003,7 +1069,7 @@ namespace equilibrium {
		auto c32_temp = graph::piecewise_2D(c32, num_cols, r_norm, z_norm);
		auto c33_temp = graph::piecewise_2D(c33, num_cols, r_norm, z_norm);

		auto psi = c00_temp
		return c00_temp
		+ c01_temp*z_norm
		+ c02_temp(z_normz_norm)
		+ c03_temp(z_normz_norm*z_norm)
		@@ -1019,38 +1085,63 @@ namespace equilibrium {
		+ c31_temp(r_normr_normr_norm)z_norm
		+ c32_temp(r_normr_normr_norm)(z_norm*z_norm)
		+ c33_temp(r_normr_normr_norm)(z_normz_normz_norm);
		}

		auto psi_norm = (psi - psimin)/dpsi;
		//------------------------------------------------------------------------------
		/// @brief Set cache values.
		///
		/// Sets the cached values if x and y do not match.
		///
		/// @param[in] x X position.
		/// @param[in] y Y position.
		/// @param[in] z Z position.
		//------------------------------------------------------------------------------
		void set_cache(graph::shared_leaf<T, SAFE_MATH> x,
		graph::shared_leaf<T, SAFE_MATH> y,
		graph::shared_leaf<T, SAFE_MATH> z) {
		if (!x->is_match(x_cache) \|\|
		!y->is_match(y_cache) \|\|
		!z->is_match(z_cache)) {
		x_cache = x;
		y_cache = y;
		z_cache = z;

		auto r = graph::sqrt(xx + yy);
		auto r_norm = (r - rmin)/dr;
		auto z_norm = (z - zmin)/dz;

		auto n0_temp = graph::piecewise_1D(ne_c0, psi_norm);
		auto n1_temp = graph::piecewise_1D(ne_c1, psi_norm);
		auto n2_temp = graph::piecewise_1D(ne_c2, psi_norm);
		auto n3_temp = graph::piecewise_1D(ne_c3, psi_norm);
		auto psi = build_psi(r_norm, z_norm);
		psi_norm_cache = (psi - psimin)/dpsi;

		auto n0_temp = graph::piecewise_1D(ne_c0, psi_norm_cache);
		auto n1_temp = graph::piecewise_1D(ne_c1, psi_norm_cache);
		auto n2_temp = graph::piecewise_1D(ne_c2, psi_norm_cache);
		auto n3_temp = graph::piecewise_1D(ne_c3, psi_norm_cache);

		ne_cache = ne_scale*(n0_temp +
		n1_temp*psi_norm +
		n2_temppsi_normpsi_norm +
		n3_temppsi_normpsi_norm*psi_norm);
		n1_temp*psi_norm_cache +
		n2_temppsi_norm_cachepsi_norm_cache +
		n3_temppsi_norm_cachepsi_norm_cache*psi_norm_cache);

		auto t0_temp = graph::piecewise_1D(te_c0, psi_norm);
		auto t1_temp = graph::piecewise_1D(te_c1, psi_norm);
		auto t2_temp = graph::piecewise_1D(te_c2, psi_norm);
		auto t3_temp = graph::piecewise_1D(te_c3, psi_norm);
		auto t0_temp = graph::piecewise_1D(te_c0, psi_norm_cache);
		auto t1_temp = graph::piecewise_1D(te_c1, psi_norm_cache);
		auto t2_temp = graph::piecewise_1D(te_c2, psi_norm_cache);
		auto t3_temp = graph::piecewise_1D(te_c3, psi_norm_cache);

		te_cache = te_scale*(t0_temp +
		t1_temp*psi_norm +
		t2_temppsi_normpsi_norm +
		t3_temppsi_normpsi_norm*psi_norm);
		t1_temp*psi_norm_cache +
		t2_temppsi_norm_cachepsi_norm_cache +
		t3_temppsi_norm_cachepsi_norm_cache*psi_norm_cache);

		auto p0_temp = graph::piecewise_1D(pres_c0, psi_norm);
		auto p1_temp = graph::piecewise_1D(pres_c1, psi_norm);
		auto p2_temp = graph::piecewise_1D(pres_c2, psi_norm);
		auto p3_temp = graph::piecewise_1D(pres_c3, psi_norm);
		auto p0_temp = graph::piecewise_1D(pres_c0, psi_norm_cache);
		auto p1_temp = graph::piecewise_1D(pres_c1, psi_norm_cache);
		auto p2_temp = graph::piecewise_1D(pres_c2, psi_norm_cache);
		auto p3_temp = graph::piecewise_1D(pres_c3, psi_norm_cache);

		auto pressure = pres_scale*(p0_temp +
		p1_temp*psi_norm +
		p2_temppsi_normpsi_norm +
		p3_temppsi_normpsi_norm*psi_norm);
		p1_temp*psi_norm_cache +
		p2_temppsi_norm_cachepsi_norm_cache +
		p3_temppsi_norm_cachepsi_norm_cache*psi_norm_cache);

		auto q = graph::constant<T, SAFE_MATH> (static_cast<T> (1.60218E-19));

		@@ -1271,6 +1362,45 @@ namespace equilibrium {
		set_cache(x, y, z);
		return b_cache;
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		auto x_axis = graph::variable<T, SAFE_MATH> (1, "x");
		auto y_axis = graph::variable<T, SAFE_MATH> (1, "y");
		auto z_axis = graph::variable<T, SAFE_MATH> (1, "z");
		x_axis->set(static_cast<T> (1.7));
		y_axis->set(static_cast<T> (0.0));
		z_axis->set(static_cast<T> (0.0));
		auto b_vec = get_magnetic_field(x_axis, y_axis, z_axis);
		auto b_mod = b_vec->length();

		graph::input_nodes<T, SAFE_MATH> inputs {
		graph::variable_cast(x_axis),
		graph::variable_cast(y_axis),
		graph::variable_cast(z_axis)
		};

		workflow::manager<T, SAFE_MATH> work(device_number);
		solver::newton(work, {
		x_axis, z_axis
		}, inputs, psi_norm_cache, static_cast<T> (1.0E-30), 1000, static_cast<T> (0.1));
		work.add_item(inputs, {b_mod}, {}, "bmod_at_axis");
		work.compile();
		work.run();

		T result;
		work.copy_to_host(b_mod, &result);

		return graph::constant<T, SAFE_MATH> (result);
		}
		};

		//------------------------------------------------------------------------------
		@@ -2014,6 +2144,23 @@ namespace equilibrium {
		return bvec_cache;
		}

		//------------------------------------------------------------------------------
		/// @brief Get the characteristic field.
		///
		/// Use the value at the y intercept.
		///
		/// @params[in] device_number Device to use.
		/// @returns The characteristic field.
		//------------------------------------------------------------------------------
		virtual graph::shared_leaf<T, SAFE_MATH>
		get_characteristic_field(const size_t device_number=0) final {
		auto s_axis = graph::zero<T, SAFE_MATH> ();
		auto u_axis = graph::zero<T, SAFE_MATH> ();
		auto v_axis = graph::zero<T, SAFE_MATH> ();
		auto b_vec = get_magnetic_field(s_axis, u_axis, v_axis);
		return b_vec->length();
		}

		//------------------------------------------------------------------------------
		/// @brief Get the x position.
		///

graph_framework/jit.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -133,7 +133,7 @@ namespace jit {

		for (auto &in : inputs) {
		if (usage.find(in.get()) == usage.end()) {
		usage[in.get()] == 0;
		usage[in.get()] = 0;
		}
		}