DumpEMFields.cpp

Go to the documentation of this file.

//
// 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"
#include <filesystem>
#include <fstream>
#include "AbstractObjects/OpalData.h"
#include "Attributes/Attributes.h"
#include "Physics/Units.h"
#include "Utilities/OpalException.h"
#include "Utilities/Util.h"
 
extern Inform* gmsg;
 
std::unordered_set<std::unique_ptr<DumpEMFields>> DumpEMFields::dumpsSet_m;
 
DumpEMFields::DumpEMFields()
    : 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.") {
    // would be nice if "steps" could be integer
    itsAttr[FILE_NAME] =
            Attributes::makeString("FILE_NAME", "Name of the file to which field data is dumped");
 
    itsAttr[COORDINATE_SYSTEM] = Attributes::makePredefinedString(
            "COORDINATE_SYSTEM", "Choose to use CARTESIAN or CYLINDRICAL coordinates",
            {"CARTESIAN", "CYLINDRICAL"}, "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 (vertical) [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");
 
    itsAttr[CYL_ORIGIN_X] = Attributes::makeReal(
            "CYL_ORIGIN_X", "(Cylindrical) The X coordinate of the origin [m]. Default=0", 0.0);
 
    itsAttr[CYL_ORIGIN_Y] = Attributes::makeReal(
            "CYL_ORIGIN_Y", "(Cylindrical) The Y coordinate of the origin [m]. Default=0", 0.0);
 
    itsAttr[CYL_ORIGIN_Z] = Attributes::makeReal(
            "CYL_ORIGIN_Z", "(Cylindrical) The Z coordinate of the origin [m]. Default=0", 0.0);
 
    registerOwnership(AttributeHandler::STATEMENT);
}
 
DumpEMFields::DumpEMFields(const std::string& name, DumpEMFields* parent) : Action(name, parent) {}
 
DumpEMFields* DumpEMFields::clone(const std::string& name) {
    auto* dumper = new DumpEMFields(name, this);
    if (grid_m != nullptr) {
        dumper->grid_m.reset(grid_m->clone());
    }
    dumper->filename_m          = filename_m;
    dumper->coordinates_m       = coordinates_m;
    dumper->cylindricalOrigin_m = cylindricalOrigin_m;
    return dumper;
}
 
void DumpEMFields::parseCoordinateSystem() {
    static const std::map<std::string, CoordinateSystem> stringCoordinateSystem_s = {
            {"CARTESIAN", CoordinateSystem::CARTESIAN},
            {"CYLINDRICAL", CoordinateSystem::CYLINDRICAL}};
    coordinates_m = stringCoordinateSystem_s.at(Attributes::getString(itsAttr[COORDINATE_SYSTEM]));
}
 
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(std::unique_ptr<DumpEMFields>(this->clone("")));
}
 
void DumpEMFields::buildGrid() {
    std::vector<double> spacing(4);
    std::vector<int> gridSize(4);
    std::vector<double> origin(4);
    parseCoordinateSystem();
    if (coordinates_m == CoordinateSystem::CARTESIAN) {
        origin[0]       = Attributes::getReal(itsAttr[X_START]);
        spacing[0]      = Attributes::getReal(itsAttr[DX]);
        const double nx = Attributes::getReal(itsAttr[X_STEPS]);
        checkInt(nx, "X_STEPS");
        gridSize[0] = nx;
 
        origin[1]       = Attributes::getReal(itsAttr[Y_START]);
        spacing[1]      = Attributes::getReal(itsAttr[DY]);
        const double ny = Attributes::getReal(itsAttr[Y_STEPS]);
        checkInt(ny, "Y_STEPS");
        gridSize[1] = ny;
 
        origin[2]       = Attributes::getReal(itsAttr[Z_START]);
        spacing[2]      = Attributes::getReal(itsAttr[DZ]);
        const double nz = Attributes::getReal(itsAttr[Z_STEPS]);
        checkInt(nz, "Z_STEPS");
        gridSize[2] = nz;
    } else /* CoordinateSystem::CYLINDRICAL */ {
        origin[0]           = Attributes::getReal(itsAttr[R_START]);
        spacing[0]          = Attributes::getReal(itsAttr[DR]);
        const double rSteps = Attributes::getReal(itsAttr[R_STEPS]);
        checkInt(rSteps, "R_STEPS");
        gridSize[0] = rSteps;
 
        origin[1]         = Attributes::getReal(itsAttr[PHI_START]);
        spacing[1]        = Attributes::getReal(itsAttr[DPHI]);
        const double nphi = Attributes::getReal(itsAttr[PHI_STEPS]);
        checkInt(nphi, "PHI_STEPS");
        gridSize[1] = nphi;
 
        origin[2]       = Attributes::getReal(itsAttr[Y_START]);
        spacing[2]      = Attributes::getReal(itsAttr[DY]);
        const double ny = Attributes::getReal(itsAttr[Y_STEPS]);
        checkInt(ny, "Y_STEPS");
        gridSize[2] = ny;
    }
 
    origin[3]       = Attributes::getReal(itsAttr[T_START]);
    spacing[3]      = Attributes::getReal(itsAttr[DT]);
    const double nt = Attributes::getReal(itsAttr[T_STEPS]);
    checkInt(nt, "T_STEPS");
    gridSize[3] = nt;
 
    grid_m = std::make_unique<interpolation::NDGrid>(4, &gridSize[0], &spacing[0], &origin[0]);
 
    filename_m             = Attributes::getString(itsAttr[FILE_NAME]);
    cylindricalOrigin_m[0] = Attributes::getReal(itsAttr[CYL_ORIGIN_X]);
    cylindricalOrigin_m[1] = Attributes::getReal(itsAttr[CYL_ORIGIN_Y]);
    cylindricalOrigin_m[2] = Attributes::getReal(itsAttr[CYL_ORIGIN_Z]);
}
 
void DumpEMFields::writeFields(const std::set<std::shared_ptr<Component>>& elements) {
    for (auto& item : dumpsSet_m) {
        item->writeFieldThis(elements);
    }
    dumpsSet_m.clear();
}
 
void DumpEMFields::checkInt(double value, const std::string& name, const double tolerance) {
    value += tolerance;  // prevent rounding error
    if (std::abs(std::floor(value) - value) > 2 * tolerance) {
        throw OpalException(
                "DumpEMFields::checkInt",
                "Value for " + name + " should be an integer but a real value was found");
    }
    if (std::floor(value) < 0.5) {
        throw OpalException("DumpEMFields::checkInt", "Value for " + name + " should be 1 or more");
    }
}
 
void DumpEMFields::writeHeader(std::ofstream& fout) const {
    fout << grid_m->end().toInteger() << "\n";
    if (coordinates_m == CoordinateSystem::CARTESIAN) {
        fout << 1 << "  x [m]\n";
        fout << 2 << "  y [m]\n";
        fout << 3 << "  z [m]\n";
        fout << 4 << "  t [ns]\n";
        fout << 5 << "  Bx [T]\n";
        fout << 6 << "  By [T]\n";
        fout << 7 << "  Bz [T]\n";
        fout << 8 << "  Ex [MV/m]\n";
        fout << 9 << "  Ey [MV/m]\n";
        fout << 10 << " Ez [MV/m]\n";
    } else /*CoordinateSystem::CYLINDRICAL*/ {
        fout << 1 << "  r [m]\n";
        fout << 2 << "  phi [deg]\n";
        fout << 3 << "  y [m]\n";
        fout << 4 << "  t [ns]\n";
        fout << 5 << "  Br   [T]\n";
        fout << 6 << "  Bphi [T]\n";
        fout << 7 << "  By   [T]\n";
        fout << 8 << "  Er   [MV/m]\n";
        fout << 9 << "  Ephi [MV/m]\n";
        fout << 10 << " Ey   [MV/m]\n";
    }
    fout << 0 << std::endl;
}
 
void DumpEMFields::writeFieldLine(
        const std::set<std::shared_ptr<Component>>& elements, const Vector_t<double, 3>& point,
        const double& time, std::ofstream& fout) const {
    // Translate the coordinate if necessary
    Vector_t<double, 3> E{};
    Vector_t<double, 3> B{};
    Vector_t<double, 3> localPoint = point;
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
        // point is (r, phi, y), but result and origin is (x, y, z) - y is vertical
        localPoint[0] = std::cos(point[1]) * point[0] + cylindricalOrigin_m[0];
        localPoint[2] = std::sin(point[1]) * point[0] + cylindricalOrigin_m[2];
        localPoint[1] = point[2] + cylindricalOrigin_m[1];
    }
    // Collect the fields
    for (auto& element : elements) {
        auto transform             = element->getCSTrafoGlobal2Local();
        Vector_t<double, 3> localR = transform.transformTo(localPoint);
        Vector_t<double, 3> localB{};
        Vector_t<double, 3> localE{};
        element->apply(localR, {}, time, localE, localB);
        transform.invert();
        localB = transform.rotateTo(localB);
        localE = transform.rotateTo(localE);
        E += localE;
        B += localB;
    }
    // Translate the results and output them
    Vector_t<double, 3> Bout = B;
    Vector_t<double, 3> Eout = E;
    if (coordinates_m == CoordinateSystem::CYLINDRICAL) {
        // point and field out is (r, phi, y), but field in is (x, y, z) - y is vertical
        Bout[0] = B[0] * std::cos(point[1]) + B[2] * std::sin(point[1]);
        Bout[1] = -B[0] * std::sin(point[1]) + B[2] * std::cos(point[1]);
        Bout[2] = B[1];
        Eout[0] = E[0] * std::cos(point[1]) + E[2] * std::sin(point[1]);
        Eout[1] = -E[0] * std::sin(point[1]) + E[2] * std::cos(point[1]);
        Eout[2] = E[1];
        fout << point[0] << " " << point[1] * Units::rad2deg << " " << point[2] << " " << time
             << " ";
    } else {
        fout << point[0] << " " << point[1] << " " << point[2] << " " << time << " ";
    }
    fout << Bout[0] << " " << Bout[1] << " " << Bout[2] << " ";
    fout << Eout[0] * Units::Vpm2MVpm << " " << Eout[1] * Units::Vpm2MVpm << " "
         << Eout[2] * Units::Vpm2MVpm << "\n";
}
 
void DumpEMFields::writeFieldThis(const std::set<std::shared_ptr<Component>>& elements) {
    if (grid_m == nullptr) {
        throw OpalException(
                "DumpEMFields::writeFieldThis",
                "The grid was nullptr; there was a problem with the "
                "DumpEMFields initialisation.");
    }
 
    *gmsg << level5 << *this << endl;
 
    std::string fname;
    if (std::filesystem::path(filename_m).is_absolute() == true) {
        fname = filename_m;
    } else {
        fname = Util::combineFilePath(
                {OpalData::getInstance()->getAuxiliaryOutputDirectory(), filename_m});
    }
 
    std::vector<double> point_std(4);
    Vector_t<double, 3> point(0., 0., 0.);
    std::ofstream fout;
    fout.open(fname.c_str(), std::ofstream::out);
    if (!fout.good()) {
        throw OpalException(
                "DumpEMFields::writeFieldThis", "Failed to open DumpEMFields file " + filename_m);
    }
    // set precision
    writeHeader(fout);
    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];
        }
        double time = point_std[3];
        writeFieldLine(elements, point, time, fout);
    }
    if (failWrite_m) {  // For unit tests only - force a write failure if required
        fout.setstate(std::ios::badbit);
    }
    if (!fout.good()) {
        throw OpalException(
                "DumpEMFields::writeFieldThis",
                "Something went wrong during writing " + filename_m);
    }
    fout.close();
}
 
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: " << Attributes::getString(itsAttr[COORDINATE_SYSTEM]) << '\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'
           << "* Z_START   = " << Attributes::getReal(itsAttr[Z_START]) << " [m]\n"
           << "* DZ        = " << Attributes::getReal(itsAttr[DZ]) << " [m]\n"
           << "* Z_STEPS   = " << Attributes::getReal(itsAttr[Z_STEPS]) << '\n';
    } else /*if (coordinates_m == CoordinateSystem::CYLINDRICAL)*/ {
        os << "* Coordinate system: " << Attributes::getString(itsAttr[COORDINATE_SYSTEM]) << '\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'
           << "* Y_START   = " << Attributes::getReal(itsAttr[Y_START]) << " [m]\n"
           << "* DY        = " << Attributes::getReal(itsAttr[DY]) << " [m]\n"
           << "* Y_STEPS   = " << Attributes::getReal(itsAttr[Y_STEPS]) << '\n';
    }
    os << "* 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;
}