Commit fed401f6 authored by frey_m's avatar frey_m
Browse files

Merge branch '496-completion-of-the-amr-interface' into 'master'

Resolve "Completion of the AMR interface"

Part of #496

See merge request !312
parents ba764f33 d5bdff2c
......@@ -8,6 +8,7 @@ SET (_EXPR_SRCS
SumErrSqRadialPeak.cpp
Parser/expression.cpp
Parser/evaluator.cpp
SeptumExpr.h
)
add_sources(${_EXPR_SRCS})
......
//
// 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
......@@ -5,6 +5,7 @@ SET (_UTIL_SRCS
PeakReader.cpp
ProbeReader.cpp
ProbeHistReader.cpp
SDDSParser.cpp
SDDSParser/array.cpp
SDDSParser/associate.cpp
......
//
// 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);
}
//
// 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
//
// 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:
......
//
// 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
......