H5PartWrapperForPT.cpp

Go to the documentation of this file.

//
// Class H5PartWrapperForPT
//   A class that manages all calls to H5Part for the Parallel-T tracker.
//
// Copyright (c) 200x-2021, Paul Scherrer Institut, Villigen PSI, Switzerland
// 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 "Structure/H5PartWrapperForPT.h"
 
#include "AbstractObjects/OpalData.h"
#include "PartBunch/PartBunch.h"
 
#include "OPALconfig.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
#include "Utilities/OpalException.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "h5core/h5_types.h"
 
#include <cmath>
#include <set>
#include <sstream>
 
H5PartWrapperForPT::H5PartWrapperForPT(const std::string& fileName, h5_int32_t flags)
    : H5PartWrapper(fileName, flags) {
}
 
H5PartWrapperForPT::H5PartWrapperForPT(
    const std::string& fileName, int restartStep, std::string sourceFile, h5_int32_t flags)
    : H5PartWrapper(fileName, restartStep, sourceFile, flags) {
    if (restartStep == -1) {
        restartStep = H5GetNumSteps(file_m) - 1;
        OpalData::getInstance()->setRestartStep(restartStep);
    }
}
 
H5PartWrapperForPT::~H5PartWrapperForPT() {
}
 
void H5PartWrapperForPT::readHeader() {
    h5_int64_t numFileAttributes = H5GetNumFileAttribs(file_m);
 
    const h5_size_t lengthAttributeName = 256;
    char attributeName[lengthAttributeName];
    h5_int64_t attributeType;
    h5_size_t numAttributeElements;
    std::set<std::string> attributeNames;
 
    for (h5_int64_t i = 0; i < numFileAttributes; ++i) {
        REPORTONERROR(H5GetFileAttribInfo(
            file_m, i, attributeName, lengthAttributeName, &attributeType, &numAttributeElements));
 
        attributeNames.insert(attributeName);
    }
 
    if (attributeNames.find("dump frequency") != attributeNames.end()) {
        h5_int64_t dumpfreq;
        READFILEATTRIB(Int64, file_m, "dump frequency", &dumpfreq);
        OpalData::getInstance()->setRestartDumpFreq(dumpfreq);
    }
 
    if (attributeNames.find("dPhiGlobal") != attributeNames.end()) {
        h5_float64_t dphi;
        READFILEATTRIB(Float64, file_m, "dPhiGlobal", &dphi);
        OpalData::getInstance()->setGlobalPhaseShift(dphi);
    }
 
    if (attributeNames.find("nAutoPhaseCavities") != attributeNames.end()) {
        auto opal = OpalData::getInstance();
        h5_float64_t phase;
        h5_int64_t numAutoPhaseCavities = 0;
        if (!H5HasFileAttrib(file_m, "nAutoPhaseCavities")
            || H5ReadFileAttribInt64(file_m, "nAutoPhaseCavities", &numAutoPhaseCavities)
                   != H5_SUCCESS) {
            numAutoPhaseCavities = 0;
        } else {
            for (long i = 0; i < numAutoPhaseCavities; ++i) {
                std::string elementName  = "Cav-" + std::to_string(i + 1) + "-name";
                std::string elementPhase = "Cav-" + std::to_string(i + 1) + "-value";
 
                char name[128];
                READFILEATTRIB(String, file_m, elementName.c_str(), name);
                READFILEATTRIB(Float64, file_m, elementPhase.c_str(), &phase);
 
                opal->setMaxPhase(name, phase);
            }
        }
    }
}
 
void H5PartWrapperForPT::readStep(
    PartBunch_t* bunch, h5_ssize_t firstParticle, h5_ssize_t lastParticle) {
    h5_ssize_t numStepsInSource = H5GetNumSteps(file_m);
    h5_ssize_t readStep         = numStepsInSource - 1;
    REPORTONERROR(H5SetStep(file_m, readStep));
 
    readStepHeader(bunch);
    readStepData(bunch, firstParticle, lastParticle);
}
 
void H5PartWrapperForPT::readStepHeader(PartBunch_t* bunch) {
    double actualT;
    READSTEPATTRIB(Float64, file_m, "TIME", &actualT);
    bunch->setT(actualT);
 
    double spos;
    READSTEPATTRIB(Float64, file_m, "SPOS", &spos);
    bunch->set_sPos(spos);
 
    h5_int64_t ltstep;
    READSTEPATTRIB(Int64, file_m, "LocalTrackStep", &ltstep);
    bunch->setLocalTrackStep((long long)(ltstep + 1));
 
    h5_int64_t gtstep;
    READSTEPATTRIB(Int64, file_m, "GlobalTrackStep", &gtstep);
    bunch->setGlobalTrackStep((long long)(gtstep + 1));
 
    Vector_t<double, 3> RefPartR;
    READSTEPATTRIB(Float64, file_m, "RefPartR", (h5_float64_t*)&RefPartR);
    bunch->RefPartR_m = RefPartR;
 
    Vector_t<double, 3> RefPartP;
    READSTEPATTRIB(Float64, file_m, "RefPartP", (h5_float64_t*)&RefPartP);
    bunch->RefPartP_m = RefPartP;
 
    Vector_t<double, 3> TaitBryant;
    READSTEPATTRIB(Float64, file_m, "TaitBryantAngles", (h5_float64_t*)&TaitBryant);
    Quaternion rotTheta(std::cos(0.5 * TaitBryant[0]), 0, std::sin(0.5 * TaitBryant[0]), 0);
    Quaternion rotPhi(std::cos(0.5 * TaitBryant[1]), std::sin(0.5 * TaitBryant[1]), 0, 0);
    Quaternion rotPsi(std::cos(0.5 * TaitBryant[2]), 0, 0, std::sin(0.5 * TaitBryant[2]));
    Quaternion rotation = rotTheta * (rotPhi * rotPsi);
    // ADA bunch->toLabTrafo_m = CoordinateSystemTrafo(-rotation.conjugate().rotate(RefPartR),
    // rotation);
}
 
void H5PartWrapperForPT::readStepData(
    PartBunch_t* bunch, h5_ssize_t firstParticle, h5_ssize_t lastParticle) {
    h5_ssize_t numParticles = getNumParticles();
    if (lastParticle >= numParticles || firstParticle > lastParticle) {
        throw OpalException(
            "H5PartWrapperForPT::readStepData",
            "the provided particle numbers don't match the number of particles in the file");
    }
 
    REPORTONERROR(H5PartSetView(file_m, firstParticle, lastParticle));
 
    numParticles = lastParticle - firstParticle + 1;
 
    std::vector<char> buffer(numParticles * sizeof(h5_float64_t));
    char* buffer_ptr        = Util::c_data(buffer);
    h5_float64_t* f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
    h5_int32_t* i32buffer   = reinterpret_cast<h5_int32_t*>(buffer_ptr);
 
    READDATA(Float64, file_m, "x", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // \todo bunch->R(n)(0) = f64buffer[n];
        // bunch->Bin(n)  = 0;
    }
 
    READDATA(Float64, file_m, "y", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->R(n)(1) = f64buffer[n];
    }
 
    READDATA(Float64, file_m, "z", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->R(n)(2) = f64buffer[n];
    }
 
    READDATA(Float64, file_m, "px", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->P(n)(0) = f64buffer[n];
    }
 
    READDATA(Float64, file_m, "py", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->P(n)(1) = f64buffer[n];
    }
 
    READDATA(Float64, file_m, "pz", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->P(n)(2) = f64buffer[n];
    }
 
    READDATA(Float64, file_m, "q", f64buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->Q(n) = f64buffer[n];
    }
 
    READDATA(Int32, file_m, "id", i32buffer);
    for (long int n = 0; n < numParticles; ++n) {
        // bunch->ID(n) = i32buffer[n];
    }
 
    REPORTONERROR(H5PartSetView(file_m, -1, -1));
}
 
void H5PartWrapperForPT::writeHeader() {
    std::stringstream OPAL_version;
    OPAL_version << OPAL_PROJECT_NAME << " " << OPAL_PROJECT_VERSION << " # git rev. "
                 << Util::getGitRevision();
    WRITESTRINGFILEATTRIB(file_m, "OPAL_version", OPAL_version.str().c_str());
 
    WRITESTRINGFILEATTRIB(file_m, "idUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "xUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "yUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "zUnit", "m");
 
    WRITESTRINGFILEATTRIB(file_m, "pxUnit", "#beta#gamma");
    WRITESTRINGFILEATTRIB(file_m, "pyUnit", "#beta#gamma");
    WRITESTRINGFILEATTRIB(file_m, "pzUnit", "#beta#gamma");
 
    WRITESTRINGFILEATTRIB(file_m, "ptypeUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "poriginUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "qUnit", "C");
 
    if (Options::ebDump) {
        WRITESTRINGFILEATTRIB(file_m, "ExUnit", "MV/m");
        WRITESTRINGFILEATTRIB(file_m, "EyUnit", "MV/m");
        WRITESTRINGFILEATTRIB(file_m, "EzUnit", "MV/m");
 
        WRITESTRINGFILEATTRIB(file_m, "BxUnit", "T");
        WRITESTRINGFILEATTRIB(file_m, "ByUnit", "T");
        WRITESTRINGFILEATTRIB(file_m, "BzUnit", "T");
    }
 
    if (Options::rhoDump) {
        WRITESTRINGFILEATTRIB(file_m, "rhoUnit", "C/m^3");
    }
 
    WRITESTRINGFILEATTRIB(file_m, "TaitBryantAnglesUnit", "rad");
 
    WRITESTRINGFILEATTRIB(file_m, "SPOSUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "TIMEUnit", "s");
    WRITESTRINGFILEATTRIB(file_m, "#gammaUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "ENERGYUnit", "MeV");
    WRITESTRINGFILEATTRIB(file_m, "#varepsilonUnit", "m rad");
    WRITESTRINGFILEATTRIB(file_m, "#varepsilon-geomUnit", "m rad");
 
    WRITESTRINGFILEATTRIB(file_m, "#sigmaUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "RMSXUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "centroidUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "RMSPUnit", "#beta#gamma");
    WRITESTRINGFILEATTRIB(file_m, "MEANPUnit", "#beta#gamma");
 
    WRITESTRINGFILEATTRIB(file_m, "minXUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "maxXUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "minPUnit", "#beta#gamma");
    WRITESTRINGFILEATTRIB(file_m, "maxPUnit", "#beta#gamma");
 
    WRITESTRINGFILEATTRIB(file_m, "dEUnit", "MeV");
 
    WRITESTRINGFILEATTRIB(file_m, "MASSUnit", "GeV");
    WRITESTRINGFILEATTRIB(file_m, "CHARGEUnit", "C");
 
    WRITESTRINGFILEATTRIB(file_m, "E-refUnit", "MV/m");
    WRITESTRINGFILEATTRIB(file_m, "B-refUnit", "T");
 
    WRITESTRINGFILEATTRIB(file_m, "StepUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "LocalTrackStepUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "GlobalTrackStepUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "NumBunchUnit", "1");
    WRITESTRINGFILEATTRIB(file_m, "NumPartUnit", "1");
 
    WRITESTRINGFILEATTRIB(file_m, "RefPartRUnit", "m");
    WRITESTRINGFILEATTRIB(file_m, "RefPartPUnit", "#beta#gamma");
    WRITESTRINGFILEATTRIB(file_m, "SteptoLastInjUnit", "1");
 
    WRITESTRINGFILEATTRIB(file_m, "dump frequencyUnit", "1")
    WRITESTRINGFILEATTRIB(file_m, "dPhiGlobalUnit", "rad")
 
    
    h5_int64_t dumpfreq = Options::psDumpFreq;
    WRITEFILEATTRIB(Int64, file_m, "dump frequency", &dumpfreq, 1);

    h5_float64_t dphi = OpalData::getInstance()->getGlobalPhaseShift();
    WRITEFILEATTRIB(Float64, file_m, "dPhiGlobal", &dphi, 1);
}
 
void H5PartWrapperForPT::writeStep(
    PartBunch_t* bunch, const std::map<std::string, double>& additionalStepAttributes) {
    
    if (bunch->getTotalNum() == 0)
        return;
 
    open(H5_O_APPENDONLY);
    bunch->calcBeamParameters();
    writeStepHeader(bunch, additionalStepAttributes);
    writeStepData(bunch);
    close();
}
 
void H5PartWrapperForPT::writeStepHeader(
    PartBunch_t* bunch, const std::map<std::string, double>& additionalStepAttributes) {
    double actPos                = bunch->get_sPos();
    double t                     = bunch->getT();
    Vector_t<double, 3> rmin     = bunch->get_origin();
    Vector_t<double, 3> rmax     = bunch->get_maxExtent();
    Vector_t<double, 3> centroid = bunch->get_centroid();
 
    Vector_t<double, 3> maxP(0.0);
    Vector_t<double, 3> minP(0.0);
 
    Vector_t<double, 3> xsigma     = bunch->get_rrms();
    Vector_t<double, 3> psigma     = bunch->get_prms();
    Vector_t<double, 3> vareps     = bunch->get_norm_emit();
    Vector_t<double, 3> geomvareps = bunch->get_emit();
    Vector_t<double, 3> RefPartR   = bunch->RefPartR_m;
    Vector_t<double, 3> RefPartP   = bunch->RefPartP_m;
    Vector_t<double, 3>
        TaitBryant;  // ADA = Util::getTaitBryantAngles(bunch->toLabTrafo_m.getRotation());
    Vector_t<double, 3> pmean = bunch->get_pmean();
 
    double meanEnergy   = bunch->get_meanKineticEnergy();
    double energySpread = bunch->getdE();
    double I_0 =
        4.0 * Physics::pi * Physics::epsilon_0 * Physics::c * bunch->getM() / bunch->getQ();
    double sigma = ((xsigma[0] * xsigma[0]) + (xsigma[1] * xsigma[1]))
                   / (2.0 * bunch->get_gamma() * I_0
                      * (geomvareps[0] * geomvareps[0] + geomvareps[1] * geomvareps[1]));
 
    h5_int64_t localTrackStep  = (h5_int64_t)bunch->getLocalTrackStep();
    h5_int64_t globalTrackStep = (h5_int64_t)bunch->getGlobalTrackStep();
 
    double mass   = Units::eV2GeV * bunch->getM();
    double charge = bunch->getCharge();
 
    h5_int64_t numBunch      = 1;
    h5_int64_t SteptoLastInj = 0;
 
    // ADA    bunch->get_PBounds(minP, maxP);
 
    /* ------------------------------------------------------------------------ */
 
    REPORTONERROR(H5SetStep(file_m, numSteps_m));
 
    char const* OPALFlavour = "opal-t";
    WRITESTRINGSTEPATTRIB(file_m, "OPAL_flavour", OPALFlavour);
    WRITESTEPATTRIB(Float64, file_m, "SPOS", &actPos, 1);
    WRITESTEPATTRIB(Float64, file_m, "RefPartR", (h5_float64_t*)&RefPartR, 3);
    WRITESTEPATTRIB(Float64, file_m, "centroid", (h5_float64_t*)&centroid, 3);
    WRITESTEPATTRIB(Float64, file_m, "RMSX", (h5_float64_t*)&xsigma, 3);
 
    WRITESTEPATTRIB(Float64, file_m, "RefPartP", (h5_float64_t*)&RefPartP, 3);
    WRITESTEPATTRIB(Float64, file_m, "MEANP", (h5_float64_t*)&pmean, 3);
    WRITESTEPATTRIB(Float64, file_m, "RMSP", (h5_float64_t*)&psigma, 3);
    WRITESTEPATTRIB(Float64, file_m, "TaitBryantAngles", (h5_float64_t*)&TaitBryant, 3);
 
    WRITESTEPATTRIB(Float64, file_m, "#varepsilon", (h5_float64_t*)&vareps, 3);
    WRITESTEPATTRIB(Float64, file_m, "#varepsilon-geom", (h5_float64_t*)&geomvareps, 3);
 
    WRITESTEPATTRIB(Float64, file_m, "minX", (h5_float64_t*)&rmin, 3);
    WRITESTEPATTRIB(Float64, file_m, "maxX", (h5_float64_t*)&rmax, 3);
 
    WRITESTEPATTRIB(Float64, file_m, "minP", (h5_float64_t*)&minP, 3);
    WRITESTEPATTRIB(Float64, file_m, "maxP", (h5_float64_t*)&maxP, 3);
 
    WRITESTEPATTRIB(Int64, file_m, "Step", &numSteps_m, 1);
    WRITESTEPATTRIB(Int64, file_m, "LocalTrackStep", &localTrackStep, 1);
    WRITESTEPATTRIB(Int64, file_m, "GlobalTrackStep", &globalTrackStep, 1);
 
    WRITESTEPATTRIB(Float64, file_m, "#sigma", &sigma, 1);
 
    WRITESTEPATTRIB(Float64, file_m, "TIME", &t, 1);
 
    WRITESTEPATTRIB(Float64, file_m, "ENERGY", &meanEnergy, 1);
    WRITESTEPATTRIB(Float64, file_m, "dE", &energySpread, 1);

    WRITESTEPATTRIB(Float64, file_m, "MASS", &mass, 1);
 
    WRITESTEPATTRIB(Float64, file_m, "CHARGE", &charge, 1);
 
    WRITESTEPATTRIB(Int64, file_m, "NumBunch", &numBunch, 1);
 
    WRITESTEPATTRIB(Int64, file_m, "SteptoLastInj", &SteptoLastInj, 1);
 
    try {
        Vector_t<double, 3> referenceB(
            additionalStepAttributes.at("B-ref_x"), additionalStepAttributes.at("B-ref_z"),
            additionalStepAttributes.at("B-ref_y"));
        Vector_t<double, 3> referenceE(
            additionalStepAttributes.at("E-ref_x"), additionalStepAttributes.at("E-ref_z"),
            additionalStepAttributes.at("E-ref_y"));
 
        WRITESTEPATTRIB(Float64, file_m, "B-ref", (h5_float64_t*)&referenceB, 3);
        WRITESTEPATTRIB(Float64, file_m, "E-ref", (h5_float64_t*)&referenceE, 3);
    } catch (std::out_of_range& m) {
        *ippl::Error << m.what() << endl;
 
        throw OpalException(
            "H5PartWrapperForPC::writeStepHeader", "some additional step attribute not found");
    }
 
    ++numSteps_m;
}
 
void H5PartWrapperForPT::writeStepData(PartBunch_t* bunch) {
 
    size_t numLocalParticles = bunch->getLocalNum();
       
    auto rViewDevice  = bunch->getParticleContainer()->R.getView();
    auto rView = Kokkos::create_mirror_view(rViewDevice);
    Kokkos::deep_copy(rView,rViewDevice);
    
    REPORTONERROR(H5PartSetNumParticles(file_m, (h5_size_t)numLocalParticles));
    std::vector<char> buffer(numLocalParticles * sizeof(h5_float64_t));
    char* buffer_ptr        = Util::c_data(buffer);
    h5_float64_t* f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
    //    h5_int64_t* i64buffer   = reinterpret_cast<h5_int64_t*>(buffer_ptr);
    h5_int32_t* i32buffer   = reinterpret_cast<h5_int32_t*>(buffer_ptr);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = rView(i)(0);
    WRITEDATA(Float64, file_m, "x", f64buffer);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = rView(i)(1);
    WRITEDATA(Float64, file_m, "y", f64buffer);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = rView(i)(1);
    WRITEDATA(Float64, file_m, "z", f64buffer);
 
    auto pViewDevice  = bunch->getParticleContainer()->P.getView();
    auto pView = Kokkos::create_mirror_view(pViewDevice);
    Kokkos::deep_copy(pView,pViewDevice);
 
    for (long unsigned i = 0; i < numLocalParticles; i++) 
        f64buffer[i] = pView(i)(0);
    WRITEDATA(Float64, file_m, "px", f64buffer);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = pView(i)(1);
    WRITEDATA(Float64, file_m, "py", f64buffer);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = pView(i)(2);
    WRITEDATA(Float64, file_m, "pz", f64buffer);
 
    auto qViewDevice  = bunch->getParticleContainer()->Q.getView();
    auto qView = Kokkos::create_mirror_view(qViewDevice);
    Kokkos::deep_copy(qView,qViewDevice);
    
    for (size_t i = 0; i < numLocalParticles; ++i)
        f64buffer[i] = qView(i);
    WRITEDATA(Float64, file_m, "q", f64buffer);
 
 
    /*   
 
    auto idViewDevice  = bunch->getParticleContainer()->ID.getView();
    auto idView = Kokkos::create_mirror_view(idViewDevice);
    Kokkos::deep_copy(idView,idViewDevice);
      
    for (long unsigned i = 0; i < numLocalParticles; i++)
        i64buffer[i] = idView(i);
    WRITEDATA(Int64, file_m, "id", i64buffer);
    */
    
    auto binViewDevice  = bunch->getParticleContainer()->Bin.getView();
    auto binView = Kokkos::create_mirror_view(binViewDevice);
    Kokkos::deep_copy(binView,binViewDevice);
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        i32buffer[i] = binView(i);
    WRITEDATA(Int32, file_m, "bin", i32buffer);
 
    auto spViewDevice  = bunch->getParticleContainer()->Sp.getView();
    auto spView = Kokkos::create_mirror_view(spViewDevice);
    Kokkos::deep_copy(spView,spViewDevice);
 
 
    for (size_t i = 0; i < numLocalParticles; ++i)
        i32buffer[i] = spView(i);
    WRITEDATA(Int32, file_m, "sp", i32buffer);
    
    if (Options::ebDump) {
        
        auto EViewDevice  = bunch->getParticleContainer()->E.getView();
        auto EView = Kokkos::create_mirror_view(EViewDevice);
        Kokkos::deep_copy(EView,EViewDevice);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = EView(i)(0);
        WRITEDATA(Float64, file_m, "Ex", f64buffer);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = EView(i)(1);
        WRITEDATA(Float64, file_m, "Ey", f64buffer);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = EView(i)(2);
        WRITEDATA(Float64, file_m, "Ez", f64buffer);
 
        auto BViewDevice  = bunch->getParticleContainer()->B.getView();
        auto BView = Kokkos::create_mirror_view(BViewDevice);
        Kokkos::deep_copy(BView,BViewDevice);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = BView(i)(0);
        WRITEDATA(Float64, file_m, "Bx", f64buffer);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = BView(i)(1);
        WRITEDATA(Float64, file_m, "By", f64buffer);
 
        for (size_t i = 0; i < numLocalParticles; ++i)
            f64buffer[i] = BView(i)(2);
        WRITEDATA(Float64, file_m, "Bz", f64buffer);
    }
    
    /*
    if (Options::rhoDump) {
        NDIndex<3> idx = bunch->getFieldLayout().getLocalNDIndex();
        NDIndex<3> elem;
        h5_err_t herr = H5Block3dSetView(
            file_m, idx[0].min(), idx[0].max(), idx[1].min(), idx[1].max(), idx[2].min(),
            idx[2].max());
        reportOnError(herr, __FILE__, __LINE__);
 
        std::unique_ptr<h5_float64_t[]> data(
            new h5_float64_t[(idx[0].max() + 1) * (idx[1].max() + 1) * (idx[2].max() + 1)]);
 
        int ii = 0;
        // h5block uses the fortran convention of storing data:
        // INTEGER, DIMENSION(2,3) :: a
        // => {a(1,1), a(2,1), a(1,2), a(2,2), a(1,3), a(2,3)}
        for (int i = idx[2].min(); i <= idx[2].max(); ++i) {
            for (int j = idx[1].min(); j <= idx[1].max(); ++j) {
                for (int k = idx[0].min(); k <= idx[0].max(); ++k) {
                    data[ii] = bunch->getRho(k, j, i);
                    ++ii;
                }
            }
        }
        herr = H5Block3dWriteScalarFieldFloat64(file_m, "rho", data.get());
        reportOnError(herr, __FILE__, __LINE__);

        herr = H5Block3dSetFieldOrigin(
            file_m, "rho", (h5_float64_t)bunch->get_origin()(0),
            (h5_float64_t)bunch->get_origin()(1), (h5_float64_t)bunch->get_origin()(2));
        reportOnError(herr, __FILE__, __LINE__);
 
        herr = H5Block3dSetFieldSpacing(
            file_m, "rho", (h5_float64_t)bunch->get_hr()(0), (h5_float64_t)bunch->get_hr()(1),
            (h5_float64_t)bunch->get_hr()(2));
        reportOnError(herr, __FILE__, __LINE__);
    }
    */
}