TimeIntegrator.cpp 4.06 KB
Newer Older
1
#include "Arnoldi.hpp"
2
#include "TimeIntegrator.hpp"
3
#include "GenericRungeKutta.hpp"
4
#include "ExponentialIntegrator.hpp"
5
#include "ImplicitEuler.hpp"
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56


bool TimeIntegrator::Method(ILinearTimeAssemble *assemble, Vector &u) {
  mout.StartBlock("TI");
  vout(1) << "Starting " << Name() << " (Order=" << GetOrder() << ")" << endl;
  double t = assemble->GetTimeSeries().FirstTStep();

  Matrix massMatrix(u);
  Matrix systemMatrix(u);

  assemble->Initialize(u);
  assemble->MassMatrix(massMatrix);
  assemble->SystemMatrix(systemMatrix);

  (*solvers.at(0))(massMatrix);

  Vector du(u);
  Matrix B(u);

  while (t < assemble->GetTimeSeries().LastTStep()) {
    t = assemble->GetTimeSeries().NextTStep(t);
    double dt = assemble->GetTimeSeries().StepSize();
    assemble->InitTimeStep();
    AssembleRHS(assemble, u);
    StageFunctions(dt, u, massMatrix, systemMatrix, du);
    u += du;
    if (NotDiverged(u)) assemble->FinishTimeStep(t, u);
    else {
      mout.EndBlock();
      return false;
    }
  }
  mout.EndBlock();
  return true;
}

bool TimeIntegrator::NotDiverged(Vector &u) {
  if (u.norm() <= 1e10) {
    return true;
  } else {
    u = 0.0;
    Warning("Time integrator diverged")
    return false;
  }
}

void TimeIntegrator::PrintInfo() const {
  mout.PrintInfo("Time Integrator", verbose,
                 PrintInfoEntry("Name", Name(), 1),
                 PrintInfoEntry("Order", GetOrder(), 1));
}
57

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
/*
 * NonLineaTimeIntegrator
 */

bool NonLinearTimeIntegrator::NotDiverged(Vector &u) {
  if (u.norm() <= 1e10) {
    return true;
  } else {
    u = 0.0;
    Warning("Time integrator diverged")
    return false;
  }
}

void NonLinearTimeIntegrator::PrintInfo() const {
  mout.PrintInfo("Time Integrator", verbose,
                 PrintInfoEntry("Name", Name(), 1),
                 PrintInfoEntry("Order", GetOrder(), 1));
}

/*
 * TimeIntegratorCreator
 */

82
TimeIntegrator *TimeIntegratorCreator::Create() {
83
  try {
84
    switch (type) {
85
      case EXPLICIT_EULER:
86
        return new GenericExplicitEuler(linearSolvers.at(0));
87
      case HEUN:
88
        return new GenericHeun(linearSolvers.at(0));
89
      case RUNGE:
90
        return new GenericRunge(linearSolvers.at(0));
91
      case RUNGE_KUTTA:
92
        return new GenericRungeKutta4(linearSolvers.at(0));
93
      case IMPLICIT_EULER: // Todo figure out how many are really needed
94
95
        return new GenericImplicitEuler(linearSolvers.at(0), linearSolvers.at(1),
                                        linearSolvers.at(2), linearSolvers.at(3));
96
      case IMPLICIT_MIDPOINT:
97
98
        return new GenericImplicitMidpointRule(linearSolvers.at(0), linearSolvers.at(1),
                                               linearSolvers.at(2), linearSolvers.at(3));
99
      case CRANK_NICOLSON:
100
101
        return new CrankNicolson(linearSolvers.at(0), linearSolvers.at(1),
                                 linearSolvers.at(2), linearSolvers.at(3));
102
      case DI_RK_ORDER2:
103
104
        return new DiagonalImplicitRungeKutta2(linearSolvers.at(0), linearSolvers.at(1),
                                               linearSolvers.at(2), linearSolvers.at(3));
105
      case DI_RK_ORDER3:
106
107
        return new DiagonalImplicitRungeKutta3(linearSolvers.at(0), linearSolvers.at(1),
                                               linearSolvers.at(2), linearSolvers.at(3));
108
      case ARNOLDI:
109
        return new Arnoldi(linearSolvers.at(0));
110
      case ARNOLDI2:
111
        return new Arnoldi2(linearSolvers.at(0), linearSolvers.at(1));
112
      case EXPONENTIAL:
113
        return new ExponentialIntegrator(linearSolvers.at(0), linearSolvers.at(1));
114
      default: Exit(to_string(type) + " TimeIntegrator not implemented")
115
    }
116
117
  }
  catch (std::out_of_range) {
118
    Exit("Number of solvers does not match selected time integrator")
119
  }
120
}
121
122
123
124
125
126
127
128
129
130
131
132
133
134

NonLinearTimeIntegrator *TimeIntegratorCreator::CreateNonLinearTimeIntegrator() {
  try {
    switch (type) {
      case IMPLICIT_EULER:
        return new ImplicitEuler(nonLinearSolvers.at(0));
      default: Exit(to_string(type) + " TimeIntegrator not implemented")
    }
  }
  catch (std::out_of_range) {
    Exit("Number of solvers does not match selected time integrator")
  }
}