DOpE
fractional_step_theta_problem.h
Go to the documentation of this file.
1 
24 #ifndef _FRACTIONALSTEPTHETAPROBLEM_H_
25 #define _FRACTIONALSTEPTHETAPROBLEM_H_
26 
27 #include "initialproblem.h"
28 #include "primal_ts_base.h"
29 
30 namespace DOpE
31 {
57  template<typename OPTPROBLEM, typename SPARSITYPATTERN, typename VECTOR,
58  int dopedim, int dealdim, template<int, int> class FE = dealii::FESystem,
59  template<int, int> class DH = dealii::DoFHandler>
60  class FractionalStepThetaProblem : public PrimalTSBase<OPTPROBLEM,
61  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>
62  {
63  public:
64  FractionalStepThetaProblem(OPTPROBLEM& OP) :
65  PrimalTSBase<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim,
66  FE, DH>(OP)
67  {
68  _fs_theta = 1.0 - std::sqrt(2.0) / 2.0;
69  _fs_theta_prime = 1.0 - 2.0 * _fs_theta;
70  _fs_alpha = (1.0 - 2.0 * _fs_theta) / (1.0 - _fs_theta);
71  _fs_beta = 1.0 - _fs_alpha;
72  _initial_problem = NULL;
73  }
74 
76  {
77  if (_initial_problem != NULL)
78  delete _initial_problem;
79  }
80 
81  /******************************************************/
82 
83  std::string
85  {
86  return "Fractional-Step-Theta";
87  }
88  /******************************************************/
89 
91  FractionalStepThetaProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR,
92  dopedim, dealdim, FE, DH>, VECTOR, dealdim>&
94  {
95  if (_initial_problem == NULL)
96  {
97  _initial_problem = new InitialProblem<
98  FractionalStepThetaProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR,
99  dopedim, dealdim, FE, DH>, VECTOR, dealdim>(*this);
100  }
101  return *_initial_problem;
102  }
103 
104  /******************************************************/
105  FractionalStepThetaProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim,
106  dealdim, FE, DH>&
108  {
109  return *this;
110  }
111  /******************************************************/
112 
113  template<typename EDC>
114  void
115  ElementEquation(const EDC& edc,
116  dealii::Vector<double> &local_vector, double scale, double /*scale_ico*/)
117  {
118  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
119  {
120  dealii::Vector<double> tmp(local_vector);
121 
122  tmp = 0.0;
123  this->GetProblem().ElementEquation(edc, tmp,
124  scale * _fs_alpha
125  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
126  scale
127  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
128  local_vector += tmp;
129 
130  tmp = 0.0;
131  this->GetProblem().ElementTimeEquation(edc, tmp, scale / (_fs_theta));
132  local_vector += tmp;
133 
134  this->GetProblem().ElementTimeEquationExplicit(edc, local_vector,
135  scale / (_fs_theta));
136  }
137  else if (this->GetPart() == "Old_for_1st_cycle")
138  {
139  dealii::Vector<double> tmp(local_vector);
140  tmp = 0.0;
141 
142  // The explicit parts with old_time_values; e.g. for fluid problems: laplace, convection, etc.
143  this->GetProblem().ElementEquation(edc, tmp,
144  scale * _fs_beta
145  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
146  0.);
147  local_vector += tmp;
148 
149  this->GetProblem().ElementTimeEquation(edc, local_vector,
150  (-1) * scale / (_fs_theta));
151  }
152  else if (this->GetPart() == "Old_for_3rd_cycle")
153  {
154  dealii::Vector<double> tmp(local_vector);
155  tmp = 0.0;
156 
157  // The explicit parts with old_time_values; e.g. for fluid problems: laplace, convection, etc.
158  this->GetProblem().ElementEquation(edc, tmp,
159  scale * _fs_beta
160  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
161  0.);
162  local_vector += tmp;
163 
164  this->GetProblem().ElementTimeEquation(edc, local_vector,
165  (-1) * scale / (_fs_theta));
166  }
167  else if (this->GetPart() == "New_for_2nd_cycle")
168  {
169  dealii::Vector<double> tmp(local_vector);
170 
171  tmp = 0.0;
172  this->GetProblem().ElementEquation(edc, tmp,
173  scale * _fs_beta
174  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
175  scale
176  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
177  local_vector += tmp;
178 
179  tmp = 0.0;
180  this->GetProblem().ElementTimeEquation(edc, tmp,
181  scale / (_fs_theta_prime));
182  local_vector += tmp;
183 
184  this->GetProblem().ElementTimeEquationExplicit(edc, local_vector,
185  scale / (_fs_theta_prime));
186  }
187  else if (this->GetPart() == "Old_for_2nd_cycle")
188  {
189  dealii::Vector<double> tmp(local_vector);
190  tmp = 0.0;
191 
192  // The explicit parts with old_time_values; e.g. for fluid problems: laplace, convection, etc.
193  this->GetProblem().ElementEquation(edc, tmp,
194  scale * _fs_alpha
195  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
196  0.);
197  local_vector += tmp;
198 
199  this->GetProblem().ElementTimeEquation(edc, local_vector,
200  (-1) * scale / (_fs_theta_prime));
201  }
202  else
203  {
204  abort();
205  }
206  }
207 
208  /******************************************************/
209 
210  template<typename EDC>
211  void
212  ElementRhs(const EDC& edc,
213  dealii::Vector<double> &local_vector, double scale)
214  {
215  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
216  {
217 
218  }
219  else if (this->GetPart() == "Old_for_1st_cycle"
220  || this->GetPart() == "Old_for_3rd_cycle")
221  {
222  this->GetProblem().ElementRhs(edc, local_vector,
223  scale
224  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
225  }
226  else if (this->GetPart() == "New_for_2nd_cycle")
227  {
228  this->GetProblem().ElementRhs(edc, local_vector,
229  scale
230  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
231  }
232  else if (this->GetPart() == "Old_for_2nd_cycle")
233  {
234  }
235  else
236  {
237  abort();
238  }
239  }
240  /******************************************************/
241 
242  void
244  const std::map<std::string, const dealii::Vector<double>*> &param_values,
245  const std::map<std::string, const VECTOR*> &domain_values,
246  VECTOR& rhs_vector, double scale)
247  {
248  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
249  {
250 
251  }
252  else if (this->GetPart() == "Old_for_1st_cycle"
253  || this->GetPart() == "Old_for_3rd_cycle")
254  {
255  this->GetProblem().PointRhs(param_values, domain_values, rhs_vector,
256  scale
257  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
258  }
259  else if (this->GetPart() == "New_for_2nd_cycle")
260  {
261  this->GetProblem().PointRhs(param_values, domain_values, rhs_vector,
262  scale
263  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
264  }
265  else if (this->GetPart() == "Old_for_2nd_cycle")
266  {
267  }
268  else
269  {
270  abort();
271  }
272  }
273 
274  /******************************************************/
275 
276  template<typename EDC>
277  void
278  ElementMatrix(const EDC& edc,
279  dealii::FullMatrix<double> &local_matrix)
280  {
281  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
282  {
283  dealii::FullMatrix<double> m(local_matrix);
284  m = 0.;
285  this->GetProblem().ElementMatrix(edc, m,
286  _fs_alpha
287  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
288  1.
289  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
290  local_matrix.add(1.0, m);
291 
292  m = 0.;
293  this->GetProblem().ElementTimeMatrix(edc, m);
294  local_matrix.add(1.0 / (_fs_theta), m);
295 
296  m = 0.;
297  this->GetProblem().ElementTimeMatrixExplicit(edc, m);
298  local_matrix.add(1.0 / (_fs_theta), m);
299  }
300  else if (this->GetPart() == "New_for_2nd_cycle")
301  {
302  dealii::FullMatrix<double> m(local_matrix);
303  m = 0.;
304  this->GetProblem().ElementMatrix(edc, local_matrix,
305  _fs_beta
306  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
307  1.
308  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
309  local_matrix.add(1.0, m);
310 
311  m = 0.;
312  this->GetProblem().ElementTimeMatrix(edc, m);
313  local_matrix.add(1.0 / (_fs_theta_prime), m);
314 
315  m = 0.;
316  this->GetProblem().ElementTimeMatrixExplicit(edc, m);
317  local_matrix.add(1.0 / (_fs_theta_prime), m);
318  }
319  }
320 
321  /******************************************************/
322 
323  template<typename FDC>
324  void
325  FaceEquation(const FDC& fdc,
326  dealii::Vector<double> &local_vector, double scale, double)
327  {
328  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
329  {
330  this->GetProblem().FaceEquation(fdc, local_vector,
331  scale * _fs_alpha
332  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
333  scale
334  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
335 
336  }
337  else if ((this->GetPart() == "Old_for_1st_cycle")
338  || (this->GetPart() == "Old_for_3rd_cycle"))
339  {
340  this->GetProblem().FaceEquation(fdc, local_vector,
341  scale * _fs_beta
342  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
343  0);
344  }
345  else if (this->GetPart() == "New_for_2nd_cycle")
346  {
347  this->GetProblem().FaceEquation(fdc, local_vector,
348  scale * _fs_beta
349  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
350  scale
351  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
352  }
353  else if (this->GetPart() == "Old_for_2nd_cycle")
354  {
355  this->GetProblem().FaceEquation(fdc, local_vector,
356  scale * _fs_alpha
357  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
358  0.);
359  }
360  else
361  {
362  abort();
363  }
364 
365  }
366 
367  /******************************************************/
368 
369  template<typename FDC>
370  void
371  InterfaceEquation(const FDC& fdc,
372  dealii::Vector<double> &local_vector, double scale,
373  double /*scale_ico*/)
374  {
375  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
376  {
377  this->GetProblem().InterfaceEquation(fdc, local_vector,
378  scale * _fs_alpha
379  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
380  scale
381  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
382 
383  }
384  else if ((this->GetPart() == "Old_for_1st_cycle")
385  || (this->GetPart() == "Old_for_3rd_cycle"))
386  {
387  this->GetProblem().InterfaceEquation(fdc, local_vector,
388  scale * _fs_beta
389  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
390  0);
391  }
392  else if (this->GetPart() == "New_for_2nd_cycle")
393  {
394  this->GetProblem().InterfaceEquation(fdc, local_vector,
395  scale * _fs_beta
396  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
397  scale
398  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
399  }
400  else if (this->GetPart() == "Old_for_2nd_cycle")
401  {
402  this->GetProblem().InterfaceEquation(fdc, local_vector,
403  scale * _fs_alpha
404  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
405  0.);
406  }
407  else
408  {
409  abort();
410  }
411 
412  }
413 
414  /******************************************************/
415 
416  template<typename FDC>
417  void
418  FaceRhs(const FDC& fdc,
419  dealii::Vector<double> &local_vector, double scale = 1.)
420  {
421  this->GetProblem().FaceRhs(fdc, local_vector,
422  scale
423  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
424  }
425 
426  /******************************************************/
427 
428  template<typename FDC>
429  void
430  FaceMatrix(const FDC& fdc,
431  dealii::FullMatrix<double> &local_matrix)
432  {
433  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
434  {
435  dealii::FullMatrix<double> m(local_matrix);
436  m = 0.;
437  this->GetProblem().FaceMatrix(fdc, m,
438  _fs_alpha
439  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
440  1.
441  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
442  local_matrix.add(1., m);
443  }
444  else if (this->GetPart() == "New_for_2nd_cycle")
445  {
446  dealii::FullMatrix<double> m(local_matrix);
447  m = 0.;
448  this->GetProblem().FaceMatrix(fdc, m,
449  _fs_beta
450  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
451  1.
452  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
453  local_matrix.add(1.0, m);
454  }
455 
456  }
457 
458  /******************************************************/
459 
460  template<typename FDC>
461  void
462  InterfaceMatrix(const FDC& fdc,
463  dealii::FullMatrix<double> &local_matrix)
464  {
465  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
466  {
467  dealii::FullMatrix<double> m(local_matrix);
468  m = 0.;
469  this->GetProblem().InterfaceMatrix(fdc, m,
470  _fs_alpha
471  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
472  1.
473  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
474  local_matrix.add(1.0, m);
475  }
476  else if (this->GetPart() == "New_for_2nd_cycle")
477  {
478  dealii::FullMatrix<double> m(local_matrix);
479  m = 0.;
480  this->GetProblem().InterfaceMatrix(fdc, m,
481  _fs_beta
482  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
483  1.
484  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
485  local_matrix.add(1.0, m);
486  }
487 
488  }
489 
490  /******************************************************/
491 
492  template<typename FDC>
493  void
494  BoundaryEquation(const FDC& fdc,
495  dealii::Vector<double> &local_vector, double scale, double)
496  {
497  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
498  {
499  this->GetProblem().BoundaryEquation(fdc, local_vector,
500  scale * _fs_alpha
501  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
502  scale
503  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
504  }
505  else if ((this->GetPart() == "Old_for_1st_cycle")
506  || (this->GetPart() == "Old_for_3rd_cycle"))
507  {
508  this->GetProblem().BoundaryEquation(fdc, local_vector,
509  scale * _fs_beta
510  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
511  0.);
512  }
513  else if (this->GetPart() == "New_for_2nd_cycle")
514  {
515  this->GetProblem().BoundaryEquation(fdc, local_vector,
516  scale * _fs_beta
517  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
518  scale
519  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
520  }
521  else if (this->GetPart() == "Old_for_2nd_cycle")
522  {
523  this->GetProblem().BoundaryEquation(fdc, local_vector,
524  scale * _fs_alpha
525  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
526  0);
527  }
528  else
529  {
530  abort();
531  }
532 
533  }
534 
535  /******************************************************/
536 
537  template<typename FDC>
538  void
539  BoundaryRhs(const FDC& fdc,
540  dealii::Vector<double> &local_vector, double scale)
541  {
542  this->GetProblem().BoundaryRhs(fdc, local_vector,
543  scale
544  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
545  }
546 
547  /******************************************************/
548 
549  template<typename FDC>
550  void
551  BoundaryMatrix(const FDC& fdc,
552  dealii::FullMatrix<double> &local_matrix)
553  {
554  if (this->GetPart() == "New_for_1st_and_3rd_cycle")
555  {
556  dealii::FullMatrix<double> m(local_matrix);
557  m = 0.;
558  this->GetProblem().BoundaryMatrix(fdc, m,
559  _fs_alpha
560  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
561  1.
562  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
563  local_matrix.add(1.0, m);
564  }
565  else if (this->GetPart() == "New_for_2nd_cycle")
566  {
567  dealii::FullMatrix<double> m(local_matrix);
568  m = 0.;
569  this->GetProblem().BoundaryMatrix(fdc, m,
570  _fs_beta
571  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize(),
572  1.
573  * this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize());
574  local_matrix.add(1.0, m);
575  }
576  }
577  private:
578  // parameters for FS scheme
579  double _fs_theta;
580  double _fs_theta_prime;
581  double _fs_alpha;
582  double _fs_beta;
583 
585  FractionalStepThetaProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR,
586  dopedim, dealdim, FE, DH>, VECTOR, dealdim> * _initial_problem;
587  };
588 }
589 
590 #endif
std::string GetName()
Definition: fractional_step_theta_problem.h:84
~FractionalStepThetaProblem()
Definition: fractional_step_theta_problem.h:75
void FaceRhs(const FDC &fdc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: fractional_step_theta_problem.h:418
FractionalStepThetaProblem< OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH > & GetBaseProblem()
Definition: fractional_step_theta_problem.h:107
void BoundaryMatrix(const FDC &fdc, dealii::FullMatrix< double > &local_matrix)
Definition: fractional_step_theta_problem.h:551
void InterfaceMatrix(const FDC &fdc, dealii::FullMatrix< double > &local_matrix)
Definition: fractional_step_theta_problem.h:462
Definition: fractional_step_theta_problem.h:60
Definition: primal_ts_base.h:49
void FaceEquation(const FDC &fdc, dealii::Vector< double > &local_vector, double scale, double)
Definition: fractional_step_theta_problem.h:325
void PointRhs(const std::map< std::string, const dealii::Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values, VECTOR &rhs_vector, double scale)
Definition: fractional_step_theta_problem.h:243
void ElementEquation(const EDC &edc, dealii::Vector< double > &local_vector, double scale, double)
Definition: fractional_step_theta_problem.h:115
void ElementRhs(const EDC &edc, dealii::Vector< double > &local_vector, double scale)
Definition: fractional_step_theta_problem.h:212
Definition: initialproblem.h:42
InitialProblem< FractionalStepThetaProblem< OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH >, VECTOR, dealdim > & GetInitialProblem()
Definition: fractional_step_theta_problem.h:93
void InterfaceEquation(const FDC &fdc, dealii::Vector< double > &local_vector, double scale, double)
Definition: fractional_step_theta_problem.h:371
const SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > * GetSpaceTimeHandler() const
Definition: ts_base.h:442
OPTPROBLEM & GetProblem()
Definition: ts_base.h:466
void BoundaryEquation(const FDC &fdc, dealii::Vector< double > &local_vector, double scale, double)
Definition: fractional_step_theta_problem.h:494
void FaceMatrix(const FDC &fdc, dealii::FullMatrix< double > &local_matrix)
Definition: fractional_step_theta_problem.h:430
FractionalStepThetaProblem(OPTPROBLEM &OP)
Definition: fractional_step_theta_problem.h:64
std::string GetPart() const
Definition: ts_base.h:479
void BoundaryRhs(const FDC &fdc, dealii::Vector< double > &local_vector, double scale)
Definition: fractional_step_theta_problem.h:539
void ElementMatrix(const EDC &edc, dealii::FullMatrix< double > &local_matrix)
Definition: fractional_step_theta_problem.h:278