AmrPartBunch.cpp 8.28 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
//
// Class AmrPartBunch
//   This class is used to represent a bunch in AMR mode.
//
// Copyright (c) 2017 - 2019, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// 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/>.
//
21 22
#include "AmrPartBunch.h"

23
#include "Utilities/OpalException.h"
24

25
AmrPartBunch::AmrPartBunch(const PartData *ref)
26 27 28 29
    : PartBunchBase<double, 3>(new AmrPartBunch::pbase_t(new AmrLayout_t()), ref)
    , amrobj_mp(nullptr)
    , amrpbase_mp(dynamic_cast<AmrPartBunch::pbase_t*>(pbase.get()))
    , fieldlayout_m(nullptr)
frey_m's avatar
frey_m committed
30
{
frey_m's avatar
AMR:  
frey_m committed
31
    amrpbase_mp->initializeAmr();
frey_m's avatar
frey_m committed
32 33
}

frey_m's avatar
frey_m committed
34 35 36 37 38 39 40 41
AmrPartBunch::AmrPartBunch(const PartData *ref, pbase_t* pbase_p)
    : PartBunchBase<double, 3>(new AmrPartBunch::pbase_t(new AmrLayout_t(&pbase_p->getAmrLayout())), ref)
    , amrobj_mp(nullptr)
    , amrpbase_mp(dynamic_cast<AmrPartBunch::pbase_t*>(pbase.get()))
    , fieldlayout_m(nullptr)
{
    amrpbase_mp->initializeAmr();
}
frey_m's avatar
frey_m committed
42 43

AmrPartBunch::~AmrPartBunch() {
frey_m's avatar
frey_m committed
44
    
frey_m's avatar
frey_m committed
45 46
}

frey_m's avatar
AMR:  
frey_m committed
47 48 49

AmrPartBunch::pbase_t *AmrPartBunch::getAmrParticleBase() {
    return amrpbase_mp;
frey_m's avatar
frey_m committed
50 51
}

52

frey_m's avatar
AMR:  
frey_m committed
53 54 55 56 57
const AmrPartBunch::pbase_t *AmrPartBunch::getAmrParticleBase() const {
    return amrpbase_mp;
}


58
void AmrPartBunch::initialize(FieldLayout_t *fLayout) {
frey_m's avatar
frey_m committed
59
//     Layout_t* layout = static_cast<Layout_t*>(&getLayout());
60 61
}

62

frey_m's avatar
frey_m committed
63
void AmrPartBunch::do_binaryRepart() {
frey_m's avatar
frey_m committed
64
    
frey_m's avatar
frey_m committed
65 66 67 68 69 70 71 72 73 74 75 76
    if ( amrobj_mp && !amrobj_mp->isRefined() ) {
        /* In the first call to this function we
         * intialize all fine levels
         */
        amrobj_mp->initFineLevels();

    } else {
        bool isForbidTransform = amrpbase_mp->isForbidTransform();

        if ( !isForbidTransform ) {
            /*
             * regrid in boosted frame
frey_m's avatar
frey_m committed
77
             */
frey_m's avatar
frey_m committed
78 79 80 81 82 83 84 85 86 87 88 89
            this->updateLorentzFactor();
            // Lorentz transform + mapping to [-1,1]
            amrpbase_mp->domainMapping();
            amrpbase_mp->setForbidTransform(true);
        }

        this->update();

        if ( !isForbidTransform ) {
            amrpbase_mp->setForbidTransform(false);
            // map particles back + undo Lorentz transform
            amrpbase_mp->domainMapping(true);
frey_m's avatar
frey_m committed
90 91
        }
    }
frey_m's avatar
frey_m committed
92 93 94
}


frey_m's avatar
frey_m committed
95 96 97 98 99 100
Vector_t AmrPartBunch::get_hr() const {
    const double& scalefactor = amrpbase_mp->getScalingFactor();
    return hr_m * scalefactor;
}


frey_m's avatar
frey_m committed
101 102 103 104 105 106 107 108 109 110 111 112 113
void AmrPartBunch::set_meshEnlargement(double dh) {
    // set dh_m = dh
    PartBunchBase<double, 3>::set_meshEnlargement(dh);
    
    // update base geometry with new mesh enlargement
    AmrLayout_t* layout_p = &amrpbase_mp->getAmrLayout();
    layout_p->setBoundingBox(dh);
    
    // if amrobj_mp != nullptr --> we need to regrid
    this->do_binaryRepart();
}


114 115 116 117
AmrPartBunch::VectorPair_t AmrPartBunch::getEExtrema() {
    return amrobj_mp->getEExtrema();
}

118
double AmrPartBunch::getRho(int x, int y, int z) {
119 120 121 122 123 124 125 126 127 128
    /* This function is called in
     * H5PartWrapperForPC::writeStepData(PartBunchBase<double, 3>* bunch)
     * and
     * H5PartWrapperForPT::writeStepData(PartBunchBase<double, 3>* bunch)
     * in case of Options::rhoDump = true.
     * 
     * Currently, we do not support writing multilevel grid data that's why
     * we average the values down to the coarsest level.
     */
    return amrobj_mp->getRho(x, y, z);
129 130
}

131

frey_m's avatar
frey_m committed
132 133
FieldLayout_t &AmrPartBunch::getFieldLayout() {
    //TODO Implement
134
    throw OpalException("AmrPartBunch::getFieldLayout() ", "Not yet Implemented.");
135
    return *fieldlayout_m;
frey_m's avatar
frey_m committed
136
}
137 138


139 140 141 142 143 144 145 146
void AmrPartBunch::boundp() {
    IpplTimings::startTimer(boundpTimer_m);
    //if(!R.isDirty() && stateOfLastBoundP_ == unit_state_) return;
    if ( !(R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_) return; //-DW

    stateOfLastBoundP_ = unit_state_;
    
    if ( amrobj_mp ) {
frey_m's avatar
frey_m committed
147 148 149 150
        /* we do an explicit domain mapping of the particles and then
         * forbid it during the regrid process, this way it's only
         * executed ones --> saves computation
         */
frey_m's avatar
frey_m committed
151
        bool isForbidTransform = amrpbase_mp->isForbidTransform();
frey_m's avatar
frey_m committed
152
            
frey_m's avatar
frey_m committed
153
        if ( !isForbidTransform ) {
154 155
            this->updateLorentzFactor();
            // Lorentz transform + mapping to [-1,1]
frey_m's avatar
frey_m committed
156 157
            amrpbase_mp->domainMapping();
            amrpbase_mp->setForbidTransform(true);
frey_m's avatar
frey_m committed
158 159 160 161 162
        }
        
        this->update();
        
        if ( !isForbidTransform ) {
frey_m's avatar
frey_m committed
163
            amrpbase_mp->setForbidTransform(false);
164
            // map particles back + undo Lorentz transform
frey_m's avatar
frey_m committed
165
            amrpbase_mp->domainMapping(true);
166
        }
frey_m's avatar
frey_m committed
167
        
168 169
    } else {
        // At this point an amrobj_mp needs already be set
frey_m's avatar
AMR:  
frey_m committed
170 171
        throw GeneralClassicException("AmrPartBunch::boundp() ",
                                      "AmrObject pointer is not set.");
172 173 174
    }
    
    R.resetDirtyFlag();
175
    
176 177 178 179
    IpplTimings::stopTimer(boundpTimer_m);
}


180
void AmrPartBunch::computeSelfFields() {
181 182
    IpplTimings::startTimer(selfFieldTimer_m);
    
183
    if ( !fs_m->hasValidSolver() )
frey_m's avatar
AMR:  
frey_m committed
184 185
        throw OpalException("AmrPartBunch::computeSelfFields() ",
                            "No field solver.");
186
    
187
    amrobj_mp->computeSelfFields();
188 189
    
    IpplTimings::stopTimer(selfFieldTimer_m);
190 191 192
}


193 194 195 196
void AmrPartBunch::computeSelfFields(int bin) {
    IpplTimings::startTimer(selfFieldTimer_m);
    amrobj_mp->computeSelfFields(bin);
    IpplTimings::stopTimer(selfFieldTimer_m);
197 198 199 200
}


void AmrPartBunch::computeSelfFields_cycl(double gamma) {
201
    IpplTimings::startTimer(selfFieldTimer_m);
202
    amrobj_mp->computeSelfFields_cycl(gamma);
203
    IpplTimings::stopTimer(selfFieldTimer_m);
204 205 206
}


207 208
void AmrPartBunch::computeSelfFields_cycl(int bin) {
    IpplTimings::startTimer(selfFieldTimer_m);
frey_m's avatar
frey_m committed
209 210 211 212 213 214 215
    
    /* make sure it is refined in multi-bunch case
     */
    if ( !amrobj_mp->isRefined() ) {
        amrobj_mp->initFineLevels();
    }
    
216
    amrobj_mp->computeSelfFields_cycl(bin);
frey_m's avatar
frey_m committed
217
    
218
    IpplTimings::stopTimer(selfFieldTimer_m);
219
}
frey_m's avatar
frey_m committed
220 221


222 223 224 225 226
void AmrPartBunch::setAmrDomainRatio(const std::vector<double>& ratio) {
    AmrLayout_t* layout_p = &amrpbase_mp->getAmrLayout();
    layout_p->setDomainRatio(ratio);
}

frey_m's avatar
frey_m committed
227
void AmrPartBunch::gatherLevelStatistics() {
frey_m's avatar
frey_m committed
228
    int nLevel = amrobj_mp->maxLevel() + 1;
frey_m's avatar
AMR:  
frey_m committed
229
    
230
    std::unique_ptr<size_t[]> partPerLevel( new size_t[nLevel] );
frey_m's avatar
frey_m committed
231
    globalPartPerLevel_m.reset( new size_t[nLevel] );
frey_m's avatar
AMR:  
frey_m committed
232
    
233
    for (int i = 0; i < nLevel; ++i)
frey_m's avatar
frey_m committed
234
        partPerLevel[i] = globalPartPerLevel_m[i] = 0.0;
frey_m's avatar
frey_m committed
235 236 237
    
    // do not modify LocalNumPerLevel in here!!!
    auto& LocalNumPerLevel = amrpbase_mp->getLocalNumPerLevel();
frey_m's avatar
AMR:  
frey_m committed
238
        
frey_m's avatar
frey_m committed
239 240
    for (size_t i = 0; i < LocalNumPerLevel.size(); ++i)
        partPerLevel[i] = LocalNumPerLevel[i];
241 242 243 244
    
    reduce(*partPerLevel.get(),
           *globalPartPerLevel_m.get(),
           nLevel, std::plus<size_t>());
frey_m's avatar
frey_m committed
245 246 247 248 249
}


const size_t& AmrPartBunch::getLevelStatistics(int l) const {
    return globalPartPerLevel_m[l];
frey_m's avatar
AMR:  
frey_m committed
250 251 252
}


frey_m's avatar
frey_m committed
253
void AmrPartBunch::updateLorentzFactor(int bin) {
254 255 256 257 258
    double gamma = this->get_gamma();

    if ( this->weHaveBins() ) {
        gamma = this->getBinGamma(bin);
    }
frey_m's avatar
frey_m committed
259 260 261 262 263 264 265 266 267 268 269 270
    
    
    /* At the beginning of simulation the Lorentz factor
     * is not yet set since PartBunchBase::calcBeamParameters
     * is not yet called. Therefore, we do this ugly workaround.
     */
    bool init = true;
    if ( gamma >= 1.0 ) {
        init = false;
    }
    
    if ( init ) {
271 272 273
        gamma = 1.0;
    }

274 275 276 277 278
    updateLorentzFactor(gamma);
}


void AmrPartBunch::updateLorentzFactor(double gamma) {
279 280 281 282 283 284 285 286 287 288 289 290 291
    // keep all 1.0, except longitudinal direction
    Vector_t lorentzFactor(1.0, 1.0, 1.0);

    if (OpalData::getInstance()->isInOPALCyclMode()) {
        lorentzFactor[1] = gamma;
    } else {
        lorentzFactor[2] = gamma;
    }

    amrpbase_mp->setLorentzFactor(lorentzFactor);
}


292 293
void AmrPartBunch::updateFieldContainers_m() {
    
frey_m's avatar
frey_m committed
294 295
}

frey_m's avatar
frey_m committed
296
void AmrPartBunch::updateDomainLength(Vektor<int, 3>& grid) {
297
    grid = amrobj_mp->getBaseLevelGridPoints();
frey_m's avatar
frey_m committed
298
}
299 300 301 302


void AmrPartBunch::updateFields(const Vector_t& hr, const Vector_t& origin) {
    //TODO regrid; called in boundp()
303
//     amrobj_mp->updateMesh();
304
}