Commit 3ff473cc authored by chrisrogers1234's avatar chrisrogers1234
Browse files

Fixes #56 - OpalRing components throw segmentation fault and #59 - Unit tests...

Fixes #56 - OpalRing components throw segmentation fault and #59 - Unit tests throw segmentation fault; partial implementation of #68 DumpEmFields
parent a0898c93
......@@ -73,6 +73,7 @@
#include "Utilities/OpalException.h"
#include "BasicActions/DumpFields.h"
#include "BasicActions/DumpEMFields.h"
#include "Structure/H5PartWrapperForPC.h"
#include "Structure/BoundaryGeometry.h"
......@@ -1145,6 +1146,7 @@ void ParallelCyclotronTracker::execute() {
extE_m = Vector_t(0.0, 0.0, 0.0);
extB_m = Vector_t(0.0, 0.0, 0.0);
DumpFields::writeFields((((*FieldDimensions.begin())->second).second));
DumpEMFields::writeFields((((*FieldDimensions.begin())->second).second));
if(timeIntegrator_m == 0) {
*gmsg << "* 4th order Runge-Kutta integrator" << endl;
......@@ -5573,4 +5575,4 @@ void ParallelCyclotronTracker::evaluateSpaceChargeField() {
localToGlobal(itsBunch->Bf, phi);
localToGlobal(itsBunch->R, phi, meanR);
}
}
\ No newline at end of file
}
......@@ -2,6 +2,7 @@ set (_SRCS
Call.cpp
Dump.cpp
DumpFields.cpp
DumpEMFields.cpp
Echo.cpp
Help.cpp
Option.cpp
......@@ -26,6 +27,7 @@ add_opal_sources(${_SRCS})
set (HDRS
Call.h
DumpFields.h
DumpEMFields.h
Dump.h
Echo.h
Help.h
......@@ -42,4 +44,4 @@ set (HDRS
What.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/BasicActions")
\ No newline at end of file
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/BasicActions")
/*
* Copyright (c) 2016, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <fstream>
#include "Fields/Interpolation/NDGrid.h" // classic
#include "AbsBeamline/Component.h" // classic
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "BasicActions/DumpEMFields.h"
// extern Inform *gmsg;
std::unordered_set<DumpEMFields*> DumpEMFields::dumpsSet_m;
std::string DumpEMFields::dumpemfields_docstring =
std::string("The \"DUMPEMFIELDS\" statement dumps a field map to a user-defined")+
std::string(" field file, for checking that fields are generated correctly.")+
std::string(" The fields are written out on a grid in Cartesian space and time.");
DumpEMFields::DumpEMFields() :
Action(13, "DUMPEMFIELDS", dumpemfields_docstring.c_str()),
grid_m(NULL), filename_m("") {
// would be nice if "steps" could be integer
itsAttr[0] = Attributes::makeReal
("X_START", "Start point in the grid in x [mm]");
itsAttr[1] = Attributes::makeReal
("DX", "Grid step size in x [mm]");
itsAttr[2] = Attributes::makeReal
("X_STEPS", "Number of steps in x");
itsAttr[3] = Attributes::makeReal
("Y_START", "Start point in the grid in y [mm]");
itsAttr[4] = Attributes::makeReal
("DY", "Grid step size in y [mm]");
itsAttr[5] = Attributes::makeReal
("Y_STEPS", "Number of steps in y");
itsAttr[6] = Attributes::makeReal
("Z_START", "Start point in the grid in z [mm]");
itsAttr[7] = Attributes::makeReal
("DZ", "Grid step size in z [mm]");
itsAttr[8] = Attributes::makeReal
("Z_STEPS", "Number of steps in z");
itsAttr[9] = Attributes::makeReal
("T_START", "Start point in the grid in time [ns]");
itsAttr[10] = Attributes::makeReal
("DT", "Grid step size in time [ns]");
itsAttr[11] = Attributes::makeReal
("T_STEPS", "Number of steps in time");
itsAttr[12] = Attributes::makeString
("FILE_NAME", "Name of the file to which field data is dumped");
}
DumpEMFields::~DumpEMFields() {
if (grid_m != NULL) {
delete grid_m;
}
if (dumpsSet_m.find(this) != dumpsSet_m.end()) {
dumpsSet_m.erase(this);
}
}
DumpEMFields* DumpEMFields::clone(const std::string &name) {
DumpEMFields* dumper = new DumpEMFields();
if (grid_m != NULL) {
dumper->grid_m = grid_m->clone();
}
dumper->filename_m = filename_m;
if (dumpsSet_m.find(this) != dumpsSet_m.end()) {
dumpsSet_m.insert(dumper);
}
return dumper;
}
void DumpEMFields::execute() {
buildGrid();
// the routine for action (OpalParser/OpalParser) calls execute and then
// deletes 'this'; so we must build a copy that lasts until the field maps
// are constructed and we are ready for tracking (which is when the field
// maps are written). Hence the clone call below.
dumpsSet_m.insert(this->clone(""));
}
void DumpEMFields::buildGrid() {
std::vector<double> spacing(4);
std::vector<double> origin(4);
std::vector<int> gridSize(4);
origin[0] = Attributes::getReal(itsAttr[0]);
spacing[0] = Attributes::getReal(itsAttr[1]);
double nx = Attributes::getReal(itsAttr[2]);
origin[1] = Attributes::getReal(itsAttr[3]);
spacing[1] = Attributes::getReal(itsAttr[4]);
double ny = Attributes::getReal(itsAttr[5]);
origin[2] = Attributes::getReal(itsAttr[6]);
spacing[2] = Attributes::getReal(itsAttr[7]);
double nz = Attributes::getReal(itsAttr[8]);
origin[3] = Attributes::getReal(itsAttr[9]);
spacing[3] = Attributes::getReal(itsAttr[10]);
double nt = Attributes::getReal(itsAttr[11]);
checkInt(nx, "X_STEPS");
checkInt(ny, "Y_STEPS");
checkInt(nz, "Z_STEPS");
checkInt(nt, "T_STEPS");
gridSize[0] = nx;
gridSize[1] = ny;
gridSize[2] = nz;
gridSize[3] = nt;
if (grid_m != NULL) {
delete grid_m;
}
grid_m = new interpolation::NDGrid(4, &gridSize[0], &spacing[0], &origin[0]);
filename_m = Attributes::getString(itsAttr[12]);
}
void DumpEMFields::writeFields(Component* field) {
typedef std::unordered_set<DumpEMFields*>::iterator dump_iter;
for (dump_iter it = dumpsSet_m.begin(); it != dumpsSet_m.end(); ++it) {
(*it)->writeFieldThis(field);
}
}
void DumpEMFields::checkInt(double real, std::string name, double tolerance) {
if (fabs(floor(real) - real) > tolerance) {
throw OpalException("DumpEMFields::checkInt",
"Value for "+name+
" should be an integer but a real value was found");
}
if (floor(real) < 0.5) {
throw OpalException("DumpEMFields::checkInt",
"Value for "+name+" should be 1 or more");
}
}
void DumpEMFields::writeFieldThis(Component* field) {
if (grid_m == NULL) {
throw OpalException("DumpEMFields::writeFieldThis",
"The grid was NULL; there was a problem with the DumpEMFields initialisation.");
}
if (field == NULL) {
throw OpalException("DumpEMFields::writeFieldThis",
"The field to be written was NULL.");
}
double time = 0.;
std::vector<double> point_std(4);
Vector_t point(0., 0., 0.);
Vector_t centroid(0., 0., 0.);
std::ofstream fout;
try {
fout.open(filename_m.c_str(), std::ofstream::out);
} catch (std::exception& exc) {
throw OpalException("DumpEMFields::writeFieldThis",
"Failed to open DumpEMFields file "+filename_m);
}
if (!fout.good()) {
throw OpalException("DumpEMFields::writeFieldThis",
"Failed to open DumpEMFields file "+filename_m);
}
// set precision
fout << grid_m->end().toInteger() << "\n";
fout << 1 << " x [mm]\n";
fout << 2 << " y [mm]\n";
fout << 3 << " z [mm]\n";
fout << 4 << " t [ns]\n";
fout << 5 << " Bx [kGauss]\n";
fout << 6 << " By [kGauss]\n";
fout << 7 << " Bz [kGauss]\n";
fout << 8 << " Ex [?]\n";
fout << 9 << " Ey [?]\n";
fout << 10 << " Ez [?]\n";
fout << 0 << std::endl;
for (interpolation::Mesh::Iterator it = grid_m->begin();
it < grid_m->end();
++it) {
Vector_t E(0., 0., 0.);
Vector_t B(0., 0., 0.);
it.getPosition(&point_std[0]);
for (size_t i = 0; i < 3; ++i) {
point[i] = point_std[i];
}
time = point_std[3];
field->apply(point, centroid, time, E, B);
fout << point[0] << " " << point[1] << " " << point[2] << " " << time << " ";
fout << B[0] << " " << B[1] << " " << B[2] << " ";
fout << E[0] << " " << E[1] << " " << E[2] << "\n";
}
if (!fout.good()) {
throw OpalException("DumpEMFields::writeFieldThis",
"Something went wrong during writing "+filename_m);
}
fout.close();
}
/*
* Copyright (c) 2017, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef OPAL_BASICACTIONS_DUMPEMFIELDS_HH
#define OPAL_BASICACTIONS_DUMPEMFIELDS_HH
#include <unordered_set>
#include <string>
#include "AbstractObjects/Action.h"
namespace interpolation {
class NDGrid;
}
class Component;
/** DumpEMFields dumps the dynamically changing fields of a Ring in a user-
* defined grid.
*
* The idea is to print out the field map across a 4D grid in space-time for
* debugging purposes. The problem is to manage the DumpEMFields object through
* three phases of program execution; initial construction, parsing and then
* actual field map writing (where we need to somehow let DumpFields know what
* the field maps are). So for each DumpFields object created, we store in a
* set. When the execute() method is called, DumpFields builds a grid using
* the parsed information.
*
* When the ParallelCyclotronTracker is about to start tracking, it calls
* writeFields method which loops over the static set of DumpFields and writes
* each one. It is not the cleanest implementation, but I can't see a better
* way.
*
* The DumpEMFields themselves operate by iterating over a NDGrid object
* and looking up the field/writing it out on each grid point.
*
*/
class DumpEMFields : public Action {
public:
/** Constructor */
DumpEMFields();
/** Destructor deletes grid_m and if in the dumps set, take it out */
virtual ~DumpEMFields();
/** Make a clone (overloadable copy-constructor).
* @param name not used
* If this is in the dumpsSet_m, so will the clone. Not sure how the
* itsAttr stuff works, so this may not get properly copied?
*/
virtual DumpEMFields *clone(const std::string &name);
/** Builds the grid but does not write the field map
*
* Builds a grid of points in x-y-z space using the NDGrid algorithm.
* Checks that X_STEPS, Y_STEPS, Z_STEPS are integers or throws
* OpalException.
*/
virtual void execute();
/** Write the fields for all defined DumpEMFields objects
* @param field borrowed reference to the Component object that holds the
* field map; caller owns the memory.
* Iterates over the DumpEMFields in the dumpsSet_m and calls writeFieldThis
* on each DumpEMFields. This writes each field map in turn. Format is:
* <number of rows>
* <column 1> <units>
* <column 2> <units>
* <column 3> <units>
* <column 4> <units>
* <column 5> <units>
* <column 6> <units>
* 0
* <field map data>
*/
static void writeFields(Component* field);
private:
virtual void writeFieldThis(Component* field);
virtual void buildGrid();
static void checkInt(double value, std::string name, double tolerance = 1e-9);
interpolation::NDGrid* grid_m;
std::string filename_m;
static std::unordered_set<DumpEMFields*> dumpsSet_m;
static std::string dumpemfields_docstring;
DumpEMFields(const DumpEMFields& dump); // disabled
DumpEMFields& operator=(const DumpEMFields& dump); // disabled
};
#endif // ifdef OPAL_DUMPFIELDS_HH
......@@ -38,9 +38,9 @@ namespace interpolation {
class Component;
/** DumpFields dumps the fields of a Ring in a user-defined grid
/** DumpFields dumps the static magnetic field of a Ring in a user-defined grid
*
* The idea is to print out the field map across a grid in space-time for
* The idea is to print out the field map across a grid in space for
* debugging purposes. The problem is to manage the DumpFields object through
* three phases of program execution; initial construction, parsing and then
* actual field map writing (where we need to somehow let DumpFields know what
......@@ -56,6 +56,8 @@ class Component;
* The DumpFields themselves operate by iterating over a ThreeDGrid object
* and looking up the field/writing it out on each grid point.
*
* In order to dump time dependent fields, for example RF, see the DumpEMFields
* action.
*/
class DumpFields : public Action {
public:
......
......@@ -17,7 +17,7 @@
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
//c
// ------------------------------------------------------------------------
#include "AbsBeamline/Component.h"
......
......@@ -105,8 +105,10 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
<< " out of the field map boundary" << endl;
lossDS_m->addParticle(refPartBunch_m->R[id], refPartBunch_m->P[id], id);
lossDS_m->save();
refPartBunch_m->Bin[id] = -1;
}
return flagNeedUpdate;
}
......
......@@ -33,6 +33,7 @@ SBend3D::SBend3D(const std::string &name)
: Component(name), map_m(NULL),
planarArcGeometry_m(1., 1.), fieldUnits_m(1.), lengthUnits_m(1.),
polyOrder_m(1), smoothOrder_m(1), dummy() {
setElType(isDrift);
}
SBend3D::SBend3D(const SBend3D &right)
......@@ -44,6 +45,7 @@ SBend3D::SBend3D(const SBend3D &right)
RefPartBunch_m = right.RefPartBunch_m;
if (right.map_m != NULL)
map_m = new SectorMagneticFieldMap(*right.map_m);
setElType(isDrift);
}
SBend3D::~SBend3D() {
......@@ -117,4 +119,4 @@ void SBend3D::setFieldMapFileName(std::string name) {
void SBend3D::accept(BeamlineVisitor& visitor) const {
visitor.visitSBend3D(*this);
}
\ No newline at end of file
}
......@@ -4,6 +4,7 @@ set (_SRCS
Mesh.cpp
MMatrix.cpp
MVector.cpp
NDGrid.cpp
PolynomialPatch.cpp
PPSolveFactory.cpp
SolveFactory.cpp
......@@ -25,6 +26,7 @@ set (HDRS
Mesh.h
MMatrix.h
MVector.h
NDGrid.cpp
PolynomialCoefficient.h
PolynomialPatch.h
PPSolveFactory.h
......@@ -35,4 +37,4 @@ set (HDRS
VectorMap.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Interpolation")
\ No newline at end of file
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Interpolation")
// MAUS WARNING: THIS IS LEGACY CODE.
#include <math.h>
#include <iomanip>
#include "Utilities/OpalException.h"
#include "Fields/Interpolation/Mesh.h"
#include "Fields/Interpolation/NDGrid.h"
namespace interpolation {
/// ////// NDGrid ///////
NDGrid::NDGrid() : coord_m(), maps_m(), constantSpacing_m(false) {
}
NDGrid::NDGrid(std::vector<int> size, std::vector<const double *> gridCoordinates)
: coord_m(), maps_m(), constantSpacing_m(false) {
for (unsigned int i=0; i<size.size(); i++) {
if (size[i] < 1) {
throw(OpalException("NDGrid::NDGrid(...)",
"ND Grid must be at least 1x1x...x1 grid"));
}
coord_m.push_back(std::vector<double>(gridCoordinates[i],
gridCoordinates[i] + size[i]));
}
setConstantSpacing();
}
NDGrid::NDGrid(int nDimensions, int* size, double* spacing, double* min)
: coord_m(nDimensions), maps_m(), constantSpacing_m(true) {
for (int i=0; i<nDimensions; i++) {
if (size[i] < 1) {
throw(OpalException("NDGrid::NDGrid(...)",
"ND Grid must be at least 1x1x...x1 grid"));
}
coord_m[i] = std::vector<double>(size[i]);
for (unsigned int j=0; j<coord_m[i].size(); j++) {
coord_m[i][j] = min[i] + j*spacing[i];
}
}
}
NDGrid::NDGrid(std::vector< std::vector<double> > gridCoordinates)
: coord_m(gridCoordinates), maps_m(), constantSpacing_m(false) {
for (unsigned int i=0; i<gridCoordinates.size(); i++) {
if (gridCoordinates[i].size() < 1) {
throw (OpalException("NDGrid::NDGrid(...)",
"ND Grid must be at least 1x1x...x1 grid"));
}
}
setConstantSpacing();
}
double* NDGrid::newCoordArray ( const int& dimension) const {
double * array = new double[coord_m[dimension].size() ];
for (unsigned int i=0; i<coord_m[dimension].size(); i++) {
array[i] = coord_m[dimension][i];
}
return array;
}
//Mesh::Iterator wraps around a std::vector<int>
//it[0] is least significant, it[max] is most signifcant
Mesh::Iterator& NDGrid::addEquals(Mesh::Iterator& lhs, int difference) const {
if (difference < 0) {
return subEquals(lhs, -difference);
}
std::vector<int> index(coord_m.size(),0);
std::vector<int> content(coord_m.size(),1);
for (int i = int(index.size()-2); i >= 0; i--) {
content[i] = content[i+1]*coord_m[i+1].size(); //content could be member variable
}
for (int i = 0; i < int(index.size()); i++) {
index[i] = difference/content[i];
difference -= index[i] * content[i];
}
for (unsigned int i=0; i<index.size(); i++) {
lhs.state_m[i] += index[i];
}
for (int i = int(index.size())-1; i > 0; i--) {
if (lhs.state_m[i] > int(coord_m[i].size())) {
lhs.state_m[i-1]++;
lhs.state_m[i] -= coord_m[i].size();
}
}
return lhs;
}
Mesh::Iterator& NDGrid::subEquals(Mesh::Iterator& lhs, int difference) const {
if (difference < 0) {
return addEquals(lhs, -difference);
}
std::vector<int> index(coord_m.size(),0);
std::vector<int> content(coord_m.size(),1);
for (int i = int(index.size()-2); i >= 0; i--) {
content[i] = content[i+1]*coord_m[i+1].size(); //content could be member variable
}
for (int i = 0; i < int(index.size()); i++) {
index[i] = difference/content[i];
difference -= index[i] * content[i];
}
for (unsigned int i = 0; i < index.size(); i++) {
lhs.state_m[i] -= index[i];
}
for (int i=int(index.size())-1; i>0; i--) {
if (lhs.state_m[i] < 1) {
lhs.state_m[i-1]--;
lhs.state_m[i] += coord_m[i].size();
}
}
return lhs;
}
Mesh::Iterator& NDGrid::addEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
return addEquals(lhs, rhs.toInteger());
}
Mesh::Iterator& NDGrid::subEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
return subEquals(lhs, rhs.toInteger());
}
Mesh::Iterator& NDGrid::addOne(Mesh::Iterator& lhs) const {
int i = coord_m.size()-1;