FieldSolver.cpp 32.9 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// ------------------------------------------------------------------------
// $RCSfile: FieldSolver.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.3.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: FieldSolver
//   The class for the OPAL FIELDSOLVER command.
//
// ------------------------------------------------------------------------
//
// $Date: 2003/08/11 22:09:00 $
// $Author: ADA $
//
// ------------------------------------------------------------------------

#include "Structure/FieldSolver.h"
20 21
#include "Solvers/FFTPoissonSolver.h"
#include "Solvers/FFTBoxPoissonSolver.h"
Andreas Adelmann's avatar
Andreas Adelmann committed
22
#include "Solvers/P3MPoissonSolver.h"
23
#ifdef HAVE_SAAMG_SOLVER
24 25
#include "Solvers/MGPoissonSolver.h"
#endif
gsell's avatar
gsell committed
26 27 28 29 30 31 32 33 34
#include "AbstractObjects/Expressions.h"
#include "AbstractObjects/OpalData.h"
#include "Attributes/Attributes.h"
#include "Expressions/SAutomatic.h"
#include "Expressions/SRefExpr.h"
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
#include "BoundaryGeometry.h"
#include "AbstractObjects/Element.h"
35
#include "Algorithms/PartBunchBase.h"
gsell's avatar
gsell committed
36

37
#ifdef ENABLE_AMR
frey_m's avatar
AMR:  
frey_m committed
38
    #include "Amr/AmrBoxLib.h"
39
#ifdef AMREX_ENABLE_FBASELIB
40
    #include "Solvers/BoxLibSolvers/FMGPoissonSolver.h"
41
#endif
42
    #include "Solvers/AMReXSolvers/MLPoissonSolver.h"
43
    #include "Amr/AmrDefs.h"
frey_m's avatar
frey_m committed
44
    #include "Algorithms/AmrPartBunch.h"
45 46 47 48

    #include <algorithm>
    #include <cctype>
    #include <functional>
49 50
#endif

51
#ifdef HAVE_AMR_MG_SOLVER
frey_m's avatar
frey_m committed
52
    #include "Solvers/AMR_MG/AmrMultiGrid.h"
53 54
#endif

gsell's avatar
gsell committed
55 56 57 58 59 60 61 62 63 64
using namespace Expressions;
using namespace Physics;

//TODO: o add a FIELD for DISCRETIZATION, MAXITERS, TOL...

// Class FieldSolver
// ------------------------------------------------------------------------

// The attributes of class FieldSolver.
namespace {
kraus's avatar
kraus committed
65 66 67 68 69 70
    namespace deprecated {
        enum {
            BCFFTT,
            SIZE
        };
    }
gsell's avatar
gsell committed
71
    enum {
kraus's avatar
kraus committed
72
        FSTYPE = deprecated::SIZE,   // The field solver name
gsell's avatar
gsell committed
73 74 75 76 77 78 79
        // FOR FFT BASED SOLVER
        MX,         // mesh sixe in x
        MY,         // mesh sixe in y
        MT,         //  mesh sixe in z
        PARFFTX,    // parallelized grind in x
        PARFFTY,    // parallelized grind in y
        PARFFTT,    // parallelized grind in z
frey_m's avatar
frey_m committed
80 81 82
        BCFFTX,     // boundary condition in x [FFT + AMR_MG only]
        BCFFTY,     // boundary condition in y [FFT + AMR_MG only]
        BCFFTZ,     // boundary condition in z [FFT + AMR_MG only]
gsell's avatar
gsell committed
83 84
        GREENSF,    // holds greensfunction to be used [FFT only]
        BBOXINCR,   // how much the boundingbox is increased
85
        GEOMETRY,   // geometry of boundary [SAAMG only]
frey_m's avatar
frey_m committed
86
        ITSOLVER,   // iterative solver [SAAMG + AMR_MG]
87 88 89 90
        INTERPL,    // interpolation used for boundary points [SAAMG only]
        TOL,        // tolerance of the SAAMG preconditioned solver [SAAMG only]
        MAXITERS,   // max number of iterations [SAAMG only]
        PRECMODE,   // preconditioner mode [SAAMG only]
Andreas Adelmann's avatar
Andreas Adelmann committed
91
        RC,         // cutoff radius for PP interactions
kraus's avatar
kraus committed
92 93
        ALPHA,      // Green’s function splitting parameter
        EPSILON,    // regularization for PP interaction
94
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
95 96 97
        AMR_MAXLEVEL,       // AMR, maximum refinement level
        AMR_REFX,           // AMR, refinement ratio in x
        AMR_REFY,           // AMR, refinement ratio in y
frey_m's avatar
frey_m committed
98 99 100 101 102 103 104
        AMR_REFZ,           // AMR, refinement ration in z
        AMR_MAXGRIDX,       // AMR, maximum grid size in x (default: 16)
        AMR_MAXGRIDY,       // AMR, maximum grid size in y (default: 16)
        AMR_MAXGRIDZ,       // AMR, maximum grid size in z (default: 16)
        AMR_BFX,            // AMR, blocking factor in x (maxgrid needs to be a multiple, default: 8)
        AMR_BFY,            // AMR, blocking factor in y (maxgrid needs to be a multiple, default: 8)
        AMR_BFZ,            // AMR, blocking factor in z (maxgrid needs to be a multiple, default: 8)
frey_m's avatar
frey_m committed
105 106 107 108 109
        AMR_TAGGING,
        AMR_DENSITY,
        AMR_MAX_NUM_PART,
        AMR_MIN_NUM_PART,
        AMR_SCALING,
110
        AMR_DOMAIN_RATIO,
111 112
#endif
#ifdef HAVE_AMR_MG_SOLVER
frey_m's avatar
frey_m committed
113 114 115 116 117 118
        // AMR_MG = Adaptive-Mesh-Refinement Multi-Grid
        AMR_MG_SMOOTHER,    // AMR, smoother for level solution
        AMR_MG_NSWEEPS,     // AMR, number of smoothing sweeps
        AMR_MG_PREC,        // AMR, preconditioner for bottom solver
        AMR_MG_INTERP,      // AMR, interpolater for solution from level l to l+1
        AMR_MG_NORM,        // AMR, norm convergence criteria
119
        AMR_MG_VERBOSE,     // AMR, enable solver info writing (SDDS file)
120
        AMR_MG_REBALANCE,   // AMR, rebalance smoothed aggregation (SA) preconditioner
frey_m's avatar
frey_m committed
121
        AMR_MG_REUSE,       // AMR, reuse type of SA (none, RP, RAP, S or full)
frey_m's avatar
frey_m committed
122
        AMR_MG_TOL,         // AMR, tolerance of solver
123
#endif
gsell's avatar
gsell committed
124 125 126 127 128 129 130 131 132 133
        // FOR XXX BASED SOLVER
        SIZE
    };
}


FieldSolver::FieldSolver():
    Definition(SIZE, "FIELDSOLVER",
               "The \"FIELDSOLVER\" statement defines data for a the field solver ") {

134 135 136 137 138 139
    itsAttr[FSTYPE] = Attributes::makeUpperCaseString("FSTYPE", "Name of the attached field solver: "
                                                      "FFT, "
                                                      "FFTPERIODIC, "
                                                      "SAAMG, "
                                                      "AMR, "
                                                      "NONE ");
gsell's avatar
gsell committed
140 141 142 143 144

    itsAttr[MX] = Attributes::makeReal("MX", "Meshsize in x");
    itsAttr[MY] = Attributes::makeReal("MY", "Meshsize in y");
    itsAttr[MT] = Attributes::makeReal("MT", "Meshsize in z(t)");

145 146 147
    itsAttr[PARFFTX] = Attributes::makeBool("PARFFTX",
                                            "True, dimension 0 i.e x is parallelized",
                                            false);
kraus's avatar
kraus committed
148

149 150 151
    itsAttr[PARFFTY] = Attributes::makeBool("PARFFTY",
                                            "True, dimension 1 i.e y is parallelized",
                                            false);
kraus's avatar
kraus committed
152

153 154 155
    itsAttr[PARFFTT] = Attributes::makeBool("PARFFTT",
                                            "True, dimension 2 i.e z(t) is parallelized",
                                            true);
gsell's avatar
gsell committed
156 157

    //FFT ONLY:
158
    itsAttr[BCFFTX] = Attributes::makeUpperCaseString("BCFFTX",
frey_m's avatar
frey_m committed
159 160
                                                      "Boundary conditions in x: open, dirichlet (box), periodic",
                                                      "OPEN");
kraus's avatar
kraus committed
161

162
    itsAttr[BCFFTY] = Attributes::makeUpperCaseString("BCFFTY",
frey_m's avatar
frey_m committed
163 164
                                                      "Boundary conditions in y: open, dirichlet (box), periodic",
                                                      "OPEN");
kraus's avatar
kraus committed
165

166
    itsAttr[BCFFTZ] = Attributes::makeUpperCaseString("BCFFTZ",
frey_m's avatar
frey_m committed
167 168
                                                      "Boundary conditions in z(t): open, periodic",
                                                      "OPEN");
kraus's avatar
kraus committed
169

170
    itsAttr[deprecated::BCFFTT] = Attributes::makeUpperCaseString("BCFFTT",
frey_m's avatar
frey_m committed
171 172
                                                                  "Boundary conditions in z(t): open, periodic",
                                                                  "OPEN");
gsell's avatar
gsell committed
173

174
    itsAttr[GREENSF]  = Attributes::makeUpperCaseString("GREENSF",
frey_m's avatar
frey_m committed
175 176
                                                        "Which Greensfunction to be used [STANDARD | INTEGRATED]",
                                                        "INTEGRATED");
kraus's avatar
kraus committed
177

178 179 180
    itsAttr[BBOXINCR] = Attributes::makeReal("BBOXINCR",
                                             "Increase of bounding box in % ",
                                             2.0);
gsell's avatar
gsell committed
181 182

    // P3M only:
183 184 185
    itsAttr[RC] = Attributes::makeReal("RC",
                                       "cutoff radius for PP interactions",
                                       0.0);
kraus's avatar
kraus committed
186

187 188 189
    itsAttr[ALPHA] = Attributes::makeReal("ALPHA",
                                          "Green’s function splitting parameter",
                                          0.0);
kraus's avatar
kraus committed
190

191 192 193
    itsAttr[EPSILON] = Attributes::makeReal("EPSILON",
                                            "regularization for PP interaction",
                                            0.0);
gsell's avatar
gsell committed
194

195
    //SAAMG and in case of FFT with dirichlet BC in x and y
196
    itsAttr[GEOMETRY] = Attributes::makeUpperCaseString("GEOMETRY",
frey_m's avatar
frey_m committed
197 198
                                                        "GEOMETRY to be used as domain boundary",
                                                        "");
kraus's avatar
kraus committed
199

200
    itsAttr[ITSOLVER] = Attributes::makeUpperCaseString("ITSOLVER",
frey_m's avatar
frey_m committed
201 202
                                                        "Type of iterative solver [CG | BiCGSTAB | GMRES]",
                                                        "CG");
kraus's avatar
kraus committed
203

204
    itsAttr[INTERPL] = Attributes::makeUpperCaseString("INTERPL",
frey_m's avatar
frey_m committed
205 206
                                                        "interpolation used for boundary points [CONSTANT | LINEAR | QUADRATIC]",
                                                        "LINEAR");
kraus's avatar
kraus committed
207

208 209 210
    itsAttr[TOL] = Attributes::makeReal("TOL",
                                        "Tolerance for iterative solver",
                                        1e-8);
kraus's avatar
kraus committed
211

212 213 214
    itsAttr[MAXITERS] = Attributes::makeReal("MAXITERS",
                                             "Maximum number of iterations of iterative solver",
                                             100);
kraus's avatar
kraus committed
215

216
    itsAttr[PRECMODE] = Attributes::makeUpperCaseString("PRECMODE",
frey_m's avatar
frey_m committed
217 218
                                                        "Preconditioner Mode [STD | HIERARCHY | REUSE]",
                                                        "HIERARCHY");
219

220
    // AMR
221
#ifdef ENABLE_AMR
222 223 224
    itsAttr[AMR_MAXLEVEL] = Attributes::makeReal("AMR_MAXLEVEL",
                                                 "Maximum number of levels in AMR",
                                                 0);
kraus's avatar
kraus committed
225

226 227 228
    itsAttr[AMR_REFX] = Attributes::makeReal("AMR_REFX",
                                             "Refinement ration in x-direction in AMR",
                                             2);
kraus's avatar
kraus committed
229

230 231 232
    itsAttr[AMR_REFY] = Attributes::makeReal("AMR_REFY",
                                             "Refinement ration in y-direction in AMR",
                                             2);
kraus's avatar
kraus committed
233

frey_m's avatar
frey_m committed
234
    itsAttr[AMR_REFZ] = Attributes::makeReal("AMR_REFZ",
235 236
                                             "Refinement ration in z-direction in AMR",
                                             2);
kraus's avatar
kraus committed
237

frey_m's avatar
frey_m committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
    itsAttr[AMR_MAXGRIDX] = Attributes::makeReal("AMR_MAXGRIDX",
                                                 "Maximum grid size in x for AMR",
                                                 16);

    itsAttr[AMR_MAXGRIDY] = Attributes::makeReal("AMR_MAXGRIDY",
                                                 "Maximum grid size in y for AMR",
                                                 16);

    itsAttr[AMR_MAXGRIDZ] = Attributes::makeReal("AMR_MAXGRIDZ",
                                                 "Maximum grid size in z for AMR",
                                                 16);

    itsAttr[AMR_BFX] = Attributes::makeReal("AMR_BFX",
                                            "Blocking factor in x for AMR (AMR_MAXGRIDX needs to be a multiple",
                                            8);

    itsAttr[AMR_BFY] = Attributes::makeReal("AMR_BFY",
                                            "Blocking factor in y for AMR (AMR_MAXGRIDY needs to be a multiple",
                                            8);

    itsAttr[AMR_BFZ] = Attributes::makeReal("AMR_BFZ",
                                            "Blocking factor in y for AMR (AMR_MAXGRIDZ needs to be a multiple",
                                            8);
kraus's avatar
kraus committed
261

262
    itsAttr[AMR_TAGGING] = Attributes::makeUpperCaseString("AMR_TAGGING",
frey_m's avatar
frey_m committed
263 264
                                                           "Refinement criteria [CHARGE_DENSITY | POTENTIAL | EFIELD]",
                                                           "CHARGE_DENSITY");
kraus's avatar
kraus committed
265

frey_m's avatar
frey_m committed
266 267 268
    itsAttr[AMR_DENSITY] = Attributes::makeReal("AMR_DENSITY",
                                               "Tagging value for charge density refinement [C / cell volume]",
                                               1.0e-14);
kraus's avatar
kraus committed
269

frey_m's avatar
frey_m committed
270
    itsAttr[AMR_MAX_NUM_PART] = Attributes::makeReal("AMR_MAX_NUM_PART",
271 272
                                                     "Tagging value for max. #particles",
                                                     1);
kraus's avatar
kraus committed
273

frey_m's avatar
frey_m committed
274
    itsAttr[AMR_MIN_NUM_PART] = Attributes::makeReal("AMR_MIN_NUM_PART",
275 276
                                                     "Tagging value for min. #particles",
                                                     1);
kraus's avatar
kraus committed
277

frey_m's avatar
frey_m committed
278 279 280 281
    itsAttr[AMR_SCALING] = Attributes::makeReal("AMR_SCALING",
                                                "Scaling value for maximum value tagging "
                                                "(only POTENTIAL / CHARGE_DENSITY / "
                                                "MOMENTA", 0.75);
282 283 284 285 286 287

    itsAttr[AMR_DOMAIN_RATIO] = Attributes::makeRealArray("AMR_DOMAIN_RATIO",
                                                         "Box ratio of AMR computation domain. Default: [-1, 1]^3");

    // default
    Attributes::setRealArray(itsAttr[AMR_DOMAIN_RATIO], {1.0, 1.0, 1.0});
288
#endif
kraus's avatar
kraus committed
289

290
#ifdef HAVE_AMR_MG_SOLVER
291
    itsAttr[AMR_MG_SMOOTHER] = Attributes::makeUpperCaseString("AMR_MG_SMOOTHER",
frey_m's avatar
frey_m committed
292
                                                               "Smoothing of level solution", "GS");
kraus's avatar
kraus committed
293

frey_m's avatar
frey_m committed
294 295 296
    itsAttr[AMR_MG_NSWEEPS] = Attributes::makeReal("AMR_MG_NSWEEPS",
                                                   "Number of relaxation steps",
                                                   8);
kraus's avatar
kraus committed
297

298
    itsAttr[AMR_MG_PREC] = Attributes::makeUpperCaseString("AMR_MG_PREC",
frey_m's avatar
frey_m committed
299 300
                                                           "Preconditioner of bottom solver",
                                                           "NONE");
kraus's avatar
kraus committed
301

302
    itsAttr[AMR_MG_INTERP] = Attributes::makeUpperCaseString("AMR_MG_INTERP",
frey_m's avatar
frey_m committed
303 304
                                                             "Interpolater between levels",
                                                             "PC");
kraus's avatar
kraus committed
305

306
    itsAttr[AMR_MG_NORM] = Attributes::makeUpperCaseString("AMR_MG_NORM",
frey_m's avatar
frey_m committed
307 308
                                                           "Norm for convergence criteria",
                                                           "LINF");
309

310 311 312
    itsAttr[AMR_MG_VERBOSE] = Attributes::makeBool("AMR_MG_VERBOSE",
                                                   "Write solver info in SDDS format (*.solver)",
                                                   false);
313 314 315 316 317

    itsAttr[AMR_MG_REBALANCE] = Attributes::makeBool("AMR_MG_REBALANCE",
                                                     "Rebalancing of Smoothed Aggregation "
                                                     "Preconditioner",
                                                     false);
318

319
    itsAttr[AMR_MG_REUSE] = Attributes::makeUpperCaseString("AMR_MG_REUSE",
frey_m's avatar
frey_m committed
320 321
                                                            "Reuse type of Smoothed Aggregation",
                                                            "RAP");
frey_m's avatar
frey_m committed
322 323 324 325

    itsAttr[AMR_MG_TOL] = Attributes::makeReal("AMR_MG_TOL",
                                               "AMR MG solver tolerance (default: 1.0e-10)",
                                               1.0e-10);
326
#endif
327

gsell's avatar
gsell committed
328 329
    mesh_m = 0;
    FL_m = 0;
frey_m's avatar
frey_m committed
330
    PL_m.reset(nullptr);
331

332
    solver_m = 0;
333 334

    registerOwnership(AttributeHandler::STATEMENT);
gsell's avatar
gsell committed
335 336 337
}


338
FieldSolver::FieldSolver(const std::string &name, FieldSolver *parent):
gsell's avatar
gsell committed
339
    Definition(name, parent)
340 341 342
{
    mesh_m = 0;
    FL_m = 0;
frey_m's avatar
frey_m committed
343
    PL_m.reset(nullptr);
344 345
    solver_m = 0;
}
gsell's avatar
gsell committed
346 347 348


FieldSolver::~FieldSolver() {
349 350 351 352 353 354 355 356 357
    if (mesh_m) {
        delete mesh_m;
        mesh_m = 0;
    }
    if (FL_m) {
        delete FL_m;
        FL_m = 0;
    }
    if (solver_m) {
358
       delete solver_m;
359 360
       solver_m = 0;
    }
gsell's avatar
gsell committed
361 362
}

363
FieldSolver *FieldSolver::clone(const std::string &name) {
gsell's avatar
gsell committed
364 365 366 367 368 369 370
    return new FieldSolver(name, this);
}

void FieldSolver::execute() {
    update();
}

371
FieldSolver *FieldSolver::find(const std::string &name) {
gsell's avatar
gsell committed
372 373 374 375 376 377 378 379
    FieldSolver *fs = dynamic_cast<FieldSolver *>(OpalData::getInstance()->find(name));

    if(fs == 0) {
        throw OpalException("FieldSolver::find()", "FieldSolver \"" + name + "\" not found.");
    }
    return fs;
}

380
std::string FieldSolver::getType() {
381
    return Attributes::getString(itsAttr[FSTYPE]);
382
}
gsell's avatar
gsell committed
383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427

double FieldSolver::getMX() const {
    return Attributes::getReal(itsAttr[MX]);
}

double FieldSolver::getMY() const {
    return Attributes::getReal(itsAttr[MY]);
}

double FieldSolver::getMT() const {
    return Attributes::getReal(itsAttr[MT]);
}

void FieldSolver::setMX(double value) {
    Attributes::setReal(itsAttr[MX], value);
}

void FieldSolver::setMY(double value) {
    Attributes::setReal(itsAttr[MY], value);
}

void FieldSolver::setMT(double value) {
    Attributes::setReal(itsAttr[MT], value);
}

void FieldSolver::update() {

}

void FieldSolver::initCartesianFields() {

    e_dim_tag decomp[3] = {SERIAL, SERIAL, SERIAL};

    NDIndex<3> domain;
    domain[0] = Index((int)getMX() + 1);
    domain[1] = Index((int)getMY() + 1);
    domain[2] = Index((int)getMT() + 1);

    if(Attributes::getBool(itsAttr[PARFFTX]))
        decomp[0] = PARALLEL;
    if(Attributes::getBool(itsAttr[PARFFTY]))
        decomp[1] = PARALLEL;
    if(Attributes::getBool(itsAttr[PARFFTT]))
        decomp[2] = PARALLEL;

428
    if(Attributes::getString(itsAttr[FSTYPE]) == "FFTPERIODIC") {
gsell's avatar
gsell committed
429 430 431 432
        decomp[0] = decomp[1] = SERIAL;
        decomp[2] = PARALLEL;
    }
    // create prototype mesh and layout objects for this problem domain
433

434 435 436 437 438 439 440 441
    if ( !isAmrSolverType() ) {
        mesh_m   = new Mesh_t(domain);
        FL_m     = new FieldLayout_t(*mesh_m, decomp);
        PL_m.reset(new Layout_t(*FL_m, *mesh_m));
        // OpalData::getInstance()->setMesh(mesh_m);
        // OpalData::getInstance()->setFieldLayout(FL_m);
        // OpalData::getInstance()->setLayout(PL_m);
    }
gsell's avatar
gsell committed
442 443
}

444
bool FieldSolver::hasPeriodicZ() {
kraus's avatar
kraus committed
445
    if (itsAttr[deprecated::BCFFTT])
446
        return (Attributes::getString(itsAttr[deprecated::BCFFTT]) == "PERIODIC");
kraus's avatar
kraus committed
447

448
    return (Attributes::getString(itsAttr[BCFFTZ]) == "PERIODIC");
adelmann's avatar
adelmann committed
449 450
}

frey_m's avatar
AMR:  
frey_m committed
451
inline bool FieldSolver::isAmrSolverType() const {
frey_m's avatar
frey_m committed
452
    return Options::amr;
453 454
}

frey_m's avatar
frey_m committed
455 456
void FieldSolver::initSolver(PartBunchBase<double, 3> *b) {
    itsBunch_m = b;
457 458 459 460 461
    fsType_m = Attributes::getString(itsAttr[FSTYPE]);
    std::string greens = Attributes::getString(itsAttr[GREENSF]);
    std::string bcx = Attributes::getString(itsAttr[BCFFTX]);
    std::string bcy = Attributes::getString(itsAttr[BCFFTY]);
    std::string bcz = Attributes::getString(itsAttr[deprecated::BCFFTT]);
kraus's avatar
kraus committed
462
    if (bcz == "") {
463
        bcz = Attributes::getString(itsAttr[BCFFTZ]);
kraus's avatar
kraus committed
464 465
    }

frey_m's avatar
frey_m committed
466 467
    if ( isAmrSolverType() ) {
        Inform m("FieldSolver::initAmrSolver");
468
        fsType_m = "AMR";
kraus's avatar
kraus committed
469

470
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
471
        initAmrObject_m();
kraus's avatar
kraus committed
472

frey_m's avatar
frey_m committed
473 474
        initAmrSolver_m();
#endif
475
    } else if(fsType_m == "FFT") {
476
        bool sinTrafo = ((bcx == "DIRICHLET") && (bcy == "DIRICHLET") && (bcz == "DIRICHLET"));
frey_m's avatar
frey_m committed
477 478 479
        if(sinTrafo) {
            std::cout << "FFTBOX ACTIVE" << std::endl;
            //we go over all geometries and add the Geometry Elements to the geometry list
480
            std::string geoms = Attributes::getString(itsAttr[GEOMETRY]);
frey_m's avatar
frey_m committed
481 482 483 484 485 486 487 488 489 490 491 492 493
            std::string tmp = "";
            //split and add all to list
            std::vector<BoundaryGeometry *> geometries;
            for(unsigned int i = 0; i <= geoms.length(); i++) {
                if(i == geoms.length() || geoms[i] == ',') {
                    BoundaryGeometry *geom = BoundaryGeometry::find(tmp);
                    if(geom != 0)
                        geometries.push_back(geom);
                    tmp.clear();
                } else
                    tmp += geoms[i];
            }
            BoundaryGeometry *ttmp = geometries[0];
494
            solver_m = new FFTBoxPoissonSolver(mesh_m, FL_m, greens, ttmp->getA());
frey_m's avatar
frey_m committed
495 496 497
            itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
            fsType_m = "FFTBOX";
        } else {
498
            solver_m = new FFTPoissonSolver(mesh_m, FL_m, greens, bcz);
frey_m's avatar
frey_m committed
499 500
            itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
        }
501
    } else if (fsType_m == "P3M") {
frey_m's avatar
frey_m committed
502 503 504 505 506 507 508 509 510
        solver_m = new P3MPoissonSolver(mesh_m,
                                        FL_m,
                                        Attributes::getReal(itsAttr[RC]),
                                        Attributes::getReal(itsAttr[ALPHA]),
                                        Attributes::getReal(itsAttr[EPSILON]));

        PL_m->setAllCacheDimensions(Attributes::getReal(itsAttr[RC]));
        PL_m->enableCaching();

511
    } else if(fsType_m == "SAAMG") {
frey_m's avatar
frey_m committed
512 513
#ifdef HAVE_SAAMG_SOLVER
        //we go over all geometries and add the Geometry Elements to the geometry list
514
        std::string geoms = Attributes::getString(itsAttr[GEOMETRY]);
frey_m's avatar
frey_m committed
515 516 517 518 519 520 521 522 523 524 525 526 527
        std::string tmp = "";
        //split and add all to list
        std::vector<BoundaryGeometry *> geometries;
        for(unsigned int i = 0; i <= geoms.length(); i++) {
            if(i == geoms.length() || geoms[i] == ',') {
                BoundaryGeometry *geom = OpalData::getInstance()->getGlobalGeometry();
                if(geom != 0) {
                    geometries.push_back(geom);
                }
                tmp.clear();
            } else
            tmp += geoms[i];
        }
528
        solver_m = new MGPoissonSolver(dynamic_cast<PartBunch*>(itsBunch_m), mesh_m, FL_m,
529
                                       geometries,
530 531
                                       Attributes::getString(itsAttr[ITSOLVER]),
                                       Attributes::getString(itsAttr[INTERPL]),
532 533
                                       Attributes::getReal(itsAttr[TOL]),
                                       Attributes::getReal(itsAttr[MAXITERS]),
534
                                       Attributes::getString(itsAttr[PRECMODE]));
frey_m's avatar
frey_m committed
535 536
        itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
#else
537 538
        throw OpalException("FieldSolver::initSolver",
                            "SAAMG Solver not enabled! Please build OPAL with -DENABLE_SAAMG_SOLVER=1");
539
#endif
frey_m's avatar
frey_m committed
540
    } else {
541 542
        solver_m = 0;
        INFOMSG("no solver attached" << endl);
543
    }
gsell's avatar
gsell committed
544 545 546 547 548 549 550
}

bool FieldSolver::hasValidSolver() {
    return (solver_m != 0);
}

Inform &FieldSolver::printInfo(Inform &os) const {
551
    std::string fsType = Attributes::getString(itsAttr[FSTYPE]);
552

gsell's avatar
gsell committed
553 554
    os << "* ************* F I E L D S O L V E R ********************************************** " << endl;
    os << "* FIELDSOLVER  " << getOpalName() << '\n'
555
       << "* TYPE         " << fsType << '\n'
gsell's avatar
gsell committed
556 557 558 559 560
       << "* N-PROCESSORS " << Ippl::getNodes() << '\n'
       << "* MX           " << Attributes::getReal(itsAttr[MX])   << '\n'
       << "* MY           " << Attributes::getReal(itsAttr[MY])   << '\n'
       << "* MT           " << Attributes::getReal(itsAttr[MT])   << '\n'
       << "* BBOXINCR     " << Attributes::getReal(itsAttr[BBOXINCR]) << endl;
561

562
    if (fsType == "P3M")
563 564 565
        os << "* RC           " << Attributes::getReal(itsAttr[RC]) << '\n'
           << "* ALPHA        " << Attributes::getReal(itsAttr[ALPHA]) << '\n'
           << "* EPSILON      " << Attributes::getReal(itsAttr[EPSILON]) << endl;
566

567

568
    if (fsType == "FFT") {
569
        os << "* GRRENSF      " << Attributes::getString(itsAttr[GREENSF]) << endl;
570
    } else if (fsType == "SAAMG") {
gsell's avatar
gsell committed
571
        os << "* GEOMETRY     " << Attributes::getString(itsAttr[GEOMETRY]) << '\n'
572 573
           << "* ITSOLVER     " << Attributes::getString(itsAttr[ITSOLVER])   << '\n'
           << "* INTERPL      " << Attributes::getString(itsAttr[INTERPL])  << '\n'
gsell's avatar
gsell committed
574 575
           << "* TOL          " << Attributes::getReal(itsAttr[TOL])        << '\n'
           << "* MAXITERS     " << Attributes::getReal(itsAttr[MAXITERS]) << '\n'
576
           << "* PRECMODE     " << Attributes::getString(itsAttr[PRECMODE])   << endl;
577
    } else if (Options::amr) {
578
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
579 580 581
        os << "* AMR_MAXLEVEL     " << Attributes::getReal(itsAttr[AMR_MAXLEVEL]) << '\n'
           << "* AMR_REFX         " << Attributes::getReal(itsAttr[AMR_REFX]) << '\n'
           << "* AMR_REFY         " << Attributes::getReal(itsAttr[AMR_REFY]) << '\n'
frey_m's avatar
frey_m committed
582 583 584 585 586 587 588
           << "* AMR_REFZ         " << Attributes::getReal(itsAttr[AMR_REFZ]) << '\n'
           << "* AMR_MAXGRIDX     " << Attributes::getReal(itsAttr[AMR_MAXGRIDX]) << '\n'
           << "* AMR_MAXGRIDY     " << Attributes::getReal(itsAttr[AMR_MAXGRIDY]) << '\n'
           << "* AMR_MAXGRIDZ     " << Attributes::getReal(itsAttr[AMR_MAXGRIDZ]) << '\n'
           << "* AMR_BFX          " << Attributes::getReal(itsAttr[AMR_BFX]) << '\n'
           << "* AMR_BFY          " << Attributes::getReal(itsAttr[AMR_BFY]) << '\n'
           << "* AMR_BFZ          " << Attributes::getReal(itsAttr[AMR_BFZ]) << '\n'
589
           << "* AMR_TAGGING      " << this->getTagging_m() <<'\n'
frey_m's avatar
frey_m committed
590 591 592 593
           << "* AMR_DENSITY      " << Attributes::getReal(itsAttr[AMR_DENSITY]) << '\n'
           << "* AMR_MAX_NUM_PART " << Attributes::getReal(itsAttr[AMR_MAX_NUM_PART]) << '\n'
           << "* AMR_MIN_NUM_PART " << Attributes::getReal(itsAttr[AMR_MIN_NUM_PART]) << '\n'
           << "* AMR_DENSITY      " << Attributes::getReal(itsAttr[AMR_DENSITY]) << '\n'
frey_m's avatar
frey_m committed
594
           << "* AMR_SCALING      " << Attributes::getReal(itsAttr[AMR_SCALING]) << endl;
595 596 597 598 599 600
        auto length = Attributes::getRealArray(itsAttr[AMR_DOMAIN_RATIO]);
        os << "* AMR_DOMAIN_RATIO ( ";
        for (auto& l : length) {
            os << l << " ";
        }
        os << ")" << endl;
601
#endif
602
    }
603

604
#ifdef HAVE_AMR_MG_SOLVER
frey_m's avatar
frey_m committed
605
    if (fsType == "AMR_MG") {
606
        os << "* ITSOLVER (AMR_MG)    "
607
           << Attributes::getString(itsAttr[ITSOLVER]) << '\n'
608
           << "* AMR_MG_PREC          "
609
           << Attributes::getString(itsAttr[AMR_MG_PREC]) << '\n'
610
           << "* AMR_MG_REBALANCE     "
frey_m's avatar
frey_m committed
611
           << Attributes::getBool(itsAttr[AMR_MG_REBALANCE]) << '\n'
frey_m's avatar
frey_m committed
612
           << "* AMR_MG_REUSE         "
613
           << Attributes::getString(itsAttr[AMR_MG_REUSE]) << '\n'
614
           << "* AMR_MG_SMOOTHER      "
615
           << Attributes::getString(itsAttr[AMR_MG_SMOOTHER]) << '\n'
616
           << "* AMR_MG_NSWEEPS       "
frey_m's avatar
frey_m committed
617
           << Attributes::getReal(itsAttr[AMR_MG_NSWEEPS]) << '\n'
618
           << "* AMR_MG_INTERP        "
619
           << Attributes::getString(itsAttr[AMR_MG_INTERP]) << '\n'
620
           << "* AMR_MG_NORM          "
621
           << Attributes::getString(itsAttr[AMR_MG_NORM]) << '\n'
frey_m's avatar
frey_m committed
622 623
           << "* AMR_MG_TOL           "
           << Attributes::getReal(itsAttr[AMR_MG_TOL]) << '\n'
624 625 626
           << "* AMR_MG_VERBOSE       "
           << Attributes::getBool(itsAttr[AMR_MG_VERBOSE]) << '\n'
           << "* BCFFTX               "
627
           << Attributes::getString(itsAttr[BCFFTX]) << '\n'
628
           << "* BCFFTY               "
629
           << Attributes::getString(itsAttr[BCFFTY]) << '\n';
kraus's avatar
kraus committed
630
        if (itsAttr[deprecated::BCFFTT]) {
631
            os << "* BCFFTT (deprec.) "
632
               << Attributes::getString(itsAttr[deprecated::BCFFTT]) << endl;
kraus's avatar
kraus committed
633 634
        } else {
            os << "* BCFFTZ           "
635
               << Attributes::getString(itsAttr[BCFFTZ]) << endl;
kraus's avatar
kraus committed
636
        }
637 638 639
    }
#endif

gsell's avatar
gsell committed
640
    if(Attributes::getBool(itsAttr[PARFFTX]))
641
        os << "* XDIM         parallel  " << endl;
gsell's avatar
gsell committed
642
    else
643
        os << "* XDIM         serial  " << endl;
gsell's avatar
gsell committed
644 645

    if(Attributes::getBool(itsAttr[PARFFTY]))
646
        os << "* YDIM         parallel  " << endl;
gsell's avatar
gsell committed
647
    else
648
        os << "* YDIM         serial  " << endl;
gsell's avatar
gsell committed
649 650

    if(Attributes::getBool(itsAttr[PARFFTT]))
651
        os << "* Z(T)DIM      parallel  " << endl;
gsell's avatar
gsell committed
652
    else
653
        os << "* Z(T)DIM      serial  " << endl;
gsell's avatar
gsell committed
654

655 656 657 658
    if ( !isAmrSolverType() ) {
        INFOMSG(level3 << *mesh_m << endl);
        INFOMSG(level3 << *PL_m << endl);
    }
659

gsell's avatar
gsell committed
660 661 662 663
    if(solver_m)
        os << *solver_m << endl;
    os << "* ********************************************************************************** " << endl;
    return os;
664
}
frey_m's avatar
frey_m committed
665

666
#ifdef ENABLE_AMR
667 668 669 670 671 672
std::string FieldSolver::getTagging_m() const {
    std::string tagging = Attributes::getString(itsAttr[AMR_TAGGING]);

    std::function<bool(const std::string&)> all_digits = [](const std::string& s) {
        // 15. Feb. 2019
        // https://stackoverflow.com/questions/19678572/how-to-validate-that-there-are-only-digits-in-a-string
kraus's avatar
kraus committed
673
        return std::all_of(s.begin(),  s.end(),
674 675 676 677 678 679 680 681 682
        [](char c) { return std::isdigit(c); });
    };

    if ( all_digits(tagging) )
        tagging = AmrObject::enum2string(std::stoi(tagging));
    return tagging;
}


frey_m's avatar
frey_m committed
683
void FieldSolver::initAmrObject_m() {
kraus's avatar
kraus committed
684

685 686 687 688
    auto domain = Attributes::getRealArray(itsAttr[AMR_DOMAIN_RATIO]);
    dynamic_cast<AmrPartBunch*>(itsBunch_m)->setAmrDomainRatio(
        Attributes::getRealArray(itsAttr[AMR_DOMAIN_RATIO])
    );
689
    itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) * 0.01);
kraus's avatar
kraus committed
690

frey_m's avatar
AMR:  
frey_m committed
691
    // setup initial info for creating the object
frey_m's avatar
frey_m committed
692 693 694 695 696 697 698 699 700 701 702 703 704 705
    AmrObject::AmrInfo info;
    info.grid[0]     = (int)this->getMX();
    info.grid[1]     = (int)this->getMY();
    info.grid[2]     = (int)this->getMT();
    info.maxgrid[0]  = Attributes::getReal(itsAttr[AMR_MAXGRIDX]);
    info.maxgrid[1]  = Attributes::getReal(itsAttr[AMR_MAXGRIDY]);
    info.maxgrid[2]  = Attributes::getReal(itsAttr[AMR_MAXGRIDZ]);
    info.bf[0]       = Attributes::getReal(itsAttr[AMR_BFX]);
    info.bf[1]       = Attributes::getReal(itsAttr[AMR_BFY]);
    info.bf[2]       = Attributes::getReal(itsAttr[AMR_BFZ]);
    info.maxlevel    = Attributes::getReal(itsAttr[AMR_MAXLEVEL]);
    info.refratio[0] = Attributes::getReal(itsAttr[AMR_REFX]);
    info.refratio[1] = Attributes::getReal(itsAttr[AMR_REFY]);
    info.refratio[2] = Attributes::getReal(itsAttr[AMR_REFZ]);
kraus's avatar
kraus committed
706

frey_m's avatar
AMR:  
frey_m committed
707
    itsAmrObject_mp = AmrBoxLib::create(info, dynamic_cast<AmrPartBunch*>(itsBunch_m));
kraus's avatar
kraus committed
708

709
    itsAmrObject_mp->setTagging( this->getTagging_m() );
kraus's avatar
kraus committed
710

frey_m's avatar
frey_m committed
711
    itsAmrObject_mp->setScalingFactor( Attributes::getReal(itsAttr[AMR_SCALING]) );
kraus's avatar
kraus committed
712

frey_m's avatar
frey_m committed
713
    itsAmrObject_mp->setChargeDensity( Attributes::getReal(itsAttr[AMR_DENSITY]) );
kraus's avatar
kraus committed
714

frey_m's avatar
frey_m committed
715 716 717
    itsAmrObject_mp->setMaxNumParticles(
        Attributes::getReal(itsAttr[AMR_MAX_NUM_PART])
    );
kraus's avatar
kraus committed
718

frey_m's avatar
frey_m committed
719 720 721
    itsAmrObject_mp->setMinNumParticles(
        Attributes::getReal(itsAttr[AMR_MIN_NUM_PART])
    );
frey_m's avatar
frey_m committed
722 723 724 725
}


void FieldSolver::initAmrSolver_m() {
726
    std::string fsType = Attributes::getString(itsAttr[FSTYPE]);
kraus's avatar
kraus committed
727

728
    if ( fsType == "ML" ) {
frey_m's avatar
frey_m committed
729
        if ( dynamic_cast<AmrBoxLib*>( itsAmrObject_mp.get() ) == 0 )
730
            throw OpalException("FieldSolver::initAmrSolver_m()",
731 732 733 734
                                "ML solver requires AMReX.");
        solver_m = new MLPoissonSolver(static_cast<AmrBoxLib*>(itsAmrObject_mp.get()));
#ifdef AMREX_ENABLE_FBASELIB
    } else if (fsType == "FMG") {
kraus's avatar
kraus committed
735

736 737
        if ( dynamic_cast<AmrBoxLib*>( itsAmrObject_mp.get() ) == 0 )
            throw OpalException("FieldSolver::initAmrSolver_m()",
738 739 740 741
                                "FMultiGrid solver requires AMReX.");

        solver_m = new FMGPoissonSolver(static_cast<AmrBoxLib*>(itsAmrObject_mp.get()));
#endif
742
    } else if (fsType == "HYPRE") {
743 744
        throw OpalException("FieldSolver::initAmrSolver_m()",
                            "HYPRE solver not yet implemented.");
745
    } else if (fsType == "HPGMG") {
746 747
        throw OpalException("FieldSolver::initAmrSolver_m()",
                            "HPGMG solver not yet implemented.");
frey_m's avatar
frey_m committed
748
    } else if (fsType == "AMR_MG") {
749
#ifdef HAVE_AMR_MG_SOLVER
750 751
        if ( dynamic_cast<AmrBoxLib*>( itsAmrObject_mp.get() ) == 0 )
            throw OpalException("FieldSolver::initAmrSolver_m()",
752
                                "FMultiGrid solver requires AMReX.");
kraus's avatar
kraus committed
753 754 755 756 757

        std::string bcz = Attributes::getString(itsAttr[deprecated::BCFFTT]);
        if (bcz == "") {
            bcz = Attributes::getString(itsAttr[BCFFTZ]);
        }
758 759
        solver_m = new AmrMultiGrid(static_cast<AmrBoxLib*>(itsAmrObject_mp.get()),
                                    Attributes::getString(itsAttr[ITSOLVER]),
frey_m's avatar
frey_m committed
760
                                    Attributes::getString(itsAttr[AMR_MG_PREC]),
761
                                    Attributes::getBool(itsAttr[AMR_MG_REBALANCE]),
frey_m's avatar
frey_m committed
762
                                    Attributes::getString(itsAttr[AMR_MG_REUSE]),
763 764
                                    Attributes::getString(itsAttr[BCFFTX]),
                                    Attributes::getString(itsAttr[BCFFTY]),
kraus's avatar
kraus committed
765
                                    bcz,
frey_m's avatar
frey_m committed
766 767 768 769
                                    Attributes::getString(itsAttr[AMR_MG_SMOOTHER]),
                                    Attributes::getReal(itsAttr[AMR_MG_NSWEEPS]),
                                    Attributes::getString(itsAttr[AMR_MG_INTERP]),
                                    Attributes::getString(itsAttr[AMR_MG_NORM]));
770

771 772
        dynamic_cast<AmrMultiGrid*>(solver_m)->setVerbose(
            Attributes::getBool(itsAttr[AMR_MG_VERBOSE]));
frey_m's avatar
frey_m committed
773 774 775

        dynamic_cast<AmrMultiGrid*>(solver_m)->setTolerance(
            Attributes::getReal(itsAttr[AMR_MG_TOL]));
776 777 778 779 780
#else
        throw OpalException("FieldSolver::initAmrSolver_m()",
                            "Multigrid solver not enabled! "
                            "Please build OPAL with -DENABLE_AMR_MG_SOLVER=1");
#endif
781
    } else
782 783
        throw OpalException("FieldSolver::initAmrSolver_m()",
                            "Unknown solver " + fsType + ".");
frey_m's avatar
frey_m committed
784
}
785
#endif