ProblemStat.hpp 15.3 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
29
30
#include <dune/amdis/StandardProblemIteration.hpp>

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

34
namespace AMDiS
35
{
36

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

43
  public: // typedefs and static constants
44

45
46
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
47
    using MeshView = typename Mesh::LeafGridView;
48

49
50
51
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

52
53
    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
54
55


56
    /// Dimension of the mesh
57
    static constexpr int dim = Traits::dim;
58

59
    /// Dimension of the world
60
    static constexpr int dow = Traits::dimworld;
61

62
    /// Number of problem components
63
    static constexpr int nComponents = Traits::nComponents;
64
65


66
67
    template <size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
68

69
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
70

71
72
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
73

74
    using OperatorType = Operator<MeshView>;
75
76
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
77

78
  public:
79
80
    /**
     * \brief Constructor. Takes the name of the problem that is used to
81
     * access values correpsonding to this püroblem in the parameter file.
82
     *
83
     * Parameters read by ProblemStatSeq, with name 'PROB'
84
     *   PROB->names: \ref componentNames
85
     **/
86
    ProblemStatSeq(std::string name)
87
88
      : name(name)
    {
89
      Parameters::get(name + "->names", componentNames);
90
      for (size_t i = componentNames.size(); i < nComponents; ++i)
91
        componentNames.push_back("solution[" + std::to_string(i) + "]");
92
    }
93
94


95
    /**
96
     * \brief Initialisation of the problem.
97
     *
98
99
100
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
101
    void initialize(Flag initFlag,
Praetorius, Simon's avatar
Praetorius, Simon committed
102
103
                    Self* adoptProblem = NULL,
                    Flag adoptFlag = INIT_NOTHING);
104

105
106

    /// Adds an operator to \ref A.
107
    /** @{ */
108
    void addMatrixOperator(OperatorType& op,
Praetorius, Simon's avatar
Praetorius, Simon committed
109
			                     int i, int j,
110
                           double* factor = NULL,
111
                           double* estFactor = NULL);
112

Praetorius, Simon's avatar
Praetorius, Simon committed
113
    void addMatrixOperator(std::shared_ptr<OperatorType> op,
114
                           int i, int j,
115
                           double* factor = NULL,
116
                           double* estFactor = NULL);
117
118

    void addMatrixOperator(std::map< std::pair<int,int>,
Praetorius, Simon's avatar
Praetorius, Simon committed
119
                                     std::shared_ptr<OperatorType> > ops);
120
    void addMatrixOperator(std::map< std::pair<int,int>,
Praetorius, Simon's avatar
Praetorius, Simon committed
121
                                     std::pair<std::shared_ptr<OperatorType>, bool> > ops);
122
    /** @} */
123
124
125


    /// Adds an operator to \ref rhs.
126
    /** @{ */
127
    void addVectorOperator(OperatorType& op,
128
                           int i,
129
                           double* factor = NULL,
130
                           double* estFactor = NULL);
131

Praetorius, Simon's avatar
Praetorius, Simon committed
132
    void addVectorOperator(std::shared_ptr<OperatorType> op,
133
                           int i,
134
                           double* factor = NULL,
135
                           double* estFactor = NULL);
136

Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
    void addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops);
    void addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops);
139
    /** @} */
140
141
142

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
143
    void addDirichletBC(Predicate const& predicate,
Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
                        int row, int col,
                        Values const& values);
146

147

148
    /// Implementation of \ref ProblemStatBase::solve
149
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
150
151
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
152
153


154
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
155
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
156
157
158
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
159

160
161

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

Praetorius, Simon's avatar
Praetorius, Simon committed
164

165
  public: // get-methods
166
167

    /// Returns nr of components \ref nComponents
168
    static constexpr int getNumComponents()
169
170
171
    {
      return nComponents;
    }
172
173


174
175
176
    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
177
178


179
180
181
    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
182

183
    /// Return the i'th \ref solution component
184
    template <size_t I = 0>
185
    decltype(auto) getSolution(const index_<I> _i = {})
186
    {
187
      return (*solution)[_i];
188
    }
189
190


191
192
193
    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
194
195


196
197
    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }
198

Praetorius, Simon's avatar
Praetorius, Simon committed
199
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
200
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
201
      linearSolver = solver;
202
203
    }

204
205
    /// Return the \ref mesh
    auto getMesh() { return mesh; }
206

207
208
209
    void setMesh(std::shared_ptr<Mesh> const& mesh_)
    {
      mesh = mesh_;
Praetorius, Simon's avatar
Praetorius, Simon committed
210
      meshView = std::make_shared<MeshView>(mesh->leafGridView());
211
212
213
214
215
216
217

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


218
    /// Return the \ref meshView
219
    auto getMeshView() { return meshView; }
220

221
222
    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }
223
224


225
    /// Return the I'th \ref feSpaces component
226
    template <size_t I = 0>
227
    FeSpace<I> const& getFeSpace(const index_<I> = {}) const
228
    {
229
      return std::get<I>(*feSpaces);
230
    }
231

232
    /// Implementation of \ref ProblemStatBase::getName
233
234
235
236
    virtual std::string getName() const override
    {
      return name;
    }
237

238
    /// Return a vector of names of the problem components
239
240
241
242
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
243

244
  protected: // initialization methods
245

246
247
    void createMesh()
    {
248
249
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
250
251
252
      test_exit(!meshName.empty(),
                "No mesh name specified for '", name, "->mesh'!");

253
      mesh = MeshCreator<Mesh>::create(meshName);
Praetorius, Simon's avatar
Praetorius, Simon committed
254
      meshView = std::make_shared<MeshView>(mesh->leafGridView());
255
256
257
258
259
260

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

263
264
    void createFeSpaces()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
265
      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
266
    }
267

268
269
    void createMatricesAndVectors()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
270
271
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
272

273
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
Praetorius, Simon's avatar
Praetorius, Simon committed
274
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
275
    }
276

277
    void createSolver()
278
    {
279
280
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
281

282
283
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
284

285
      linearSolver = solverCreator->create(name + "->solver");
286
    }
287

288
289
    void createFileWriter()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
290
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", *meshView, componentNames);
291
    }
292

Praetorius, Simon's avatar
Praetorius, Simon committed
293

294
  protected: // sub-methods to assemble DOFMatrix
295

296
297
    template <class LhsData, class RhsData, class Elements>
    void assemble(LhsData lhs, RhsData rhs, Elements const& elements);
298

299
    template <class RowView, class ColView>
300
    bool getElementMatrix(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
301
302
303
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
                          std::list<std::shared_ptr<OperatorType>>& operators,
304
                          std::list<double*> const& factors);
305

306
    template <class RowView>
307
    bool getElementVector(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
308
309
                          ElementVector& elementVector,
                          std::list<std::shared_ptr<OperatorType>>& operators,
310
                          std::list<double*> const& factors);
311
312


313
314
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
Praetorius, Simon's avatar
Praetorius, Simon committed
315
316
317
                          RowIndexSet const& rowIndexSet,
                          ColIndexSet const& colIndexSet,
                          ElementMatrix const& elementMatrix);
318

319
320
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
Praetorius, Simon's avatar
Praetorius, Simon committed
321
322
                          IndexSet const& indexSet,
                          ElementVector const& elementvector);
323
324


325
  private:
326
327
328
329
330
331
    /// 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;
332

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

336
    /// Name of the mesh
337
    std::string meshName;
338

339
    /// A gridView object
Praetorius, Simon's avatar
Praetorius, Simon committed
340
    std::shared_ptr<MeshView> meshView;
341

342
343
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
344

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

348
    /// A FileWriter object
Praetorius, Simon's avatar
Praetorius, Simon committed
349
    std::shared_ptr<FileWriter<Traits>> filewriter;
350

351
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
352
    std::shared_ptr<LinearSolverType> linearSolver;
353

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

357
    /// A block-vector with the solution components
Praetorius, Simon's avatar
Praetorius, Simon committed
358
    std::shared_ptr<SystemVectorType> solution;
359
360

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


365
366
    template <class T>
    using MatrixEntries = std::map< std::pair<int,int>, std::list<T> >;
367

368
369
    template <class T>
    using VectorEntries = std::map< int, std::list<T> >;
370
371
372


    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
373
    /// each matrix block
Praetorius, Simon's avatar
Praetorius, Simon committed
374
    MatrixEntries<std::shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
375
376

    /// A map (i,j) -> list<OperatorType> string the differential operators for
377
    /// each matrix block
Praetorius, Simon's avatar
Praetorius, Simon committed
378
    MatrixEntries<std::shared_ptr<OperatorType>> matrixOperators;
379
    MatrixEntries<double*>                  matrixFactors;
380
381
    std::map< std::pair<int,int>, bool >    matrixAssembled; // if false, do reassemble
    std::map< std::pair<int,int>, bool >    matrixChanging;  // if true, or vectorAssembled false, do reassemble
382
383

    /// A map (i) -> list<OperatorType> string the differential operators for
384
    /// each rhs block
Praetorius, Simon's avatar
Praetorius, Simon committed
385
    VectorEntries<std::shared_ptr<OperatorType>> vectorOperators;
386
    VectorEntries<double*>                  vectorFactors;
387
388
    std::map< int, bool >                   vectorAssembled; // if false, do reassemble
    std::map< int, bool >                   vectorChanging;  // if true, or vectorAssembled false, do reassemble
389

Praetorius, Simon's avatar
Praetorius, Simon committed
390
    template <int R, int C>
391
    struct MatrixData
392
    {
393
      using DOFMatrixType =
Praetorius, Simon's avatar
Praetorius, Simon committed
394
        std::tuple_element_t<C, std::tuple_element_t<R, typename SystemMatrixType::DOFMatrices>>;
395

396
      DOFMatrixType&                       matrix;
Praetorius, Simon's avatar
Praetorius, Simon committed
397
      std::list<std::shared_ptr<OperatorType>>& operators;
398
399
      std::list<double*> const&            factors;
      bool                                 assemble;
400

Praetorius, Simon's avatar
Praetorius, Simon committed
401
      std::pair<int,int> row_col = {R, C};
402
    };
403

Praetorius, Simon's avatar
Praetorius, Simon committed
404
    template <int I>
405
406
    struct VectorData
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
407
      using DOFVectorType = std::tuple_element_t<I, typename SystemVectorType::DOFVectors>;
408

409
      DOFVectorType&                       vector;
Praetorius, Simon's avatar
Praetorius, Simon committed
410
      std::list<std::shared_ptr<OperatorType>>& operators;
411
412
      std::list<double*> const&            factors;
      bool                                 assemble;
413

Praetorius, Simon's avatar
Praetorius, Simon committed
414
      int row = I;
415
    };
416
  };
417
418


419
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
420
421
422
423
//   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>>;
424
#endif
425

426
  namespace Impl
427
428
  {
    template <class ProblemStatType>
Praetorius, Simon's avatar
Praetorius, Simon committed
429
430
431
    struct ProblemStat
        : public ProblemStatType
        , public StandardProblemIteration
432
433
434
435
436
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
437
438
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
439
440
      {}

441
      /**
442
443
444
445
446
447
       * \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.
       *
448
449
       * Implementation of \ref StandardProblemIteration::oneIteration.
       **/
450
451
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
452
        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
453
454
455
456
457
458
459
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

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

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

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

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

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

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

476
      /// Implementation of \ref ProblemStatBase::markElements.
477
      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
478
    };
479

480
  } // end namespace Impl
481

482
  template <class Traits>
483
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
484

485
} // end namespace AMDiS
486
487
488

#include "ProblemStat.inc.hpp"