diff --git a/src/Algorithms/CavityAutophaser.cpp b/src/Algorithms/CavityAutophaser.cpp
index 9fa6b29ecacc0e848295ab56a1fcd8a27e7dc419..4c791f71400bc77f3267baced2efad76627cac83 100644
--- a/src/Algorithms/CavityAutophaser.cpp
+++ b/src/Algorithms/CavityAutophaser.cpp
@@ -72,7 +72,7 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
<< std::left << std::setw(68) << std::setfill('*') << ss.str()
<< std::setfill(' ') << endl);
if (!apVeto) {
- double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
+ double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * 1e-6;
double AstraPhase = 0.0;
double designEnergy = element->getDesignEnergy();
@@ -191,7 +191,7 @@ double CavityAutophaser::guessCavityPhase(double t) {
return orig_phi;
}
- Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
+ Phimax = element->getAutoPhaseEstimate(Util::getKineticEnergy(refP, itsReference_m.getM()) * 1e-6,
t,
itsReference_m.getQ(),
itsReference_m.getM() * 1e-6);
@@ -288,7 +288,7 @@ double CavityAutophaser::track(Vector_t /*R*/,
out);
rfc->setPhasem(initialPhase);
- double finalKineticEnergy = Util::getEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
+ double finalKineticEnergy = Util::getKineticEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
return finalKineticEnergy;
}
\ No newline at end of file
diff --git a/src/Algorithms/CavityAutophaser.h b/src/Algorithms/CavityAutophaser.h
index 29db048b24349a39a7e9b93b5518d18aafc5571a..a4e49c8da834b0c8cd8d8010aea35439c4315403 100644
--- a/src/Algorithms/CavityAutophaser.h
+++ b/src/Algorithms/CavityAutophaser.h
@@ -28,8 +28,6 @@ private:
const double dt,
const double phase,
std::ofstream *out = NULL) const;
- double getEnergyMeV(const Vector_t &P);
- double getMomentum(double kineticEnergyMeV);
const PartData &itsReference_m;
std::shared_ptr<Component> itsCavity_m;
@@ -39,14 +37,4 @@ private:
};
-inline
-double CavityAutophaser::getEnergyMeV(const Vector_t &P) {
- return itsReference_m.getM() * 1e-6 * (sqrt(dot(P,P) + 1) - 1);
-}
-
-inline
-double CavityAutophaser::getMomentum(double kineticEnergyMeV) {
- return sqrt(std::pow(kineticEnergyMeV / (itsReference_m.getM() * 1e-6) + 1, 2) - 1);
-}
-
#endif
\ No newline at end of file
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index f38da15e80bd5e36a201c9a9637a7aeaee16c12c..23e2a6c2fdfb97e3d3d8cec97a629ab3756b0452 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -1162,7 +1162,7 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
selectDT();
if ((back_track && itsOpalBeamline_m.containsSource()) ||
- Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
+ Util::getKineticEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
double gamma = 0.1 / itsBunch_m->getM() + 1.0;
double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index b8ebb79459043ae993e9c3b3b05712ad5183e541..33cecafa90dc8f8fee454843f7be5d0b6e1894e3 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -494,22 +494,22 @@ ElementBase::ElementType RFCavity::getType() const {
double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
const double dt = 1e-13;
- const double p0 = Util::getP(E0, mass);
+ const double p0 = Util::getBetaGamma(E0, mass);
const double origPhase =getPhasem();
double dphi = Physics::pi / 18;
double phi = 0.0;
setPhasem(phi);
- std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
+ std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
double phimax = 0.0;
- double Emax = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+ double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
phi += dphi;
for (unsigned int j = 0; j < 2; ++ j) {
for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
setPhasem(phi);
ret = trackOnAxisParticle(p0, t0, dt, q, mass);
- double Ekin = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+ double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
if (Ekin > Emax) {
Emax = Ekin;
phimax = phi;
@@ -635,7 +635,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
double cosine_part = 0.0, sine_part = 0.0;
- double p0 = std::sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
+ double p0 = Util::getBetaGamma(E0, mass);
cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
@@ -675,11 +675,11 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
const double zend = length_m + startField_m;
Vector_t z(0.0, 0.0, zbegin);
- double dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
+ double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
Vector_t Ef(0.0), Bf(0.0);
if (out) *out << std::setw(18) << z[2]
- << std::setw(18) << Util::getEnergy(p, mass)
+ << std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
while (z(2) + dz < zend && z(2) + dz > zbegin) {
z /= cdt;
@@ -700,7 +700,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
t += dt;
if (out) *out << std::setw(18) << z[2]
- << std::setw(18) << Util::getEnergy(p, mass)
+ << std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
}
diff --git a/src/Classic/Algorithms/PartBunch.cpp b/src/Classic/Algorithms/PartBunch.cpp
index b475a680dd1b5f7604acef6ec212cf60b14f0d23..2be8fe9af747162e8cb6e4722b4b1cb65c3e7aa8 100644
--- a/src/Classic/Algorithms/PartBunch.cpp
+++ b/src/Classic/Algorithms/PartBunch.cpp
@@ -123,8 +123,7 @@ void PartBunch::computeSelfFields(int binNumber) {
if(fs_m->hasValidSolver()) {
/// Mesh the whole domain
- if(fs_m->getFieldSolverType() == "SAAMG")
- resizeMesh();
+ resizeMesh();
/// Scatter charge onto space charge grid.
this->Q *= this->dt;
@@ -328,12 +327,14 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
+ if (fs_m->getFieldSolverType() != "SAAMG") {
+ return;
+ }
+
double xmin = fs_m->solver_m->getXRangeMin();
double xmax = fs_m->solver_m->getXRangeMax();
double ymin = fs_m->solver_m->getYRangeMin();
double ymax = fs_m->solver_m->getYRangeMax();
- double zmin = fs_m->solver_m->getZRangeMin();
- double zmax = fs_m->solver_m->getZRangeMax();
if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
ymin > rmin_m[1] || ymax < rmax_m[1]) {
@@ -354,14 +355,14 @@ void PartBunch::resizeMesh() {
boundp();
get_bounds(rmin_m, rmax_m);
}
- Vector_t mymin = Vector_t(xmin, ymin , zmin);
- Vector_t mymax = Vector_t(xmax, ymax , zmax);
- for (int i = 0; i < 3; i++)
- hr_m[i] = (mymax[i] - mymin[i])/nr_m[i];
+ Vector_t origin = Vector_t(0.0, 0.0, 0.0);
+
+ // update the mesh origin and mesh spacing hr_m
+ fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
getMesh().set_meshSpacing(&(hr_m[0]));
- getMesh().set_origin(mymin);
+ getMesh().set_origin(origin);
rho_m.initialize(getMesh(),
getFieldLayout(),
@@ -384,8 +385,7 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
//mesh the whole domain
- if(fs_m->getFieldSolverType() == "SAAMG")
- resizeMesh();
+ resizeMesh();
//scatter charges onto grid
this->Q *= this->dt;
@@ -510,8 +510,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
- if(fs_m->getFieldSolverType() == "SAAMG")
- resizeMesh();
+ resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
@@ -639,8 +638,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
- if(fs_m->getFieldSolverType() == "SAAMG")
- resizeMesh();
+ resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h
index 76807acfa037f60c17faccd72b4f51a3a7b43870..47b865bd6587ccaacf28de35295d892cb421313a 100644
--- a/src/Classic/Algorithms/PartBunchBase.h
+++ b/src/Classic/Algorithms/PartBunchBase.h
@@ -411,6 +411,8 @@ public:
// virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
virtual FieldLayout_t &getFieldLayout() = 0;
+ virtual void resizeMesh() { };
+
/*
* Wrapped member functions of IpplParticleBase
*/
diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h
index edf7178517c527c010c2bfcfd183278465f0fd1a..42479de4b024055ced8e0e5cc28a80a702a933e2 100644
--- a/src/Classic/Utilities/Util.h
+++ b/src/Classic/Utilities/Util.h
@@ -6,6 +6,7 @@
#include <string>
#include <cstring>
+#include <limits>
#include <sstream>
#include <type_traits>
#include <functional>
@@ -28,16 +29,23 @@ namespace Util {
}
inline
- double getEnergy(Vector_t p, double mass) {
+ double getKineticEnergy(Vector_t p, double mass) {
return (getGamma(p) - 1.0) * mass;
}
inline
- double getP(double E, double mass) {
- double gamma = E / mass + 1;
- return std::sqrt(std::pow(gamma, 2.0) - 1.0);
+ double getBetaGamma(double Ekin, double mass) {
+ double value = std::sqrt(std::pow(Ekin / mass + 1.0, 2) - 1.0);
+ if (value < std::numeric_limits<double>::epsilon())
+ value = std::sqrt(2 * Ekin / mass);
+ return value;
}
+ inline
+ double convertMomentumeVToBetaGamma(double p, double mass) {
+ return p / mass;
+ }
+
inline
std::string getTimeString(double time, unsigned int precision = 3) {
std::string timeUnit(" [ps]");
@@ -162,7 +170,7 @@ namespace Util {
}
Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = "");
-
+
std::string toUpper(const std::string &str);
std::string combineFilePath(std::initializer_list<std::string>);
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index eddb76da941932e356999563f1aaa32b001ce506..744df64db31d2ac6f4f04aeedaea6926d611cce6 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -17,7 +17,6 @@
#include <string>
#include <vector>
#include <numeric>
-#include <limits>
// IPPL
#include "DataSource/DataConnect.h"
@@ -875,14 +874,6 @@ void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUn
}
-double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
- double value = std::copysign(std::sqrt(std::pow(std::abs(valueIneV) / massIneV + 1.0, 2) - 1.0), valueIneV);
- if (std::abs(value) < std::numeric_limits<double>::epsilon())
- value = std::copysign(std::sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
-
- return value;
-}
-
void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
setDistParametersBinomial(massIneV);
@@ -1069,9 +1060,9 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
saveProcessor = 0;
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
- P(0) = converteVToBetaGamma(P(0), massIneV);
- P(1) = converteVToBetaGamma(P(1), massIneV);
- P(2) = converteVToBetaGamma(P(2), massIneV);
+ P(0) = Util::convertMomentumeVToBetaGamma(P(0), massIneV);
+ P(1) = Util::convertMomentumeVToBetaGamma(P(1), massIneV);
+ P(2) = Util::convertMomentumeVToBetaGamma(P(2), massIneV);
}
pmean_m += P;
@@ -1422,7 +1413,7 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
// This is PC from BEAM
double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
- deltaP = converteVToBetaGamma(deltaP, beam->getM());
+ deltaP = Util::convertMomentumeVToBetaGamma(deltaP, beam->getM());
}
avrgpz_m = beam->getP()/beam->getM() + deltaP;
@@ -3650,9 +3641,9 @@ void Distribution::setSigmaP_m(double massIneV) {
// Check what input units we are using for momentum.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
- sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
- sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
- sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
+ sigmaP_m[0] = Util::convertMomentumeVToBetaGamma(sigmaP_m[0], massIneV);
+ sigmaP_m[1] = Util::convertMomentumeVToBetaGamma(sigmaP_m[1], massIneV);
+ sigmaP_m[2] = Util::convertMomentumeVToBetaGamma(sigmaP_m[2], massIneV);
}
}
@@ -3985,14 +3976,14 @@ void Distribution::setupEmissionModel(PartBunchBase<double, 3> *beam) {
void Distribution::setupEmissionModelAstra(PartBunchBase<double, 3> *beam) {
double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
- pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+ pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
pmean_m = Vector_t(0.0, 0.0, 0.5 * pTotThermal_m);
}
void Distribution::setupEmissionModelNone(PartBunchBase<double, 3> *beam) {
double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
- pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+ pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
double avgPz = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
size_t numParticles = pzDist_m.size();
reduce(avgPz, avgPz, OpAddAssign());
@@ -4146,9 +4137,9 @@ void Distribution::shiftDistCoordinates(double massIneV) {
// Check input momentum units.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
- deltaPx = converteVToBetaGamma(deltaPx, massIneV);
- deltaPy = converteVToBetaGamma(deltaPy, massIneV);
- deltaPz = converteVToBetaGamma(deltaPz, massIneV);
+ deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+ deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
+ deltaPz = Util::convertMomentumeVToBetaGamma(deltaPz, massIneV);
}
size_t endIdx = startIdx + particlesPerDist_m[i];
@@ -4429,8 +4420,8 @@ void Distribution::adjustPhaseSpace(double massIneV) {
double deltaPy = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPY]);
// Check input momentum units.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
- deltaPx = converteVToBetaGamma(deltaPx, massIneV);
- deltaPy = converteVToBetaGamma(deltaPy, massIneV);
+ deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+ deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
}
double avrg[6];
diff --git a/src/Distribution/Distribution.h b/src/Distribution/Distribution.h
index 15e728fa666d756889d65ad9ddfa49ef0f5a030e..31817e3c968886ab2d141f20d0e3c6876441a383 100644
--- a/src/Distribution/Distribution.h
+++ b/src/Distribution/Distribution.h
@@ -292,7 +292,6 @@ private:
void checkIfEmitted();
void checkParticleNumber(size_t &numberOfParticles);
void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits);
- double converteVToBetaGamma(double valueIneV, double massIneV);
size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
class BinomialBehaviorSplitter {
diff --git a/src/Elements/OpalElement.cpp b/src/Elements/OpalElement.cpp
index 63bf8b771bf9200523104f97bbb7684f2ca46d29..684fc18f62354414ef97a6bbe5bb05471f603e96 100644
--- a/src/Elements/OpalElement.cpp
+++ b/src/Elements/OpalElement.cpp
@@ -568,7 +568,8 @@ void OpalElement::update() {
rotation = rotTheta * (rotPhi * rotPsi);
} else {
if (itsAttr[ORIENTATION]) {
- throw OpalException("Line::parse","Parameter orientation is array of 3 values (theta, phi, psi);\n" +
+ throw OpalException("OpalElement::update",
+ "Parameter orientation is array of 3 values (theta, phi, psi);\n" +
std::to_string(dir.size()) + " values provided");
}
}
@@ -577,7 +578,8 @@ void OpalElement::update() {
origin = Vector_t(ori[0], ori[1], ori[2]);
} else {
if (itsAttr[ORIGIN]) {
- throw OpalException("Line::parse","Parameter origin is array of 3 values (x, y, z);\n" +
+ throw OpalException("OpalElement::update",
+ "Parameter origin is array of 3 values (x, y, z);\n" +
std::to_string(ori.size()) + " values provided");
}
}
diff --git a/src/Elements/OpalRBend.cpp b/src/Elements/OpalRBend.cpp
index dec75dc88fd2a27239cc8538a2487d68a35b513b..0f1b9e0be8c28d849f5a493ee21c4a8839024680 100644
--- a/src/Elements/OpalRBend.cpp
+++ b/src/Elements/OpalRBend.cpp
@@ -24,7 +24,6 @@
#include "Structure/OpalWake.h"
#include "Structure/ParticleMatterInteraction.h"
#include "Utilities/OpalException.h"
-#include <cmath>
OpalRBend::OpalRBend():
OpalBend("RBEND",
@@ -47,10 +46,8 @@ OpalRBend::OpalRBend(const std::string &name, OpalRBend *parent):
OpalRBend::~OpalRBend() {
- if(owk_m)
- delete owk_m;
- if(parmatint_m)
- delete parmatint_m;
+ delete owk_m;
+ delete parmatint_m;
}
@@ -131,7 +128,7 @@ void OpalRBend::update() {
double k0s = itsAttr[K0S] ? Attributes::getReal(itsAttr[K0S]) : 0.0;
//JMJ 4/10/2000: above line replaced
// length ? angle / length : angle;
- // to avoid closed orbit created by RBEND with defalt K0.
+ // to avoid closed orbit created by RBEND with default K0.
field.setNormalComponent(1, factor * k0);
field.setSkewComponent(1, factor * Attributes::getReal(itsAttr[K0S]));
field.setNormalComponent(2, factor * Attributes::getReal(itsAttr[K1]));
@@ -177,8 +174,11 @@ void OpalRBend::update() {
}
// Energy in eV.
- if(itsAttr[DESIGNENERGY]) {
+ if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+ } else if (bend->getName() != "RBEND") {
+ throw OpalException("OpalRBend::update",
+ "RBend requires non-zero DESIGNENERGY");
}
bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
diff --git a/src/Elements/OpalRBend3D.cpp b/src/Elements/OpalRBend3D.cpp
index f64500bfddc8945e4520fc85fbce6fe40bdd9ad1..e1175d63fcd64bce123bf92242e2bd0dfbb925b3 100644
--- a/src/Elements/OpalRBend3D.cpp
+++ b/src/Elements/OpalRBend3D.cpp
@@ -72,10 +72,8 @@ OpalRBend3D::OpalRBend3D(const std::string &name, OpalRBend3D *parent):
}
OpalRBend3D::~OpalRBend3D() {
- if(owk_m)
- delete owk_m;
- if(parmatint_m)
- delete parmatint_m;
+ delete owk_m;
+ delete parmatint_m;
}
OpalRBend3D *OpalRBend3D::clone(const std::string &name) {
@@ -130,8 +128,11 @@ void OpalRBend3D::update() {
bend->setEntranceAngle(e1);
// Energy in eV.
- if(itsAttr[DESIGNENERGY]) {
+ if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+ } else if (bend->getName() != "RBEND3D") {
+ throw OpalException("OpalRBend3D::update",
+ "RBend3D requires non-zero DESIGNENERGY");
}
bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
diff --git a/src/Elements/OpalSBend.cpp b/src/Elements/OpalSBend.cpp
index 785c1bf3d17951b734bc898323be340f58d297e2..2ec776a807aa5568bbd0e2005e294229f88bbb4a 100644
--- a/src/Elements/OpalSBend.cpp
+++ b/src/Elements/OpalSBend.cpp
@@ -24,7 +24,6 @@
#include "Structure/OpalWake.h"
#include "Structure/ParticleMatterInteraction.h"
#include "Utilities/OpalException.h"
-#include <cmath>
OpalSBend::OpalSBend():
OpalBend("SBEND",
@@ -47,10 +46,8 @@ OpalSBend::OpalSBend(const std::string &name, OpalSBend *parent):
OpalSBend::~OpalSBend() {
- if(owk_m)
- delete owk_m;
- if(parmatint_m)
- delete parmatint_m;
+ delete owk_m;
+ delete parmatint_m;
}
@@ -131,7 +128,7 @@ void OpalSBend::update() {
double k0s = itsAttr[K0S] ? Attributes::getReal(itsAttr[K0S]) : 0.0;
//JMJ 4/10/2000: above line replaced
// length ? angle / length : angle;
- // to avoid closed orbit created by SBEND with defalt K0.
+ // to avoid closed orbit created by SBEND with default K0.
field.setNormalComponent(1, factor * k0);
field.setSkewComponent(1, factor * Attributes::getReal(itsAttr[K0S]));
field.setNormalComponent(2, factor * Attributes::getReal(itsAttr[K1]));
@@ -163,11 +160,11 @@ void OpalSBend::update() {
if(itsAttr[GREATERTHANPI])
throw OpalException("OpalSBend::update",
- "GREATERTHANPI not supportet any more");
+ "GREATERTHANPI not supported anymore");
if(itsAttr[ROTATION])
throw OpalException("OpalSBend::update",
- "ROTATION not supportet any more; use PSI instead");
+ "ROTATION not supported anymore; use PSI instead");
if(itsAttr[FMAPFN])
bend->setFieldMapFN(Attributes::getString(itsAttr[FMAPFN]));
@@ -183,15 +180,18 @@ void OpalSBend::update() {
bend->setExitAngle(e2);
// Units are eV.
- if(itsAttr[DESIGNENERGY]) {
+ if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+ } else if (bend->getName() != "SBEND") {
+ throw OpalException("OpalSBend::update",
+ "SBend requires non-zero DESIGNENERGY");
}
bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
if(itsAttr[APERT])
- throw OpalException("OpalRBend::fillRegisteredAttributes",
- "APERTURE in RBEND not supported; use GAP and HAPERT instead");
+ throw OpalException("OpalSBend::update",
+ "APERTURE in SBEND not supported; use GAP and HAPERT instead");
if(itsAttr[HAPERT]) {
double hapert = Attributes::getReal(itsAttr[HAPERT]);
diff --git a/src/Solvers/EllipticDomain.cpp b/src/Solvers/EllipticDomain.cpp
index 28b6e8f2780ca407ad43de1898c5edfac73781a6..d3365cd53028bd7542eceaadb4ad370342b7e714 100644
--- a/src/Solvers/EllipticDomain.cpp
+++ b/src/Solvers/EllipticDomain.cpp
@@ -1,6 +1,8 @@
//
// Class EllipticDomain
-// :FIXME: add brief description
+// This class provides an elliptic beam pipe. The mesh adapts to the bunch size
+// in the longitudinal direction. At the intersection of the mesh with the
+// beam pipe, three stencil interpolation methods are available.
//
// Copyright (c) 2008, Yves Ineichen, ETH Zürich,
// 2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
@@ -38,52 +40,36 @@ EllipticDomain::EllipticDomain(Vector_t nr, Vector_t hr) {
setHr(hr);
}
-EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl) {
- SemiMajor = semimajor;
- SemiMinor = semiminor;
+EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
+ Vector_t hr, std::string interpl)
+ : semiMajor_m(semimajor)
+ , semiMinor_m(semiminor)
+{
setNr(nr);
setHr(hr);
- if(interpl == "CONSTANT")
- interpolationMethod = CONSTANT;
- else if(interpl == "LINEAR")
- interpolationMethod = LINEAR;
- else if(interpl == "QUADRATIC")
- interpolationMethod = QUADRATIC;
+ if (interpl == "CONSTANT")
+ interpolationMethod_m = CONSTANT;
+ else if (interpl == "LINEAR")
+ interpolationMethod_m = LINEAR;
+ else if (interpl == "QUADRATIC")
+ interpolationMethod_m = QUADRATIC;
}
-EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl) {
- SemiMajor = bgeom->getA();
- SemiMinor = bgeom->getB();
+EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
+ std::string interpl)
+ : EllipticDomain(bgeom->getA(), bgeom->getB(), nr, hr, interpl)
+{
setMinMaxZ(bgeom->getS(), bgeom->getLength());
- setNr(nr);
- setHr(hr);
-
- if(interpl == "CONSTANT")
- interpolationMethod = CONSTANT;
- else if(interpl == "LINEAR")
- interpolationMethod = LINEAR;
- else if(interpl == "QUADRATIC")
- interpolationMethod = QUADRATIC;
}
-EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl) {
- SemiMajor = bgeom->getA();
- SemiMinor = bgeom->getB();
- setMinMaxZ(bgeom->getS(), bgeom->getLength());
- Vector_t hr_m;
- hr_m[0] = (getXRangeMax()-getXRangeMin())/nr[0];
- hr_m[1] = (getYRangeMax()-getYRangeMin())/nr[1];
- hr_m[2] = (getZRangeMax()-getZRangeMin())/nr[2];
- setHr(hr_m);
-
- if(interpl == "CONSTANT")
- interpolationMethod = CONSTANT;
- else if(interpl == "LINEAR")
- interpolationMethod = LINEAR;
- else if(interpl == "QUADRATIC")
- interpolationMethod = QUADRATIC;
-}
+EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl)
+ : EllipticDomain(bgeom, nr,
+ Vector_t(2.0 * bgeom->getA() / nr[0],
+ 2.0 * bgeom->getB() / nr[1],
+ (bgeom->getLength() - bgeom->getS()) / nr[2]),
+ interpl)
+{ }
EllipticDomain::~EllipticDomain() {
//nothing so far
@@ -93,164 +79,115 @@ EllipticDomain::~EllipticDomain() {
// for this geometry we only have to calculate the intersection on
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
-void EllipticDomain::compute(Vector_t hr){
- //there is nothing to be done if the mesh spacings have not changed
- if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
+void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
+ // there is nothing to be done if the mesh spacings in x and y have not changed
+ if (hr[0] == getHr()[0] &&
+ hr[1] == getHr()[1])
+ {
hasGeometryChanged_m = false;
return;
}
+
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
- IdxMap.clear();
- CoordMap.clear();
+ idxMap_m.clear();
+ coordMap_m.clear();
//clear previous intersection points
- IntersectYDir.clear();
- IntersectXDir.clear();
+ intersectYDir_m.clear();
+ intersectXDir_m.clear();
// build a index and coordinate map
int idx = 0;
int x, y;
- for(x = 0; x < nr[0]; x++) {
- for(y = 0; y < nr[1]; y++) {
- if(isInside(x, y, 1)) {
- IdxMap[toCoordIdx(x, y)] = idx;
- CoordMap[idx++] = toCoordIdx(x, y);
+ /* we need to scan the full x-y-plane on all cores
+ * in order to figure out the number of valid
+ * grid points per plane --> otherwise we might
+ * get not unique global indices in the Tpetra::CrsMatrix
+ */
+ for (x = 0; x < nr[0]; ++x) {
+ for (y = 0; y < nr[1]; ++y) {
+ if (isInside(x, y, 1)) {
+ idxMap_m[toCoordIdx(x, y)] = idx;
+ coordMap_m[idx++] = toCoordIdx(x, y);
nxy_m++;
}
-
}
}
- switch(interpolationMethod) {
+ switch (interpolationMethod_m) {
case CONSTANT:
break;
case LINEAR:
case QUADRATIC:
- double smajsq = SemiMajor * SemiMajor;
- double sminsq = SemiMinor * SemiMinor;
+ double smajsq = semiMajor_m * semiMajor_m;
+ double sminsq = semiMinor_m * semiMinor_m;
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
- double mx = (nr[0] - 1) * hr[0] / 2.0;
- double my = (nr[1] - 1) * hr[1] / 2.0;
- //calculate intersection with the ellipse
- for(x = 0; x < nr[0]; x++) {
- pos = x * hr[0] - mx;
- if (pos <= -SemiMajor || pos >= SemiMajor)
+ // calculate intersection with the ellipse
+ for (x = localId[0].first(); x <= localId[0].last(); x++) {
+ pos = - semiMajor_m + hr[0] * (x + 0.5);
+ if (std::abs(pos) >= semiMajor_m)
{
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- }else{
- yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
- IntersectYDir.insert(std::pair<int, double>(x, yd));
- IntersectYDir.insert(std::pair<int, double>(x, -yd));
+ intersectYDir_m.insert(std::pair<int, double>(x, 0));
+ intersectYDir_m.insert(std::pair<int, double>(x, 0));
+ } else {
+ yd = std::abs(std::sqrt(sminsq - sminsq * pos * pos / smajsq));
+ intersectYDir_m.insert(std::pair<int, double>(x, yd));
+ intersectYDir_m.insert(std::pair<int, double>(x, -yd));
}
}
- for(y = 0; y < nr[1]; y++) {
- pos = y * hr[1] - my;
- if (pos <= -SemiMinor || pos >= SemiMinor)
+ for (y = localId[0].first(); y < localId[1].last(); y++) {
+ pos = - semiMinor_m + hr[1] * (y + 0.5);
+ if (std::abs(pos) >= semiMinor_m)
{
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- }else{
- xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
- IntersectXDir.insert(std::pair<int, double>(y, xd));
- IntersectXDir.insert(std::pair<int, double>(y, -xd));
+ intersectXDir_m.insert(std::pair<int, double>(y, 0));
+ intersectXDir_m.insert(std::pair<int, double>(y, 0));
+ } else {
+ xd = std::abs(std::sqrt(smajsq - smajsq * pos * pos / sminsq));
+ intersectXDir_m.insert(std::pair<int, double>(y, xd));
+ intersectXDir_m.insert(std::pair<int, double>(y, -xd));
}
}
}
}
-void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
- //there is nothing to be done if the mesh spacings have not changed
- if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
- hasGeometryChanged_m = false;
- return;
- }
- setHr(hr);
- hasGeometryChanged_m = true;
- //reset number of points inside domain
- nxy_m = 0;
-
- // clear previous coordinate maps
- IdxMap.clear();
- CoordMap.clear();
- //clear previous intersection points
- IntersectYDir.clear();
- IntersectXDir.clear();
-
- // build a index and coordinate map
- int idx = 0;
- int x, y;
-
- for(x = localId[0].first(); x<= localId[0].last(); x++) {
- for(y = localId[1].first(); y <= localId[1].last(); y++) {
- if(isInside(x, y, 1)) {
- IdxMap[toCoordIdx(x, y)] = idx;
- CoordMap[idx++] = toCoordIdx(x, y);
- nxy_m++;
- }
- }
- }
-
- switch(interpolationMethod) {
- case CONSTANT:
- break;
- case LINEAR:
- case QUADRATIC:
+void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+ const Vector_t& rmax, double dh)
+{
+ Vector_t mymax = Vector_t(0.0, 0.0, 0.0);
- double smajsq = SemiMajor * SemiMajor;
- double sminsq = SemiMinor * SemiMinor;
- double yd = 0.0;
- double xd = 0.0;
- double pos = 0.0;
- double mx = (nr[0] - 1) * hr[0] / 2.0;
- double my = (nr[1] - 1) * hr[1] / 2.0;
+ // apply bounding box increment dh, i.e., "BBOXINCR" input argument
+ double zsize = rmax[2] - rmin[2];
- //calculate intersection with the ellipse
- for(x = localId[0].first(); x <= localId[0].last(); x++) {
- pos = x * hr[0] - mx;
- if (pos <= -SemiMajor || pos >= SemiMajor)
- {
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- }else{
- yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
- IntersectYDir.insert(std::pair<int, double>(x, yd));
- IntersectYDir.insert(std::pair<int, double>(x, -yd));
- }
+ setMinMaxZ(rmin[2] - zsize * (1.0 + dh),
+ rmax[2] + zsize * (1.0 + dh));
- }
+ origin = Vector_t(-semiMajor_m, -semiMinor_m, getMinZ());
+ mymax = Vector_t( semiMajor_m, semiMinor_m, getMaxZ());
- for(y = localId[0].first(); y < localId[1].last(); y++) {
- pos = y * hr[1] - my;
- if (pos <= -SemiMinor || pos >= SemiMinor)
- {
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- }else{
- xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
- IntersectXDir.insert(std::pair<int, double>(y, xd));
- IntersectXDir.insert(std::pair<int, double>(y, -xd));
- }
- }
- }
+ for (int i = 0; i < 3; ++i)
+ hr[i] = (mymax[i] - origin[i]) / nr[i];
}
-void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
+void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- // determine which interpolation method we use for points near the boundary
- switch(interpolationMethod) {
+ // determine which interpolation method we use for points near the boundary
+ switch (interpolationMethod_m) {
case CONSTANT:
constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
break;
@@ -266,9 +203,10 @@ void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &
assert(C > 0);
}
-void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
-
-
+void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E,
+ double &S, double &N, double &F,
+ double &B, double &C, double &scaleFactor)
+{
int x = 0, y = 0, z = 0;
getCoord(idx, x, y, z);
getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
@@ -283,113 +221,99 @@ void EllipticDomain::getNeighbours(int idx, int &W, int &E, int &S, int &N, int
}
-void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) {
-
- if(x > 0)
+void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E,
+ int &S, int &N, int &F, int &B)
+{
+ if (x > 0)
W = getIdx(x - 1, y, z);
else
W = -1;
- if(x < nr[0] - 1)
+
+ if (x < nr[0] - 1)
E = getIdx(x + 1, y, z);
else
E = -1;
- if(y < nr[1] - 1)
+ if (y < nr[1] - 1)
N = getIdx(x, y + 1, z);
else
N = -1;
- if(y > 0)
+
+ if (y > 0)
S = getIdx(x, y - 1, z);
else
S = -1;
- if(z > 0)
+ if (z > 0)
F = getIdx(x, y, z - 1);
else
F = -1;
- if(z < nr[2] - 1)
+
+ if (z < nr[2] - 1)
B = getIdx(x, y, z + 1);
else
B = -1;
}
-void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV, double &EV, double &SV, double &NV, double &FV, double &BV, double &CV, double &scaleFactor) {
-
+void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV,
+ double &EV, double &SV, double &NV,
+ double &FV, double &BV, double &CV,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- WV = -1/(hr[0]*hr[0]);
- EV = -1/(hr[0]*hr[0]);
- NV = -1/(hr[1]*hr[1]);
- SV = -1/(hr[1]*hr[1]);
- FV = -1/(hr[2]*hr[2]);
- BV = -1/(hr[2]*hr[2]);
- CV = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
+ WV = -1.0 / (hr[0] * hr[0]);
+ EV = -1.0 / (hr[0] * hr[0]);
+ NV = -1.0 / (hr[1] * hr[1]);
+ SV = -1.0 / (hr[1] * hr[1]);
+ FV = -1.0 / (hr[2] * hr[2]);
+ BV = -1.0 / (hr[2] * hr[2]);
+
+ CV = 2.0 / (hr[0] * hr[0])
+ + 2.0 / (hr[1] * hr[1])
+ + 2.0 / (hr[2] * hr[2]);
// we are a right boundary point
- if(!isInside(x + 1, y, z))
- EV= 0.0;
+ if (!isInside(x + 1, y, z))
+ EV = 0.0;
// we are a left boundary point
- if(!isInside(x - 1, y, z))
+ if (!isInside(x - 1, y, z))
WV = 0.0;
// we are a upper boundary point
- if(!isInside(x, y + 1, z))
+ if (!isInside(x, y + 1, z))
NV = 0.0;
// we are a lower boundary point
- if(!isInside(x, y - 1, z))
+ if (!isInside(x, y - 1, z))
SV = 0.0;
- if(z == 1 || z == nr[2] - 2) {
-
- // case where we are on the Robin BC in Z-direction
- // where we distinguish two cases
- // IFF: this values should not matter because they
- // never make it into the discretization matrix
- if(z == 1)
- FV = 0.0;
- else
- BV = 0.0;
-
- // add contribution of Robin discretization to center point
- // d the distance between the center of the bunch and the boundary
- //double cx = (x-(nr[0]-1)/2)*hr[0];
- //double cy = (y-(nr[1]-1)/2)*hr[1];
- //double cz = hr[2]*(nr[2]-1);
- //double d = sqrt(cx*cx+cy*cy+cz*cz);
- double d = hr[2] * (nr[2] - 1) / 2;
- CV += 2 / (d * hr[2]);
- //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
-
- // scale all stencil-points in z-plane with 0.5 (Robin discretization)
- WV /= 2.0;
- EV /= 2.0;
- NV /= 2.0;
- SV /= 2.0;
- CV /= 2.0;
- scaleFactor *= 0.5;
- }
-
+ // handle boundary condition in z direction
+ robinBoundaryStencil(z, FV, BV, CV);
}
-void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
-
+void EllipticDomain::linearInterpolation(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
- double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
+ double cx = - semiMajor_m + hr[0] * (x + 0.5);
+ double cy = - semiMinor_m + hr[1] * (y + 0.5);
double dx = 0.0;
- std::multimap<int, double>::iterator it = IntersectXDir.find(y);
- if(cx < 0)
+ std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
+
+ if (cx < 0)
++it;
dx = it->second;
double dy = 0.0;
- it = IntersectYDir.find(x);
- if(cy < 0)
+ it = intersectYDir_m.find(x);
+ if (cy < 0)
++it;
dy = it->second;
@@ -400,73 +324,55 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double
double ds = hr[1];
C = 0.0;
- //we are a right boundary point
- if(!isInside(x + 1, y, z)) {
- C += 1 / ((dx - cx) * de);
+ // we are a right boundary point
+ if (!isInside(x + 1, y, z)) {
+ C += 1.0 / ((dx - cx) * de);
E = 0.0;
} else {
- C += 1 / (de * de);
- E = -1 / (de * de);
+ C += 1.0 / (de * de);
+ E = -1.0 / (de * de);
}
- //we are a left boundary point
- if(!isInside(x - 1, y, z)) {
- C += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
+ // we are a left boundary point
+ if (!isInside(x - 1, y, z)) {
+ C += 1.0 / ((std::abs(dx) - std::abs(cx)) * dw);
W = 0.0;
} else {
- C += 1 / (dw * dw);
- W = -1 / (dw * dw);
+ C += 1.0 / (dw * dw);
+ W = -1.0 / (dw * dw);
}
- //we are a upper boundary point
- if(!isInside(x, y + 1, z)) {
- C += 1 / ((dy - cy) * dn);
+ // we are a upper boundary point
+ if (!isInside(x, y + 1, z)) {
+ C += 1.0 / ((dy - cy) * dn);
N = 0.0;
} else {
- C += 1 / (dn * dn);
- N = -1 / (dn * dn);
+ C += 1.0 / (dn * dn);
+ N = -1.0 / (dn * dn);
}
- //we are a lower boundary point
- if(!isInside(x, y - 1, z)) {
- C += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
+ // we are a lower boundary point
+ if (!isInside(x, y - 1, z)) {
+ C += 1.0 / ((std::abs(dy) - std::abs(cy)) * ds);
S = 0.0;
} else {
- C += 1 / (ds * ds);
- S = -1 / (ds * ds);
+ C += 1.0 / (ds * ds);
+ S = -1.0 / (ds * ds);
}
- F = -1 / (hr[2] * hr[2]);
- B = -1 / (hr[2] * hr[2]);
- C += 2 / (hr[2] * hr[2]);
-
// handle boundary condition in z direction
-/*
- if(z == 0 || z == nr[2] - 1) {
-
- // case where we are on the NEUMAN BC in Z-direction
- // where we distinguish two cases
- if(z == 0)
- F = 0.0;
- else
- B = 0.0;
-
- //hr[2]*(nr2[2]-1)/2 = radius
- double d = hr[2] * (nr[2] - 1) / 2;
- C += 2 / (d * hr[2]);
-
- W /= 2.0;
- E /= 2.0;
- N /= 2.0;
- S /= 2.0;
- C /= 2.0;
- scaleFactor *= 0.5;
-
- }
-*/
+ F = -1.0 / (hr[2] * hr[2]);
+ B = -1.0 / (hr[2] * hr[2]);
+ C += 2.0 / (hr[2] * hr[2]);
+ robinBoundaryStencil(z, F, B, C);
}
-void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
+void EllipticDomain::quadraticInterpolation(int x, int y, int z,
+ double &W, double &E, double &S,
+ double &N, double &F, double &B, double &C,
+ double &scaleFactor)
+{
+ scaleFactor = 1.0;
double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
@@ -474,14 +380,14 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, doub
// since every vector for elliptic domains has ALWAYS size 2 we
// can catch all cases manually
double dx = 0.0;
- std::multimap<int, double>::iterator it = IntersectXDir.find(y);
- if(cx < 0)
+ std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
+ if (cx < 0)
++it;
dx = it->second;
double dy = 0.0;
- it = IntersectYDir.find(x);
- if(cy < 0)
+ it = intersectYDir_m.find(x);
+ if (cy < 0)
++it;
dy = it->second;
@@ -489,197 +395,80 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, doub
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
- W = 1.0;
- E = 1.0;
- N = 1.0;
- S = 1.0;
- F = 1.0;
- B = 1.0;
+
+ W = 0.0;
+ E = 0.0;
+ S = 0.0;
+ N = 0.0;
C = 0.0;
- //TODO: = cx+hr[0] > dx && cx > 0
- //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
- ////we are a right boundary point
- ////if(!isInside(x+1,y,z)) {
- //de = dx-cx;
- //}
-
- //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
- ////we are a left boundary point
- ////if(!isInside(x-1,y,z)) {
- //dw = std::abs(dx)-std::abs(cx);
- //}
-
- //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
- ////we are a upper boundary point
- ////if(!isInside(x,y+1,z)) {
- //dn = dy-cy;
- //}
-
- //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
- ////we are a lower boundary point
- ////if(!isInside(x,y-1,z)) {
- //ds = std::abs(dy)-std::abs(cy);
- //}
-
- //TODO: = cx+hr[0] > dx && cx > 0
- //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
- //we are a right boundary point
- if(!isInside(x + 1, y, z)) {
- de = dx - cx;
+ // we are a right boundary point
+ if (!isInside(x + 1, y, z)) {
+ double s = dx - cx;
+ C -= 2.0 * (s - 1.0) / (s * de);
E = 0.0;
+ W += (s - 1.0) / ((s + 1.0) * de);
+ } else {
+ C += 1.0 / (de * de);
+ E = -1.0 / (de * de);
}
- //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
- //we are a left boundary point
- if(!isInside(x - 1, y, z)) {
- dw = std::abs(dx) - std::abs(cx);
+ // we are a left boundary point
+ if (!isInside(x - 1, y, z)) {
+ double s = std::abs(dx) - std::abs(cx);
+ C -= 2.0 * (s - 1.0) / (s * de);
W = 0.0;
+ E += (s - 1.0) / ((s + 1.0) * de);
+ } else {
+ C += 1.0 / (dw * dw);
+ W = -1.0 / (dw * dw);
}
- //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
- //we are a upper boundary point
- if(!isInside(x, y + 1, z)) {
- dn = dy - cy;
+ // we are a upper boundary point
+ if (!isInside(x, y + 1, z)) {
+ double s = dy - cy;
+ C -= 2.0 * (s - 1.0) / (s * dn);
N = 0.0;
+ S += (s - 1.0) / ((s + 1.0) * dn);
+ } else {
+ C += 1.0 / (dn * dn);
+ N = -1.0 / (dn * dn);
}
- //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
- //we are a lower boundary point
- if(!isInside(x, y - 1, z)) {
- ds = std::abs(dy) - std::abs(cy);
+ // we are a lower boundary point
+ if (!isInside(x, y - 1, z)) {
+ double s = std::abs(dy) - std::abs(cy);
+ C -= 2.0 * (s - 1.0) / (s * dn);
S = 0.0;
+ N += (s - 1.0) / ((s + 1.0) * dn);
+ } else {
+ C += 1.0 / (ds * ds);
+ S = -1.0 / (ds * ds);
}
- //2/dw*(dw_de)
- W *= -1.0 / (dw * (dw + de));
- E *= -1.0 / (de * (dw + de));
- N *= -1.0 / (dn * (dn + ds));
- S *= -1.0 / (ds * (dn + ds));
- F = -1 / (hr[2] * (hr[2] + hr[2]));
- B = -1 / (hr[2] * (hr[2] + hr[2]));
-
- //TODO: problem when de,dw,dn,ds == 0
- //is NOT a regular BOUND PT
- C += 1 / de * 1 / dw;
- C += 1 / dn * 1 / ds;
- C += 1 / hr[2] * 1 / hr[2];
- scaleFactor = 0.5;
-
-
- //for regular gridpoints no problem with symmetry, just boundary
- //z direction is right
- //implement isLastInside(dir)
- //we have LOCAL x,y coordinates!
-
- /*
- if(dw != 0 && !wIsB)
- W = -1/dw * (dn+ds) * 2*hr[2];
- else
- W = 0;
- if(de != 0 && !eIsB)
- E = -1/de * (dn+ds) * 2*hr[2];
- else
- E = 0;
- if(dn != 0 && !nIsB)
- N = -1/dn * (dw+de) * 2*hr[2];
- else
- N = 0;
- if(ds != 0 && !sIsB)
- S = -1/ds * (dw+de) * 2*hr[2];
- else
- S = 0;
- F = -(dw+de)*(dn+ds)/hr[2];
- B = -(dw+de)*(dn+ds)/hr[2];
- */
-
- //if(dw != 0)
- //W = -2*hr[2]*(dn+ds)/dw;
- //else
- //W = 0;
- //if(de != 0)
- //E = -2*hr[2]*(dn+ds)/de;
- //else
- //E = 0;
- //if(dn != 0)
- //N = -2*hr[2]*(dw+de)/dn;
- //else
- //N = 0;
- //if(ds != 0)
- //S = -2*hr[2]*(dw+de)/ds;
- //else
- //S = 0;
- //F = -(dw+de)*(dn+ds)/hr[2];
- //B = -(dw+de)*(dn+ds)/hr[2];
-
- //// RHS scaleFactor for current 3D index
- //// Factor 0.5 results from the SW/quadratic extrapolation
- //scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr[2]);
-
- // catch the case where a point lies on the boundary
- //FIXME: do this more elegant!
- //double m1 = dw*de;
- //double m2 = dn*ds;
- //if(de == 0)
- //m1 = dw;
- //if(dw == 0)
- //m1 = de;
- //if(dn == 0)
- //m2 = ds;
- //if(ds == 0)
- //m2 = dn;
- ////XXX: dn+ds || dw+de can be 0
- ////C = 2*(dn+ds)*(dw+de)/hr[2];
- //C = 2/hr[2];
- //if(dw != 0 || de != 0)
- //C *= (dw+de);
- //if(dn != 0 || ds != 0)
- //C *= (dn+ds);
- //if(dw != 0 || de != 0)
- //C += (2*hr[2])*(dn+ds)*(dw+de)/m1;
- //if(dn != 0 || ds != 0)
- //C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
-
- //handle Neumann case
- //if(z == 0 || z == nr[2]-1) {
-
- //if(z == 0)
- //F = 0.0;
- //else
- //B = 0.0;
-
- ////neumann stuff
- //W = W/2.0;
- //E = E/2.0;
- //N = N/2.0;
- //S = S/2.0;
- //C /= 2.0;
-
- //scaleFactor /= 2.0;
- //}
-
// handle boundary condition in z direction
- if(z == 0 || z == nr[2] - 1) {
+ F = -1.0 / (hr[2] * hr[2]);
+ B = -1.0 / (hr[2] * hr[2]);
+ C += 2.0 / (hr[2] * hr[2]);
+ robinBoundaryStencil(z, F, B, C);
+}
+
+void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) {
+ if (z == 0 || z == nr[2] - 1) {
- // case where we are on the NEUMAN BC in Z-direction
+ // case where we are on the Robin BC in Z-direction
// where we distinguish two cases
- if(z == 0)
+ // IFF: this values should not matter because they
+ // never make it into the discretization matrix
+ if (z == 0)
F = 0.0;
else
B = 0.0;
- //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
- //hr[2]*(nr2[2]-1)/2 = radius
- double d = hr[2] * (nr[2] - 1) / 2;
- C += 2 / (d * hr[2]);
-
- W /= 2.0;
- E /= 2.0;
- N /= 2.0;
- S /= 2.0;
- C /= 2.0;
- scaleFactor /= 2.0;
-
+ // add contribution of Robin discretization to center point
+ // d the distance between the center of the bunch and the boundary
+ double d = 0.5 * hr[2] * (nr[2] - 1);
+ C += 2.0 / (d * hr[2]);
}
}
diff --git a/src/Solvers/EllipticDomain.h b/src/Solvers/EllipticDomain.h
index 691e82cffb148c0180103eef7fd127c843681fb2..e5478b5ae1fe0cdd90739b84cffa1c9232a5449b 100644
--- a/src/Solvers/EllipticDomain.h
+++ b/src/Solvers/EllipticDomain.h
@@ -1,6 +1,8 @@
//
// Class EllipticDomain
-// :FIXME: add brief description
+// This class provides an elliptic beam pipe. The mesh adapts to the bunch size
+// in the longitudinal direction. At the intersection of the mesh with the
+// beam pipe, three stencil interpolation methods are available.
//
// Copyright (c) 2008, Yves Ineichen, ETH Zürich,
// 2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
@@ -32,100 +34,127 @@
#include <cmath>
#include "IrregularDomain.h"
#include "Structure/BoundaryGeometry.h"
+#include "Utilities/OpalException.h"
class EllipticDomain : public IrregularDomain {
public:
EllipticDomain(Vector_t nr, Vector_t hr);
- EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl);
- EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
+
+ EllipticDomain(double semimajor, double semiminor, Vector_t nr,
+ Vector_t hr, std::string interpl);
+
+ EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
+ Vector_t hr, std::string interpl);
+
EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
~EllipticDomain();
/// returns discretization at (x,y,z)
- void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+ void getBoundaryStencil(int x, int y, int z,
+ double &W, double &E, double &S,
+ double &N, double &F, double &B,
+ double &C, double &scaleFactor);
+
/// returns discretization at 3D index
- void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+ void getBoundaryStencil(int idx, double &W, double &E,
+ double &S, double &N, double &F,
+ double &B, double &C, double &scaleFactor);
+
/// returns index of neighbours at (x,y,z)
- void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B);
+ void getNeighbours(int x, int y, int z, int &W, int &E,
+ int &S, int &N, int &F, int &B);
+
/// returns index of neighbours at 3D index
void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
+
/// returns type of boundary condition
std::string getType() {return "Elliptic";}
+
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z) {
- double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
- double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
- return ((xx * xx / (SemiMajor * SemiMajor) + yy * yy / (SemiMinor * SemiMinor) < 1) && z != 0 && z != nr[2] - 1);
+ double xx = - semiMajor_m + hr[0] * (x + 0.5);
+ double yy = - semiMinor_m + hr[1] * (y + 0.5);
+
+ bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
+ yy * yy / (semiMinor_m * semiMinor_m) < 1);
+
+ return (isInsideEllipse && z > 0 && z < nr[2] - 1);
}
int getNumXY(int /*z*/) { return nxy_m; }
- /// set semi-minor
- void setSemiMinor(double sm) {SemiMinor = sm;}
- /// set semi-major
- void setSemiMajor(double sm) {SemiMajor = sm;}
- /// calculates intersection
- void compute(Vector_t hr);
+
+ /// calculates intersection
+ void compute(Vector_t /*hr*/) {
+ throw OpalException("EllipticDomain::compute()", "This function is not available.");
+ }
+
void compute(Vector_t hr, NDIndex<3> localId);
- double getXRangeMin() { return -SemiMajor; }
- double getXRangeMax() { return SemiMajor; }
- double getYRangeMin() { return -SemiMinor; }
- double getYRangeMax() { return SemiMinor; }
+ double getXRangeMin() { return -semiMajor_m; }
+ double getXRangeMax() { return semiMajor_m; }
+ double getYRangeMin() { return -semiMinor_m; }
+ double getYRangeMax() { return semiMinor_m; }
double getZRangeMin() { return zMin_m; }
double getZRangeMax() { return zMax_m; }
+ bool hasGeometryChanged() { return hasGeometryChanged_m; }
- //TODO: ?
- int getStartIdx() {return 0;}
- bool hasGeometryChanged() { return hasGeometryChanged_m; }
+ void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+ const Vector_t& rmax, double dh);
private:
/// Map from a single coordinate (x or y) to a list of intersection values with
/// boundary.
- typedef std::multimap<int, double> EllipticPointList;
+ typedef std::multimap<int, double> EllipticPointList_t;
/// all intersection points with grid lines in X direction
- EllipticPointList IntersectXDir;
+ EllipticPointList_t intersectXDir_m;
/// all intersection points with grid lines in Y direction
- EllipticPointList IntersectYDir;
+ EllipticPointList_t intersectYDir_m;
/// mapping (x,y,z) -> idx
- std::map<int, int> IdxMap;
+ std::map<int, int> idxMap_m;
/// mapping idx -> (x,y,z)
- std::map<int, int> CoordMap;
+ std::map<int, int> coordMap_m;
/// semi-major of the ellipse
- double SemiMajor;
+ double semiMajor_m;
+
/// semi-minor of the ellipse
- double SemiMinor;
+ double semiMinor_m;
+
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
+
/// interpolation type
- int interpolationMethod;
+ int interpolationMethod_m;
+
/// flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
/// conversion from (x,y) to index in xy plane
inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
+
/// conversion from (x,y,z) to index on the 3D grid
inline int getIdx(int x, int y, int z) {
- if(isInside(x, y, z) && x >= 0 && y >= 0 && z >= 0)
- return IdxMap[toCoordIdx(x, y)] + (z - 1) * nxy_m;
+ if (isInside(x, y, z))
+ return idxMap_m[toCoordIdx(x, y)] + (z - 1) * nxy_m;
else
return -1;
}
+
/// conversion from a 3D index to (x,y,z)
inline void getCoord(int idx, int &x, int &y, int &z) {
int ixy = idx % nxy_m;
- int xy = CoordMap[ixy];
+ int xy = coordMap_m[ixy];
int inr = (int)nr[0];
x = xy % inr;
y = (xy - x) / nr[0];
@@ -133,10 +162,23 @@ private:
}
/// different interpolation methods for boundary points
- void constantInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
- void linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
- void quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
-
+ void constantInterpolation(int x, int y, int z,
+ double &W, double &E, double &S,
+ double &N, double &F, double &B,
+ double &C, double &scaleFactor);
+
+ void linearInterpolation(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor);
+
+ void quadraticInterpolation(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor);
+
+ /// function to handle the open boundary condition in longitudinal direction
+ void robinBoundaryStencil(int z, double &F, double &B, double &C);
};
#endif //#ifdef ELLIPTICAL_DOMAIN_H
diff --git a/src/Solvers/IrregularDomain.h b/src/Solvers/IrregularDomain.h
index 4569f663d9cfbd9a3bdb78ce6875b5b69a438fd0..eb00435e1772f8dca9f019ebc0848650c29e0047 100644
--- a/src/Solvers/IrregularDomain.h
+++ b/src/Solvers/IrregularDomain.h
@@ -129,6 +129,22 @@ public:
virtual bool hasGeometryChanged() = 0;
virtual ~IrregularDomain() {};
+ virtual void resizeMesh(Vector_t& origin, Vector_t& hr,
+ const Vector_t& /*rmin*/, const Vector_t& /*rmax*/, double /*dh*/) {
+ double xmin = getXRangeMin();
+ double xmax = getXRangeMax();
+ double ymin = getYRangeMin();
+ double ymax = getYRangeMax();
+ double zmin = getZRangeMin();
+ double zmax = getZRangeMax();
+
+ origin = Vector_t(xmin, ymin , zmin);
+ Vector_t mymax = Vector_t(xmax, ymax , zmax);
+
+ for (int i = 0; i < 3; i++)
+ hr[i] = (mymax[i] - origin[i]) / nr[i];
+ };
+
protected:
// a irregular domain is always defined on a grid
diff --git a/src/Solvers/MGPoissonSolver.cpp b/src/Solvers/MGPoissonSolver.cpp
index 57836bfd6ea0634b7c4b8ebb34c2525aefff80b8..ddcf129c7ebe3b38e3ae1374d1de41320d266817 100644
--- a/src/Solvers/MGPoissonSolver.cpp
+++ b/src/Solvers/MGPoissonSolver.cpp
@@ -353,8 +353,10 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
+ NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
if (bp_m->isInside(idx, idy, idz))
- RHS->getDataNonConst()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
+ RHS->replaceGlobalValue(bp_m->getIdx(idx, idy, idz),
+ 4.0 * M_PI * rho.localElement(l) / scaleFactor);
}
}
}
@@ -508,6 +510,7 @@ void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
}
}
}
+
int indexbase = 0;
map_p = Teuchos::rcp(new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
&MyGlobalElements[0], NumMyElements, indexbase, comm_mp));
diff --git a/src/Solvers/MGPoissonSolver.h b/src/Solvers/MGPoissonSolver.h
index 44bf05aca7449937087cb2a7c415721e9972b1fe..c2e1de9a9fe87dea819218f1b629574b29cd4234 100644
--- a/src/Solvers/MGPoissonSolver.h
+++ b/src/Solvers/MGPoissonSolver.h
@@ -125,8 +125,16 @@ public:
void extrapolateLHS();
+ void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+ const Vector_t& rmax, double dh)
+ {
+ bp_m->resizeMesh(origin, hr, rmin, rmax, dh);
+ }
+
Inform &print(Inform &os) const;
+
+
private:
bool isMatrixfilled_m;
@@ -171,6 +179,7 @@ private:
/// Map holding the processor distribution of data
Teuchos::RCP<TpetraMap_t> map_p;
+
/// communicator used by Trilinos
Teuchos::RCP<const Comm_t> comm_mp;
diff --git a/src/Solvers/PoissonSolver.h b/src/Solvers/PoissonSolver.h
index e2dd5a350b3ec2c1f8a9c2f65e07f8c9592891d6..84aa754cee851194d88d2da0b47501c89b172ed4 100644
--- a/src/Solvers/PoissonSolver.h
+++ b/src/Solvers/PoissonSolver.h
@@ -58,6 +58,11 @@ public:
virtual void test(PartBunchBase<double, 3> *bunch) = 0 ;
virtual ~PoissonSolver(){};
+ virtual void resizeMesh(Vector_t& /*origin*/, Vector_t& /*hr*/,
+ const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
+ double /*dh*/)
+ { };
+
};
inline Inform &operator<<(Inform &os, const PoissonSolver &/*fs*/) {
diff --git a/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp b/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
index e6b4a98ddd711a46cac0ae49fd5612973f4bba94..b839e4d8009c6adc1e2a0670f7400b738bfe66c8 100644
--- a/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
+++ b/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
@@ -25,10 +25,10 @@ TEST(GreenWakeFunctionTest, TestApply)
bool const_length = true;
std::string fname = "";
- std::vector<double> finalWakeValues = {1.61757e+10, 2.78525e+19};
- std::vector<double> finalEnergyValues = {-0.00125303, -2.15739e+06};
+ std::vector<double> finalWakeValues = { 1.61757e+10, 2.78525e+19};
+ std::vector<double> finalEnergyValues = {-2.446072e-5, -42118.09};
std::vector<double> relativeErrorWake = {1e+6, 1e+15};
- std::vector<double> relativeErrorEnergy = {1e-7, 1e+2};
+ std::vector<double> relativeErrorEnergy = {1e-9, 1};
for (int acmode : acmodes) {
GreenWakeFunction gwf("opal", nullptr, filters, nbin, Z0, radius, sigma, acmode, tau, 0, const_length, fname);
@@ -37,7 +37,7 @@ TEST(GreenWakeFunctionTest, TestApply)
// determine K and charge
double charge = 0.8e-9; // nC
double K = 0.20536314319923724e-9; //K normalizes nC data in lambda.h?
- gwf.NBin_m = 294;
+ gwf.NBin_m = 10;
std::cout << "# Z0 = " << gwf.Z0_m << std::endl
<< "# radius = " << gwf.radius_m << std::endl