From 93ede306494274e1a93d8e4508365b01a86073de Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 10 Oct 2024 13:49:54 -0700 Subject: [PATCH] Use `amrex::SmallMatrix` Replace current small matrix logic and math with the new `amrex::SmallMatrix` class. --- src/particles/ReferenceParticle.H | 2 +- .../diagnostics/CovarianceMatrixMath.H | 54 +------- .../diagnostics/EmittanceInvariants.H | 14 +- .../diagnostics/EmittanceInvariants.cpp | 56 ++++---- .../ReducedBeamCharacteristics.cpp | 82 +++++------ src/particles/elements/RFCavity.H | 91 ++++++------ src/particles/elements/SoftQuad.H | 88 ++++++------ src/particles/elements/SoftSol.H | 130 +++++++++--------- 8 files changed, 234 insertions(+), 283 deletions(-) diff --git a/src/particles/ReferenceParticle.H b/src/particles/ReferenceParticle.H index 226025b74..71df2b734 100644 --- a/src/particles/ReferenceParticle.H +++ b/src/particles/ReferenceParticle.H @@ -41,7 +41,7 @@ namespace impactx amrex::ParticleReal charge = 0.0; ///< reference charge, in C amrex::ParticleReal sedge = 0.0; ///< value of s at entrance of the current beamline element - amrex::Array2D map; ///< linearized map + amrex::SmallMatrix map; ///< linearized map (zero-indexed) /** Get reference particle relativistic gamma * diff --git a/src/particles/diagnostics/CovarianceMatrixMath.H b/src/particles/diagnostics/CovarianceMatrixMath.H index a3b0099cf..6ea712ba1 100644 --- a/src/particles/diagnostics/CovarianceMatrixMath.H +++ b/src/particles/diagnostics/CovarianceMatrixMath.H @@ -13,9 +13,9 @@ #include #include -#include #include #include +#include #include #include @@ -23,7 +23,6 @@ namespace impactx::diagnostics { - /** Function to return the roots of a cubic polynomial ax^3 + bx^2 + cx + d. * The trigonometric form of Cardano's formula is used. * This implementation expects three real roots, which is verified @@ -163,57 +162,6 @@ namespace impactx::diagnostics return roots; } - - /** Function to take the trace of a square 6x6 matrix. - * - * @param[in] A a square matrix - * @returns the trace of A - */ - amrex::ParticleReal - TraceMat ( - amrex::Array2D const & A - ) - { - int const dim = 6; - amrex::ParticleReal trA = 0.0; - - for (int i = 1; i < dim+1; i++) { - trA += A(i,i); - } - return trA; - } - - - /** Function to multiply two square matrices of dimension 6. - * - * @param[in] A a square matrix - * @param[in] B square matrix - * @returns the matrix C = AB - */ - amrex::Array2D - MultiplyMat ( - amrex::Array2D const & A, - amrex::Array2D const & B - ) - { - amrex::Array2D C; - int const dim = 6; - - for (int i = 1; i < dim+1; i++) { - for (int j = 1; j < dim+1; j++) { - C(i,j) = 0; - - for (int k = 1; k < dim+1; k++) { - C(i,j) += A(i,k) * B(k,j); - } - - } - - } - return C; - } - - } // namespace impactx::diagnostics #endif // COVARIANCE_MATRIX_MATH_H diff --git a/src/particles/diagnostics/EmittanceInvariants.H b/src/particles/diagnostics/EmittanceInvariants.H index 177371953..b27571ecc 100644 --- a/src/particles/diagnostics/EmittanceInvariants.H +++ b/src/particles/diagnostics/EmittanceInvariants.H @@ -12,10 +12,11 @@ #include "particles/ImpactXParticleContainer.H" -#include #include +#include #include +#include @@ -44,7 +45,7 @@ namespace impactx::diagnostics amrex::ParticleReal > KineticInvariants ( - amrex::Array2D const & Sigma + amrex::SmallMatrix const & Sigma ); /** Returns the three eigenemittances @@ -65,11 +66,12 @@ namespace impactx::diagnostics * @returns tuple containing eigenemittances e1, e2, and e3 */ std::tuple< - amrex::ParticleReal, - amrex::ParticleReal, - amrex::ParticleReal> + amrex::ParticleReal, + amrex::ParticleReal, + amrex::ParticleReal + > Eigenemittances ( - amrex::Array2D const & Sigma + amrex::SmallMatrix const & Sigma ); } // namespace impactx::diagnostics diff --git a/src/particles/diagnostics/EmittanceInvariants.cpp b/src/particles/diagnostics/EmittanceInvariants.cpp index 20c6df32a..52c3339a7 100644 --- a/src/particles/diagnostics/EmittanceInvariants.cpp +++ b/src/particles/diagnostics/EmittanceInvariants.cpp @@ -30,7 +30,7 @@ namespace impactx::diagnostics * calculation of the three eigenemittances. * * input - Sigma symmetric 6x6 covariance matrix - * returns - tuple containing invarants I2, I4, and I6 + * returns - tuple containing invariants I2, I4, and I6 */ std::tuple< amrex::ParticleReal, @@ -38,46 +38,45 @@ namespace impactx::diagnostics amrex::ParticleReal > KineticInvariants ( - amrex::Array2D const & Sigma + amrex::SmallMatrix const & Sigma ) { using namespace amrex::literals; - std::tuple invariants; + std::tuple invariants; amrex::ParticleReal I2 = 0.0_prt; amrex::ParticleReal I4 = 0.0_prt; amrex::ParticleReal I6 = 0.0_prt; // Intermediate matrices used for storage. - amrex::Array2D S1; - amrex::Array2D S2; - amrex::Array2D S4; - amrex::Array2D S6; + amrex::SmallMatrix S1; + amrex::SmallMatrix S2; + amrex::SmallMatrix S4; + amrex::SmallMatrix S6; // Construct the matrix S1 = Sigma*J. This is a // permutation of the columns of Sigma with // a change of sign. - for (int i = 1; i < 7; i++) { - for (int j = 1; j < 7; j++) { + for (int i = 0; i < 6; i++) { + for (int j = 0; j < 6; j++) { if (j % 2 != 0) { - S1(i,j) = -Sigma(i,j+1); // if j is odd + S1(i,j) = +Sigma(i,j-1); // if j is odd } else { - S1(i,j) = +Sigma(i,j-1); // if j is even + S1(i,j) = -Sigma(i,j+1); // if j is even } } } // Carry out necessary matrix multiplications (3 are needed). - S2 = impactx::diagnostics::MultiplyMat(S1,S1); - S4 = impactx::diagnostics::MultiplyMat(S2,S2); - S6 = impactx::diagnostics::MultiplyMat(S2,S4); + S2 = S1 * S1; + S4 = S2 * S2; + S6 = S2 * S4; // Define the three kinematic invariants (should be nonnegative). - I2 = -impactx::diagnostics::TraceMat(S2)/2.0_prt; - I4 = +impactx::diagnostics::TraceMat(S4)/2.0_prt; - I6 = -impactx::diagnostics::TraceMat(S6)/2.0_prt; - + I2 = -S2.trace() / 2.0_prt; + I4 = +S4.trace() / 2.0_prt; + I6 = -S6.trace() / 2.0_prt; invariants = std::make_tuple(I2,I4,I6); return invariants; @@ -95,11 +94,12 @@ namespace impactx::diagnostics * returns - tuple containing eigenemittances e1, e2, and e3 */ std::tuple< - amrex::ParticleReal, - amrex::ParticleReal, - amrex::ParticleReal> + amrex::ParticleReal, + amrex::ParticleReal, + amrex::ParticleReal + > Eigenemittances ( - amrex::Array2D const & Sigma + amrex::SmallMatrix const & Sigma ) { BL_PROFILE("impactx::diagnostics::Eigenemittances"); @@ -123,8 +123,8 @@ namespace impactx::diagnostics // doi:10.48550/arXiv.1305.1532. amrex::ParticleReal a = 1.0_prt; amrex::ParticleReal b = -I2; - amrex::ParticleReal c = (pow(I2,2)-I4)/2.0_prt; - amrex::ParticleReal d = -pow(I2,3)/6.0_prt + I2*I4/2.0_prt - I6/3.0_prt; + amrex::ParticleReal c = (std::pow(I2, 2) - I4) / 2.0_prt; + amrex::ParticleReal d = -std::pow(I2, 3) / 6.0_prt + I2 * I4 / 2.0_prt - I6 / 3.0_prt; // Return the cubic coefficients //std::cout << "Return a,b,c,d " << a << " " << b << " " << c << " " << d << "\n"; @@ -133,10 +133,10 @@ namespace impactx::diagnostics // Caution: The order of e1,e2,e3 should be consistent with the // order ex,ey,et in the limit of uncoupled transport. // The order below is important and has been checked. - roots = CubicRootsTrig(a,b,c,d); - amrex::ParticleReal e1 = sqrt(std::abs(std::get<1>(roots))); - amrex::ParticleReal e2 = sqrt(std::abs(std::get<2>(roots))); - amrex::ParticleReal e3 = sqrt(std::abs(std::get<0>(roots))); + roots = CubicRootsTrig(a, b, c, d); + amrex::ParticleReal e1 = std::sqrt(std::abs(std::get<1>(roots))); + amrex::ParticleReal e2 = std::sqrt(std::abs(std::get<2>(roots))); + amrex::ParticleReal e3 = std::sqrt(std::abs(std::get<0>(roots))); emittances = std::make_tuple(e1,e2,e3); return emittances; diff --git a/src/particles/diagnostics/ReducedBeamCharacteristics.cpp b/src/particles/diagnostics/ReducedBeamCharacteristics.cpp index 99a11560d..1cbdcef2b 100644 --- a/src/particles/diagnostics/ReducedBeamCharacteristics.cpp +++ b/src/particles/diagnostics/ReducedBeamCharacteristics.cpp @@ -16,10 +16,11 @@ #include // for TinyProfiler #include // for AMREX_GPU_DEVICE -#include // for ParticleReal -#include // for ReduceOps #include // for ParallelDescriptor #include // for ParticleReduce +#include // for ParticleReal +#include // for ReduceOps +#include // for SmallMatrix #include // for TypeMultiplier @@ -320,45 +321,46 @@ namespace impactx::diagnostics if (compute_eigenemittances) { // Store the covariance matrix in dynamical variables: - amrex::Array2D Sigma; - Sigma(1,1) = x_ms; - Sigma(1,2) = xpx * bg; - Sigma(1,3) = xy; - Sigma(1,4) = xpy * bg; - Sigma(1,5) = xt; - Sigma(1,6) = xpt * bg; - Sigma(2,1) = xpx * bg; - Sigma(2,2) = px_ms * bg2; - Sigma(2,3) = pxy * bg; - Sigma(2,4) = pxpy * bg2; - Sigma(2,5) = pxt * bg; - Sigma(2,6) = pxpt * bg2; - Sigma(3,1) = xy; - Sigma(3,2) = pxy * bg; - Sigma(3,3) = y_ms; - Sigma(3,4) = ypy * bg; - Sigma(3,5) = yt; - Sigma(3,6) = ypt * bg; - Sigma(4,1) = xpy * bg; - Sigma(4,2) = pxpy * bg2; - Sigma(4,3) = ypy * bg; - Sigma(4,4) = py_ms * bg2; - Sigma(4,5) = pyt * bg; - Sigma(4,6) = pypt * bg2; - Sigma(5,1) = xt; - Sigma(5,2) = pxt * bg; - Sigma(5,3) = yt; - Sigma(5,4) = pyt * bg; - Sigma(5,5) = t_ms; - Sigma(5,6) = tpt * bg; - Sigma(6,1) = xpt * bg; - Sigma(6,2) = pxpt * bg2; - Sigma(6,3) = ypt * bg; - Sigma(6,4) = pypt * bg2; - Sigma(6,5) = tpt * bg; - Sigma(6,6) = pt_ms * bg2; + // note: Sigma is zero-indexed + amrex::SmallMatrix Sigma; + Sigma(0 ,0) = x_ms; + Sigma(0, 1) = xpx * bg; + Sigma(0, 2) = xy; + Sigma(0, 3) = xpy * bg; + Sigma(0, 4) = xt; + Sigma(0, 5) = xpt * bg; + Sigma(1, 0) = xpx * bg; + Sigma(1, 1) = px_ms * bg2; + Sigma(1, 2) = pxy * bg; + Sigma(1, 3) = pxpy * bg2; + Sigma(1, 4) = pxt * bg; + Sigma(1, 5) = pxpt * bg2; + Sigma(2, 0) = xy; + Sigma(2, 1) = pxy * bg; + Sigma(2, 2) = y_ms; + Sigma(2, 3) = ypy * bg; + Sigma(2, 4) = yt; + Sigma(2, 5) = ypt * bg; + Sigma(3, 0) = xpy * bg; + Sigma(3, 1) = pxpy * bg2; + Sigma(3, 2) = ypy * bg; + Sigma(3, 3) = py_ms * bg2; + Sigma(3, 4) = pyt * bg; + Sigma(3, 5) = pypt * bg2; + Sigma(4, 0) = xt; + Sigma(4, 1) = pxt * bg; + Sigma(4, 2) = yt; + Sigma(4, 3) = pyt * bg; + Sigma(4, 4) = t_ms; + Sigma(4, 5) = tpt * bg; + Sigma(5, 0) = xpt * bg; + Sigma(5, 1) = pxpt * bg2; + Sigma(5, 2) = ypt * bg; + Sigma(5, 3) = pypt * bg2; + Sigma(5, 4) = tpt * bg; + Sigma(5, 5) = pt_ms * bg2; // Calculate eigenemittances - std::tuple emittances = Eigenemittances(Sigma); + std::tuple emittances = Eigenemittances(Sigma); emittance_1 = std::get<0>(emittances); emittance_2 = std::get<1>(emittances); emittance_3 = std::get<2>(emittances); diff --git a/src/particles/elements/RFCavity.H b/src/particles/elements/RFCavity.H index e4aae613e..d1204d284 100644 --- a/src/particles/elements/RFCavity.H +++ b/src/particles/elements/RFCavity.H @@ -20,9 +20,9 @@ #include #include -#include #include #include +#include #include #include @@ -218,8 +218,8 @@ namespace RFCavityData amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; - // get the linear map - amrex::Array2D const R = refpart.map; + // get the linear map (note: zero-indexed) + amrex::SmallMatrix const R = refpart.map; // symplectic linear map for the RF cavity is computed using the // Hamiltonian formalism as described in: @@ -228,18 +228,21 @@ namespace RFCavityData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - xout = R(1,1)*x + R(1,2)*px + R(1,3)*y - + R(1,4)*py + R(1,5)*t + R(1,6)*pt; - pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y - + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - yout = R(3,1)*x + R(3,2)*px + R(3,3)*y - + R(3,4)*py + R(3,5)*t + R(3,6)*pt; - pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y - + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - tout = R(5,1)*x + R(5,2)*px + R(5,3)*y - + R(5,4)*py + R(5,5)*t + R(5,6)*pt; - ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y - + R(6,4)*py + R(6,5)*t + R(6,6)*pt; + // note: R is zero-indexed + // clang-format off + xout = R(0, 0) * x + R(0, 1) * px + R(0, 2) * y + + R(0, 3) * py + R(0, 4) * t + R(0, 5) * pt; + pxout = R(1, 0) * x + R(1, 1) * px + R(1, 2) * y + + R(1, 3) * py + R(1, 4) * t + R(1, 5) * pt; + yout = R(2, 0) * x + R(2, 1) * px + R(2, 2) * y + + R(2, 3) * py + R(2, 4) * t + R(2, 5) * pt; + pyout = R(3, 0) * x + R(3, 1) * px + R(3, 3) * y + + R(3, 3) * py + R(3, 4) * t + R(3, 5) * pt; + tout = R(4, 0) * x + R(4, 1) * px + R(4, 2) * y + + R(4, 3) * py + R(4, 4) * t + R(4, 5) * pt; + ptout = R(5, 0) * x + R(5, 1) * px + R(5, 2) * y + + R(5, 3) * py + R(5, 4) * t + R(5, 5) * pt; + // clang-format on // assign updated values x = xout; @@ -274,14 +277,7 @@ namespace RFCavityData amrex::ParticleReal const sedge = refpart.sedge; // initialize linear map (deviation) values - for (int i=1; i<7; i++) { - for (int j=1; j<7; j++) { - if (i == j) - refpart.map(i, j) = 1.0_prt; - else - refpart.map(i, j) = 0.0_prt; - } - } + refpart.map = decltype(refpart.map)::Identity(); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -314,13 +310,13 @@ namespace RFCavityData amrex::ParticleReal scale_in = 1.0_prt; amrex::ParticleReal scale_fin = 1.0_prt; - for (int i=1; i<7; i++) { - for (int j=1; j<7; j++) { - if( i % 2 == 0) + for (int i=0; i<6; i++) { + for (int j=0; j<6; j++) { + if( i % 2 != 0) scale_fin = bgf; else scale_fin = 1.0_prt; - if( j % 2 == 0) + if( j % 2 != 0) scale_in = bgi; else scale_in = 1.0_prt; @@ -422,11 +418,12 @@ namespace RFCavityData } // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const betgam = refpart.beta_gamma(); - refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,3); - refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,3); + // note: map and R are zero-indexed + refpart.map(4, 4) = R(4, 4) + tau * R(5, 4) / std::pow(betgam, 3); + refpart.map(4, 5) = R(4, 5) + tau * R(5, 5) / std::pow(betgam, 3); } /** This pushes the reference particle and the linear map matrix @@ -462,19 +459,20 @@ namespace RFCavityData refpart.pt = pt; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const s = tau/refpart.beta_gamma(); amrex::ParticleReal const L = E0*ezp * std::sin(k*t+phi)/(2.0_prt*k); - refpart.map(1,1) = (1.0_prt-s*L)*R(1,1) + s*R(2,1); - refpart.map(1,2) = (1.0_prt-s*L)*R(1,2) + s*R(2,2); - refpart.map(2,1) = -s * std::pow(L,2)*R(1,1) + (1.0_prt+s*L)*R(2,1); - refpart.map(2,2) = -s * std::pow(L,2)*R(1,2) + (1.0_prt+s*L)*R(2,2); + // note: map and R are zero-indexed + refpart.map(0, 0) = (1.0_prt - s * L) * R(0, 0) + s * R(1, 0); + refpart.map(0, 1) = (1.0_prt - s * L) * R(0, 1) + s * R(1, 1); + refpart.map(1, 0) = -s * std::pow(L, 2) * R(0, 0) + (1.0_prt + s * L) * R(1, 0); + refpart.map(1, 1) = -s * std::pow(L, 2) * R(0, 1) + (1.0_prt + s * L) * R(1, 1); - refpart.map(3,3) = (1.0_prt-s*L)*R(3,3) + s*R(4,3); - refpart.map(3,4) = (1.0_prt-s*L)*R(3,4) + s*R(4,4); - refpart.map(4,3) = -s * std::pow(L,2)*R(3,3) + (1.0_prt+s*L)*R(4,3); - refpart.map(4,4) = -s * std::pow(L,2)*R(3,4) + (1.0_prt+s*L)*R(4,4); + refpart.map(2, 2) = (1.0_prt - s * L) * R(2, 2) + s * R(3, 2); + refpart.map(2, 3) = (1.0_prt - s * L) * R(2, 3) + s * R(3, 3); + refpart.map(3, 2) = -s * std::pow(L, 2) * R(2, 2) + (1.0_prt + s * L) * R(3, 2); + refpart.map(3, 3) = -s * std::pow(L, 2) * R(2, 3) + (1.0_prt + s * L) * R(3, 3); } /** This pushes the reference particle and the linear map matrix @@ -514,18 +512,19 @@ namespace RFCavityData refpart.pt = pt - E0*(ezintf-ezint) * std::cos(k*t+phi); // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const M = E0*(ezintf-ezint)*k * std::sin(k*t+phi); amrex::ParticleReal const L = E0*(ezpf-ezp) * std::sin(k*t+phi)/(2.0_prt*k)+M/2.0_prt; - refpart.map(2,1) = L*R(1,1) + R(2,1); - refpart.map(2,2) = L*R(1,2) + R(2,2); + // note: map and R are zero-indexed + refpart.map(1, 0) = L * R(0, 0) + R(1, 0); + refpart.map(1, 1) = L * R(0, 1) + R(1, 1); - refpart.map(4,3) = L*R(3,3) + R(4,3); - refpart.map(4,4) = L*R(3,4) + R(4,4); + refpart.map(3, 2) = L * R(2, 2) + R(3, 2); + refpart.map(3, 3) = L * R(2, 3) + R(3, 3); - refpart.map(6,5) = M*R(5,5) + R(6,5); - refpart.map(6,6) = M*R(5,6) + R(6,6); + refpart.map(5, 4) = M * R(4, 4) + R(5, 4); + refpart.map(5, 5) = M * R(4, 5) + R(5, 5); } /** Close and deallocate all data and handles. diff --git a/src/particles/elements/SoftQuad.H b/src/particles/elements/SoftQuad.H index 2637d638e..9134e51e2 100644 --- a/src/particles/elements/SoftQuad.H +++ b/src/particles/elements/SoftQuad.H @@ -20,9 +20,9 @@ #include #include -#include #include #include +#include #include #include @@ -223,8 +223,8 @@ namespace SoftQuadrupoleData amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; - // get the linear map - amrex::Array2D const R = refpart.map; + // get the linear map (note: zero-indexed) + amrex::SmallMatrix const R = refpart.map; // symplectic linear map for a quadrupole is computed using the // Hamiltonian formalism as described in: @@ -233,18 +233,21 @@ namespace SoftQuadrupoleData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - xout = R(1,1)*x + R(1,2)*px + R(1,3)*y - + R(1,4)*py + R(1,5)*t + R(1,6)*pt; - pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y - + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - yout = R(3,1)*x + R(3,2)*px + R(3,3)*y - + R(3,4)*py + R(3,5)*t + R(3,6)*pt; - pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y - + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - tout = R(5,1)*x + R(5,2)*px + R(5,3)*y - + R(5,4)*py + R(5,5)*t + R(5,6)*pt; - ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y - + R(6,4)*py + R(6,5)*t + R(6,6)*pt; + // note: R is zero-indexed + // clang-format off + xout = R(0, 0) * x + R(0, 1) * px + R(0, 2) * y + + R(0, 3) * py + R(0, 4) * t + R(0, 5) * pt; + pxout = R(1, 0) * x + R(1, 1) * px + R(1, 2) * y + + R(1, 3) * py + R(1, 4) * t + R(1, 5) * pt; + yout = R(2, 0) * x + R(2, 1) * px + R(2, 2) * y + + R(2, 3) * py + R(2, 4) * t + R(2, 5) * pt; + pyout = R(4, 0) * x + R(3, 1) * px + R(3, 2) * y + + R(3, 3) * py + R(3, 4) * t + R(3, 5) * pt; + tout = R(4, 0) * x + R(4, 1) * px + R(4, 2) * y + + R(4, 3) * py + R(4, 4) * t + R(4, 5) * pt; + ptout = R(5, 0) * x + R(5, 1) * px + R(5, 2) * y + + R(5, 3) * py + R(5, 4) * t + R(5, 5) * pt; + // clang-format on // assign updated values x = xout; @@ -279,12 +282,7 @@ namespace SoftQuadrupoleData amrex::ParticleReal const sedge = refpart.sedge; // initialize linear map (deviation) values - for (int i=1; i<7; i++) { - for (int j=1; j<7; j++) { - auto const default_value = (i == j) ? 1.0_prt : 0.0_prt; - refpart.map(i, j) = default_value; - } - } + refpart.map = decltype(refpart.map)::Identity(); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -302,10 +300,10 @@ namespace SoftQuadrupoleData /* // print computed linear map: - for(int i=1; i<7; ++i){ - for(int j=1; j<7; ++j){ + for(int i=0; i<6; ++i){ + for(int j=0; j<6; ++j){ amrex::PrintToFile("QuadMap.txt") << i << " " << - j << " " << refpart.map(i,j) << "\n"; + j << " " << refpart.map(i, j) << "\n"; } } // @@ -410,21 +408,22 @@ namespace SoftQuadrupoleData zeval = z + tau; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const betgam = refpart.beta_gamma(); - refpart.map(1,1) = R(1,1) + tau*R(2,1); - refpart.map(1,2) = R(1,2) + tau*R(2,2); - refpart.map(1,3) = R(1,3) + tau*R(2,3); - refpart.map(1,4) = R(1,4) + tau*R(2,4); + // note: map and R are zero-indexed + refpart.map(0, 0) = R(0, 0) + tau * R(1, 0); + refpart.map(0, 1) = R(0, 1) + tau * R(1, 1); + refpart.map(0, 2) = R(0, 2) + tau * R(1, 2); + refpart.map(0, 3) = R(0, 3) + tau * R(1, 3); - refpart.map(3,1) = R(3,1) + tau*R(4,1); - refpart.map(3,2) = R(3,2) + tau*R(4,2); - refpart.map(3,3) = R(3,3) + tau*R(4,3); - refpart.map(3,4) = R(3,4) + tau*R(4,4); + refpart.map(2, 0) = R(2, 0) + tau * R(3, 0); + refpart.map(2, 1) = R(2, 1) + tau * R(3, 1); + refpart.map(2, 2) = R(2, 2) + tau * R(3, 2); + refpart.map(2, 3) = R(2, 3) + tau * R(3, 3); - refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,2); - refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,2); + refpart.map(4, 4) = R(4, 4) + tau * R(5, 4) / std::pow(betgam, 2); + refpart.map(4, 5) = R(4, 5) + tau * R(5, 5) / std::pow(betgam, 2); } @@ -457,18 +456,19 @@ namespace SoftQuadrupoleData refpart.pt = pt; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const alpha = G0*bz; - refpart.map(2,1) = R(2,1) - tau*alpha*R(1,1); - refpart.map(2,2) = R(2,2) - tau*alpha*R(1,2); - refpart.map(2,3) = R(2,3) - tau*alpha*R(1,3); - refpart.map(2,4) = R(2,4) - tau*alpha*R(1,4); + // note: map and R are zero-indexed + refpart.map(1, 0) = R(1, 0) - tau * alpha * R(0, 0); + refpart.map(1, 1) = R(1, 1) - tau * alpha * R(0, 1); + refpart.map(1, 2) = R(1, 2) - tau * alpha * R(0, 2); + refpart.map(1, 3) = R(1, 3) - tau * alpha * R(0, 3); - refpart.map(4,1) = R(4,1) + tau*alpha*R(3,1); - refpart.map(4,2) = R(4,2) + tau*alpha*R(3,2); - refpart.map(4,3) = R(4,3) + tau*alpha*R(3,3); - refpart.map(4,4) = R(4,4) + tau*alpha*R(3,4); + refpart.map(3, 0) = R(3, 0) + tau * alpha * R(2, 0); + refpart.map(3, 1) = R(3, 1) + tau * alpha * R(2, 1); + refpart.map(3, 2) = R(3, 2) + tau * alpha * R(2, 2); + refpart.map(3, 3) = R(3, 3) + tau * alpha * R(2, 3); } diff --git a/src/particles/elements/SoftSol.H b/src/particles/elements/SoftSol.H index f5263fe20..749821fe4 100644 --- a/src/particles/elements/SoftSol.H +++ b/src/particles/elements/SoftSol.H @@ -20,9 +20,9 @@ #include #include -#include #include #include +#include #include #include @@ -234,8 +234,8 @@ namespace SoftSolenoidData amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; - // get the linear map - amrex::Array2D const R = refpart.map; + // get the linear map (note: zero-indexed) + amrex::SmallMatrix const R = refpart.map; // symplectic linear map for a solenoid is computed using the // Hamiltonian formalism as described in: @@ -244,18 +244,21 @@ namespace SoftSolenoidData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - xout = R(1,1)*x + R(1,2)*px + R(1,3)*y - + R(1,4)*py + R(1,5)*t + R(1,6)*pt; - pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y - + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - yout = R(3,1)*x + R(3,2)*px + R(3,3)*y - + R(3,4)*py + R(3,5)*t + R(3,6)*pt; - pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y - + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - tout = R(5,1)*x + R(5,2)*px + R(5,3)*y - + R(5,4)*py + R(5,5)*t + R(5,6)*pt; - ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y - + R(6,4)*py + R(6,5)*t + R(6,6)*pt; + // note: R is zero-indexed + // clang-format off + xout = R(0, 0) * x + R(0, 1) * px + R(0, 2) * y + + R(0, 3) * py + R(0, 4) * t + R(0, 5) * pt; + pxout = R(1, 0) * x + R(1, 1) * px + R(1, 2) * y + + R(1, 3) * py + R(1, 4) * t + R(1, 5) * pt; + yout = R(2, 0) * x + R(2, 1) * px + R(2, 2) * y + + R(2, 3) * py + R(2, 4) * t + R(2, 5) * pt; + pyout = R(3, 0) * x + R(3, 1) * px + R(3, 2) * y + + R(3, 3) * py + R(3, 4) * t + R(3, 5) * pt; + tout = R(4, 0) * x + R(4, 1) * px + R(4, 2) * y + + R(4, 3) * py + R(4, 4) * t + R(4, 5) * pt; + ptout = R(5, 0) * x + R(5, 1) * px + R(5, 2) * y + + R(5, 3) * py + R(5, 4) * t + R(5, 5) * pt; + // clang-format on // assign updated values x = xout; @@ -290,12 +293,7 @@ namespace SoftSolenoidData amrex::ParticleReal const sedge = refpart.sedge; // initialize linear map (deviation) values - for (int i=1; i<7; i++) { - for (int j=1; j<7; j++) { - auto const default_value = (i == j) ? 1.0_prt : 0.0_prt; - refpart.map(i, j) = default_value; - } - } + refpart.map = decltype(refpart.map)::Identity(); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -312,10 +310,10 @@ namespace SoftSolenoidData amrex::ParticleReal const ptf = refpart.pt; /* print computed linear map: - for(int i=1; i<7; ++i){ - for(int j=1; j<7; ++j){ + for(int i=0; i<6; ++i){ + for(int j=0; j<6; ++j){ amrex::PrintToFile("SolMap.txt") << i << " " << - j << " " << refpart.map(i,j) << "\n"; + j << " " << refpart.map(i, j) << "\n"; } } */ @@ -419,22 +417,22 @@ namespace SoftSolenoidData zeval = z + tau; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const betgam = refpart.beta_gamma(); - refpart.map(1,1) = R(1,1) + tau*R(2,1); - refpart.map(1,2) = R(1,2) + tau*R(2,2); - refpart.map(1,3) = R(1,3) + tau*R(2,3); - refpart.map(1,4) = R(1,4) + tau*R(2,4); - - refpart.map(3,1) = R(3,1) + tau*R(4,1); - refpart.map(3,2) = R(3,2) + tau*R(4,2); - refpart.map(3,3) = R(3,3) + tau*R(4,3); - refpart.map(3,4) = R(3,4) + tau*R(4,4); + // note: map and R are zero-indexed + refpart.map(0, 0) = R(0, 0) + tau * R(2,1); + refpart.map(0, 1) = R(0, 1) + tau * R(2,2); + refpart.map(0, 2) = R(0, 2) + tau * R(2,3); + refpart.map(0, 3) = R(0, 3) + tau * R(2,4); - refpart.map(5,5) = R(5,5) + tau*R(6,5)/pow(betgam,2); - refpart.map(5,6) = R(5,6) + tau*R(6,6)/pow(betgam,2); + refpart.map(2, 0) = R(2, 0) + tau * R(3, 0); + refpart.map(2, 1) = R(2, 1) + tau * R(3, 1); + refpart.map(2, 2) = R(2, 2) + tau * R(3, 2); + refpart.map(2, 3) = R(2, 3) + tau * R(3, 3); + refpart.map(4, 4) = R(4, 4) + tau * R(5, 4) / std::pow(betgam, 2); + refpart.map(4, 5) = R(4, 5) + tau * R(5, 5) / std::pow(betgam, 2); } /** This pushes the reference particle and the linear map matrix @@ -469,19 +467,20 @@ namespace SoftSolenoidData refpart.pt = pt; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const alpha = B0*bz/2.0_prt; amrex::ParticleReal const alpha2 = std::pow(alpha,2); - refpart.map(2,1) = R(2,1) - tau*alpha2*R(1,1); - refpart.map(2,2) = R(2,2) - tau*alpha2*R(1,2); - refpart.map(2,3) = R(2,3) - tau*alpha2*R(1,3); - refpart.map(2,4) = R(2,4) - tau*alpha2*R(1,4); + // note: map and R are zero-indexed + refpart.map(1, 0) = R(1, 0) - tau * alpha2 * R(0, 0); + refpart.map(1, 1) = R(1, 1) - tau * alpha2 * R(0, 1); + refpart.map(1, 2) = R(1, 2) - tau * alpha2 * R(0, 2); + refpart.map(1, 3) = R(1, 3) - tau * alpha2 * R(0, 3); - refpart.map(4,1) = R(4,1) - tau*alpha2*R(3,1); - refpart.map(4,2) = R(4,2) - tau*alpha2*R(3,2); - refpart.map(4,3) = R(4,3) - tau*alpha2*R(3,3); - refpart.map(4,4) = R(4,4) - tau*alpha2*R(3,4); + refpart.map(3, 0) = R(3, 0) - tau * alpha2 * R(2, 0); + refpart.map(3, 1) = R(3, 1) - tau * alpha2 * R(2, 1); + refpart.map(3, 2) = R(3, 2) - tau * alpha2 * R(2, 2); + refpart.map(3, 3) = R(3, 3) - tau * alpha2 * R(2, 3); } @@ -518,30 +517,31 @@ namespace SoftSolenoidData refpart.pt = pt; // push the linear map equations - amrex::Array2D const R = refpart.map; + amrex::SmallMatrix const R = refpart.map; amrex::ParticleReal const theta = tau*B0*bz/2.0_prt; amrex::ParticleReal const cs = std::cos(theta); amrex::ParticleReal const sn = std::sin(theta); - refpart.map(1,1) = R(1,1)*cs + R(3,1)*sn; - refpart.map(1,2) = R(1,2)*cs + R(3,2)*sn; - refpart.map(1,3) = R(1,3)*cs + R(3,3)*sn; - refpart.map(1,4) = R(1,4)*cs + R(3,4)*sn; - - refpart.map(2,1) = R(2,1)*cs + R(4,1)*sn; - refpart.map(2,2) = R(2,2)*cs + R(4,2)*sn; - refpart.map(2,3) = R(2,3)*cs + R(4,3)*sn; - refpart.map(2,4) = R(2,4)*cs + R(4,4)*sn; - - refpart.map(3,1) = R(3,1)*cs - R(1,1)*sn; - refpart.map(3,2) = R(3,2)*cs - R(1,2)*sn; - refpart.map(3,3) = R(3,3)*cs - R(1,3)*sn; - refpart.map(3,4) = R(3,4)*cs - R(1,4)*sn; - - refpart.map(4,1) = R(4,1)*cs - R(2,1)*sn; - refpart.map(4,2) = R(4,2)*cs - R(2,2)*sn; - refpart.map(4,3) = R(4,3)*cs - R(2,3)*sn; - refpart.map(4,4) = R(4,4)*cs - R(2,4)*sn; + // note: map and R are zero-indexed + refpart.map(0, 0) = R(0, 0) * cs + R(2, 0) * sn; + refpart.map(0, 1) = R(0, 1) * cs + R(2, 1) * sn; + refpart.map(0, 2) = R(0, 2) * cs + R(2, 2) * sn; + refpart.map(0, 3) = R(0, 3) * cs + R(2, 3) * sn; + + refpart.map(1, 0) = R(1, 0) * cs + R(3, 0) * sn; + refpart.map(1, 1) = R(1, 1) * cs + R(3, 1) * sn; + refpart.map(1, 2) = R(1, 2) * cs + R(3, 2) * sn; + refpart.map(1, 3) = R(1, 3) * cs + R(3, 3) * sn; + + refpart.map(2, 0) = R(2, 0) * cs - R(0, 0) * sn; + refpart.map(2, 1) = R(2, 1) * cs - R(0, 1) * sn; + refpart.map(2, 2) = R(2, 2) * cs - R(0, 2) * sn; + refpart.map(2, 3) = R(2, 3) * cs - R(0, 3) * sn; + + refpart.map(3, 0) = R(3, 0) * cs - R(1, 0) * sn; + refpart.map(3, 1) = R(3, 1) * cs - R(1, 1) * sn; + refpart.map(3, 2) = R(3, 2) * cs - R(1, 2) * sn; + refpart.map(3, 3) = R(3, 3) * cs - R(1, 3) * sn; }