Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit bf938d0e authored by ext-rogers_c's avatar ext-rogers_c
Browse files

More ScalingFFAGMagnetTest testing and minor tweaks to field calculation

parent 789c8812
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/Classic/AbsBeamline/ScalingFFAGMagnet.cpp
+ 22
− 17
View file @ bf938d0e
...@@ -33,7 +33,7 @@ ...@@ -33,7 +33,7 @@
ScalingFFAGMagnet::ScalingFFAGMagnet(const std::string &name) ScalingFFAGMagnet::ScalingFFAGMagnet(const std::string &name)
: Component(name), : Component(name),
planarArcGeometry_m(1., 1.), dummy() { planarArcGeometry_m(1., 1.), dummy(), endField_m(NULL) {
setElType(isDrift); setElType(isDrift);
} }
...@@ -46,7 +46,9 @@ ScalingFFAGMagnet::ScalingFFAGMagnet(const ScalingFFAGMagnet &right) ...@@ -46,7 +46,9 @@ ScalingFFAGMagnet::ScalingFFAGMagnet(const ScalingFFAGMagnet &right)
phiEnd_m(right.phiEnd_m), azimuthalExtent_m(right.azimuthalExtent_m), phiEnd_m(right.phiEnd_m), azimuthalExtent_m(right.azimuthalExtent_m),
verticalExtent_m(right.verticalExtent_m), centre_m(right.centre_m), verticalExtent_m(right.verticalExtent_m), centre_m(right.centre_m),
dfCoefficients_m(right.dfCoefficients_m) { dfCoefficients_m(right.dfCoefficients_m) {
delete endField_m; if (endField_m != NULL) {
delete endField_m;
}
endField_m = right.endField_m->clone(); endField_m = right.endField_m->clone();
RefPartBunch_m = right.RefPartBunch_m; RefPartBunch_m = right.RefPartBunch_m;
setElType(isDrift); setElType(isDrift);
...@@ -55,7 +57,9 @@ ScalingFFAGMagnet::ScalingFFAGMagnet(const ScalingFFAGMagnet &right) ...@@ -55,7 +57,9 @@ ScalingFFAGMagnet::ScalingFFAGMagnet(const ScalingFFAGMagnet &right)
} }
ScalingFFAGMagnet::~ScalingFFAGMagnet() { ScalingFFAGMagnet::~ScalingFFAGMagnet() {
delete endField_m; if (endField_m != NULL) {
delete endField_m;
}
} }
ElementBase* ScalingFFAGMagnet::clone() const { ElementBase* ScalingFFAGMagnet::clone() const {
...@@ -144,7 +148,7 @@ bool ScalingFFAGMagnet::getFieldValueCylindrical(const Vector_t &pos, Vector_t & ...@@ -144,7 +148,7 @@ bool ScalingFFAGMagnet::getFieldValueCylindrical(const Vector_t &pos, Vector_t &
if (z < -verticalExtent_m || z > verticalExtent_m) { if (z < -verticalExtent_m || z > verticalExtent_m) {
return true; return true;
} }
std::vector<double> fringeDerivatives(maxOrder_m, 0.); std::vector<double> fringeDerivatives(maxOrder_m+1, 0.);
for (size_t i = 0; i < fringeDerivatives.size(); ++i) { for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
fringeDerivatives[i] = endField_m->function(phiSpiral, i); // d^i_phi f fringeDerivatives[i] = endField_m->function(phiSpiral, i); // d^i_phi f
} }
...@@ -154,13 +158,15 @@ bool ScalingFFAGMagnet::getFieldValueCylindrical(const Vector_t &pos, Vector_t & ...@@ -154,13 +158,15 @@ bool ScalingFFAGMagnet::getFieldValueCylindrical(const Vector_t &pos, Vector_t &
for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) { for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) {
f2n += dfCoefficients_m[n][i]*fringeDerivatives[i]; f2n += dfCoefficients_m[n][i]*fringeDerivatives[i];
} }
double f2nplus1 = 0;
for (size_t i = 0; i < dfCoefficients_m[n+1].size() && n+1 < dfCoefficients_m.size(); ++i) {
f2nplus1 += dfCoefficients_m[n+1][i]*fringeDerivatives[i];
}
deltaB[1] = f2n*h*pow(z/r, n); // Bz = sum(f_2n * h * (z/r)^2n deltaB[1] = f2n*h*pow(z/r, n); // Bz = sum(f_2n * h * (z/r)^2n
deltaB[2] = f2nplus1*h*pow(z/r, n+1); // Bphi = sum(f_2n+1 * h * (z/r)^2n+1 if (maxOrder_m >= n+1) {
deltaB[0] = (f2n*(k_m-n)/(n+1) - tanDelta_m*f2nplus1)*h*pow(z/r, n+1); double f2nplus1 = 0;
for (size_t i = 0; i < dfCoefficients_m[n+1].size() && n+1 < dfCoefficients_m.size(); ++i) {
f2nplus1 += dfCoefficients_m[n+1][i]*fringeDerivatives[i];
}
deltaB[2] = f2nplus1*h*pow(z/r, n+1); // Bphi = sum(f_2n+1 * h * (z/r)^2n+1
deltaB[0] = (f2n*(k_m-n)/(n+1) - tanDelta_m*f2nplus1)*h*pow(z/r, n+1);
}
B += deltaB; B += deltaB;
} }
return false; return false;
...@@ -173,17 +179,14 @@ bool ScalingFFAGMagnet::apply(const Vector_t &R, const Vector_t &P, ...@@ -173,17 +179,14 @@ bool ScalingFFAGMagnet::apply(const Vector_t &R, const Vector_t &P,
} }
void ScalingFFAGMagnet::calculateDfCoefficients() { void ScalingFFAGMagnet::calculateDfCoefficients() {
dfCoefficients_m = std::vector<std::vector<double> >(maxOrder_m); dfCoefficients_m = std::vector<std::vector<double> >(maxOrder_m+1);
if (maxOrder_m == 0) {
return;
}
dfCoefficients_m[0] = std::vector<double>(1, 1.); // f_0 = 1.*0th derivative dfCoefficients_m[0] = std::vector<double>(1, 1.); // f_0 = 1.*0th derivative
for (size_t n = 0; n+1 < maxOrder_m; n += 2) { // n indexes the power in z for (size_t n = 0; n < maxOrder_m; n += 2) { // n indexes the power in z
dfCoefficients_m[n+1] = std::vector<double>(dfCoefficients_m[n].size()+1, 0); dfCoefficients_m[n+1] = std::vector<double>(dfCoefficients_m[n].size()+1, 0);
for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) { // i indexes the derivative for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) { // i indexes the derivative
dfCoefficients_m[n+1][i+1] = dfCoefficients_m[n][i]/(n+1); dfCoefficients_m[n+1][i+1] = dfCoefficients_m[n][i]/(n+1);
} }
if (n+2 == maxOrder_m) { if (n+1 == maxOrder_m) {
break; break;
} }
dfCoefficients_m[n+2] = std::vector<double>(dfCoefficients_m[n].size()+2, 0); dfCoefficients_m[n+2] = std::vector<double>(dfCoefficients_m[n].size()+2, 0);
...@@ -199,7 +202,9 @@ void ScalingFFAGMagnet::calculateDfCoefficients() { ...@@ -199,7 +202,9 @@ void ScalingFFAGMagnet::calculateDfCoefficients() {
} }
void ScalingFFAGMagnet::setEndField(endfieldmodel::EndFieldModel* endField) { void ScalingFFAGMagnet::setEndField(endfieldmodel::EndFieldModel* endField) {
delete endField_m; if (endField_m != NULL) {
delete endField_m;
}
endField_m = endField; endField_m = endField;
} }
tests/classic_src/AbsBeamline/ScalingFFAGMagnetTest.cpp
+ 87
− 10
View file @ bf938d0e
...@@ -30,6 +30,7 @@ ...@@ -30,6 +30,7 @@
#include <sstream> #include <sstream>
#include "gtest/gtest.h" #include "gtest/gtest.h"
#include "opal_test_utilities/SilenceTest.h"
#include "Classic/AbsBeamline/EndFieldModel/Tanh.h" #include "Classic/AbsBeamline/EndFieldModel/Tanh.h"
#include "Classic/AbsBeamline/ScalingFFAGMagnet.h" #include "Classic/AbsBeamline/ScalingFFAGMagnet.h"
...@@ -151,12 +152,12 @@ public: ...@@ -151,12 +152,12 @@ public:
void printCoefficients() { void printCoefficients() {
std::vector<std::vector<double> > coeff = sector_m->getDfCoefficients(); std::vector<std::vector<double> > coeff = sector_m->getDfCoefficients();
std::cerr << "Coefficients" << std::endl; std::cout << "Coefficients" << std::endl;
for (size_t n = 0; n < coeff.size(); ++n) { for (size_t n = 0; n < coeff.size(); ++n) {
for (size_t i = 0; i < coeff[n].size(); ++i) { for (size_t i = 0; i < coeff[n].size(); ++i) {
std::cerr << coeff[n][i] << " "; std::cout << coeff[n][i] << " ";
} }
std::cerr << std::endl; std::cout << std::endl;
} }
} }
...@@ -283,12 +284,74 @@ TEST_F(ScalingFFAGMagnetTest, PlacementTest) { ...@@ -283,12 +284,74 @@ TEST_F(ScalingFFAGMagnetTest, PlacementTest) {
} }
} }
TEST_F(ScalingFFAGMagnetTest, DFCoefficientsTest) {
sector_m->setTanDelta(0.0);
sector_m->setMaxOrder(5);
sector_m->setFieldIndex(5);
sector_m->initialise();
// hard coded - calculated by hand
double ref[5][5] = {
{1., -999., -999., -999., -999.}, // n = 0
{0., 1., -999., -999., -999.}, // n = 1
{-25./2., 0., -1./2., -999., -999.}, // n = 2
{0., -25./6., 0., -1./6., -999.}, // n = 3
{+25./2.*3./4., 0., +1./2.*3./4.+25./6./4., 0., +1./6./4.}, // n = 4
};
std::vector< std::vector<double> > coeffs = sector_m->getDfCoefficients();
ASSERT_GE(coeffs.size(), 5);
for (size_t n = 0; n < 5; ++n) {
ASSERT_EQ(coeffs[n].size(), n+1);
for (size_t i = 0; i < coeffs[n].size(); ++i) {
EXPECT_NEAR(coeffs[n][i], ref[n][i], 1e-9) << " n: " << n << " i: " << i;
}
}
}
TEST_F(ScalingFFAGMagnetTest, DFCoefficientsTanDeltaTest) {
sector_m->setTanDelta(2.0);
sector_m->setMaxOrder(4); // BUG - max order is 1 to high
sector_m->setFieldIndex(5);
sector_m->initialise();
// hard coded - calculated by hand
double ref[5][5] = {
{1., -999., -999., -999., -999.}, // n = 0
{0., 1., -999., -999., -999.}, // n = 1
{-25./2., 10.*2./2., -5./2., -999., -999.}, // n = 2
{0., -25./6., +10./3., -5./6., -999.}, // n = 3
// i gave up at 4 - head exploding
{+25./2.*3./4., -25./6./4.*2.*3.*2.-10.*3/4., -999., -999., -999.}, // n = 4
};
std::vector< std::vector<double> > coeffs = sector_m->getDfCoefficients();
ASSERT_GE(coeffs.size(), 4);
for (size_t n = 0; n < 4; ++n) {
ASSERT_EQ(coeffs[n].size(), n+1);
for (size_t i = 0; i < coeffs[n].size(); ++i) {
EXPECT_NEAR(coeffs[n][i], ref[n][i], 1e-9) << " n: " << n << " i: " << i;
}
}
}
TEST_F(ScalingFFAGMagnetTest, TanhTest) {
OpalTestUtilities::SilenceTest silence;
double numericalDerivative = sector_m->getEndField()->function(-psi0_m, 0);
for (size_t order = 0; order < 5; ++order) {
double analyticalDerivative = sector_m->getEndField()->function(-psi0_m, order);
if (fabs(numericalDerivative)+fabs(analyticalDerivative) > 1e-3) {
EXPECT_NEAR(analyticalDerivative, numericalDerivative, fabs(analyticalDerivative)*1e-3);
}
std::cout << order << " " << analyticalDerivative << " " << numericalDerivative << std::endl;
numericalDerivative = sector_m->getEndField()->function(-psi0_m*0.9999, order)-
sector_m->getEndField()->function(-psi0_m*1.0001, order);
numericalDerivative /= -psi0_m*0.9999 + psi0_m*1.0001;
}
}
TEST_F(ScalingFFAGMagnetTest, BTwoDTest) { TEST_F(ScalingFFAGMagnetTest, BTwoDTest) {
std::ofstream fout("/tmp/b_twod.out"); std::ofstream fout("/tmp/b_twod.out");
bool passtest = true; bool passtest = true;
for (double y = 0.; y < 0.025; y += 0.015) { for (double y = 0.; y < 0.025; y += 0.015) {
for (double r = r0_m; r < r0_m+1; r += 0.01) { for (double r = r0_m; r < r0_m+1; r += 0.02) {
for (double psi = -2.*psi0_m; psi < 4.00001*psi0_m; psi += psi0_m/10.) { for (double psi = -2.*psi0_m; psi < 4.00001*psi0_m; psi += psi0_m/5.) {
passtest &= printLine(Vector_t(r, y, psi), 0., fout, 1e-1); passtest &= printLine(Vector_t(r, y, psi), 0., fout, 1e-1);
} }
} }
...@@ -309,21 +372,35 @@ TEST_F(ScalingFFAGMagnetTest, ConvergenceYTest) { ...@@ -309,21 +372,35 @@ TEST_F(ScalingFFAGMagnetTest, ConvergenceYTest) {
} }
TEST_F(ScalingFFAGMagnetTest, ConvergenceOrderTest) { TEST_F(ScalingFFAGMagnetTest, ConvergenceOrderTest) {
for (double y = 0.05; y > 0.0002; y /= 10.) { OpalTestUtilities::SilenceTest silence;
std::vector<double> divBVec(3); for (double y = 0.5; y > 0.2; y /= 10.) { // 50 cm off midplane
std::cout << "order y divB |curlB| curlB" << std::endl;
std::vector<double> divBVec(13);
std::vector<double> curlBVec(13);
double delta = y/100.; double delta = y/100.;
for (size_t i = 0; i < divBVec.size(); ++i) { for (size_t i = 0; i < divBVec.size(); ++i) {
sector_m->setMaxOrder((i+1)*2); sector_m->setMaxOrder(i);
sector_m->initialise(); sector_m->initialise();
Vector_t pos(r0_m, y, psi0_m*2); Vector_t pos(r0_m, y, psi0_m*2);
double divB = getDivBCylindrical(pos, Vector_t(delta, delta, delta/r0_m)); double divB = getDivBCylindrical(pos, Vector_t(delta, delta, delta/r0_m));
Vector_t curlB = getCurlBCyl(pos, Vector_t(delta, delta, delta/r0_m));
double curlBMag =
sqrt(curlB[0]*curlB[0] + curlB[1]*curlB[1] + curlB[2]*curlB[2]);
divB = fabs(divB); divB = fabs(divB);
divBVec[i] = divB; divBVec[i] = divB;
if (i > 0) { curlBVec[i] = curlBMag;
EXPECT_LT(divBVec[i], divBVec[i-1]) << " with i " std::cout << i << " " << y << " " << divB << " "
<< curlBMag << " " << curlB << std::endl;
if (i > 1 && i % 2 == 1) {
EXPECT_LT(divBVec[i], divBVec[i-2]) << " with i "
<< i << std::endl;
}
if (i > 1 && i % 2 == 0) {
EXPECT_LT(curlBVec[i], curlBVec[i-2]) << " with i "
<< i << std::endl; << i << std::endl;
} }
} }
std::cout << std::endl;
} }
sector_m->setMaxOrder(10); sector_m->setMaxOrder(10);
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment