Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

ProblemStat.hpp 13.3 KB
Newer Older
1 2 3 4
#pragma once

#include <list>
#include <map>
5 6 7 8 9
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
10

11 12
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
13 14
#include <dune/grid/common/grid.hh>

15 16 17 18 19 20 21
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
22
#include <amdis/LocalAssemblerList.hpp>
23
#include <amdis/Marker.hpp>
24 25 26 27 28 29 30 31 32 33 34 35 36
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

37
#include <amdis/utility/TreeData.hpp>
38
#include <amdis/utility/TreePath.hpp>
39

40
namespace AMDiS
41
{
42 43 44 45
  // forward declaration
  template <class Traits>
  class ProblemInstat;

46
  template <class Traits>
47
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
48
      : public ProblemStatBase
49
      , public StandardProblemIteration
50
  {
51
    using Self = ProblemStat;
52

53 54
    friend class ProblemInstat<Traits>;

55
  public: // typedefs and static constants
56

57
    using GlobalBasis = typename Traits::GlobalBasis;
58 59 60
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
61
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
62

63
    /// Dimension of the grid
64
    static constexpr int dim = Grid::dimension;
65

66
    /// Dimension of the world
67
    static constexpr int dow = Grid::dimensionworld;
68

69
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
70
    using SystemVector = DOFVector<GlobalBasis, double>;
71

72
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
73

74
  public:
75 76
    /**
     * \brief Constructor. Takes the name of the problem that is used to
77
     * access values corresponding to this problem in the parameter file.
78
     **/
79
    explicit ProblemStat(std::string name)
80
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
81
      , name_(std::move(name))
82
    {}
83

84 85
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
86
    ProblemStat(std::string name, Grid& grid)
87
      : ProblemStat(std::move(name))
88
    {
89 90 91 92 93 94 95 96 97 98 99
      grid_ = Dune::stackobject_to_shared_ptr(grid);
      Parameters::get(name_ + "->mesh", gridName_);
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(std::move(name), grid)
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
100
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
101

102

103
    /**
104
     * \brief Initialisation of the problem.
105
     *
106 107 108
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
109
    void initialize(Flag initFlag,
110
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
111
                    Flag adoptFlag = INIT_NOTHING);
112

113 114

    /// Adds an operator to \ref A.
115
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
116 117
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
118

119
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
120 121
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
122
    /** @} */
123 124 125


    /// Adds an operator to \ref rhs.
126
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
127 128
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(Operator const& op, TreePath = {});
129

130
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
131 132
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
133
    /** @} */
134

135

136
    /// Adds a Dirichlet boundary condition
137
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
138
    void addDirichletBC(Predicate const& predicate,
139
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
140
                        Values const& values);
141

142
  public:
143

144
    /// Implementation of \ref ProblemStatBase::solve
145
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
146 147
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
148

149
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
150
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
151 152 153
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
154

155
    /// Writes output files.
156
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
157

Praetorius, Simon's avatar
Praetorius, Simon committed
158

159
  public: // get-methods
160

161 162 163 164 165 166 167 168 169 170 171
    /// Returns a reference to system-matrix, \ref systemMatrix_
    SystemMatrix&       getSystemMatrix()       { return *systemMatrix_; }
    SystemMatrix const& getSystemMatrix() const { return *systemMatrix_; }

    /// Returns a reference to the solution vector, \ref solution_
    SystemVector&       getSolutionVector()       { return *solution_; }
    SystemVector const& getSolutionVector() const { return *solution_; }

    /// Return a reference to the rhs system-vector, \ref rhs
    SystemVector&       getRhsVector()       { return *rhs_; }
    SystemVector const& getRhsVector() const { return *rhs_; }
172

173 174 175 176 177 178

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {})
    {
      auto&& tp = makeTreePath(path);
179
      return makeDOFVectorView(*solution_, tp);
180
    }
181

182 183 184
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
185
    {
186
      auto&& tp = makeTreePath(path);
187
      return makeDOFVectorView(*solution_, tp);
188
    }
189 190


191 192 193
    /// Return a reference to the linear solver, \ref linearSolver
    LinearSolverType&       getSolver()       { return *linearSolver_; }
    LinearSolverType const& getSolver() const { return *linearSolver_; }
194

195
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
196
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
197
    {
198
      linearSolver_ = solver;
199 200
    }

201 202 203
    /// Return a reference to the grid, \ref grid
    Grid&       grid()       { return *grid_; }
    Grid const& grid() const { return *grid_; }
204

205
    /// Set the grid. Stores pointer and initializes feSpaces
206
    /// matrices and vectors, as well as the file-writer.
207
    void setGrid(std::shared_ptr<Grid> const& grid)
208
    {
209
      grid_ = grid;
210

211
      createGlobalBasis();
212 213 214 215
      createMatricesAndVectors();
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
216
    /// Return the gridView of the leaf-level
217
    GridView const& gridView() { return globalBasis_->gridView(); }
218

219
    /// Return the \ref feSpaces
220
    GlobalBasis const& globalBasis() { return *globalBasis_; }
221 222


223
    /// Implementation of \ref ProblemStatBase::getName
224 225
    virtual std::string getName() const override
    {
226
      return name_;
227
    }
228

229
  protected: // initialization methods
230

231 232 233 234 235 236 237 238 239 240 241
    template <class T, class GV>
    using HasCreate = decltype(T::create(std::declval<GV>()));

    void createGlobalBasis()
    {
      createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
      initGlobalBasis(*globalBasis_);
    }

    void createGlobalBasis(std::true_type)
    {
242
      assert( bool(grid_) );
243
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
244
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
245 246 247 248
    }

    void createGlobalBasis(std::false_type)
    {
249
      error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
250 251
    }

252
    void createGrid()
253
    {
254 255
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
256
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
257

258
      grid_ = MeshCreator<Grid>::create(gridName_);
259

260
      msg("Create grid:");
261 262 263
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
264
      msg("");
265
    }
266

267 268
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
269 270 271 272
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis_->localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
273
    }
274

275 276
    void createMatricesAndVectors()
    {
277 278 279
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
      solution_ = std::make_shared<SystemVector>(*globalBasis_);
      rhs_ = std::make_shared<SystemVector>(*globalBasis_);
280

281
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
282 283
      {
        std::string i = to_string(treePath);
284 285 286
        estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates_[i].size(); j++)
          estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
287
      });
288
    }
289

290
    void createSolver()
291
    {
292
      std::string solverName = "cg";
293
      Parameters::get(name_ + "->solver->name", solverName);
294

295
      auto solverCreator
296
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
297

298
      linearSolver_ = solverCreator->create(name_ + "->solver");
299
    }
300

301
    void createMarker();
302

303
    void createFileWriter();
304

Praetorius, Simon's avatar
Praetorius, Simon committed
305

306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
  public: // implementation of iteration interface methods

    /**
      * \brief Determines the execution order of the single adaption steps.
      *
      * If adapt is true, mesh adaption will be performed. This allows to avoid
      * mesh adaption, e.g. in timestep adaption loops of timestep adaptive
      * strategies.
      *
      * Implementation of \ref StandardProblemIteration::oneIteration.
      **/
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    /// Implementation of \ref ProblemStatBase::buildBeforeRefine.
    virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
    virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::coarsenMesh.
    virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::markElements.
338
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
339 340


341
  private:
342
    /// Name of this problem.
343
    std::string name_;
344

345
    /// Grid of this problem.
346
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
347

348 349 350
    /// Number of grids
    int nGrids = 1;

351
    /// Name of the grid
352
    std::string gridName_ = "none";
353

354
    /// FE spaces of this problem.
355 356
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
357

358
    /// A FileWriter object
359
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
360

361
    /// Pointer to the adaptation markers
362
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
363

364
    /// An object of the linearSolver Interface
365
    std::shared_ptr<LinearSolverType> linearSolver_;
366

367
    /// A block-matrix that is filled during assembling
368
    std::shared_ptr<SystemMatrix> systemMatrix_;
369

370
    /// A block-vector with the solution components
371
    std::shared_ptr<SystemVector> solution_;
372

373
    /// A vector with the local element error estimates
374 375
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double> > estimates_;
376

377
    /// A block-vector (load-vector) corresponding to the right.hand side
378
    /// of the equation, filled during assembling
379
    std::shared_ptr<SystemVector> rhs_;
380 381


382 383
  private: // some internal data-structures

384 385 386
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
387
  };
388

389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
  ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
  makeProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
  {
    return {std::move(name), grid, globalBasis};
  }

404

405
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
406
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
407
#endif
408

409
} // end namespace AMDiS
410 411 412

#include "ProblemStat.inc.hpp"