From d5bdff2c29edca97b39ed66cfd9e1b11908522f1 Mon Sep 17 00:00:00 2001 From: Matthias Frey <matthias.frey@psi.ch> Date: Wed, 29 Apr 2020 09:03:05 +0200 Subject: [PATCH] Merge AMR fork in OPAL master. --- optimizer/Expression/CMakeLists.txt | 1 + optimizer/Expression/SeptumExpr.h | 80 +++++++++ optimizer/Util/CMakeLists.txt | 1 + optimizer/Util/ProbeHistReader.cpp | 128 +++++++++++++ optimizer/Util/ProbeHistReader.h | 70 ++++++++ src/Amr/AbstractAmrWriter.h | 24 ++- src/Amr/AmrBoxLib.cpp | 170 +++++++++++------- src/Amr/AmrBoxLib.h | 41 ++++- src/Amr/AmrDefs.h | 20 +++ src/Amr/AmrObject.cpp | 23 +++ src/Amr/AmrObject.h | 35 ++-- src/Amr/AmrYtWriter.cpp | 26 ++- src/Amr/AmrYtWriter.h | 34 ++-- src/Amr/BoxLibLayout.h | 53 ++++-- src/Amr/BoxLibLayout.hpp | 47 ++++- src/Amr/BoxLibParticle.h | 33 +++- src/Amr/BoxLibParticle.hpp | 29 ++- src/Classic/Algorithms/AmrPartBunch.cpp | 4 +- src/Classic/Algorithms/AmrPartBunch.h | 4 +- src/Classic/Algorithms/PartBunchBase.hpp | 2 +- src/Classic/Algorithms/Tracker.h | 4 +- src/Distribution/ClosedOrbitFinder.h | 2 + src/Optimize/OptimizeCmd.cpp | 5 + src/Sample/SampleCmd.cpp | 5 + src/Solvers/AMR_MG/AmrMultiGrid.cpp | 5 + src/Solvers/AMR_MG/AmrMultiGrid.h | 8 +- src/Solvers/AMReXSolvers/MLPoissonSolver.cpp | 22 ++- src/Solvers/AMReXSolvers/MLPoissonSolver.h | 20 +++ src/Solvers/AmrPoissonSolver.h | 6 +- .../BoxLibSolvers/FMGPoissonSolver.cpp | 24 ++- src/Solvers/BoxLibSolvers/FMGPoissonSolver.h | 22 +++ src/Solvers/PoissonSolver.h | 12 +- src/Structure/DataSink.cpp | 8 +- src/Structure/FieldSolver.cpp | 10 ++ src/Structure/GridLBalWriter.cpp | 4 +- src/Structure/GridLBalWriter.h | 2 +- src/Structure/H5PartWrapperForPC.cpp | 2 + src/Structure/LBalWriter.cpp | 6 +- src/Structure/LBalWriter.h | 4 + src/opal.cpp | 7 +- 40 files changed, 846 insertions(+), 157 deletions(-) create mode 100644 optimizer/Expression/SeptumExpr.h create mode 100644 optimizer/Util/ProbeHistReader.cpp create mode 100644 optimizer/Util/ProbeHistReader.h diff --git a/optimizer/Expression/CMakeLists.txt b/optimizer/Expression/CMakeLists.txt index b1756543b..7c93c9d62 100644 --- a/optimizer/Expression/CMakeLists.txt +++ b/optimizer/Expression/CMakeLists.txt @@ -8,6 +8,7 @@ SET (_EXPR_SRCS SumErrSqRadialPeak.cpp Parser/expression.cpp Parser/evaluator.cpp + SeptumExpr.h ) add_sources(${_EXPR_SRCS}) diff --git a/optimizer/Expression/SeptumExpr.h b/optimizer/Expression/SeptumExpr.h new file mode 100644 index 000000000..2db55b002 --- /dev/null +++ b/optimizer/Expression/SeptumExpr.h @@ -0,0 +1,80 @@ +// +// Struct SeptumExpr +// Objective to obtain a nice septum in cyclotron simulations. +// +// Copyright (c) 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/>. +// +#ifndef __SEPTUM_EXPRESSION_H__ +#define __SEPTUM_EXPRESSION_H__ + +#include <string> +#include <iostream> + +#include "boost/variant/get.hpp" +#include "boost/variant/variant.hpp" +#include "boost/smart_ptr.hpp" + +#include "Util/Types.h" +#include "Util/PeakReader.h" +#include "Expression/Parser/function.hpp" + +#include "Util/ProbeHistReader.h" + +struct SeptumExpr { + + static const std::string name; + + Expressions::Result_t operator()(client::function::arguments_t args) { + if (args.size() != 1) { + throw OptPilotException("SeptumExpr::operator()", + "SeptumExpr expects 1 arguments, " + std::to_string(args.size()) + " given"); + } + + std::string probe = boost::get<std::string>(args[0]); + + bool is_valid = true; + + double result = 0.0; + + try { + boost::scoped_ptr<PeakReader> sim_peaks(new PeakReader(probe + std::string(".peaks"))); + sim_peaks->parseFile(); + + boost::scoped_ptr<ProbeHistReader> sim_hist(new ProbeHistReader(probe + std::string(".hist"))); + sim_hist->parseFile(); + + double upperBound = 0.0; + double lowerBound = 0.0; + + size_t nTurns = sim_peaks->getNumberOfPeaks(); + sim_peaks->getPeak(nTurns, upperBound); + sim_peaks->getPeak(nTurns - 1, lowerBound); + + result = sim_hist->minimum(lowerBound, upperBound); + + } catch (OptPilotException &ex) { + std::cout << "Exception while getting septum value " + << ex.what() + << std::endl; + is_valid = false; + } + + return boost::make_tuple(result, is_valid); + } +}; + +#endif diff --git a/optimizer/Util/CMakeLists.txt b/optimizer/Util/CMakeLists.txt index fa8990203..d318cff56 100644 --- a/optimizer/Util/CMakeLists.txt +++ b/optimizer/Util/CMakeLists.txt @@ -5,6 +5,7 @@ SET (_UTIL_SRCS PeakReader.cpp ProbeReader.cpp + ProbeHistReader.cpp SDDSParser.cpp SDDSParser/array.cpp SDDSParser/associate.cpp diff --git a/optimizer/Util/ProbeHistReader.cpp b/optimizer/Util/ProbeHistReader.cpp new file mode 100644 index 000000000..9f5974e2e --- /dev/null +++ b/optimizer/Util/ProbeHistReader.cpp @@ -0,0 +1,128 @@ +// +// Class ProbeHistReader +// Implements a parser and value extractor for hist files (*.hist). +// It is for example used together with the septum objective. +// A histogram file is generated by the OPAL probe element. +// +// Copyright (c) 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/>. +// +#include <iterator> +#include <regex> + +#include "ProbeHistReader.h" +#include "Util/OptPilotException.h" + +ProbeHistReader::ProbeHistReader(std::string filename) + : filename_m(filename) + , rmin_m(0.0) + , binwidth_m(0.0) + , bincount_m(0) +{ } + + +void ProbeHistReader::parseFile() { + + std::ifstream histfile; + + histfile.open(filename_m.c_str(), std::ios::in); + + if (!histfile) { + throw OptPilotException("ProbeHistReader::parseFile()", + "Error opening file " + filename_m); + } + + parseHeader(histfile); + + std::istream_iterator<size_t> it(histfile); + while ( it != std::istream_iterator<size_t>() ) { + bincount_m.push_back(*it); + ++it; + } + + histfile.close(); + + if (histfile.is_open()) { + throw OptPilotException("ProbeHistReader::parseFile()", + "Error closing file " + filename_m); + } +} + + +size_t ProbeHistReader::minimum(double lower, double upper) { + size_t lidx = (lower - rmin_m) / binwidth_m; + size_t uidx = (upper - rmin_m) / binwidth_m; + + if (lidx >= uidx) { + throw OptPilotException("ProbeHistReader::minimum()", + "Lower index >= upper index: " + std::to_string(lidx) + + " >= " + std::to_string(uidx)); + } + + if (uidx >= bincount_m.size()) { + throw OptPilotException("ProbeHistReader::minimum()", + "Index >= number of bins: " + std::to_string(uidx) + + " >= " + std::to_string(bincount_m.size())); + } + + container_t::iterator beg = std::begin(bincount_m); + std::advance(beg, lidx); + + container_t::iterator end = std::begin(bincount_m); + std::advance(end, uidx); + + container_t::iterator result = std::min_element(beg, end); + + if (result == bincount_m.end()) { + throw OptPilotException("ProbeHistReader::minimum()", + "No minimum between " + std::to_string(lower) + + " and " + std::to_string(upper) + " found."); + } + + return double(*result); +} + + +void ProbeHistReader::parseHeader(std::ifstream& ifs) { + std::string header; + std::getline(ifs, header); + + if (header.find("# Histogram bin counts") == std::string::npos) { + throw OptPilotException("ProbeHistReader::parseHeader()", + "Error reading file " + filename_m); + } + + const std::regex re("\\.*\\) (.*) mm (.*) mm (.*) (.*) mm"); + std::smatch match; + std::regex_search(header, match, re); + + if ( match.size() != 5 ) { + throw OptPilotException("ProbeHistReader::parseHeader()", + "Error parsing header of file " + filename_m); + } + + rmin_m = getValue<double>(match[1].str()); + double rmax = getValue<double>(match[2].str()); + size_t nbins = getValue<size_t>(match[3].str()); + binwidth_m = getValue<double>(match[4].str()); + + if ( rmin_m >= rmax ) { + throw OptPilotException("ProbeHistReader::parseHeader()", + "Not a valid histogram file " + filename_m); + } + + bincount_m.reserve(nbins); +} diff --git a/optimizer/Util/ProbeHistReader.h b/optimizer/Util/ProbeHistReader.h new file mode 100644 index 000000000..187f9f4cf --- /dev/null +++ b/optimizer/Util/ProbeHistReader.h @@ -0,0 +1,70 @@ +// +// Class ProbeHistReader +// Implements a parser and value extractor for hist files (*.hist). +// It is for example used together with the septum objective. +// A histogram file is generated by the OPAL probe element. +// +// Copyright (c) 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/>. +// +#ifndef __PROBE_HIST_READER_H__ +#define __PROBE_HIST_READER_H__ + +#include <fstream> +#include <sstream> +#include <string> +#include <vector> + +class ProbeHistReader { + +public: + typedef std::vector<size_t> container_t; + + ProbeHistReader(std::string filename); + + void parseFile(); + + /** + * @param lower radius [mm] + * @param upper radius [mm] + */ + size_t minimum(double lower, double upper); + +private: + void parseHeader(std::ifstream& ifs); + + template <typename T> + T getValue(const std::string& s); + + /// Histogram file + std::string filename_m; + + double rmin_m; // start radius of probe in mm + double binwidth_m; // size of each bin in mm + + container_t bincount_m; +}; + + +template <typename T> +T ProbeHistReader::getValue(const std::string& s) { + std::istringstream ss(s); + T res; + ss >> res; + return res; +} + +#endif diff --git a/src/Amr/AbstractAmrWriter.h b/src/Amr/AbstractAmrWriter.h index e592d1655..ae2e3a38e 100644 --- a/src/Amr/AbstractAmrWriter.h +++ b/src/Amr/AbstractAmrWriter.h @@ -1,3 +1,23 @@ +// +// Class AbstractAmrWriter +// Abstract base class for writing AMR data to output files. +// +// Copyright (c) 2016 - 2020, 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/>. +// #ifndef ABSTRACT_AMR_WRITER_H #define ABSTRACT_AMR_WRITER_H @@ -5,10 +25,6 @@ #include "Amr/AmrDefs.h" #include "Algorithms/AmrPartBunch.h" -/*! - * Abstract base class for writing AMR data to files in - * order to plot them. - */ class AbstractAmrWriter { public: diff --git a/src/Amr/AmrBoxLib.cpp b/src/Amr/AmrBoxLib.cpp index f49ae94c0..f5a985c73 100644 --- a/src/Amr/AmrBoxLib.cpp +++ b/src/Amr/AmrBoxLib.cpp @@ -1,3 +1,27 @@ +// +// Class AmrBoxLib +// Concrete AMR object. It is based on the AMReX library +// (cf. https://amrex-codes.github.io/ or https://ccse.lbl.gov/AMReX/). +// AMReX is the successor of BoxLib. This class represents the interface +// to AMReX and the AMR framework in OPAL. It implements the functions of +// the AmrObject class. +// +// Copyright (c) 2016 - 2020, 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/>. +// #include "AmrBoxLib.h" #include "Algorithms/AmrPartBunch.h" @@ -83,6 +107,8 @@ void AmrBoxLib::regrid(double time) { << "* Old finest level: " << finest_level << endl; + this->preRegrid_m(); + /* ATTENTION: The bunch might be updated during * the regrid process! * We regrid from base level 0 up to the finest level. @@ -189,7 +215,7 @@ AmrBoxLib::VectorPair_t AmrBoxLib::getEExtrema() { } -double AmrBoxLib::getRho(int x, int y, int z) { +double AmrBoxLib::getRho(int /*x*/, int /*y*/, int /*z*/) { //TODO throw OpalException("AmrBoxLib::getRho(x, y, z)", "Not yet Implemented."); return 0.0; @@ -202,7 +228,7 @@ void AmrBoxLib::computeSelfFields() { } -void AmrBoxLib::computeSelfFields(int bin) { +void AmrBoxLib::computeSelfFields(int /*bin*/) { //TODO throw OpalException("AmrBoxLib::computeSelfFields(int bin)", "Not yet Implemented."); } @@ -232,51 +258,11 @@ void AmrBoxLib::computeSelfFields_cycl(double gamma) { amrpbase_p->update(); amrpbase_p->setForbidTransform(false); - /// from charge (C) to charge density (C/m^3). - amrpbase_p->scatter(bunch_mp->Q, this->rho_m, bunch_mp->R, 0, finest_level, bunch_mp->Bin); - - int baseLevel = 0; - int nLevel = finest_level + 1; double invGamma = 1.0 / gamma; - - - // mesh scaling for solver - meshScaling_m = Vector_t(1.0, 1.0, 1.0); - - // charge density is in rho_m - // calculate Possion equation (with coefficient: -1/(eps)) - for (int i = 0; i <= finest_level; ++i) { - if ( this->rho_m[i]->contains_nan(false) ) - throw OpalException("AmrBoxLib::computeSelfFields_cycl(double gamma) ", - "NANs at level " + std::to_string(i) + "."); - this->rho_m[i]->mult(-1.0 / Physics::epsilon_0, 0, 1); - } - - // find maximum and normalize each level (faster convergence) - double l0norm = 1.0; - for (int i = 0; i <= finest_level; ++i) - l0norm = std::max(l0norm, this->rho_m[i]->norm0(0)); - - for (int i = 0; i <= finest_level; ++i) { - this->rho_m[i]->mult(1.0 / l0norm, 0, 1); - } + int nLevel = finest_level + 1; - PoissonSolver *solver = bunch_mp->getFieldSolver(); + double l0norm = this->solvePoisson_m(); - IpplTimings::startTimer(this->amrSolveTimer_m); - solver->solve(rho_m, phi_m, efield_m, baseLevel, finest_level); - IpplTimings::stopTimer(this->amrSolveTimer_m); - - // make sure ghost cells are filled - for (int i = 0; i <= finest_level; ++i) { - phi_m[i]->FillBoundary(geom[i].periodicity()); - for (int j = 0; j < AMREX_SPACEDIM; ++j) { - efield_m[i][j]->FillBoundary(geom[i].periodicity()); - } - } - - this->fillPhysbc_m(*(this->phi_m[0]), 0); - /* apply scale of electric-field in order to undo the transformation * + undo normalization */ @@ -286,7 +272,7 @@ void AmrBoxLib::computeSelfFields_cycl(double gamma) { this->efield_m[i][j]->mult(scalefactor * scalefactor * l0norm, 0, 1); } } - + for (int i = 0; i <= finest_level; ++i) { for (int j = 0; j < AMREX_SPACEDIM; ++j) { if ( this->efield_m[i][j]->contains_nan(false) ) @@ -294,8 +280,7 @@ void AmrBoxLib::computeSelfFields_cycl(double gamma) { "Ef: NANs at level " + std::to_string(i) + "."); } } - - + amrpbase_p->gather(bunch_mp->Ef, this->efield_m, bunch_mp->R, 0, finest_level); // undo domain change + undo Lorentz transform @@ -539,7 +524,7 @@ inline double AmrBoxLib::getT() const { } -void AmrBoxLib::redistributeGrids(int how) { +void AmrBoxLib::redistributeGrids(int /*how*/) { // // // copied + modified version of AMReX_Amr.cpp // AmrProcMap_t::InitProximityMap(); @@ -600,7 +585,7 @@ void AmrBoxLib::redistributeGrids(int how) { } -void AmrBoxLib::RemakeLevel (int lev, AmrReal_t time, +void AmrBoxLib::RemakeLevel (int lev, AmrReal_t /*time*/, const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap) { @@ -633,7 +618,7 @@ void AmrBoxLib::RemakeLevel (int lev, AmrReal_t time, } -void AmrBoxLib::MakeNewLevel (int lev, AmrReal_t time, +void AmrBoxLib::MakeNewLevel (int lev, AmrReal_t /*time*/, const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap) { @@ -714,17 +699,17 @@ void AmrBoxLib::ErrorEst(int lev, TagBoxArray_t& tags, } -void AmrBoxLib::MakeNewLevelFromScratch(int lev, AmrReal_t time, - const AmrGrid_t& ba, - const AmrProcMap_t& dm) +void AmrBoxLib::MakeNewLevelFromScratch(int /*lev*/, AmrReal_t /*time*/, + const AmrGrid_t& /*ba*/, + const AmrProcMap_t& /*dm*/) { throw OpalException("AmrBoxLib::MakeNewLevelFromScratch()", "Shouldn't be called."); } -void AmrBoxLib::MakeNewLevelFromCoarse (int lev, AmrReal_t time, - const AmrGrid_t& ba, - const AmrProcMap_t& dm) +void AmrBoxLib::MakeNewLevelFromCoarse (int /*lev*/, AmrReal_t /*time*/, + const AmrGrid_t& /*ba*/, + const AmrProcMap_t& /*dm*/) { throw OpalException("AmrBoxLib::MakeNewLevelFromCoarse()", "Shouldn't be called."); } @@ -766,6 +751,22 @@ void AmrBoxLib::doRegrid_m(int lbase, double time) { } +void AmrBoxLib::preRegrid_m() { + /* In case of E-field or potential tagging + * in combination of binning we need to make + * sure we do not lose accuracy when tagging since + * the grid data of the potential or e-field are only + * non-zero for the last bin causing the tagging to refine + * only there. + * So, we need to solve the Poisson problem first assuming + * a single bin only + */ + if ( tagging_m == POTENTIAL || tagging_m == EFIELD ) { + this->solvePoisson_m(); + } +} + + void AmrBoxLib::postRegrid_m(int old_finest) { /* ATTENTION: In this call the bunch has to be updated * since each particle needs to be assigned to its latest @@ -805,8 +806,57 @@ void AmrBoxLib::postRegrid_m(int old_finest) { } +double AmrBoxLib::solvePoisson_m() { + AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase(); + + /// from charge (C) to charge density (C/m^3). + amrpbase_p->scatter(bunch_mp->Q, this->rho_m, bunch_mp->R, 0, finest_level, bunch_mp->Bin); + + int baseLevel = 0; + + // mesh scaling for solver + meshScaling_m = Vector_t(1.0, 1.0, 1.0); + + // charge density is in rho_m + // calculate Possion equation (with coefficient: -1/(eps)) + for (int i = 0; i <= finest_level; ++i) { + if ( this->rho_m[i]->contains_nan(false) ) + throw OpalException("AmrBoxLib::solvePoisson_m() ", + "NANs at level " + std::to_string(i) + "."); + this->rho_m[i]->mult(-1.0 / Physics::epsilon_0, 0, 1); + } + + // find maximum and normalize each level (faster convergence) + double l0norm = 1.0; + for (int i = 0; i <= finest_level; ++i) + l0norm = std::max(l0norm, this->rho_m[i]->norm0(0)); + + for (int i = 0; i <= finest_level; ++i) { + this->rho_m[i]->mult(1.0 / l0norm, 0, 1); + } + + PoissonSolver *solver = bunch_mp->getFieldSolver(); + + IpplTimings::startTimer(this->amrSolveTimer_m); + solver->solve(rho_m, phi_m, efield_m, baseLevel, finest_level); + IpplTimings::stopTimer(this->amrSolveTimer_m); + + // make sure ghost cells are filled + for (int i = 0; i <= finest_level; ++i) { + phi_m[i]->FillBoundary(geom[i].periodicity()); + for (int j = 0; j < AMREX_SPACEDIM; ++j) { + efield_m[i][j]->FillBoundary(geom[i].periodicity()); + } + } + + this->fillPhysbc_m(*(this->phi_m[0]), 0); + + return l0norm; +} + + void AmrBoxLib::tagForChargeDensity_m(int lev, TagBoxArray_t& tags, - AmrReal_t time, int ngrow) + AmrReal_t /*time*/, int /*ngrow*/) { // we need to update first AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase(); @@ -978,7 +1028,7 @@ void AmrBoxLib::tagForEfield_m(int lev, TagBoxArray_t& tags, void AmrBoxLib::tagForMomenta_m(int lev, TagBoxArray_t& tags, - AmrReal_t time, int ngrow) + AmrReal_t /*time*/, int /*ngrow*/) { AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase(); // we need to update first @@ -1037,7 +1087,7 @@ void AmrBoxLib::tagForMomenta_m(int lev, TagBoxArray_t& tags, void AmrBoxLib::tagForMaxNumParticles_m(int lev, TagBoxArray_t& tags, - AmrReal_t time, int ngrow) + AmrReal_t /*time*/, int /*ngrow*/) { AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase(); // we need to update first @@ -1087,7 +1137,7 @@ void AmrBoxLib::tagForMaxNumParticles_m(int lev, TagBoxArray_t& tags, void AmrBoxLib::tagForMinNumParticles_m(int lev, TagBoxArray_t& tags, - AmrReal_t time, int ngrow) + AmrReal_t /*time*/, int /*ngrow*/) { AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase(); // we need to update first diff --git a/src/Amr/AmrBoxLib.h b/src/Amr/AmrBoxLib.h index d4b295234..8742ad72a 100644 --- a/src/Amr/AmrBoxLib.h +++ b/src/Amr/AmrBoxLib.h @@ -1,3 +1,27 @@ +// +// Class AmrBoxLib +// Concrete AMR object. It is based on the AMReX library +// (cf. https://amrex-codes.github.io/ or https://ccse.lbl.gov/AMReX/). +// AMReX is the successor of BoxLib. This class represents the interface +// to AMReX and the AMR framework in OPAL. It implements the functions of +// the AmrObject class. +// +// Copyright (c) 2016 - 2020, 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/>. +// #ifndef AMR_BOXLIB_H #define AMR_BOXLIB_H @@ -9,12 +33,6 @@ class AmrPartBunch; #include <AMReX_AmrMesh.H> #include <AMReX.H> -/*! - * Concrete AMR object. It is based on the - * <a href="https://ccse.lbl.gov/AMReX/">AMReX</a> library - * developed at LBNL. This library is the successor of - * <a href="https://ccse.lbl.gov/BoxLib/">BoxLib</a>. - */ class AmrBoxLib : public AmrObject, public amrex::AmrMesh { @@ -200,6 +218,14 @@ private: */ void doRegrid_m(int lbase, double time); + + /*! + * Called within doRegrid_m(). Especially used + * for potential and e-field tagging in combination + * with binning. + */ + void preRegrid_m(); + /*! * Called within doRegrid_m(). Redistribute * particles and delete levels on particle side. @@ -207,6 +233,9 @@ private: */ void postRegrid_m(int old_finest); + + double solvePoisson_m(); + /* ATTENTION * The tagging routines assume the particles to be in the diff --git a/src/Amr/AmrDefs.h b/src/Amr/AmrDefs.h index 5a024beea..abf736fe0 100644 --- a/src/Amr/AmrDefs.h +++ b/src/Amr/AmrDefs.h @@ -1,3 +1,23 @@ +// +// File AmrDefs +// AMR namespace with typedefs of AMReX classes. +// +// Copyright (c) 2016 - 2020, 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/>. +// #ifndef AMR_DEFS_H #define AMR_DEFS_H diff --git a/src/Amr/AmrObject.cpp b/src/Amr/AmrObject.cpp index e6a759bb8..6697c124d 100644 --- a/src/Amr/AmrObject.cpp +++ b/src/Amr/AmrObject.cpp @@ -1,3 +1,26 @@ +// +// Class AmrObject +// The AMR interface to OPAL. A new AMR library needs +// to inherit from this class in order to work properly +// with OPAL. Among other things it specifies the refinement +// strategies. +// +// Copyright (c) 2016 - 2020, 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/>. +// #include "AmrObject.h" AmrObject::AmrObject() diff --git a/src/Amr/AmrObject.h b/src/Amr/AmrObject.h index 937ee5327..1fa293413 100644 --- a/src/Amr/AmrObject.h +++ b/src/Amr/AmrObject.h @@ -1,3 +1,26 @@ +// +// Class AmrObject +// The AMR interface to OPAL. A new AMR library needs +// to inherit from this class in order to work properly +// with OPAL. Among other things it specifies the refinement +// strategies. +// +// Copyright (c) 2016 - 2020, 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/>. +// #ifndef AMR_OBECT_H #define AMR_OBECT_H @@ -5,16 +28,6 @@ #include "Algorithms/PBunchDefs.h" -// #include "Algorithms/AmrPartBunch.h" - -// class AmrPartBunch; - -/** - * The AMR interface to OPAL. A new AMR library needs - * to inherit from this class in order to work properly - * with OPAL. Among other things it specifies the refinement - * strategies. - */ class AmrObject { public: @@ -147,7 +160,7 @@ public: * Rebalance the grids among the * cores */ - virtual void redistributeGrids(int how) { } + virtual void redistributeGrids(int /*how*/) { } /*! * Used in AmrPartBunch to check if we need to refine diff --git a/src/Amr/AmrYtWriter.cpp b/src/Amr/AmrYtWriter.cpp index dbcd12666..b3a0ec8e6 100644 --- a/src/Amr/AmrYtWriter.cpp +++ b/src/Amr/AmrYtWriter.cpp @@ -1,3 +1,27 @@ +// +// Class AmrYtWriter +// This concrete class writes output files that are readable by yt +// (cf. http://yt-project.org/). We have a fork of yt in +// the repository at https://gitlab.psi.ch/frey_m/yt. +// The functions of this class are copied from AMReX and modified to fit +// our needs. +// +// Copyright (c) 2016 - 2020, 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/>. +// #include "AmrYtWriter.h" #include <AMReX_Utility.H> @@ -271,7 +295,7 @@ void AmrYtWriter::writeFields(const amr::AmrScalarFieldContainer_t& rho, void AmrYtWriter::writeBunch(const AmrPartBunch* bunch_p, - const double& time, + const double& /*time*/, const double& gamma) { /* According to diff --git a/src/Amr/AmrYtWriter.h b/src/Amr/AmrYtWriter.h index 2d6e80834..73572ccec 100644 --- a/src/Amr/AmrYtWriter.h +++ b/src/Amr/AmrYtWriter.h @@ -1,3 +1,27 @@ +// +// Class AmrYtWriter +// This concrete class writes output files that are readable by yt +// (cf. http://yt-project.org/). We have a fork of yt in +// the repository at https://gitlab.psi.ch/frey_m/yt. +// The functions of this class are copied from AMReX and modified to fit +// our needs. +// +// Copyright (c) 2016 - 2020, 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/>. +// #ifndef AMR_YT_WRITER_H #define AMR_YT_WRITER_H @@ -7,16 +31,6 @@ #include <vector> -/*! - * This concrete class writes output files - * that are readable by - * <a href="http://yt-project.org/">yt</a>. - * The format can be read by the fork - * <a href="https://Yavin_4@bitbucket.org/Yavin_4/opal-yt">opal-yt</a> - * that uses accelerator units instead of astrophysics units. - * The functions are copied from AMReX and modified to fit - * our needs. - */ class AmrYtWriter : public AbstractAmrWriter { public: diff --git a/src/Amr/BoxLibLayout.h b/src/Amr/BoxLibLayout.h index 064f1a13a..29bab5306 100644 --- a/src/Amr/BoxLibLayout.h +++ b/src/Amr/BoxLibLayout.h @@ -1,3 +1,38 @@ +// +// Class BoxLibLayout +// In contrast to AMReX, OPAL is optimized for the +// distribution of particles to cores. In AMReX the ParGDB object +// is responsible for the particle to core distribution. This +// layout is derived from this object and does all important +// bunch updates. It is the interface for AMReX and Ippl. +// +// In AMReX, the geometry, i.e. physical domain, is fixed +// during the whole computation. Particles leaving the domain +// would be deleted. In order to prevent this we map the particles +// onto the domain [-1, 1]^3. Furthermore, it makes sure +// that we have enougth grid points to represent the bunch +// when its charges are scattered on the grid for the self-field +// computation. +// +// The self-field computation and the particle-to-core update +// are performed in the particle mapped domain. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, 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/>. +// #ifndef BOXLIB_LAYOUT_H #define BOXLIB_LAYOUT_H @@ -8,24 +43,6 @@ #include <AMReX_ParGDB.H> -/*! - * In contrast to AMReX, OPAL is optimized for - * distribution particles to cores. In AMReX the ParGDB object - * is responsible for the particle to core distribution. This - * layout is derived from this object and does all important - * bunch updates. It is the interface for AMReX and Ippl. - * - * In AMReX the geometry, i.e. physical domain, is fixed - * during the whole computation. Particles leaving the domain - * would be deleted. In order to prevent this we map the particles - * onto the domain \f$[-1, 1]^3\f$. Furthermore, it makes sure - * that we have enougth grid points to represent the bunch - * when its charges are scattered on the grid for the self-field - * computation. - * - * The self-field computation and the particle-to-core update - * are performed in the particle mapped domain. - */ template<class T, unsigned Dim> class BoxLibLayout : public ParticleAmrLayout<T, Dim>, public amrex::ParGDB diff --git a/src/Amr/BoxLibLayout.hpp b/src/Amr/BoxLibLayout.hpp index 9057bc2ce..c898753d4 100644 --- a/src/Amr/BoxLibLayout.hpp +++ b/src/Amr/BoxLibLayout.hpp @@ -1,3 +1,38 @@ +// +// Class BoxLibLayout +// In contrast to AMReX, OPAL is optimized for the +// distribution of particles to cores. In AMReX the ParGDB object +// is responsible for the particle to core distribution. This +// layout is derived from this object and does all important +// bunch updates. It is the interface for AMReX and Ippl. +// +// In AMReX, the geometry, i.e. physical domain, is fixed +// during the whole computation. Particles leaving the domain +// would be deleted. In order to prevent this we map the particles +// onto the domain [-1, 1]^3. Furthermore, it makes sure +// that we have enougth grid points to represent the bunch +// when its charges are scattered on the grid for the self-field +// computation. +// +// The self-field computation and the particle-to-core update +// are performed in the particle mapped domain. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, 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/>. +// #ifndef BoxLibLayout_HPP #define BoxLibLayout_HPP @@ -128,8 +163,8 @@ void BoxLibLayout<T, Dim>::setDomainRatio(const std::vector<double>& ratio) { template<class T, unsigned Dim> -void BoxLibLayout<T, Dim>::update(IpplParticleBase< BoxLibLayout<T,Dim> >& PData, - const ParticleAttrib<char>* canSwap) +void BoxLibLayout<T, Dim>::update(IpplParticleBase< BoxLibLayout<T,Dim> >& /*PData*/, + const ParticleAttrib<char>* /*canSwap*/) { /* Exit since we need AmrParticleBase with grids and levels for particles for this layout * if IpplParticleBase is used something went wrong @@ -179,10 +214,8 @@ void BoxLibLayout<T, Dim>::update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData, std::multimap<unsigned, unsigned> p2n; //node ID, particle - int *msgsend = new int[N]; - std::fill(msgsend, msgsend+N, 0); - int *msgrecv = new int[N]; - std::fill(msgrecv, msgrecv+N, 0); + std::vector<int> msgsend(N); + std::vector<int> msgrecv(N); unsigned sent = 0; size_t lBegin = LocalNumPerLevel.begin(lev_min); @@ -321,8 +354,6 @@ void BoxLibLayout<T, Dim>::update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData, delete buffers[j]; } - delete[] msgsend; - delete[] msgrecv; delete format; // there is extra work to do if there are multipple nodes, to distribute diff --git a/src/Amr/BoxLibParticle.h b/src/Amr/BoxLibParticle.h index 421838f68..51b69e8b2 100644 --- a/src/Amr/BoxLibParticle.h +++ b/src/Amr/BoxLibParticle.h @@ -1,3 +1,26 @@ +// +// Class BoxLibParticle +// Particle class for AMReX. It works together with BoxLibLayout. +// The class does the scatter and gather operations of attributes +// to and from the grid. Ippl implements the same functionality in the +// attribute class. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, 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/>. +// #ifndef BOXLIB_PARTICLE_H #define BOXLIB_PARTICLE_H @@ -10,15 +33,7 @@ #include <AMReX_Geometry.H> #include <AMReX_RealBox.H> -/*! - * Particle class for AMReX. It works together with BoxLibLayout. - * - * Particle class that does the scatter and gather operations of attributes - * to and from grid. In contrast to Ippl where it is implemented in the - * attribute class, i.e. src/ippl/src/Particle/ParticleAttrib.h we do it in - * the particle class. This way Ippl does not need to be modified if another - * AMR framework is used. - */ + template<class PLayout> class BoxLibParticle : public virtual AmrParticleBase<PLayout> { diff --git a/src/Amr/BoxLibParticle.hpp b/src/Amr/BoxLibParticle.hpp index d5a5a312b..7b33099a6 100644 --- a/src/Amr/BoxLibParticle.hpp +++ b/src/Amr/BoxLibParticle.hpp @@ -1,3 +1,26 @@ +// +// Class BoxLibParticle +// Particle class for AMReX. It works together with BoxLibLayout. +// The class does the scatter and gather operations of attributes +// to and from the grid. Ippl implements the same functionality in the +// attribute class. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, 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/>. +// #ifndef BOXLIB_PARTICLE_HPP #define BOXLIB_PARTICLE_HPP @@ -59,7 +82,7 @@ void BoxLibParticle<PLayout>::scatter(ParticleAttrib<FT>& attrib, AmrScalarField template<class PLayout> template <class FT, unsigned Dim, class PT> void BoxLibParticle<PLayout>::scatter(ParticleAttrib<FT>& attrib, AmrField_t& f, - ParticleAttrib<Vektor<PT, Dim> >& pp, + ParticleAttrib<Vektor<PT, Dim> >& /*pp*/, const ParticleAttrib<int>& pbin, int bin, int level) { @@ -80,7 +103,7 @@ void BoxLibParticle<PLayout>::scatter(ParticleAttrib<FT>& attrib, AmrField_t& f, template<class PLayout> template<class FT, unsigned Dim, class PT> void BoxLibParticle<PLayout>::gather(ParticleAttrib<FT>& attrib, AmrVectorFieldContainer_t& f, - ParticleAttrib<Vektor<PT, Dim> >& pp, + ParticleAttrib<Vektor<PT, Dim> >& /*pp*/, int lbase, int lfine) { this->InterpolateFort(attrib, f, lbase, lfine); @@ -159,7 +182,7 @@ void BoxLibParticle<PLayout>::AssignCellDensitySingleLevelFort(ParticleAttrib<AT const ParticleAttrib<int>& pbin, int bin, int ncomp, - int particle_lvl_offset) const + int /*particle_lvl_offset*/) const { // BL_PROFILE("ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt>::AssignCellDensitySingleLevelFort()"); diff --git a/src/Classic/Algorithms/AmrPartBunch.cpp b/src/Classic/Algorithms/AmrPartBunch.cpp index b3630abff..68121145e 100644 --- a/src/Classic/Algorithms/AmrPartBunch.cpp +++ b/src/Classic/Algorithms/AmrPartBunch.cpp @@ -55,7 +55,7 @@ const AmrPartBunch::pbase_t *AmrPartBunch::getAmrParticleBase() const { } -void AmrPartBunch::initialize(FieldLayout_t *fLayout) { +void AmrPartBunch::initialize(FieldLayout_t */*fLayout*/) { // Layout_t* layout = static_cast<Layout_t*>(&getLayout()); } @@ -298,7 +298,7 @@ void AmrPartBunch::updateDomainLength(Vektor<int, 3>& grid) { } -void AmrPartBunch::updateFields(const Vector_t& hr, const Vector_t& origin) { +void AmrPartBunch::updateFields(const Vector_t& /*hr*/, const Vector_t& /*origin*/) { //TODO regrid; called in boundp() // amrobj_mp->updateMesh(); } diff --git a/src/Classic/Algorithms/AmrPartBunch.h b/src/Classic/Algorithms/AmrPartBunch.h index a280b959a..88845b507 100644 --- a/src/Classic/Algorithms/AmrPartBunch.h +++ b/src/Classic/Algorithms/AmrPartBunch.h @@ -72,8 +72,8 @@ public: this->amrobj_mp = fs->getAmrObject(); } - virtual void setBinCharge(int bin, double q) { }; - virtual void setBinCharge(int bin) { }; + virtual void setBinCharge(int /*bin*/, double /*q*/) { }; + virtual void setBinCharge(int /*bin*/) { }; /* * AmrPartBunch only diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp index f44a1e0ad..d2bb03b97 100644 --- a/src/Classic/Algorithms/PartBunchBase.hpp +++ b/src/Classic/Algorithms/PartBunchBase.hpp @@ -1911,7 +1911,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() { std::vector<double> loc_moments(4 * Dim + Dim * ( 2 * Dim + 1 )); long int totalNum = this->getTotalNum(); - if (OpalData::getInstance()->isInOPALCyclMode()) { + if (!Options::amr && OpalData::getInstance()->isInOPALCyclMode()) { /* FIXME After issue 287 is resolved this shouldn't be necessary * anymore */ diff --git a/src/Classic/Algorithms/Tracker.h b/src/Classic/Algorithms/Tracker.h index 52b5a5be5..4feb91235 100644 --- a/src/Classic/Algorithms/Tracker.h +++ b/src/Classic/Algorithms/Tracker.h @@ -131,10 +131,10 @@ public: virtual void visitMapIntegrator(const MapIntegrator &); /// set total number of tracked bunches - virtual void setNumBunch(int) {}; + virtual void setNumBunch(short) {}; /// get total number of tracked bunches - virtual int getNumBunch() { return 0; } + virtual short getNumBunch() { return 0; } // standing wave structures FieldList cavities_m; diff --git a/src/Distribution/ClosedOrbitFinder.h b/src/Distribution/ClosedOrbitFinder.h index b5201f659..2c7491539 100644 --- a/src/Distribution/ClosedOrbitFinder.h +++ b/src/Distribution/ClosedOrbitFinder.h @@ -464,6 +464,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc out << std::left << std::setw(15) << "# energy[MeV]" << std::setw(15) << "radius_ini[m]" + << std::setw(15) << "momentum_ini[Beta Gamma]" << std::setw(15) << "radius_avg[m]" << std::setw(15) << "nu_r" << std::setw(15) << "nu_z" @@ -555,6 +556,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc out << std::left << std::setw(15) << E << std::setw(15) << reo + << std::setw(15) << peo << std::setw(15) << ravg_m << std::setw(15) << tunes.first << std::setw(15) << tunes.second << std::endl; diff --git a/src/Optimize/OptimizeCmd.cpp b/src/Optimize/OptimizeCmd.cpp index a35f3cf9d..9810b6aed 100644 --- a/src/Optimize/OptimizeCmd.cpp +++ b/src/Optimize/OptimizeCmd.cpp @@ -37,6 +37,7 @@ #include "Expression/NumberOfPeaks.h" #include "Expression/SumErrSqRadialPeak.h" #include "Expression/ProbeVariable.h" +#include "Expression/SeptumExpr.h" #include <boost/filesystem.hpp> @@ -235,6 +236,10 @@ void OptimizeCmd::execute() { funcs.insert(std::pair<std::string, client::function::type> ("statVariableAt", ff)); + ff = SeptumExpr(); + funcs.insert(std::pair<std::string, client::function::type> + ("septum", ff)); + ////////////////////////////////////////////////////////////////////////// std::vector<std::string> arguments(opal->getArguments()); diff --git a/src/Sample/SampleCmd.cpp b/src/Sample/SampleCmd.cpp index cbcdb62e6..de55db880 100644 --- a/src/Sample/SampleCmd.cpp +++ b/src/Sample/SampleCmd.cpp @@ -54,6 +54,7 @@ #include "Expression/NumberOfPeaks.h" #include "Expression/SumErrSqRadialPeak.h" #include "Expression/ProbeVariable.h" +#include "Expression/SeptumExpr.h" #include <boost/filesystem.hpp> @@ -238,6 +239,10 @@ void SampleCmd::execute() { funcs.insert(std::pair<std::string, client::function::type> ("statVariableAt", ff)); + ff = SeptumExpr(); + funcs.insert(std::pair<std::string, client::function::type> + ("septum", ff)); + ////////////////////////////////////////////////////////////////////////// std::set<std::string> objExpressions; // check if all unique objective expressions diff --git a/src/Solvers/AMR_MG/AmrMultiGrid.cpp b/src/Solvers/AMR_MG/AmrMultiGrid.cpp index 46c86684c..8f6cb641a 100644 --- a/src/Solvers/AMR_MG/AmrMultiGrid.cpp +++ b/src/Solvers/AMR_MG/AmrMultiGrid.cpp @@ -196,6 +196,11 @@ void AmrMultiGrid::setVerbose(bool verbose) { } +void AmrMultiGrid::setTolerance(const scalar_t& eps) { + eps_m = eps; +} + + void AmrMultiGrid::initPhysicalBoundary_m(const Boundary* bc) { // make sure it's reset diff --git a/src/Solvers/AMR_MG/AmrMultiGrid.h b/src/Solvers/AMR_MG/AmrMultiGrid.h index 47c1c412c..2eaa39822 100644 --- a/src/Solvers/AMR_MG/AmrMultiGrid.h +++ b/src/Solvers/AMR_MG/AmrMultiGrid.h @@ -220,6 +220,12 @@ public: * Enable solver info dumping into SDDS file */ void setVerbose(bool verbose); + + + /*! + * Change accuracy of solver + */ + void setTolerance(const scalar_t& eps); double getXRangeMin(unsigned short level = 0); double getXRangeMax(unsigned short level = 0); @@ -695,7 +701,7 @@ private: Norm norm_m; ///< norm for convergence criteria (l1, l2, linf) std::string snorm_m; ///< norm for convergence criteria - const scalar_t eps_m; ///< rhs scale for convergence + scalar_t eps_m; ///< rhs scale for convergence bool verbose_m; ///< If true, a SDDS file is written std::string fname_m; ///< SDDS filename diff --git a/src/Solvers/AMReXSolvers/MLPoissonSolver.cpp b/src/Solvers/AMReXSolvers/MLPoissonSolver.cpp index f0e364190..592a2ca84 100644 --- a/src/Solvers/AMReXSolvers/MLPoissonSolver.cpp +++ b/src/Solvers/AMReXSolvers/MLPoissonSolver.cpp @@ -1,3 +1,23 @@ +// +// Class MLPoissonSolver +// Interface to the C++ based AMR Poisson multigrid solver of AMReX. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, 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/>. +// #include "MLPoissonSolver.h" #include "Utilities/OpalException.h" @@ -19,7 +39,7 @@ void MLPoissonSolver::solve(AmrScalarFieldContainer_t& rho, AmrVectorFieldContainer_t& efield, unsigned short baseLevel, unsigned short finestLevel, - bool prevAsGuess) + bool /*prevAsGuess*/) { for (int i = 0; i <= finestLevel; ++i) { phi[i]->setVal(0.0, 1); diff --git a/src/Solvers/AMReXSolvers/MLPoissonSolver.h b/src/Solvers/AMReXSolvers/MLPoissonSolver.h index 8f506f649..1a02d8412 100644 --- a/src/Solvers/AMReXSolvers/MLPoissonSolver.h +++ b/src/Solvers/AMReXSolvers/MLPoissonSolver.h @@ -1,3 +1,23 @@ +// +// Class MLPoissonSolver +// Interface to the C++ based AMR Poisson multigrid solver of AMReX. +// +// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, 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/>. +// #ifndef ML_POISSON_SOLVER_H_ #define ML_POISSON_SOLVER_H_ diff --git a/src/Solvers/AmrPoissonSolver.h b/src/Solvers/AmrPoissonSolver.h index dd6e02117..03bc784e3 100644 --- a/src/Solvers/AmrPoissonSolver.h +++ b/src/Solvers/AmrPoissonSolver.h @@ -39,17 +39,17 @@ public: virtual ~AmrPoissonSolver() {} - void computePotential(Field_t &rho, Vector_t hr) { + void computePotential(Field_t &/*rho*/, Vector_t /*hr*/) { throw OpalException("AmrPoissonSolver::computePotential(Field_t, Vector_t)", "Not implemented."); } - void computePotential(Field_t &rho, Vector_t hr, double zshift) { + void computePotential(Field_t &/*rho*/, Vector_t /*hr*/, double /*zshift*/) { throw OpalException("AmrPoissonSolver::computePotential(Field_t, Vector_t, double)", "Not implemented."); } - void test(PartBunchBase<double, 3> *bunch) { + void test(PartBunchBase<double, 3> */*bunch*/) { throw OpalException("AmrPoissonSolver::test(PartBunchBase<double, 3>)", "Not implemented."); } diff --git a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp index b7075392c..8e7cda8ac 100644 --- a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp +++ b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp @@ -1,3 +1,25 @@ +// +// Class FMGPoissonSolver +// Interface to the Fortran based AMR Poisson multigrid solver of BoxLib. +// The Fortran solver is part of AMReX till version 18.07 +// (Link: https://github.com/AMReX-Codes/amrex/archive/18.07.tar.gz) +// +// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, 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/>. +// #include "FMGPoissonSolver.h" #include "Utilities/OpalException.h" @@ -26,7 +48,7 @@ void FMGPoissonSolver::solve(AmrScalarFieldContainer_t& rho, AmrVectorFieldContainer_t& efield, unsigned short baseLevel, unsigned short finestLevel, - bool prevAsGuess) + bool /*prevAsGuess*/) { const GeomContainer_t& geom = itsAmrObject_mp->Geom(); diff --git a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.h b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.h index 544564237..a14435537 100644 --- a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.h +++ b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.h @@ -1,3 +1,25 @@ +// +// Class FMGPoissonSolver +// Interface to the Fortran based AMR Poisson multigrid solver of BoxLib. +// The Fortran solver is part of AMReX till version 18.07 +// (Link: https://github.com/AMReX-Codes/amrex/archive/18.07.tar.gz) +// +// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, 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/>. +// #ifndef FMG_POISSON_SOLVER_H_ #define FMG_POISSON_SOLVER_H_ diff --git a/src/Solvers/PoissonSolver.h b/src/Solvers/PoissonSolver.h index cfd04b632..e2dd5a350 100644 --- a/src/Solvers/PoissonSolver.h +++ b/src/Solvers/PoissonSolver.h @@ -29,12 +29,12 @@ public: * @param finestLevel for solve * @param prevAsGuess use of previous solution as initial guess */ - virtual void solve(AmrScalarFieldContainer_t &rho, - AmrScalarFieldContainer_t &phi, - AmrVectorFieldContainer_t &efield, - unsigned short baseLevel, - unsigned short finestLevel, - bool prevAsGuess = true) + virtual void solve(AmrScalarFieldContainer_t &/*rho*/, + AmrScalarFieldContainer_t &/*phi*/, + AmrVectorFieldContainer_t &/*efield*/, + unsigned short /*baseLevel*/, + unsigned short /*finestLevel*/, + bool /*prevAsGuess*/ = true) { throw OpalException("PoissonSolver::solve()", "Not supported for non-AMR code."); }; diff --git a/src/Structure/DataSink.cpp b/src/Structure/DataSink.cpp index f8ae774fd..73d3992ad 100644 --- a/src/Structure/DataSink.cpp +++ b/src/Structure/DataSink.cpp @@ -309,7 +309,13 @@ void DataSink::writeMultiBunchStatistics(PartBunchBase<double, 3> *beam, void DataSink::setMultiBunchInitialPathLengh(MultiBunchHandler* mbhandler_p) { for (short b = 0; b < mbhandler_p->getNumBunch(); ++b) { MultiBunchHandler::beaminfo_t& binfo = mbhandler_p->getBunchInfo(b); - binfo.pathlength = mbWriter_m[b]->getLastValue("s"); + if (mbWriter_m[b]->exists()) { + binfo.pathlength = mbWriter_m[b]->getLastValue("s"); + } else if (statWriter_m->exists()) { + binfo.pathlength = statWriter_m->getLastValue("s"); + } else { + binfo.pathlength = 0.0; + } } } diff --git a/src/Structure/FieldSolver.cpp b/src/Structure/FieldSolver.cpp index e89ce46cd..53a9058b2 100644 --- a/src/Structure/FieldSolver.cpp +++ b/src/Structure/FieldSolver.cpp @@ -119,6 +119,7 @@ namespace { AMR_MG_VERBOSE, // AMR, enable solver info writing (SDDS file) AMR_MG_REBALANCE, // AMR, rebalance smoothed aggregation (SA) preconditioner AMR_MG_REUSE, // AMR, reuse type of SA (none, RP, RAP, S or full) + AMR_MG_TOL, // AMR, tolerance of solver #endif // FOR XXX BASED SOLVER SIZE @@ -318,6 +319,10 @@ FieldSolver::FieldSolver(): itsAttr[AMR_MG_REUSE] = Attributes::makeUpperCaseString("AMR_MG_REUSE", "Reuse type of Smoothed Aggregation", "RAP"); + + itsAttr[AMR_MG_TOL] = Attributes::makeReal("AMR_MG_TOL", + "AMR MG solver tolerance (default: 1.0e-10)", + 1.0e-10); #endif mesh_m = 0; @@ -624,6 +629,8 @@ Inform &FieldSolver::printInfo(Inform &os) const { << Attributes::getString(itsAttr[AMR_MG_INTERP]) << '\n' << "* AMR_MG_NORM " << Attributes::getString(itsAttr[AMR_MG_NORM]) << '\n' + << "* AMR_MG_TOL " + << Attributes::getReal(itsAttr[AMR_MG_TOL]) << '\n' << "* AMR_MG_VERBOSE " << Attributes::getBool(itsAttr[AMR_MG_VERBOSE]) << '\n' << "* BCFFTX " @@ -777,6 +784,9 @@ void FieldSolver::initAmrSolver_m() { dynamic_cast<AmrMultiGrid*>(solver_m)->setVerbose( Attributes::getBool(itsAttr[AMR_MG_VERBOSE])); + + dynamic_cast<AmrMultiGrid*>(solver_m)->setTolerance( + Attributes::getReal(itsAttr[AMR_MG_TOL])); #else throw OpalException("FieldSolver::initAmrSolver_m()", "Multigrid solver not enabled! " diff --git a/src/Structure/GridLBalWriter.cpp b/src/Structure/GridLBalWriter.cpp index ec26fcbd2..baee5bbbd 100644 --- a/src/Structure/GridLBalWriter.cpp +++ b/src/Structure/GridLBalWriter.cpp @@ -39,7 +39,7 @@ void GridLBalWriter::fillHeader(const PartBunchBase<double, 3> *beam) { columns_m.addColumn("t", "double", "ns", "Time"); - AmrPartBunch* amrbeam = dynamic_cast<AmrPartBunch*>(beam); + const AmrPartBunch* amrbeam = dynamic_cast<const AmrPartBunch*>(beam); if ( !amrbeam ) throw OpalException("GridLBalWriter::fillHeader()", @@ -87,7 +87,7 @@ void GridLBalWriter::fillHeader(const PartBunchBase<double, 3> *beam) { } -void GridLBalWriter::write(const PartBunchBase<double, 3> *beam) { +void GridLBalWriter::write(PartBunchBase<double, 3> *beam) { AmrPartBunch* amrbeam = dynamic_cast<AmrPartBunch*>(beam); if ( !amrbeam ) diff --git a/src/Structure/GridLBalWriter.h b/src/Structure/GridLBalWriter.h index df0d7619f..1b8c53f34 100644 --- a/src/Structure/GridLBalWriter.h +++ b/src/Structure/GridLBalWriter.h @@ -28,7 +28,7 @@ class GridLBalWriter : public SDDSWriter { public: GridLBalWriter(const std::string& fname, bool restart); - void write(const PartBunchBase<double, 3> *beam) override; + void write(PartBunchBase<double, 3> *beam); private: void fillHeader(const PartBunchBase<double, 3> *beam); diff --git a/src/Structure/H5PartWrapperForPC.cpp b/src/Structure/H5PartWrapperForPC.cpp index 36e9311eb..09fc25938 100644 --- a/src/Structure/H5PartWrapperForPC.cpp +++ b/src/Structure/H5PartWrapperForPC.cpp @@ -475,12 +475,14 @@ void H5PartWrapperForPC::writeStepData(PartBunchBase<double, 3>* bunch) { size_t IDZero = bunch->getLocalNum(); bool found = false; +#ifndef ENABLE_AMR for(size_t k = 0; k < IDZero; ++ k) { if (bunch->ID[k] == 0) { found = true; IDZero = k; } } +#endif const size_t numLocalParticles = (found? bunch->getLocalNum() - 1: bunch->getLocalNum()); const size_t skipID = IDZero; diff --git a/src/Structure/LBalWriter.cpp b/src/Structure/LBalWriter.cpp index b405cbb45..dab467926 100644 --- a/src/Structure/LBalWriter.cpp +++ b/src/Structure/LBalWriter.cpp @@ -92,12 +92,14 @@ void LBalWriter::fillHeader() { } -void LBalWriter::write(const PartBunchBase<double, 3> *beam) { - #ifdef ENABLE_AMR +void LBalWriter::write(PartBunchBase<double, 3> *beam) { + if ( AmrPartBunch* amrbeam = dynamic_cast<AmrPartBunch*>(beam) ) { amrbeam->gatherLevelStatistics(); } +#else +void LBalWriter::write(const PartBunchBase<double, 3> *beam) { #endif if ( Ippl::myNode() != 0 ) diff --git a/src/Structure/LBalWriter.h b/src/Structure/LBalWriter.h index e9970d53c..ef769c519 100644 --- a/src/Structure/LBalWriter.h +++ b/src/Structure/LBalWriter.h @@ -26,7 +26,11 @@ class LBalWriter : public SDDSWriter { public: LBalWriter(const std::string& fname, bool restart); +#ifdef ENABLE_AMR + void write(PartBunchBase<double, 3> *beam); +#else void write(const PartBunchBase<double, 3> *beam) override; +#endif private: #ifdef ENABLE_AMR diff --git a/src/opal.cpp b/src/opal.cpp index d04e64f29..092401b6e 100644 --- a/src/opal.cpp +++ b/src/opal.cpp @@ -11,6 +11,7 @@ extern Inform *gmsg; #include "Utilities/OpalException.h" #include "Fields/Fieldmap.h" #include "Structure/IpplInfoWrapper.h" +#include "Utilities/Options.h" #include "OPALconfig.h" @@ -35,7 +36,8 @@ int run_opal(char */*args*/[], std::string inputfile, int restartStep, IpplInfo::Warn->setDestination(output); #ifdef ENABLE_AMR - amrex::Initialize(comm); + if ( Options::amr ) + amrex::Initialize(comm); #endif OpalData *opal = OpalData::getInstance(); @@ -77,7 +79,8 @@ int run_opal(char */*args*/[], std::string inputfile, int restartStep, delete gmsg; #ifdef ENABLE_AMR - amrex::Finalize(true); + if ( Options::amr ) + amrex::Finalize(true); #endif //FIXME: strange side effects -- GitLab