ProblemStat.hpp 15.2 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
38
39
40
  template <class Traits>
  class ProblemStatSeq : public ProblemStatBase
  {
    using Self = ProblemStatSeq;
41

42
  public: // typedefs and static constants
43

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

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

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


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

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

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


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

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

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

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

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


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

104
105

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

112
    void addMatrixOperator(shared_ptr<OperatorType> op,
113
                           int i, int j,
114
                           double* factor = NULL,
115
                           double* estFactor = NULL);
116
117

    void addMatrixOperator(std::map< std::pair<int,int>,
118
                                     shared_ptr<OperatorType> > ops);
119
    void addMatrixOperator(std::map< std::pair<int,int>,
120
121
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
122
123
124


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

    void addVectorOperator(shared_ptr<OperatorType> op,
132
                           int i,
133
                           double* factor = NULL,
134
                           double* estFactor = NULL);
135
136

    void addVectorOperator(std::map< int,
137
                                     shared_ptr<OperatorType> > ops);
138
    void addVectorOperator(std::map< int,
139
140
                                     std::pair<shared_ptr<OperatorType>, bool> > ops);
    /** @} */
141
142
143

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
144
145
    void addDirichletBC(Predicate const& predicate,
			int row, int col,
146
147
			Values const& values);

148

149
    /// Implementation of \ref ProblemStatBase::solve
150
151
    virtual void solve(AdaptInfo& adaptInfo,
		       bool createMatrixData = true,
152
		       bool storeMatrixData = false) override;
153
154


155
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
156
157
158
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
				   Flag flag,
				   bool asmMatrix = true,
159
				   bool asmVector = true) override;
160

161
162

    /// Writes output files.
163
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
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

199
200
201
202
203
    void setSolver(std::shared_ptr<LinearSolverType> const& solver_)
    {
      linearSolver = solver_;
    }

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

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

      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);
254
      meshView = 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()
    {
265
      feSpaces = make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
266
    }
267

268
269
    void createMatricesAndVectors()
    {
270
271
      systemMatrix = make_shared<SystemMatrixType>(*feSpaces);
      solution = make_shared<SystemVectorType>(*feSpaces, componentNames);
272

273
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
274
      rhs = 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()
    {
290
291
      filewriter = make_shared<FileWriter<Traits>>(name + "->output",
                                                        *meshView,
292
                                                        componentNames);
293
    }
294

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

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

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

307
    template <class RowView>
308
    bool getElementVector(RowView const& rowView,
309
			  ElementVector& elementVector,
310
			  std::list<shared_ptr<OperatorType>>& operators,
311
                          std::list<double*> const& factors);
312
313


314
315
316
317
318
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
			  RowIndexSet const& rowIndexSet,
			  ColIndexSet const& colIndexSet,
			  ElementMatrix const& elementMatrix);
319

320
321
322
323
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
			  IndexSet const& indexSet,
			  ElementVector const& elementvector);
324

325
  private: // some internal methods
326

327
    template <size_t I = 0>
328
    typename FeSpace<I>::NodeFactory&
329
330
    getNodeFactory(const index_<I> = index_<I>())
    {
331
      return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
332
    }
333
334


335
  private:
336
337
338
339
340
341
    /// 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;
342

343
    /// Mesh of this problem.
344
    shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
345

346
    /// Name of the mesh
347
    std::string meshName;
348

349
    /// A gridView object
350
    shared_ptr<MeshView> meshView;
351

352
353
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
354

355
    /// FE spaces of this problem.
356
    shared_ptr<FeSpaces> feSpaces; // eventuell const
357

358
    /// A FileWriter object
359
    shared_ptr<FileWriter<Traits>> filewriter;
360

361
    /// An object of the linearSolver Interface
362
    shared_ptr<LinearSolverType> linearSolver;
363

364
    /// A block-matrix that is filled during assembling
365
    shared_ptr<SystemMatrixType> systemMatrix;
366

367
    /// A block-vector with the solution components
368
    shared_ptr<SystemVectorType> solution;
369
370

    /// A block-vector (load-vector) corresponding to the right.hand side
371
    /// of the equation, filled during assembling
372
    shared_ptr<SystemVectorType> rhs;
373
374


375
376
    template <class T>
    using MatrixEntries = std::map< std::pair<int,int>, std::list<T> >;
377

378
379
    template <class T>
    using VectorEntries = std::map< int, std::list<T> >;
380
381
382


    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
383
    /// each matrix block
384
    MatrixEntries<shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
385
386

    /// A map (i,j) -> list<OperatorType> string the differential operators for
387
    /// each matrix block
388
389
    MatrixEntries<shared_ptr<OperatorType>> matrixOperators;
    MatrixEntries<double*>                  matrixFactors;
390
391
    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
392
393

    /// A map (i) -> list<OperatorType> string the differential operators for
394
    /// each rhs block
395
396
    VectorEntries<shared_ptr<OperatorType>> vectorOperators;
    VectorEntries<double*>                  vectorFactors;
397
398
    std::map< int, bool >                   vectorAssembled; // if false, do reassemble
    std::map< int, bool >                   vectorChanging;  // if true, or vectorAssembled false, do reassemble
399

400
    template <int r, int c>
401
    struct MatrixData
402
    {
403
      using DOFMatrixType =
404
        std::tuple_element_t<c, std::tuple_element_t<r, typename SystemMatrixType::DOFMatrices>>;
405

406
407
408
409
      DOFMatrixType&                       matrix;
      std::list<shared_ptr<OperatorType>>& operators;
      std::list<double*> const&            factors;
      bool                                 assemble;
410

411
412
      std::pair<int,int> row_col = {r, c};
    };
413

414
415
416
    template <int r>
    struct VectorData
    {
417
      using DOFVectorType =
418
        std::tuple_element_t<r, typename SystemVectorType::DOFVectors>;
419

420
421
422
423
      DOFVectorType&                       vector;
      std::list<shared_ptr<OperatorType>>& operators;
      std::list<double*> const&            factors;
      bool                                 assemble;
424

425
426
      int row = r;
    };
427
  };
428
429


430
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
431
432
433
434
//   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>>;
435
#endif
436

437
  namespace Impl
438
439
440
441
442
443
444
445
446
  {
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType,
                         public StandardProblemIteration
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
447
448
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
449
450
      {}

451
      /**
452
453
454
455
456
457
       * \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.
       *
458
459
       * Implementation of \ref StandardProblemIteration::oneIteration.
       **/
460
461
      virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
      {
462
        for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
463
464
465
466
467
468
469
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

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

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

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

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

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

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

486
      /// Implementation of \ref ProblemStatBase::markElements.
487
      virtual Flag markElements(AdaptInfo& adaptInfo) override { return {0}; }
488
    };
489

490
  } // end namespace Impl
491

492
  template <class Traits>
493
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
494

495
} // end namespace AMDiS
496
497
498

#include "ProblemStat.inc.hpp"