ProblemStat.hpp 17 KB
Newer Older
1
2
#pragma once

3
4
5
6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7
8
9
10
11
12
#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

13
#include <dune/istl/matrix.hh>
14
15
#include <dune/grid/common/grid.hh>

16
17
18
19
20
21
22
23
24
25
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/FileWriter.hpp>
#include <dune/amdis/Flag.hpp>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
26
#include <dune/amdis/ProblemStatTraits.hpp>
27
28
#include <dune/amdis/StandardProblemIteration.hpp>

29
#include <dune/amdis/common/OptionalDelete.hpp>
30
31
#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
32
#include <dune/amdis/common/TypeDefs.hpp>
33
#include <dune/amdis/common/Utility.hpp>
34

35
namespace AMDiS
36
{
37
  struct BoundaryType { int b; };
38

39
  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
40
41
  class ProblemStatSeq
      : public ProblemStatBase
42
43
  {
    using Self = ProblemStatSeq;
44

45
  public: // typedefs and static constants
46

47
48
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
49
    using MeshView = typename Mesh::LeafGridView;
50

51
52
53
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

54
55
    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
56
57


58
    /// Dimension of the mesh
59
    static constexpr int dim = Traits::dim;
60

61
    /// Dimension of the world
62
    static constexpr int dow = Traits::dimworld;
63

64
    /// Number of problem components
65
    static constexpr int nComponents = Traits::nComponents;
66
67


68
    template <std::size_t I>
69
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
70

71
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
72

73
74
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
75

76
77
    using ElementOperator = Operator<MeshView, Element>;
    using IntersectionOperator = Operator<MeshView, typename MeshView::Intersection>;
78

79
80
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
81

82
83
    template <class OperatorType>
    struct Scaled
84
    {
85
      Scaled(std::shared_ptr<OperatorType> const& op, double* factor = nullptr, double* estFactor = nullptr, BoundaryType b = {0})
86
87
88
89
        : op(op)
        , factor(factor)
        , estFactor(estFactor)
      {}
90

91
92
93
94
95
96
      std::shared_ptr<OperatorType> op;
      double* factor;
      double* estFactor;
      BoundaryType b;
    };

97
  public:
98
99
    /**
     * \brief Constructor. Takes the name of the problem that is used to
100
     * access values correpsonding to this püroblem in the parameter file.
101
     *
102
     * Parameters read by ProblemStatSeq, with name 'PROB'
103
     *   PROB->names: \ref componentNames
104
     **/
105
    ProblemStatSeq(std::string name)
106
107
      : name(name)
    {
108
      Parameters::get(name + "->names", componentNames);
109
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
110
        componentNames.push_back("solution[" + std::to_string(i) + "]");
111
    }
112

113
114
115
116
    ProblemStatSeq(std::string name, Mesh& mesh)
      : ProblemStatSeq(name)
    {
      this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
117
118
119
120
      componentMeshes.resize(nComponents, this->mesh.get());

      meshName = "";
      Parameters::get(name + "->mesh", meshName);
121
122
    }

123

124
    /**
125
     * \brief Initialisation of the problem.
126
     *
127
128
129
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
130
    void initialize(Flag initFlag,
131
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
132
                    Flag adoptFlag = INIT_NOTHING);
133

134
135

    /// Adds an operator to \ref A.
136
    /** @{ */
137
    template <class OperatorType>
138
139
    void addMatrixOperator(OperatorType const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);
140

141
    // operator evaluated on the boundary
142
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
143
                           double* factor = nullptr, double* estFactor = nullptr);
144
145

    void addMatrixOperator(std::map< std::pair<int,int>,
146
                                     std::shared_ptr<ElementOperator> > ops);
147
    void addMatrixOperator(std::map< std::pair<int,int>,
148
                                     std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
149
    /** @} */
150
151
152


    /// Adds an operator to \ref rhs.
153
    /** @{ */
154
    template <class OperatorType>
155
156
    void addVectorOperator(OperatorType const& op, int i,
                           double* factor = nullptr, double* estFactor = nullptr);
157

158
    // operator evaluated on the boundary
159
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
160
                           double* factor = nullptr, double* estFactor = nullptr);
161

162
163
    void addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops);
    void addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
164
    /** @} */
165
166
167

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
168
    void addDirichletBC(Predicate const& predicate,
Praetorius, Simon's avatar
Praetorius, Simon committed
169
170
                        int row, int col,
                        Values const& values);
171

172

173
    /// Implementation of \ref ProblemStatBase::solve
174
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
177
178


179
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
180
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
181
182
183
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
184

185
186

    /// Writes output files.
187
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
188

Praetorius, Simon's avatar
Praetorius, Simon committed
189

190
  public: // get-methods
191
192

    /// Returns nr of components \ref nComponents
193
    static constexpr int getNumComponents()
194
195
196
    {
      return nComponents;
    }
197
198


199
200
201
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
202
203


204
205
206
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
207

208
    /// Return the i'th \ref solution component
209
210
    template <std::size_t I = 0>
    decltype(auto) getSolution(const index_t<I> _i = {})
211
    {
212
      return (*solution)[_i];
213
    }
214
215


216
217
218
    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
219
220


221
222
    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }
223

Praetorius, Simon's avatar
Praetorius, Simon committed
224
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
225
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
226
      linearSolver = solver;
227
228
    }

229
    /// Return the \ref mesh
230
    auto const& getMesh() { return mesh; }
231

232
    void setMesh(Mesh& mesh_)
233
    {
234
      mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
235
      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());
236
237
238
239
240
241
242

      createFeSpaces();
      createMatricesAndVectors();
      createFileWriter();
    }


243
244
245
246
247
    /// Return the gridView of the leaf-level
    auto leafGridView() { return mesh->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return mesh->levelGridView(level); }
248

249
250
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
251
252


253
    /// Return the I'th \ref feSpaces component
254
255
256
257
258
259
260
261
    template <std::size_t I = 0>
    FeSpace<I> const& getFeSpace(const index_t<I> = {}) const
    {
      return std::get<I>(*feSpaces);
    }

    template <std::size_t I = 0>
    FeSpace<I>& getFeSpace(const index_t<I> = {})
262
    {
263
      return std::get<I>(*feSpaces);
264
    }
265

266
    /// Implementation of \ref ProblemStatBase::getName
267
268
269
270
    virtual std::string getName() const override
    {
      return name;
    }
271

272
    /// Return a vector of names of the problem components
273
274
275
276
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
277

278
  protected: // initialization methods
279

280
281
    void createMesh()
    {
282
283
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
284
285
286
      test_exit(!meshName.empty(),
                "No mesh name specified for '", name, "->mesh'!");

287
      mesh = MeshCreator<Mesh>::create(meshName);
288
289
290
291
292
293

      msg("Create mesh:");
      msg("#elements = "   , mesh->size(0));
      msg("#faces/edges = ", mesh->size(1));
      msg("#vertices = "   , mesh->size(dim));
      msg("");
294
    }
295

296
297
    void createFeSpaces()
    {
298
      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
299
    }
300

301
302
    void createMatricesAndVectors()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
303
304
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
305

306
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
Praetorius, Simon's avatar
Praetorius, Simon committed
307
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
308
    }
309

310
    void createSolver()
311
    {
312
313
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
314

315
316
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
317

318
      linearSolver = solverCreator->create(name + "->solver");
319
    }
320

321
322
    void createFileWriter()
    {
323
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
324
    }
325

Praetorius, Simon's avatar
Praetorius, Simon committed
326

327
  protected: // sub-methods to assemble DOFMatrix
328

329
330
    template <class LhsData, class RhsData, class GV>
    void assemble(LhsData lhs, RhsData rhs, GV const& gridView);
331

332
    template <class RowView, class ColView>
333
    bool getElementMatrix(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
334
335
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
336
337
338
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);
339

340
    template <class RowView>
341
    bool getElementVector(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
342
                          ElementVector& elementVector,
343
344
345
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);
346
347


348
349
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
Praetorius, Simon's avatar
Praetorius, Simon committed
350
351
352
                          RowIndexSet const& rowIndexSet,
                          ColIndexSet const& colIndexSet,
                          ElementMatrix const& elementMatrix);
353

354
355
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
Praetorius, Simon's avatar
Praetorius, Simon committed
356
357
                          IndexSet const& indexSet,
                          ElementVector const& elementvector);
358
359


360
  private:
361
362
363
364
365
366
    /// Name of this problem.
    std::string name;

    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
    std::vector<std::string> componentNames;
367

368
    /// Mesh of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
369
    std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
370

371
    /// Name of the mesh
372
    std::string meshName = "none";
373

374
375
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
376

377
    /// FE spaces of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
378
    std::shared_ptr<FeSpaces> feSpaces; // eventuell const
379

380
    /// A FileWriter object
Praetorius, Simon's avatar
Praetorius, Simon committed
381
    std::shared_ptr<FileWriter<Traits>> filewriter;
382

383
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
384
    std::shared_ptr<LinearSolverType> linearSolver;
385

386
    /// A block-matrix that is filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
387
    std::shared_ptr<SystemMatrixType> systemMatrix;
388

389
    /// A block-vector with the solution components
Praetorius, Simon's avatar
Praetorius, Simon committed
390
    std::shared_ptr<SystemVectorType> solution;
391
392

    /// A block-vector (load-vector) corresponding to the right.hand side
393
    /// of the equation, filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
394
    std::shared_ptr<SystemVectorType> rhs;
395
396


397
    template <class T>
398
    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
399

400
    template <class T>
401
    using VectorEntries = Dune::FieldVector<T, nComponents>;
402

403
    using BoundaryCondition = std::list<std::shared_ptr<DirichletBC<WorldVector>>>;
404
405

    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
406
    /// each matrix block
407
    MatrixEntries<BoundaryCondition> dirichletBc;
408

409
    /// A map (i,j) -> list<ElementOperator> string the differential operators for
410
    /// each matrix block
411
412
413
    MatrixEntries<std::list<Scaled<ElementOperator>>> matrixOperators;
    MatrixEntries<std::list<Scaled<IntersectionOperator>>> matrixBoundaryOperators;
    MatrixEntries<std::list<Scaled<IntersectionOperator>>> matrixIntersectionOperators;
414
415
    MatrixEntries<bool> matrixAssembled; // if false, do reassemble
    MatrixEntries<bool> matrixChanging;  // if true, or vectorAssembled false, do reassemble
416

417
    /// A map (i) -> list<ElementOperator> string the differential operators for
418
    /// each rhs block
419
420
421
    VectorEntries<std::list<Scaled<ElementOperator>>> vectorOperators;
    VectorEntries<std::list<Scaled<IntersectionOperator>>> vectorBoundaryOperators;
    VectorEntries<std::list<Scaled<IntersectionOperator>>> vectorIntersectionOperators;
422
423
    VectorEntries<bool> vectorAssembled; // if false, do reassemble
    VectorEntries<bool> vectorChanging;  // if true, or vectorAssembled false, do reassemble
424

Praetorius, Simon's avatar
Praetorius, Simon committed
425
    template <int R, int C>
426
    struct MatrixData
427
    {
428
      using DOFMatrixType =
Praetorius, Simon's avatar
Praetorius, Simon committed
429
        std::tuple_element_t<C, std::tuple_element_t<R, typename SystemMatrixType::DOFMatrices>>;
430

431
      DOFMatrixType& matrix;
432
433
434
      std::list<Scaled<ElementOperator>>& operators;
      std::list<Scaled<IntersectionOperator>>& boundary_operators;
      std::list<Scaled<IntersectionOperator>>& intersection_operators;
435
      bool assemble = true;
436
    };
437

Praetorius, Simon's avatar
Praetorius, Simon committed
438
    template <int I>
439
440
    struct VectorData
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
441
      using DOFVectorType = std::tuple_element_t<I, typename SystemVectorType::DOFVectors>;
442

443
      DOFVectorType& vector;
444
445
446
      std::list<Scaled<ElementOperator>>& operators;
      std::list<Scaled<IntersectionOperator>>& boundary_operators;
      std::list<Scaled<IntersectionOperator>>& intersection_operators;
447
      bool assemble = true;
448
    };
449
  };
450
451


452
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
453
454
455
456
//   extern template class ProblemStatSeq<ProblemStatTraits<1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,2>>;
//   extern template class ProblemStatSeq<ProblemStatTraits<3>>;
457
#endif
458

459
  namespace Impl
460
461
  {
    template <class ProblemStatType>
Praetorius, Simon's avatar
Praetorius, Simon committed
462
463
464
    struct ProblemStat
        : public ProblemStatType
        , public StandardProblemIteration
465
466
467
468
469
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
470
471
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
472
473
      {}

474
475
476
477
478
479
      template <class Grid>
      ProblemStat(std::string nameStr, Grid& grid)
        : ProblemStatType(nameStr, grid)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      {}

480
      /**
481
482
483
484
485
486
       * \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.
       *
487
488
       * Implementation of \ref StandardProblemIteration::oneIteration.
       **/
489
490
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
491
        for (std::size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
492
493
494
495
496
497
498
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

        return StandardProblemIteration::oneIteration(adaptInfo, toDo);
      }
499

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

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

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

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

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

515
      /// Implementation of \ref ProblemStatBase::markElements.
516
      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
517
    };
518

519
  } // end namespace Impl
520

521
  template <class Traits>
522
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
523

524
} // end namespace AMDiS
525
526
527

#include "ProblemStat.inc.hpp"