FieldSolver.cpp 12.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"
22
#ifdef HAVE_SAAMG_SOLVER
23 24
#include "Solvers/MGPoissonSolver.h"
#endif
gsell's avatar
gsell committed
25 26 27 28 29 30 31 32 33
#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"
34
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

using namespace Expressions;
using namespace Physics;

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

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

// The attributes of class FieldSolver.
namespace {
    enum {
        FSTYPE,   // The field solver name
        // 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
        BCFFTX,     // boundary condition in x [FFT only]
        BCFFTY,     // boundary condition in y [FFT only]
        BCFFTT,     // boundary condition in z [FFT only]
        GREENSF,    // holds greensfunction to be used [FFT only]
        BBOXINCR,   // how much the boundingbox is increased
60 61 62 63 64 65
        GEOMETRY,   // geometry of boundary [SAAMG only]
        ITSOLVER,   // iterative solver [SAAMG only]
        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]
gsell's avatar
gsell committed
66 67 68 69 70 71 72 73 74 75 76
        RPP,        // defines in units of the meshsize where the PP interactions takes place [P3M only]
        // FOR XXX BASED SOLVER
        SIZE
    };
}


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

77
    itsAttr[FSTYPE] = Attributes::makeString("FSTYPE", "Name of the attached field solver: FFT, FFTPERIODIC, SAAMG, AMR, and NONE ");
gsell's avatar
gsell committed
78 79 80 81 82 83 84 85 86 87 88 89

    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)");

    itsAttr[PARFFTX] = Attributes::makeBool("PARFFTX", "True, dimension 0 i.e x is parallelized", false);
    itsAttr[PARFFTY] = Attributes::makeBool("PARFFTY", "True, dimension 1 i.e y is parallelized", false);
    itsAttr[PARFFTT] = Attributes::makeBool("PARFFTT", "True, dimension 2 i.e z(t) is parallelized", true);

    //FFT ONLY:
    itsAttr[BCFFTX] = Attributes::makeString("BCFFTX", "Boundary conditions in x: open, dirichlet (box) ");
    itsAttr[BCFFTY] = Attributes::makeString("BCFFTY", "Boundary conditions in y: open, dirichlet (box) ");
90
    itsAttr[BCFFTT] = Attributes::makeString("BCFFTT", "Boundary conditions in z(t): open, periodoc");
gsell's avatar
gsell committed
91 92

    itsAttr[GREENSF]  = Attributes::makeString("GREENSF", "Which Greensfunction to be used [STANDARD | INTEGRATED]", "INTEGRATED");
93
    itsAttr[BBOXINCR] = Attributes::makeReal("BBOXINCR", "Increase of bounding box in % ", 2.0);
gsell's avatar
gsell committed
94 95 96 97

    // P3M only:
    itsAttr[RPP]  = Attributes::makeReal("RPP", "Defines in units of the meshsize where the PP interactions takes place ", 1);

98
    //SAAMG and in case of FFT with dirichlet BC in x and y
gsell's avatar
gsell committed
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
    itsAttr[GEOMETRY] = Attributes::makeString("GEOMETRY", "GEOMETRY to be used as domain boundary", "");
    itsAttr[ITSOLVER]  = Attributes::makeString("ITSOLVER", "Type of iterative solver [CG | BiCGSTAB | GMRES]", "CG");
    itsAttr[INTERPL]  = Attributes::makeString("INTERPL", "interpolation used for boundary points [CONSTANT | LINEAR | QUADRATIC]", "LINEAR");
    itsAttr[TOL] = Attributes::makeReal("TOL", "Tolerance for iterative solver", 1e-8);
    itsAttr[MAXITERS] = Attributes::makeReal("MAXITERS", "Maximum number of iterations of iterative solver", 100);
    itsAttr[PRECMODE]  = Attributes::makeString("PRECMODE", "Preconditioner Mode [STD | HIERARCHY | REUSE]", "HIERARCHY");

    mesh_m = 0;
    FL_m = 0;
    PL_m = 0;
}


FieldSolver::FieldSolver(const string &name, FieldSolver *parent):
    Definition(name, parent)
{}


FieldSolver::~FieldSolver() {
118

gsell's avatar
gsell committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
}

FieldSolver *FieldSolver::clone(const string &name) {
    return new FieldSolver(name, this);
}

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

FieldSolver *FieldSolver::find(const string &name) {
    FieldSolver *fs = dynamic_cast<FieldSolver *>(OpalData::getInstance()->find(name));

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


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;

    if(Attributes::getString(itsAttr[FSTYPE]) == "FFTPERIODIC") {
        decomp[0] = decomp[1] = SERIAL;
        decomp[2] = PARALLEL;
    }
    // create prototype mesh and layout objects for this problem domain
    mesh_m   = new Mesh_t(domain);
    FL_m     = new FieldLayout_t(*mesh_m, decomp);
    PL_m     = new Layout_t(*FL_m, *mesh_m);
191 192 193
    // OpalData::getInstance()->setMesh(mesh_m);
    // OpalData::getInstance()->setFieldLayout(FL_m);
    // OpalData::getInstance()->setLayout(PL_m);
gsell's avatar
gsell committed
194 195
}

adelmann's avatar
adelmann committed
196 197 198 199
bool FieldSolver::hasPeriodicZ() { 
  return Attributes::getString(itsAttr[BCFFTT])==std::string("PERIODIC");
}

200 201 202 203 204
bool FieldSolver::isAMRSolver() { 
  return Attributes::getString(itsAttr[FSTYPE])==std::string("AMR");
}


gsell's avatar
gsell committed
205 206
void FieldSolver::initSolver(PartBunch &b) {
    itsBunch_m = &b;
207 208 209
     string bcx = Attributes::getString(itsAttr[BCFFTX]);
     string bcy = Attributes::getString(itsAttr[BCFFTY]);
     string bcz = Attributes::getString(itsAttr[BCFFTT]);
gsell's avatar
gsell committed
210

211
     if(Attributes::getString(itsAttr[FSTYPE]) == "FFT" || Attributes::getString(itsAttr[FSTYPE]) == "P3M") {
gsell's avatar
gsell committed
212

213
	bool sinTrafo = ((bcx == string("DIRICHLET")) && (bcy == string("DIRICHLET")) && (bcz == string("DIRICHLET")));
gsell's avatar
gsell committed
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
        if(sinTrafo) {
            std::cout << "FFTBOX ACTIVE" << std::endl;
            //we go over all geometries and add the Geometry Elements to the geometry list
            std::string geoms = Attributes::getString(itsAttr[GEOMETRY]);
            std::string tmp = "";
            //split and add all to list
            std::vector<BoundaryGeometry *> geometries;
            for(unsigned int i = 0; i <= geoms.length(); i++) {
                if(geoms[i] == ',' || i == geoms.length()) {
                    BoundaryGeometry *geom = BoundaryGeometry::find(tmp);
                    if(geom != 0)
                        geometries.push_back(geom);
                    tmp.clear();
                } else
                    tmp += geoms[i];
            }
            BoundaryGeometry *ttmp = geometries[0];
            solver_m = new FFTBoxPoissonSolver(mesh_m, FL_m, Attributes::getString(itsAttr[GREENSF]), ttmp->getA());
            itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
            fsType_m = "FFTBOX";
        } else {
adelmann's avatar
adelmann committed
235
            solver_m = new FFTPoissonSolver(mesh_m, FL_m, Attributes::getString(itsAttr[GREENSF]), bcz);
gsell's avatar
gsell committed
236 237 238 239
            itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
            fsType_m = "FFT";
        }

240
    } else if(Attributes::getString(itsAttr[FSTYPE]) == "SAAMG") {
241
#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
242 243 244 245 246
        //we go over all geometries and add the Geometry Elements to the geometry list
        std::string geoms = Attributes::getString(itsAttr[GEOMETRY]);
        std::string tmp = "";
        //split and add all to list
        std::vector<BoundaryGeometry *> geometries;
247
        for(unsigned int i = 0; i <= geoms.length(); i++) {
gsell's avatar
gsell committed
248
            if(geoms[i] == ',' || i == geoms.length()) {
249
                BoundaryGeometry *geom = BoundaryGeometry::find(Attributes::getString(itsAttr[GEOMETRY]))->clone(getOpalName() + string("_geometry"));
gsell's avatar
gsell committed
250 251 252 253 254 255 256 257 258
                if(geom != 0) {
                    geometries.push_back(geom);
                }
                tmp.clear();
            } else
                tmp += geoms[i];
        }
        solver_m = new MGPoissonSolver(b, mesh_m, FL_m, geometries, Attributes::getString(itsAttr[ITSOLVER]), Attributes::getString(itsAttr[INTERPL]), Attributes::getReal(itsAttr[TOL]), (int)Attributes::getReal(itsAttr[MAXITERS]), Attributes::getString(itsAttr[PRECMODE]));
        itsBunch_m->set_meshEnlargement(Attributes::getReal(itsAttr[BBOXINCR]) / 100.0);
259
        fsType_m = "SAAMG";
gsell's avatar
gsell committed
260
#else
261
        INFOMSG("SAAMG Solver not enabled! Please build OPAL with SAAMG_SOLVER enabled" << endl);
gsell's avatar
gsell committed
262
        INFOMSG("switching to FFT solver..." << endl);
263
        solver_m = new FFTPoissonSolver(mesh_m, FL_m, Attributes::getString(itsAttr[GREENSF]),bcz);
gsell's avatar
gsell committed
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
        fsType_m = "FFT";
#endif
    } else if(Attributes::getString(itsAttr[FSTYPE]) == "P3M") {

        fsType_m = "P3M";
    } else {
        solver_m = 0;
        INFOMSG("no solver attached" << endl);
    }

}


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

Inform &FieldSolver::printInfo(Inform &os) const {
282 283 284 285 286
  std::string fsType;
  if (Attributes::getString(itsAttr[BCFFTT])==std::string("PERIODIC"))
	fsType = Attributes::getString(itsAttr[FSTYPE])+"-zPeriodic";
      else
	fsType = Attributes::getString(itsAttr[FSTYPE]);
gsell's avatar
gsell committed
287 288
    os << "* ************* F I E L D S O L V E R ********************************************** " << endl;
    os << "* FIELDSOLVER  " << getOpalName() << '\n'
289
       << "* TYPE         " << fsType << '\n'
gsell's avatar
gsell committed
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
       << "* 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;
    if(Attributes::getString(itsAttr[FSTYPE]) == "P3M")
        os << "* RPP          " << Attributes::getReal(itsAttr[RPP]) << endl;
    if(Attributes::getString(itsAttr[FSTYPE]) == "FFT" || Attributes::getString(itsAttr[FSTYPE]) == "P3M")
        os << "* GRRENSF      " << Attributes::getString(itsAttr[GREENSF]) << endl;
    else
        os << "* GEOMETRY     " << Attributes::getString(itsAttr[GEOMETRY]) << '\n'
           << "* ITSOLVER     " << Attributes::getString(itsAttr[ITSOLVER])   << '\n'
           << "* INTERPL      " << Attributes::getString(itsAttr[INTERPL])  << '\n'
           << "* TOL          " << Attributes::getReal(itsAttr[TOL])        << '\n'
           << "* MAXITERS     " << Attributes::getReal(itsAttr[MAXITERS]) << '\n'
           << "* PRECMODE     " << Attributes::getString(itsAttr[PRECMODE])   << endl;
    if(Attributes::getBool(itsAttr[PARFFTX]))
        os << "* XDIM is parallel  " << endl;
    else
        os << "* XDIM is serial  " << endl;

    if(Attributes::getBool(itsAttr[PARFFTY]))
        os << "* YDIM is parallel  " << endl;
    else
        os << "* YDIM is serial  " << endl;

    if(Attributes::getBool(itsAttr[PARFFTT]))
        os << "* Z(T)DIM is parallel  " << endl;
    else
        os << "* Z(T)DIM is serial  " << endl;

    INFOMSG(*mesh_m << endl);
    INFOMSG(*PL_m << endl);
    if(solver_m)
        os << *solver_m << endl;
    os << "* ********************************************************************************** " << endl;
    return os;
}