RosenbrockStationary.h 5.99 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * 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 Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27
28

/** \file RosenbrockStationary.h */

#ifndef AMDIS_ROSENBROCKSTATIONARY_H
#define AMDIS_ROSENBROCKSTATIONARY_H

#include "AMDiS_fwd.h"
29
#include "ProblemStat.h"
30
#include "SystemVector.h"
31
#include "time/RosenbrockMethod.h"
32

33
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
34
#include "parallel/ParallelProblemStat.h"
35
#endif
36

37
38
namespace AMDiS {

39
40
  struct RosenbrockBoundary 
  {
41
42
43
44
45
46
    RosenbrockBoundary(AbstractFunction<double, WorldVector<double> >* fct_,
		       AbstractFunction<double, WorldVector<double> >* fctDt_,
		       DOFVector<double> *vec_,
		       int row_, int col_)
      : fct(fct_), fctDt(fctDt_), vec(vec_), row(row_), col(col_) {}
    
47
    AbstractFunction<double, WorldVector<double> > *fct;
48
    AbstractFunction<double, WorldVector<double> > *fctDt; // dt[f](t,x)
49
50
51
52
53
54
55
    
    DOFVector<double> *vec;

    int row;
    int col;
  };

56
57
58
59
60
61
62
63
64
65
66
67
  /// Realization of a Rosenbrock time-discretization of M*d_t(X) = F[x]
  /** 
   * 1/(tau*gamma) M*Y_i^k  - J(t^k, X^k)[Y_i^k]  =  F[t_i^k, X_i^k] + sum_{j=1}^{i-1} c_ij / tau * M * Y_j^k
   *
   * with stageSolution[i]:
   * X_i^k = X^k + sum_{j=1}^{i-1} a_ij * Y_j^k
   *
   * oldTime: t^k, stageTime: t_i^k
   * 
   * and new solution
   * X^{k+1} = X^k + sum_{j=1}^{s} m_j * Y_j^k
   **/
68
  class RosenbrockStationary : public ProblemStat
69
70
  {
  public:
71
    RosenbrockStationary(std::string name, int componentShift_ = 0)
72
      : ProblemStat(name),
73
	first(true),
74
	componentShift(componentShift_),
75
	minusOne(-1.0),
76
77
	stageTime(0.0),
	oldTime(0.0),
78
79
80
81
82
	tauPtr(NULL),
	tauGamma(NULL),
	minusTauGamma(NULL),	
	invTauGamma(NULL),
	minusInvTauGamma(NULL)
83
    {}
84
    
85

86
87
88
89
90
    Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) override;
    
    virtual Flag stageIteration(AdaptInfo *adaptInfo, Flag flag, 
				bool asmMatrix, bool asmVector);
    
91
    virtual void estimateTimeError(AdaptInfo* adaptInfo);
92
93
94
95
96
97
    
    /// update solution vector and oldTime value
    void acceptTimestep(AdaptInfo* adaptInfo);
    
    /// Add operators of function F
    virtual void addOperator(Operator &op, int row, int col, 
98
			     double *factor = NULL, double *estFactor = NULL);
99
100
101

    /// Add operators of jacobian J = d_X(F)
    virtual void addJacobianOperator(Operator &op, int row, int col, 
102
				     double *factor = NULL, double *estFactor = NULL);
103
104
105
106
107
108

    virtual void addTimeOperator(int i, int j);

    // getting methods
    // _________________________________________________________________________
    
109
110
111
112
113
    DOFVector<double>* getUnVec(int i)
    {
      return unVec->getDOFVector(i);
    }

114
    DOFVector<double>* getStageSolution(int i)
115
    {
116
      return stageSolution->getDOFVector(i);
117
118
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
119
120
121
122
123
    DOFVector<double>* getTimeRhsVec(int i)
    {
      return timeRhsVec->getDOFVector(i);
    }

124
125
126
127
    double* getTauGamma()
    {
      return tauGamma;
    }
128

129
    double* getTau()
130
    {
131
      return tauPtr;
132
133
    }

134
     double* getMinusTauGamma()
Thomas Witkowski's avatar
Thomas Witkowski committed
135
    {
136
      return minusTauGamma;
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
    }

139
140
141
142
    double* getInvTauGamma()
    {
      return invTauGamma;
    }
143

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    double* getMinusInvTauGamma()
    {
      return minusInvTauGamma;
    }
    
    double* getStageTime()
    {
      return &stageTime;
    }
    
    double* getOldTime()
    {
      return &oldTime;
    }
    
    double* getTauGammaI()
    {
      return &tauGammaI;
    }
163

164
165
166
167
168
169
170
171
    // setting methods
    // _________________________________________________________________________
    
    void setRosenbrockMethod(RosenbrockMethod *method)
    {
      rm = method;
      init();
    }
172
173
174
175
176
177

    void setTau(double *ptr)
    {
      tauPtr = ptr;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
181
182
183
184
    void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3)
    {
      tauGamma = ptr0;
      minusTauGamma = ptr1;
      invTauGamma = ptr2;
      minusInvTauGamma = ptr3;
    }    
185
186
    
    void setOldTime(double t)
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    {
188
      oldTime = t;
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    }
190
191
    
    void setStageTime(double t)
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    {
193
      stageTime = t;
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
196
    }


197
198
    // boundary conditions
    // _________________________________________________________________________
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
   /// Adds a Dirichlet boundary condition, where the rhs is given by an 
201
202
    /// abstract function.
    void addDirichletBC(BoundaryType type, int row, int col,
203
			AbstractFunction<double, WorldVector<double> > *fct) override;
204

205
206
207
208
    void addDirichletBC(BoundaryType type, int row, int col,
			AbstractFunction<double, WorldVector<double> > *fct,
			AbstractFunction<double, WorldVector<double> > *fctDt);
			
209
210
211
    /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
    /// vector.
    void addDirichletBC(BoundaryType type, int row, int col,
212
			DOFVector<double> *vec) override
213
214
215
216
217
218
    {
      FUNCNAME("RosenbrockStationary::addDirichletBC()");

      ERROR_EXIT("Not yet supported!\n");
    }

219
220
  protected:
    void init();
221

222
    void reset()
223
    {
224
225
226
227
228
      first = true;
      stageSolution->set(0.0);
      unVec->set(0.0);
      for (int i = 0; i < rm->getStages(); i++)
	stageSolutions[i]->set(0.0);
229
    }
230

231
  protected:    
232
    RosenbrockMethod *rm;
233

234
    SystemVector *stageSolution, *unVec, *timeRhsVec, *newUn, *tmp, *lowSol;
235
236
237

    std::vector<SystemVector*> stageSolutions;

238
239
    bool first;

240
241
    int componentShift;
    
242
    double minusOne;
243
244
245
    double stageTime; // t_n + alpha_i*tau
    double oldTime;
    double tauGammaI; // tau*gamma_i
246

247
    double *tauPtr;
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    double *tauGamma, *minusTauGamma, *invTauGamma, *minusInvTauGamma;
249
250

    std::vector<RosenbrockBoundary> boundaries;
251
252
253
254
255
  };

}

#endif // AMDIS_ROSENBROCKSTATIONARY_H