NavierStokes_Chorin.h 6.58 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
15
 *
16
 ******************************************************************************/
17
18
19
20

#ifndef NAVIER_STOKES_CHORIN_H
#define NAVIER_STOKES_CHORIN_H

21
#include "AMDiS.h"
22
23
24
25
26
27
28
29
30
31
32
33
#include "Helpers.h"
#include "POperators.h"
#include "BoundaryFunctions.h"
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>


#define LINKS 1
#define RECHTS 2
#define UNTEN 3
#define OBEN 4

34
namespace AMDiS { namespace base_problems {
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
60
61

/** \ingroup NavierStokes_Chorin
 * \brief
 * Navier-Stokes problem, using projection method {
 *       a.) Burgers equation (predictortProb)
 *       b.) Pressure Poisson problem
 *       c.) Burgers equation including pressure (predictorProb)
 *       d.) Pressure Poisson problem
 *       e.) Velocity update (correctorProb)
 * }
 */
class NavierStokes_Chorin : public ProblemIterationInterface,
                            public ProblemInstatBase
{
public:

  NavierStokes_Chorin(const std::string &name_, WorldVector<AbstractFunction<double, WorldVector<double> >* > initialVelocityFcts_);
  ~NavierStokes_Chorin();

  /// Initialisation of the problem.
  virtual void initialize(Flag initFlag,
                  ProblemStatType *adoptProblem = NULL,
                  Flag adoptFlag = INIT_NOTHING);

  virtual void initTimeInterface();
  virtual void solveInitialProblem(AdaptInfo *adaptInfo);

62
  virtual void transferInitialSolution(AdaptInfo *adaptInfo);
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

  /// Implementation of ProblemTimeInterface::initTimestep().
  virtual void initTimestep(AdaptInfo *adaptInfo);

  /// Implementation of ProblemTimeInterface::closeTimestep().
  virtual void closeTimestep(AdaptInfo *adaptInfo);

  virtual void beginIteration(AdaptInfo *adaptInfo);
  virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo= FULL_ITERATION);
  virtual void endIteration(AdaptInfo *adaptInfo);

  void clean();

  // === getting methods ===

  SystemVector *getSolution()
  {
    return correctorProb->getSolution();
  }

  double* getPressureStep() { return &pressureStep; }

85
  inline Mesh* getMesh(int comp = 0)
86
87
  { FUNCNAME("NavierStokes_Chorin::getMesh()");

88
    return predictorProb->getMesh(comp);
89
90
  }

91
  std::string getName() const { return name; };
92

93
  int getNumProblems() const { return 3; };
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
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

  ProblemStatType *getProblem(int number= 0)
  { FUNCNAME("NavierStokes_Chorin::getProblem()");

    if (number == 0)
      return predictorProb;
    else if (number == 1)
      return pressureProb;
    else if (number == 2)
      return correctorProb;
    else
      throw(std::runtime_error("problem with given number does not exist"));
  };

  ProblemStatType *getProblem(std::string name)
  { FUNCNAME("NavierStokes_Chorin::getProblem()");

    if (name == "correctorProb")
      return correctorProb;
    else if (name == "predictorProb")
      return predictorProb;
    else if (name == "pressureProb")
      return pressureProb;
    else
      throw(std::runtime_error("problem with given name '" + name + "' does not exist"));
  };

  /// setting methods

  void setAssembleMatrixOnlyOnce_butTimestepChange(int i, int j)
  { FUNCNAME("NavierStokes_Chorin::setAssembleMatrixOnlyOnce_butTimestepChange()");

    fixedMatrixTimestep.push_back(std::make_pair(i,j));
  }

  void setNumberOfTimesteps(int nTimesteps_)
  { FUNCNAME("NavierStokes_Chorin::setNumberOfTimesteps()");

    nTimesteps= nTimesteps_;
  }

  void writeFiles(AdaptInfo *adaptInfo, bool force);
  void serialize(std::ostream&) {};
  void deserialize(std::istream&) {};

protected:

  virtual Flag solvePredictorProb(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION, bool fixedMatrix = false);
  virtual Flag solvePressureProb(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION, bool fixedMatrix = false);
  virtual Flag solveCorrectorProb(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION, bool fixedMatrix = false);

  virtual void fillOperators();
  virtual void fillBoundaryConditions();
147

148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
  virtual void fillPredictorProb();
  virtual void fillPressureProb();
  virtual void fillCorrectorProb();
  virtual void addLaplaceTerm(int i);

private:

  void normalizePressure(DOFVector<double> *p, double mean);

  template <typename Matrix>
  void applySingularPertubation(Matrix& m)
  {
	  using namespace mtl;
	  typedef typename mtl::Collection<Matrix>::value_type value_type;
	  int nnz= m.nnz_local(m.num_rows() - 1);
	  matrix::inserter<Matrix, update_times<value_type> > ins(m, nnz);
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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
	  ins[m.num_rows() - 1][m.num_rows() - 1] << 1.0 + 1.e-5;
  }


  template <typename Matrix, typename Vector>
  void applyDbc(Matrix& m, Vector& vec, DegreeOfFreedom idx, bool setDiagonal, double value)
  {
    using namespace mtl;
    typename Matrix::size_type idx_= idx;

    typename traits::row<Matrix>::type r(m);
    typename traits::col<Matrix>::type c(m);
    typename traits::value<Matrix>::type v(m);
    typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
    typedef typename traits::range_generator<tag::nz, c_type>::type  ic_type;

    for (c_type cursor(begin<tag::row>(m)+idx_), cend(begin<tag::row>(m) + idx_ + 1); cursor != cend; ++cursor) {
        for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
          v(*icursor, (r(*icursor) == c(*icursor) && setDiagonal ? 1.0 : 0.0));
        }
        break;
    }

    if (setDiagonal)
      vec[idx]= value;
  }

protected:

  ProblemStatType *predictorProb;
  ProblemStatType *pressureProb;
  ProblemStatType *correctorProb;

  ProblemStatType *nsProb;

  DOFVector<double> *pressure;
  WorldVector<DOFVector<double>*> boundaryDOF;
  AbstractFunction<double, WorldVector<double> > *initialX, *initialY;
  WorldVector<AbstractFunction<double, WorldVector<double> >* > initialVelocityFcts;
  WorldVector<DOFVector<double>*> initialVel;

  bool forceDBC;
  bool calcPressure;
  bool initialVelocityIsSet;

  unsigned dim;

  int simpleAlg;
  int nTimesteps;
  int oldMeshChangeIdx;
  int poissonPertubation;
  int laplaceType;
  int nonLinTerm;

  double pressureStep;
  double oldTimestep;
  double viscosity;
  double c;

  WorldVector<double> force;

  std::vector<std::pair<int,int> > fixedMatrixTimestep;

private:

};
231
232
233

} }

234
#endif // NAVIER_STOKES_CHORIN_H