DumpEMFields.cpp 13.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
//
// Class DumpEMFields
//   DumpEMFields dumps the dynamically changing fields of a Ring in a user-
//   defined grid.
//
// Copyright (c) 2017, Chris Rogers
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "BasicActions/DumpEMFields.h"
20

21
#include "AbstractObjects/OpalData.h"
22
#include "AbsBeamline/Component.h"
23 24 25 26 27
#include "Attributes/Attributes.h"
#include "Fields/Interpolation/NDGrid.h"
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
#include "Utilities/Util.h"
28

29 30
#include <boost/filesystem.hpp>

31
#include <fstream>
32
#include <cmath>
33

34
extern Inform* gmsg;
35 36 37 38

std::unordered_set<DumpEMFields*> DumpEMFields::dumpsSet_m;

DumpEMFields::DumpEMFields() :
39 40 41 42 43
    Action(SIZE, "DUMPEMFIELDS",
           "The \"DUMPEMFIELDS\" statement dumps a field map to a user-defined "
           "field file, for checking that fields are generated correctly. "
           "The fields are written out on a grid in space and time."),
    grid_m(NULL),
44
    filename_m("") {
45

46
    // would be nice if "steps" could be integer
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
    itsAttr[FILE_NAME] = Attributes::makeString
        ("FILE_NAME", "Name of the file to which field data is dumped");

    itsAttr[COORDINATE_SYSTEM] = Attributes::makeUpperCaseString
        ("COORDINATE_SYSTEM", "Choose to use CARTESIAN (default) or CYLINDRICAL coordinates", "CARTESIAN");

    itsAttr[X_START] = Attributes::makeReal
        ("X_START", "(Cartesian) Start point in the grid in x [m]");

    itsAttr[DX] = Attributes::makeReal
        ("DX", "(Cartesian) Grid step size in x [m]");

    itsAttr[X_STEPS] = Attributes::makeReal
        ("X_STEPS", "(Cartesian) Number of steps in x");

    itsAttr[Y_START] = Attributes::makeReal
        ("Y_START", "(Cartesian) Start point in the grid in y [m]");

    itsAttr[DY] = Attributes::makeReal
        ("DY", "(Cartesian) Grid step size in y [m]");

    itsAttr[Y_STEPS] = Attributes::makeReal
        ("Y_STEPS", "(Cartesian) Number of steps in y");

    itsAttr[Z_START] = Attributes::makeReal
        ("Z_START", "Start point in the grid in z [m]");

    itsAttr[DZ] = Attributes::makeReal
        ("DZ", "Grid step size in z [m]");

    itsAttr[Z_STEPS] = Attributes::makeReal
        ("Z_STEPS", "Number of steps in z");

    itsAttr[T_START] = Attributes::makeReal
        ("T_START", "Start point in the grid in time [ns]");

    itsAttr[DT] = Attributes::makeReal
        ("DT", "Grid step size in time [ns]");

    itsAttr[T_STEPS] = Attributes::makeReal
        ("T_STEPS", "Number of steps in time");

    itsAttr[R_START] = Attributes::makeReal
        ("R_START", "(Cylindrical) Start point in the grid in radius [m]");

    itsAttr[DR] = Attributes::makeReal
        ("DR", "(Cylindrical) Grid step size in radius [m]");

    itsAttr[R_STEPS] = Attributes::makeReal
        ("R_STEPS", "(Cylindrical) Number of steps in radius");

    itsAttr[PHI_START] = Attributes::makeReal
        ("PHI_START", "(Cylindrical) Start point in the grid in phi [rad]");

    itsAttr[DPHI] = Attributes::makeReal
        ("DPHI", "(Cylindrical) Grid step size in phi [rad]");

    itsAttr[PHI_STEPS] = Attributes::makeReal
        ("PHI_STEPS", "(Cylindrical) Number of steps in phi");

    registerOwnership(AttributeHandler::STATEMENT);
108 109
}

110
DumpEMFields::DumpEMFields(const std::string& name, DumpEMFields* parent):
111
    Action(name, parent), grid_m(NULL)
112 113
{}

114
DumpEMFields::~DumpEMFields() {
snuverink_j's avatar
snuverink_j committed
115 116
    delete grid_m;
    dumpsSet_m.erase(this);
117 118
}

119 120
DumpEMFields* DumpEMFields::clone(const std::string& name) {
    DumpEMFields* dumper = new DumpEMFields(name, this);
121 122 123 124
    if (grid_m != NULL) {
        dumper->grid_m = grid_m->clone();
    }
    dumper->filename_m = filename_m;
125
    dumper->coordinates_m = coordinates_m;
126 127 128 129 130 131
    if (dumpsSet_m.find(this) != dumpsSet_m.end()) {
        dumpsSet_m.insert(dumper);
    }
    return dumper;
}

132
void DumpEMFields::parseCoordinateSystem() {
133 134 135 136 137 138

    std::string coordStr = Attributes::getString(itsAttr[COORDINATE_SYSTEM]);
    if (coordStr == "CYLINDRICAL") {
        coordinates_m = CoordinateSystem::CYLINDRICAL;
    } else if (coordStr == "CARTESIAN") {
        coordinates_m = CoordinateSystem::CARTESIAN;
139 140
    } else {
        throw OpalException("DumpEMFields::parseCoordinateSystem",
141 142
                            "Failed to parse '" + coordStr +
                            "'. OPAL expected either CYLINDRICAL or CARTESIAN.");
143 144 145
    }
}

146 147 148 149 150 151 152 153 154 155 156 157 158
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); 
159
    parseCoordinateSystem();
160

161 162 163 164
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
        origin[0] = Attributes::getReal(itsAttr[R_START]);
        spacing[0] = Attributes::getReal(itsAttr[DR]);
        double nr = Attributes::getReal(itsAttr[R_STEPS]);
165 166
        checkInt(nr, "R_STEPS");
        gridSize[0] = nr;
167

168 169 170
        origin[1] = Attributes::getReal(itsAttr[PHI_START]);
        spacing[1] = Attributes::getReal(itsAttr[DPHI]);
        double nphi = Attributes::getReal(itsAttr[PHI_STEPS]);
171 172 173
        checkInt(nphi, "PHI_STEPS");
        gridSize[1] = nphi;
    } else {
174 175 176
        origin[0] = Attributes::getReal(itsAttr[X_START]);
        spacing[0] = Attributes::getReal(itsAttr[DX]);
        double nx = Attributes::getReal(itsAttr[X_STEPS]);
177 178 179
        checkInt(nx, "X_STEPS");
        gridSize[0] = nx;

180 181 182
        origin[1] = Attributes::getReal(itsAttr[Y_START]);
        spacing[1] = Attributes::getReal(itsAttr[DY]);
        double ny = Attributes::getReal(itsAttr[Y_STEPS]);
183 184 185
        checkInt(ny, "Y_STEPS");
        gridSize[1] = ny;
    }
186 187 188
    origin[2] = Attributes::getReal(itsAttr[Z_START]);
    spacing[2] = Attributes::getReal(itsAttr[DZ]);
    double nz = Attributes::getReal(itsAttr[Z_STEPS]);
189 190
    checkInt(nz, "Z_STEPS");
    gridSize[2] = nz;
191

192 193 194
    origin[3] = Attributes::getReal(itsAttr[T_START]);
    spacing[3] = Attributes::getReal(itsAttr[DT]);
    double nt = Attributes::getReal(itsAttr[T_STEPS]);
195 196 197 198 199 200 201 202 203
    checkInt(nt, "T_STEPS");
    gridSize[3] = nt;

    if (grid_m != NULL) {
        delete grid_m;
    }

    grid_m = new interpolation::NDGrid(4, &gridSize[0], &spacing[0], &origin[0]);

204
    filename_m = Attributes::getString(itsAttr[FILE_NAME]);
205 206 207 208 209 210 211 212 213 214
}

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) {
215
    real += tolerance; // prevent rounding error
216
    if (std::abs(std::floor(real) - real) > 2*tolerance) {
217
        throw OpalException("DumpEMFields::checkInt",
218
                            "Value for " + name +
219 220
                            " should be an integer but a real value was found");
    }
221
    if (std::floor(real) < 0.5) {
222
        throw OpalException("DumpEMFields::checkInt",
223
                            "Value for " + name + " should be 1 or more");
224 225 226
    }
}

227 228
void DumpEMFields::writeHeader(std::ofstream& fout) const {
    fout << grid_m->end().toInteger() << "\n";
229
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
230 231
        fout << 1 << "  r [mm]\n";
        fout << 2 << "  phi [degree]\n";
232
    } else {
233 234
        fout << 1 << "  x [mm]\n";
        fout << 2 << "  y [mm]\n";
235
    }
236 237
    fout << 3 << "  z [mm]\n";
    fout << 4 << "  t [ns]\n";
238
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
239 240 241 242 243 244
        fout << 5 << "  Br [kGauss]\n";
        fout << 6 << "  Bphi [kGauss]\n";
        fout << 7 << "  Bz [kGauss]\n";
        fout << 8 << "  Er [MV/m]\n";
        fout << 9 << "  Ephi [MV/m]\n";
        fout << 10 << " Ez [MV/m]\n";
245
    } else {
246 247 248 249 250 251
        fout << 5 << "  Bx [kGauss]\n";
        fout << 6 << "  By [kGauss]\n";
        fout << 7 << "  Bz [kGauss]\n";
        fout << 8 << "  Ex [MV/m]\n";
        fout << 9 << "  Ey [MV/m]\n";
        fout << 10 << " Ez [MV/m]\n";
252 253 254 255 256 257 258 259 260 261 262 263
    }
    fout << 0 << std::endl;
}

void DumpEMFields::writeFieldLine(Component* field,
                                  const Vector_t& pointIn,
                                  const double& time,
                                  std::ofstream& fout) const {
    Vector_t centroid(0., 0., 0.);
    Vector_t E(0., 0., 0.);
    Vector_t B(0., 0., 0.);
    Vector_t point = pointIn;
264
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
265
        // pointIn is r, phi, z 
266 267
        point[0] = std::cos(pointIn[1])*pointIn[0];
        point[1] = std::sin(pointIn[1])*pointIn[0];
268 269 270 271 272
    }

    field->apply(point, centroid, time, E, B);
    Vector_t Bout = B;
    Vector_t Eout = E;
273
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
274
        // pointIn is r, phi, z 
275
        Bout[0] =  B[0]*std::cos(pointIn[1])+B[1]*std::sin(pointIn[1]);
276
        Bout[1] = -B[0]*std::sin(pointIn[1])+B[1]*std::cos(pointIn[1]);
277
        Eout[0] =  E[0]*std::cos(pointIn[1])+E[1]*std::sin(pointIn[1]);
278
        Eout[1] = -E[0]*std::sin(pointIn[1])+E[1]*std::cos(pointIn[1]);
279
        fout << pointIn[0] << " " << pointIn[1]*Physics::rad2deg << " " << pointIn[2] << " " << time << " ";
280 281 282 283 284 285 286 287
    } else {
        fout << pointIn[0] << " " << pointIn[1] << " " << pointIn[2] << " " << time << " ";
    }

    fout << Bout[0] << " " << Bout[1] << " " << Bout[2] << " ";
    fout << Eout[0] << " " << Eout[1] << " " << Eout[2] << "\n";
}

288 289 290 291 292 293 294 295 296
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.");
    }
297 298 299

    *gmsg << *this << endl;

300 301 302 303 304 305 306 307 308 309
    std::string fname;
    if (boost::filesystem::path(filename_m).is_absolute() == true) {
        fname = filename_m;
    } else {
        fname = Util::combineFilePath({
            OpalData::getInstance()->getAuxiliaryOutputDirectory(),
            filename_m
        });
    }

310 311 312 313
    std::vector<double> point_std(4);
    Vector_t point(0., 0., 0.);
    std::ofstream fout;
    try {
314
        fout.open(fname.c_str(), std::ofstream::out);
315 316
    } catch (std::exception& exc) {
        throw OpalException("DumpEMFields::writeFieldThis",
317
                            "Failed to open DumpEMFields file " + filename_m);
318 319 320
    }
    if (!fout.good()) {
        throw OpalException("DumpEMFields::writeFieldThis",
321
                            "Failed to open DumpEMFields file " + filename_m);
322 323
    }
    // set precision
324
    writeHeader(fout);
325 326 327 328 329 330 331
    for (interpolation::Mesh::Iterator it = grid_m->begin();
         it < grid_m->end();
         ++it) {
        it.getPosition(&point_std[0]);
        for (size_t i = 0; i < 3; ++i) {
            point[i] = point_std[i];
        }
snuverink_j's avatar
snuverink_j committed
332
        double time = point_std[3];
333
        writeFieldLine(field, point, time, fout);
334 335 336
    }
    if (!fout.good()) {
        throw OpalException("DumpEMFields::writeFieldThis",
337
                            "Something went wrong during writing " + filename_m);
338 339
    }
    fout.close();
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
}

void DumpEMFields::print(std::ostream& os) const {
    os << "* ************* D U M P  E M  F I E L D S ****************************************** " << std::endl;
    os << "* File name: " << filename_m << '\n';
    if (coordinates_m == CoordinateSystem::CARTESIAN) {
        os << "* Coordinate system: " << "Cartesian" << '\n'
           << "* X_START   = " << Attributes::getReal(itsAttr[X_START]) << " [m]\n"
           << "* DX        = " << Attributes::getReal(itsAttr[DX])      << " [m]\n"
           << "* X_STEPS   = " << Attributes::getReal(itsAttr[X_STEPS]) << '\n'
           << "* Y_START   = " << Attributes::getReal(itsAttr[Y_START]) << " [m]\n"
           << "* DY        = " << Attributes::getReal(itsAttr[DY])      << " [m]\n"
           << "* Y_STEPS   = " << Attributes::getReal(itsAttr[Y_STEPS]) << '\n';
    } else if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
        os << "* Coordinate system: " << "Cylindrical" << '\n'
           << "* R_START   = " << Attributes::getReal(itsAttr[R_START])   << " [m]\n"
           << "* DR        = " << Attributes::getReal(itsAttr[DR])        << " [m]\n"
           << "* R_STEPS   = " << Attributes::getReal(itsAttr[R_STEPS])   << '\n'
           << "* PHI_START = " << Attributes::getReal(itsAttr[PHI_START]) << " [rad]\n"
           << "* DPHI      = " << Attributes::getReal(itsAttr[DPHI])      << " [rad]\n"
           << "* PHI_STEPS = " << Attributes::getReal(itsAttr[PHI_STEPS]) << '\n';
    }
    os << "* Z_START   = " << Attributes::getReal(itsAttr[Z_START]) << " [m]\n"
       << "* DZ        = " << Attributes::getReal(itsAttr[DZ])      << " [m]\n"
       << "* Z_STEPS   = " << Attributes::getReal(itsAttr[Z_STEPS]) << '\n'
       << "* T_START   = " << Attributes::getReal(itsAttr[T_START]) << " [ns]\n"
       << "* DT        = " << Attributes::getReal(itsAttr[DT])      << " [ns]\n"
       << "* T_STEPS   = " << Attributes::getReal(itsAttr[T_STEPS]) << '\n';
    os << "* ********************************************************************************** " << std::endl;
}