Commit 9ea224b2 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

updated ref on mpp

parent e00c774d
//#ifndef _TIME_INTGRATROR2_H_
//#define _TIME_INTGRATROR2_H_
//
//#include "assemble/DGTAssemble.h"
//
//
//class TimeIntegrator {
// Solver &S;
// Solver &S2;
// Solver &S3;
// Solver &IM;
// int plot_tStep;
// int rkorder;
// int Kmin;
// int Kmax;
// double Keps;
// Scalar gamma;
// Scalar gamma0;
// int SSum;
// int MSum;
// int KSum;
// int KPre;
// Time TPre;
// Time TComp;
//public:
// TimeIntegrator(Solver &__S, Solver &__S2, Solver &__S3, Solver &__IM) :
// S(__S), S2(__S2), S3(__S3), IM(__IM), plot_tStep(1), rkorder(1),
// Kmin(2), Kmax(10), Keps(0.001), gamma(1),
// SSum(0), MSum(0), KSum(0), KPre(0) {
// ReadConfig(Settings, "plot_tStep", plot_tStep);
// ReadConfig(Settings, "rkorder", rkorder);
// ReadConfig(Settings, "Kmin", Kmin);
// ReadConfig(Settings, "Kmax", Kmax);
// ReadConfig(Settings, "Keps", Keps);
// double g;
// ReadConfig(Settings, "gamma", g);
// gamma0 = gamma = g;
// }
//
// void Euler(double dt, const Vector &u, const Matrix &A, Vector &du) {
// Vector Au(u);
// Au = A * u;
// du = IM * Au;
// du *= dt;
// MSum += 1;
// }
//
// void Heun(double dt, const Vector &u, const Matrix &A, Vector &du) {
// Vector u1(u);
// Vector u2(u);
// Vector Au(u);
// Au = A * u;
// u1 = IM * Au;
// u2 = dt * u1;
// u2 += u;
// Au = A * u2;
// u2 = IM * Au;
// u2 += u1;
// du = (0.5 * dt) * u2;
// MSum += 2;
// }
//
// void Runge(double dt, const Vector &u, const Matrix &A, Vector &du) {
// const double w3 = 1.0 / 6.0;
// Vector u1(u);
// Vector u2(u);
// Vector u3(u);
// Vector Au(u);
// Au = A * u;
// u1 = IM * Au;
//
// u2 = 0.5 * dt * u1;
// u2 += u;
// Au = A * u2;
// u2 = IM * Au;
//
// u3 = 2.0 * dt * u2;
// u3 -= dt * u1;
// u3 += u;
// Au = A * u3;
// u3 = IM * Au;
//
// u3 += 4.0 * u2;
// u3 += u1;
//
// du = (w3 * dt) * u3;
// MSum += 3;
// }
//
// void Kutta(double dt, const Vector &u, const Matrix &A, Vector &du) {
// const double w4 = 1.0 / 6.0;
// Vector u1(u);
// Vector u2(u);
// Vector u3(u);
// Vector u4(u);
// Vector Au(u);
//
// Au = A * u;
// u1 = IM * Au;
//
// u2 = 0.5 * dt * u1;
// u2 += u;
// Au = A * u2;
// u2 = IM * Au;
//
// u3 = 0.5 * dt * u2;
// u3 += u;
// Au = A * u3;
// u3 = IM * Au;
//
// u4 = dt * u3;
// u4 += u;
// Au = A * u4;
// u4 = IM * Au;
//
// u4 += 2.0 * u3;
// u4 += 2.0 * u2;
// u4 += u1;
// du = (w4 * dt) * u4;
//
// MSum += 4;
// }
//
// void
// Kutta(double dt, const Vector &u, const Matrix &A, Vector &du, const Vector &RHS1,
// const Vector &RHS2,
// const Vector &RHS3) {
// const double w4 = 1.0 / 6.0;
// Vector u1(u);
// Vector u2(u);
// Vector u3(u);
// Vector u4(u);
// Vector Au(u);
//
// Au = A * u;
// Au += RHS1;
// u1 = IM * Au;
//
// u2 = 0.5 * dt * u1;
// u2 += u;
// Au = A * u2;
// Au += RHS2;
// u2 = IM * Au;
//
// u3 = 0.5 * dt * u2;
// u3 += u;
// Au = A * u3;
// Au += RHS2;
// u3 = IM * Au;
//
// u4 = dt * u3;
// u4 += u;
// Au = A * u4;
// Au += RHS3;
// u4 = IM * Au;
//
// u4 += 2.0 * u3;
// u4 += 2.0 * u2;
// u4 += u1;
// du = (w4 * dt) * u4;
//
// MSum += 4;
// }
//
// void Pade2(double dt, const Vector &u, const Matrix &A, Vector &du) {
// Vector Au(u);
// Au = A * u;
// du = S * Au;
// KPre += S.Iter();
// du *= dt;
// ++KSum;
// }
//
// void
// Pade2(double dt, const Vector &u, const Matrix &A, Vector &du, const Vector &RHS) {
// Vector Au(u);
// Au = A * u;
// Au += RHS;
// du = S * Au;
// KPre += S.Iter();
// du *= dt;
// ++KSum;
// }
//
// void Pade3(double dt, const Vector &u,
// const Matrix &M, const Matrix &A, Vector &du) {
// Vector Au(u);
// Vector Mu(u);
// Au = A * u;
// du = S * Au;
// KPre += S.Iter();
// Mu = M * du;
// Au = A * du;
// Au *= (dt * iUnit / sqrt(60));
// Mu -= Au;
// du = S2 * Mu;
// KPre += S2.Iter();
// Mu = M * du;
// Au = A * du;
// Au *= (dt * iUnit / sqrt(60));
// Mu += Au;
// du = S3 * Mu;
// KPre += S3.Iter();
// du *= dt;
// KSum += 3;
// }
//
// void ExponentialIntegrator(double dt, const Vector &u,
// const Matrix &M, const Matrix &A, Vector &du) {
// vector<Vector> v(Kmax + 1, u);
// vector<Vector> w(Kmax + 1, u);
// Vector Sw(u);
// SmallMatrix h(Kmax + 1, Kmax);
// SmallVector s(Kmax);
// SmallVector s_old(Kmax);
// s_old = Scalar(0.0);
// h = Scalar(0.0);
// v[0] = u;
// w[0] = M * u;
// double beta = sqrt(std::real(v[0] * w[0]));
// v[0] *= 1 / beta;
// w[0] *= 1 / beta;
// s_old[0] = beta;
// S.SetReduction(Keps);
// int k = 0;
// for (; k < Kmax; ++k) {
// Sw = S * w[k];
// KPre += S.Iter();
// Sw *= gamma / dt;
// v[k + 1] = Sw;
// for (int j = 0; j <= k; ++j) {
// h[j][k] = Sw * w[j];
// v[k + 1] -= h[j][k] * v[j];
// }
// w[k + 1] = M * v[k + 1];
// h[k + 1][k] = sqrt(std::real(v[k + 1] * w[k + 1]));
// v[k + 1] *= (1.0 / h[k + 1][k]);
// w[k + 1] *= (1.0 / h[k + 1][k]);
// SmallMatrix S_k(k + 1, k + 1);
// for (int j = 0; j <= k; ++j)
// for (int l = 0; l <= k; ++l)
// S_k[j][l] = h[j][l];
// S_k.invert();
// for (int j = 0; j <= k; ++j)
// S_k[j][j] -= dt / gamma;
// S_k *= Scalar(-1.0);
// S_k.Exp();
// for (int j = 0; j <= k; ++j)
// s[j] = beta * S_k[j][0];
// for (int j = 0; j <= k; ++j)
// s_old[j] -= s[j];
// double norm_u = s.norm();
// double delta = s_old.norm() / norm_u;
// for (int j = 0; j <= k; ++j)
// s_old[j] = s[j];
// if (k == 0) continue;
// double eta = 1 + norm_u;
// if (delta < 1) eta = min(eta, delta * norm_u / (1 - delta));
// mout << "eta(" << k << ") = " << eta << endl;
// S.SetReduction(0.1 * Keps / (Keps + eta));
// if (eta < Keps) {
// KSum += k;
// break;
// }
// }
// du = 0.0;
// for (int j = 0; j <= k; ++j)
// du += s[j] * v[j];
// du -= u;
// if (k == Kmax) du[0] = 2e10;
// }
//
// void Arnoldi(double dt, const Vector &u,
// const Matrix &M, const Matrix &A, Vector &du) {
// vector<Vector> v(Kmax + 1, u);
// Vector Mv(u);
// Vector Av(u);
// SmallMatrix h(Kmax + 1, Kmax);
// SmallVector s(Kmax);
// SmallVector s_old(Kmax);
// s_old = Scalar(0.0);
// h = Scalar(0.0);
// v[0] = u;
// Mv = M * v[0];
// double beta = sqrt(std::real(v[0] * Mv));
// if (beta < Eps) {
// du = 0;
// return;
// }
// v[0] *= 1 / beta;
// s_old[0] = beta;
// int k = 0;
// for (; k < Kmax; ++k) {
// Av = A * v[k];
// v[k + 1] = IM * Av;
// for (int j = 0; j <= k; ++j) {
// h[j][k] = Av * v[j];
// v[k + 1] -= h[j][k] * v[j];
// }
// Mv = M * v[k + 1];
// h[k + 1][k] = sqrt(std::real(v[k + 1] * Mv));
// v[k + 1] *= (1.0 / h[k + 1][k]);
// SmallMatrix H_k(k + 1, k + 1);
// for (int j = 0; j <= k; ++j)
// for (int l = 0; l <= k; ++l)
// H_k[j][l] = h[j][l] * dt;
// H_k.Exp();
// for (int j = 0; j <= k; ++j)
// s[j] = beta * H_k[j][0];
// for (int j = 0; j <= k; ++j)
// s_old[j] -= s[j];
// double norm_u = s.norm();
// double delta = s_old.norm() / norm_u;
// for (int j = 0; j <= k; ++j)
// s_old[j] = s[j];
// if (k == 0) continue;
// double eta = 1 + norm_u;
// if (delta < 1) eta = min(eta, delta * norm_u / (1 - delta));
// mout << "eta(" << k << ") = " << eta << endl;
// if (eta < Keps) break;
// }
// du = 0.0;
// for (int j = 0; j <= k; ++j)
// du += s[j] * v[j];
// MSum += k + 1;
// du -= u;
// if (k == Kmax) du[0] = 2e10;
// }
//
// void Arnoldi2(double dt, const Vector &u,
// const Matrix &M, const Matrix &A, Vector &du) {
// vector<Vector> v(Kmax + 1, u);
// vector<Vector> w(Kmax + 1, u);
// Vector Sw(u);
// SmallMatrix h(Kmax + 1, Kmax);
// SmallVector s(Kmax);
// SmallVector s_old(Kmax);
// s_old = Scalar(0.0);
// h = Scalar(0.0);
// v[0] = u;
// w[0] = M * u;
// double beta = sqrt(std::real(v[0] * w[0]));
// v[0] *= 1 / beta;
// w[0] *= 1 / beta;
// s_old[0] = beta;
// S.SetReduction(Keps);
// int k = 0;
// for (; k < Kmax; ++k) {
// Sw = S * w[k];
// KPre += S.Iter();
// Sw *= gamma / dt;
// v[k + 1] = Sw;
// for (int j = 0; j <= k; ++j) {
// h[j][k] = Sw * w[j];
// v[k + 1] -= h[j][k] * v[j];
// }
// w[k + 1] = M * v[k + 1];
// h[k + 1][k] = sqrt(std::real(v[k + 1] * w[k + 1]));
// v[k + 1] *= (1.0 / h[k + 1][k]);
// w[k + 1] *= (1.0 / h[k + 1][k]);
// SmallMatrix S_k(k + 1, k + 1);
// for (int j = 0; j <= k; ++j)
// for (int l = 0; l <= k; ++l)
// S_k[j][l] = h[j][l];
// S_k.invert();
// for (int j = 0; j <= k; ++j)
// S_k[j][j] -= dt / gamma;
// S_k *= Scalar(-1.0);
// S_k.Exp();
// for (int j = 0; j <= k; ++j)
// s[j] = beta * S_k[j][0];
// for (int j = 0; j <= k; ++j)
// s_old[j] -= s[j];
// double norm_u = s.norm();
// double delta = s_old.norm() / norm_u;
// for (int j = 0; j <= k; ++j)
// s_old[j] = s[j];
// if (k == 0) continue;
// double eta = 1 + norm_u;
// if (delta < 1) eta = min(eta, delta * norm_u / (1 - delta));
// mout << "eta(" << k << ") = " << eta << endl;
// S.SetReduction(Keps / (Keps + eta));
// if (eta < Keps) {
// KSum += k;
// break;
// }
// }
// du = 0.0;
// for (int j = 0; j <= k; ++j)
// du += s[j] * v[j];
// du -= u;
// if (k == Kmax) du[0] = 2e10;
// }
//
// void Method(DGTAssemble &assemble, TimeSeries &timeSeries, Vector &u) {
// Date start;
//
// Matrix massMatrix(u);
// assemble.MassMatrix(massMatrix);
//
// Matrix fluxMatrix(u);
// assemble.FluxMatrix(fluxMatrix);
//
// double lastTStep = timeSeries.LastTStep();
// double t = timeSeries.FirstTStep();
// int step = 0;
//
// assemble.Initialize(t, u);
//
// double oldTStep = 0.0;
// double tau = 0.0;
// while (t < lastTStep) {
// oldTStep = t;
// t = timeSeries.NextTStep(t);
// tau = t - oldTStep;
//
// if (true) {
// Euler(dt, u, A, du);
// Heun(dt, u, A, du);
// Runge(dt, u, A, du);
//
// assemble.RHS(RHS, t);
// assemble.RHS(RHS2, t + 0.5 * dt);
// assemble.RHS(RHS3, t + dt);
// Kutta(dt, u, A, du, RHS, RHS2, RHS3);
// }
//
// }
//
// Matrix B(u);
// Matrix B2(u);
// Matrix B3(u);
// IM(massMatrix);
//
// Vector RHS(u);
// Vector RHS2(u);
// Vector RHS3(u);
// RHS = 0.;
// RHS2 = 0.;
// RHS3 = 0.;
//
// double dt_old = 0;
// Vector du(u);
// SSum = MSum = KSum = KPre = 0;
// while (t < T) {
// switch(rkorder) {
// case 0:
// if (dt_old == 0) {
// B = M;
// gamma = gamma0 * dt;
// B += (-gamma) * A;
// mout << "Exponential Integrator gamma = " << gamma << endl;
// S(B, 1);
// }
// ExponentialIntegrator(dt, u, M, A, du);
// break;
// case 1000:
// if (dt_old == 0)
// mout << "Arnoldi " << endl;
// Arnoldi(dt, u, M, A, du);
// // Arnoldi(dt,u,M,A,du,RHS); // Hier muss der Algorithmus noch für die rechte Seite angepasst werden
// break;
// case 2000:
// if (dt_old == 0) {
// B = M;
// gamma = gamma0 * dt;
// B += (-gamma) * A;
// S(B, 1);
// mout << "Arnoldi2 gamma = " << gamma << endl;
// }
// Arnoldi2(dt, u, M, A, du);
// break;
// case -2:
// if (abs(dt - dt_old) > 1e-13) {
// B = M;
// gamma = 0.5 * dt;
// assemble.SetGamma(gamma);
// B += (-gamma) * A;
// mout << "Implicit Midpoint Rule: gamma = " << gamma << endl;
// S(B, 1);
// }
//
// assemble.RHS(RHS, t + 0.5 * dt);
//
// Pade2(dt, u, A, du, RHS);
// break;
// case -3:
// if (abs(dt - dt_old) > 1e-13) {
// const Scalar gamma =
// 4.0 - 2 * pow(2.0 / (1.0 + sqrt(5.0)), 1.0 / 3.0)
// + pow(2.0, 2.0 / 3.0) * pow(1.0 + sqrt(5), 1.0 / 3.0);
// const Scalar alpha = 4.0 + (1.0 - iUnit * sqrt(3.0)) *
// pow(2.0 / (1.0 + sqrt(5)), 1.0 / 3.0)
// - (1.0 + iUnit * sqrt(3.0)) *
// pow(0.5 * (1.0 + sqrt(5)), 1.0 / 3.0);
// B = M;
// B += Scalar(-dt / gamma) * A;
// B2 = M;
// B2 += Scalar(-dt / alpha) * A;
// B3 = M;
// B3 += Scalar(-dt / std::conj(alpha)) * A;
// mout << "Gauss (s=3)" << endl;
// S(B, 1);
// S2(B2, 1);
// S3(B3, 1);
// }
// Pade3(dt, u, M, A, du);
// break;
// }
// u += du;
// ++SSum;
// if (u.norm() > 1e10) {
// u = Scalar(0.0);
// mout << "divergence at t = " << t << endl;
// return;
// }
// assemble.FinishTimeStep(t, u, ++step);
// dt_old = dt;
// }
// }
//
// int S_Sum() const { return SSum; }
//
// int M_Sum() const { return MSum; }
//
// int K_Sum() const { return KSum; }
//
// int K_Pre() const { return KPre; }
//
// Time T_Comp() const { return TComp; }
//
// Time T_Pre() const { return TPre; }
//};
//
//#endif
Subproject commit 737aca3f1b7ee3ac459d844f8114c1278d61739f
Subproject commit cab7faecdff8efcef7d54fd8bf573a4bb2c538dc
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment