From 9f7b8f971588ff1152b61e634bb18cb517a2b9f6 Mon Sep 17 00:00:00 2001 From: Alan Kaptanoglu Date: Fri, 16 Aug 2024 09:18:01 -0700 Subject: [PATCH] Did an initial hooking up the c++ code for exact magnet calculations. --- CMakeLists.txt | 1 + src/simsoptpp/Bcube_nonVec.cpp | 222 ++++++++++++++++----------------- src/simsoptpp/Bcube_nonVec.h | 6 + src/simsoptpp/python.cpp | 4 + 4 files changed, 120 insertions(+), 113 deletions(-) create mode 100644 src/simsoptpp/Bcube_nonVec.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 92663e5d7..db0246e8d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -131,6 +131,7 @@ pybind11_add_module(${PROJECT_NAME} src/simsoptpp/dommaschk.cpp src/simsoptpp/reiman.cpp src/simsoptpp/tracing.cpp src/simsoptpp/magneticfield_biotsavart.cpp src/simsoptpp/python_boozermagneticfield.cpp src/simsoptpp/boozerradialinterpolant.cpp + src/simsoptpp/Bcube_nonVec.cpp ) set_target_properties(${PROJECT_NAME} diff --git a/src/simsoptpp/Bcube_nonVec.cpp b/src/simsoptpp/Bcube_nonVec.cpp index cd02ae6d0..fe8a94e2d 100644 --- a/src/simsoptpp/Bcube_nonVec.cpp +++ b/src/simsoptpp/Bcube_nonVec.cpp @@ -1,150 +1,146 @@ -#include -#include - +#include "Bcube_nonVec.h" double mu0 = 4 * M_PI * 1e-7; -double heaviside(double& x1, double& x2) { +double heaviside(double x1, double x2) { if (x1 < 0) { return 0; - }; - else if (x1 > 0) { + } else if (x1 > 0) { return 1; - }; - else { + } else { return x2; }; } -double Pd(double& phi, double& theta) { - const int rows = 3; - const int cols = 3; +// double Pd(double& phi, double& theta) { +// const int rows = 3; +// const int cols = 3; - double cosT = std::cos(theta); - double sinT = std::sin(theta); - double cosP = std::cos(phi); - double sinP = std::sin(phi); +// double cosT = std::cos(theta); +// double sinT = std::sin(theta); +// double cosP = std::cos(phi); +// double sinP = std::sin(phi); - double* P[rows][cols] = { - {cosT * cosP, -cosT * sinP, sinT}, - {sinP, cosP, 0.0}, - {-sinT * cosP, sinT * sinP, cosT} - }; - return P; -} - - -double iterate_over_corners(tuple& corner, double& x, double& y, double& z) { - int i = corner[0], j = corner[1], k = corner[2]; - int summa = std::pow(-1, i+j+k); - double rijk = std::sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]); - - double atan_xy = std::atan2(y[j]*x[i],z[k]*rijk); - atan_xz = std::atan2(z[k]*x[i],y[j]*rijk); - atan_yz = std::atan2(y[j]*z[k],x[i]*rijk); - log_x = std::log(x[i] + rijk); - log_y = std::log(y[j] + rijk); - log_z = std::log(z[k] + rijk); +// double* P[rows][cols] = { +// {cosT * cosP, -cosT * sinP, sinT}, +// {sinP, cosP, 0.0}, +// {-sinT * cosP, sinT * sinP, cosT} +// }; +// return P; +// } + + +// double iterate_over_corners(tuple& corner, double& x, double& y, double& z) { +// int i = corner[0], j = corner[1], k = corner[2]; +// int summa = std::pow(-1, i+j+k); +// double rijk = std::sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]); + +// double atan_xy = std::atan2(y[j]*x[i],z[k]*rijk); +// atan_xz = std::atan2(z[k]*x[i],y[j]*rijk); +// atan_yz = std::atan2(y[j]*z[k],x[i]*rijk); +// log_x = std::log(x[i] + rijk); +// log_y = std::log(y[j] + rijk); +// log_z = std::log(z[k] + rijk); - const int rows = 3; - const int cols = 3; - double h[rows][cols] = { - {atan_xy + atan_xz, log_z, log_y}, - {log_z, atan_xy + atan_yz, log_x}, - {log_y, log_x, atan_xz + atan_yz} - }; +// const int rows = 3; +// const int cols = 3; +// double h[rows][cols] = { +// {atan_xy + atan_xz, log_z, log_y}, +// {log_z, atan_xy + atan_yz, log_x}, +// {log_y, log_x, atan_xz + atan_yz} +// }; - return h; -} +// return h; +// } -double Hd_i_prime(double& r, double& dims) { - const int rows = 3; - const int cols = 3; +// double Hd_i_prime(double& r, double& dims) { +// const int rows = 3; +// const int cols = 3; - double H[rows][cols] = {0}; +// double H[rows][cols] = {0}; - double xp = r[0], yp = r[1], zp = r[2]; +// double xp = r[0], yp = r[1], zp = r[2]; - double x[2] = {xp + dims[0]/2, xp - dims[0]/2}; - double y[2] = {yp + dims[1]/2, yp - dims[1]/2}; - double z[2] = {zp + dims[2]/2, zp - dims[2]/2}; -} +// double x[2] = {xp + dims[0]/2, xp - dims[0]/2}; +// double y[2] = {yp + dims[1]/2, yp - dims[1]/2}; +// double z[2] = {zp + dims[2]/2, zp - dims[2]/2}; +// } -double B_direct(double& points, double& magPos, double& M, double& dims, double& phiThetas) { - const int N = sizeof(points) / sizeof(points[0]); - const int D = sizeof(M) / sizeof(M[0]); +// double B_direct(double& points, double& magPos, double& M, double& dims, double& phiThetas) { +// const int N = sizeof(points) / sizeof(points[0]); +// const int D = sizeof(M) / sizeof(M[0]); - double B[N][3] = {0}; //find out how to do dot products, how to use multiple indeces - for (n = 0; n < N; ++n) { - for (d = 0; d < D; ++d) { - double P = Pd(phiThetas[d][0], phiThetas[d][1]); - double r_loc = P DOT (points[n] - magPos[d]); - - double tx = heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5); - double ty = heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5); - double tz = heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5); - double tm = 2*tx*ty*tz; +// double B[N][3] = {0}; //find out how to do dot products, how to use multiple indeces +// for (n = 0; n < N; ++n) { +// for (d = 0; d < D; ++d) { +// double P = Pd(phiThetas[d][0], phiThetas[d][1]); +// double r_loc = P DOT (points[n] - magPos[d]); + +// double tx = heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5); +// double ty = heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5); +// double tz = heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5); +// double tm = 2*tx*ty*tz; - B[n] += mu0 * P.T DOT (Hd_i_prime(r_loc,dims) @ (P @ M[d]) + tm*P@M[d]); //how does transpose work? - }; - }; - return B; -} +// B[n] += mu0 * P.T DOT (Hd_i_prime(r_loc,dims) @ (P @ M[d]) + tm*P@M[d]); //how does transpose work? +// }; +// }; +// return B; +// } -double Bn_direct(double& points, double& magPos, double& M, double& norms, double& dims, double& phiThetas) { - const int N = sizeof(points) / sizeof(points[0]); +// double Bn_direct(double& points, double& magPos, double& M, double& norms, double& dims, double& phiThetas) { +// const int N = sizeof(points) / sizeof(points[0]); - double B = B_direct(points, magPos, M, dims, phiThetas); - double Bn[N] = {0}; +// double B = B_direct(points, magPos, M, dims, phiThetas); +// double Bn[N] = {0}; - for (n = 0; n < N; ++n) { - Bn[n] = B[n] DOT norms[n]; - }; - return Bn; -} +// for (n = 0; n < N; ++n) { +// Bn[n] = B[n] DOT norms[n]; +// }; +// return Bn; +// } -double gd_i(r_loc, n_i_loc, dims) { - double tx = heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5); - double ty = heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5); - double tz = heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5); - double tm = 2*tx*ty*tz; +// double gd_i(r_loc, n_i_loc, dims) { +// double tx = heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5); +// double ty = heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5); +// double tz = heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5); +// double tm = 2*tx*ty*tz; - return mu0 * (Hd_i_prime(r_loc,dims).T + tm*np.eye(3)) DOT n_i_loc; //shortcut for identity matrix? -} +// return mu0 * (Hd_i_prime(r_loc,dims).T + tm*np.eye(3)) DOT n_i_loc; //shortcut for identity matrix? +// } -double Acube(double& points, double& magPos, double& norms, double& dims, double& phiThetas) { - const int N = sizeof(points) / sizeof(points[0]); - const int D = sizeof(magPos) / sizeof(magPos[0]); +// double Acube(double& points, double& magPos, double& norms, double& dims, double& phiThetas) { +// const int N = sizeof(points) / sizeof(points[0]); +// const int D = sizeof(magPos) / sizeof(magPos[0]); - double A[N][3*D] = {0}; - for (n = 0; n < N; ++n) { - for (d = 0; d < D; ++d) { - double P = Pd(phiThetas[d][0], phiThetas[d][1]); - double r_loc = P DOT (points[n] - magPos[d]); - double n_loc = P DOT norms[n]; +// double A[N][3*D] = {0}; +// for (n = 0; n < N; ++n) { +// for (d = 0; d < D; ++d) { +// double P = Pd(phiThetas[d][0], phiThetas[d][1]); +// double r_loc = P DOT (points[n] - magPos[d]); +// double n_loc = P DOT norms[n]; - double g = P.T DOT gd_i(r_loc,n_loc,dims); - A[n][3*d : 3*d + 3] = g; //can I still index this way? - } - }; - }; - - if ((N,3*D) != A.shape) - throw std::runtime_error("A shape altered during Acube"); //is there an assert function? - return A; -} - - -double Bn_fromMat(double& points, double& magPos, double& M, double& norms, double& dims, double& phiThetas) { - double A = Acube(points, magPos, norms, dims, phiThetas); - Ms = np.concatenate(M); //equiv for concatenate? - return A DOT Ms; -} +// double g = P.T DOT gd_i(r_loc,n_loc,dims); +// A[n][3*d : 3*d + 3] = g; //can I still index this way? +// } +// }; +// }; + +// if ((N,3*D) != A.shape) +// throw std::runtime_error("A shape altered during Acube"); //is there an assert function? +// return A; +// } + + +// double Bn_fromMat(double& points, double& magPos, double& M, double& norms, double& dims, double& phiThetas) { +// double A = Acube(points, magPos, norms, dims, phiThetas); +// Ms = np.concatenate(M); //equiv for concatenate? +// return A DOT Ms; +// } diff --git a/src/simsoptpp/Bcube_nonVec.h b/src/simsoptpp/Bcube_nonVec.h new file mode 100644 index 000000000..f30746803 --- /dev/null +++ b/src/simsoptpp/Bcube_nonVec.h @@ -0,0 +1,6 @@ +#include +#include +#include "xtensor-python/pyarray.hpp" // Numpy bindings +typedef xt::pyarray Array; + +double heaviside(double x1, double x2); \ No newline at end of file diff --git a/src/simsoptpp/python.cpp b/src/simsoptpp/python.cpp index 3082265fe..02d1b5228 100644 --- a/src/simsoptpp/python.cpp +++ b/src/simsoptpp/python.cpp @@ -22,6 +22,7 @@ typedef xt::pytensor PyTensor; #include "reiman.h" #include "simdhelpers.h" #include "boozerresidual_py.h" +#include "Bcube_nonVec.h" namespace py = pybind11; @@ -67,6 +68,9 @@ PYBIND11_MODULE(simsoptpp, m) { m.def("dipole_field_Bn" , &dipole_field_Bn, py::arg("points"), py::arg("m_points"), py::arg("unitnormal"), py::arg("nfp"), py::arg("stellsym"), py::arg("b"), py::arg("coordinate_flag") = "cartesian", py::arg("R0") = 0.0); m.def("define_a_uniform_cartesian_grid_between_two_toroidal_surfaces" , &define_a_uniform_cartesian_grid_between_two_toroidal_surfaces); + // Functions below are implemented for exact rectangular cube magnets + m.def("heaviside", &heaviside); + // Permanent magnet optimization algorithms have many default arguments m.def("MwPGP_algorithm", &MwPGP_algorithm, py::arg("A_obj"), py::arg("b_obj"), py::arg("ATb"), py::arg("m_proxy"), py::arg("m0"), py::arg("m_maxima"), py::arg("alpha"), py::arg("nu") = 1.0e100, py::arg("epsilon") = 1.0e-3, py::arg("reg_l0") = 0.0, py::arg("reg_l1") = 0.0, py::arg("reg_l2") = 0.0, py::arg("max_iter") = 500, py::arg("min_fb") = 1.0e-20, py::arg("verbose") = false); // variants of GPMO algorithm