DGTransportAssemble.cpp 12.2 KB
Newer Older
1
#include "DGTransportAssemble.hpp"
niklas.baumgarten's avatar
piles    
niklas.baumgarten committed
2
#include "Plotting.hpp"
3
4


niklas.baumgarten's avatar
niklas.baumgarten committed
5
void DGLinearTransportAssemble::MassMatrix(Matrix &massMatrix) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
6
7
  massMatrix = 0;
  for (cell c = massMatrix.cells(); c != massMatrix.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
8
    DGElement elem(massMatrix, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
9
10
11
12
13
14
15
16
    DGRowEntries M_c(massMatrix, c, c);
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      for (int i = 0; i < elem.NodalPoints(); ++i) {
        Scalar Phi_i = elem.Value(q, i);
        for (int j = 0; j < elem.NodalPoints(); ++j) {
          Scalar Phi_j = elem.Value(q, j);
          M_c(i, j) += w * Phi_i * Phi_j;
17
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
18
      }
19
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
20
  }
21
22
}

niklas.baumgarten's avatar
niklas.baumgarten committed
23
void DGLinearTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
24
25
  systemMatrix = 0;
  for (cell c = systemMatrix.cells(); c != systemMatrix.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
26
    DGElement elem(systemMatrix, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
27
28
29
30
31
32
33
34
35
36
37
38
39
    DGRowEntries A_c(systemMatrix, c, c);
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      VectorField B = problem->CellFlux(c, elem.QPoint(q));
      for (int i = 0; i < elem.NodalPoints(); ++i) {
        Scalar Phi_i = elem.Value(q, i);
        for (int j = 0; j < elem.NodalPoints(); ++j) {
          VectorField gradPhi_j = elem.Derivative(q, j);
          A_c(i, j) -= w * (gradPhi_j * B * Phi_i);
        }
      }
    }
    for (int f = 0; f < c.Faces(); ++f) {
niklas.baumgarten's avatar
niklas.baumgarten committed
40
      DGFaceElement faceElem(systemMatrix, c, f);
niklas.baumgarten's avatar
niklas.baumgarten committed
41
42
43
44
45
46
47
48
49
50
51
52
      if (systemMatrix.GetMesh().onBndDG(c, f)) {
        for (int q = 0; q < faceElem.nQ(); ++q) {
          const Point &Qf_c = faceElem.QPoint(q);
          VectorField Nq = faceElem.QNormal(q);
          Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
          if (BN > 0) continue;
          double w = faceElem.QWeight(q);
          for (int i = 0; i < faceElem.NodalPoints(); ++i) {
            Scalar phi_i = faceElem.Value(q, i);
            for (int j = 0; j < faceElem.NodalPoints(); ++j) {
              Scalar phi_j = faceElem.Value(q, j);
              A_c(i, j) += w * BN * phi_j * phi_i;
53
            }
niklas.baumgarten's avatar
niklas.baumgarten committed
54
          }
55
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
56
57
58
      } else {
        cell cf = systemMatrix.GetMesh().find_neighbour_cell(c, f);
        int f1 = systemMatrix.GetMesh().find_neighbour_face_id(c.Face(f), cf);
niklas.baumgarten's avatar
niklas.baumgarten committed
59
        DGFaceElement felem_1(systemMatrix, cf, f1);
niklas.baumgarten's avatar
niklas.baumgarten committed
60
61
62
63
64
65
66
67
68
69
70
71
72
        DGRowEntries A_cf(systemMatrix, c, cf);
        for (int q = 0; q < faceElem.nQ(); ++q) {
          const Point &Qf_c = faceElem.QPoint(q);
          VectorField Nq = faceElem.QNormal(q);
          Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
          if (BN > 0) continue;
          int q1 = felem_1.findQPointID(faceElem, Qf_c);
          double w = faceElem.QWeight(q);
          for (int i = 0; i < faceElem.NodalPoints(); ++i) {
            Scalar phi_i = faceElem.Value(q, i);
            for (int j = 0; j < faceElem.NodalPoints(); ++j) {
              Scalar phi_j = faceElem.Value(q, j);
              A_c(i, j) += w * BN * phi_j * phi_i;
73
            }
niklas.baumgarten's avatar
niklas.baumgarten committed
74
75
76
77
78
            for (int j = 0; j < felem_1.NodalPoints(); ++j) {
              Scalar phi_j = felem_1.Value(q1, j);
              A_cf(i, j) -= w * BN * phi_j * phi_i;
            }
          }
79
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
        if (cf() < c()) continue;
        DGRowEntries A_fc(systemMatrix, cf, c);
        DGRowEntries A_ff(systemMatrix, cf, cf);
        for (int q = 0; q < faceElem.nQ(); ++q) {
          double w = faceElem.QWeight(q);
          const Point &z = faceElem.QPoint(q);
          const Point &N = faceElem.QNormal(q);
          const Point &Qf_c = faceElem.QPoint(q);
          int q1 = felem_1.findQPointID(faceElem, Qf_c);
          double s = 1;
          for (int i = 0; i < faceElem.NodalPoints(); ++i) {
            Scalar phi_i = faceElem.Value(q, i);
            Scalar NDphi_i = diffusion *
              (faceElem.Derivative(q, i) * N);
            for (int j = 0; j < faceElem.NodalPoints(); ++j) {
              Scalar phi_j = faceElem.Value(q, j);
              Scalar NDphi_j =
                diffusion * (faceElem.Derivative(q, j) * N);
              A_c(i, j) -= w * (-0.5 * NDphi_i * phi_j
                - 0.5 * phi_i * NDphi_j
                + s * phi_i * phi_j);
101
            }
niklas.baumgarten's avatar
niklas.baumgarten committed
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
            for (int j = 0; j < felem_1.NodalPoints(); ++j) {
              Scalar phi_j = felem_1.Value(q1, j);
              Scalar NDphi_j =
                diffusion * (felem_1.Derivative(q1, j) * N);
              A_cf(i, j) -= w * (0.5 * NDphi_i * phi_j
                - phi_i * 0.5 * NDphi_j
                - s * phi_i * phi_j);
              A_fc(j, i) -= w * (0.5 * NDphi_i * phi_j
                - phi_i * 0.5 * NDphi_j
                - s * phi_i * phi_j);
            }
          }
          for (int i = 0; i < felem_1.NodalPoints(); ++i) {
            Scalar phi_i = felem_1.Value(q1, i);
            Scalar NDphi_i = diffusion
              * (felem_1.Derivative(q1, i) * N);
            for (int j = 0; j < felem_1.NodalPoints(); ++j) {
              Scalar phi_j = felem_1.Value(q1, j);
              Scalar NDphi_j = diffusion
                * (felem_1.Derivative(q1, j) * N);
              A_ff(i, j) -= w * (0.5 * NDphi_i * phi_j
                + phi_i * 0.5 * NDphi_j
                + s * phi_i * phi_j);
            }
          }
127
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
128
      }
129
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
130
  }
131
132
}

niklas.baumgarten's avatar
niklas.baumgarten committed
133
134
135
136
137
138
139
void DGLinearTransportAssemble::RHS(double t, Vector &rhs) const {
  rhs = 0;
  if (!problem->RHS()) return;
  for (cell c = rhs.cells(); c != rhs.cells_end(); ++c) {
    row r = rhs.find_row(c());
    for (int f = 0; f < c.Faces(); ++f) {
      if (!rhs.GetMesh().onBndDG(c, f)) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
140
      DGFaceElement felem(rhs, c, f);
niklas.baumgarten's avatar
niklas.baumgarten committed
141
142
143
144
145
146
147
148
149
150
      for (int q = 0; q < felem.nQ(); ++q) {
        const Point &z = felem.QPoint(q);
        double w = felem.QWeight(q);
        VectorField N = felem.QNormal(q);
        Scalar BN = problem->FaceNormalFlux(c, f, N, z);
        if (BN > 0) continue;
        Scalar U = problem->Solution(t, z);
        for (int i = 0; i < felem.NodalPoints(); ++i) {
          Scalar Phi_i = felem.Value(q, i);
          rhs(r, i) -= w * U * BN * Phi_i;
151
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
152
      }
153
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
154
  }
155
156
}

niklas.baumgarten's avatar
niklas.baumgarten committed
157
158
159
double DGLinearTransportAssemble::Energy(const Vector &u) const {
  double energy = 0;
  for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
160
    DGElement elem(u, c);
161
    for (int q = 0; q < elem.nQ(); ++q) {
niklas.baumgarten's avatar
niklas.baumgarten committed
162
163
164
      double w = elem.QWeight(q);
      Scalar U = elem.Value(q, u);
      energy += w * (U * U);
165
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
166
  }
167
  return 0.5 * PPM->SumOnCommSplit(energy, u.GetMesh().CommSplit());
niklas.baumgarten's avatar
niklas.baumgarten committed
168
169
170
171
172
}

void DGNonLinearTransportAssemble::Energy(const cell &c,
                                          const Vector &u,
                                          double &energy) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
173
  DGElement elem(u, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
174
175
176
177
178
  for (int q = 0; q < elem.nQ(); ++q) {
    double w = elem.QWeight(q);
    Scalar U = elem.Value(q, u);
    energy += w * (U * U);
  }
179
180
}

niklas.baumgarten's avatar
niklas.baumgarten committed
181
double DGLinearTransportAssemble::Error(double t, const Vector &u) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
182
183
  double err = 0.0;
  for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
184
    DGElement elem(u, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
185
186
187
188
189
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      Scalar U = elem.Value(q, u);
      Scalar Sol = problem->Solution(t, elem.QPoint(q));
      err += w * (U - Sol) * (U - Sol);
190
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
191
  }
192
  return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
193
194
}

niklas.baumgarten's avatar
niklas.baumgarten committed
195
196
197
198
199
200
std::pair<double, double> DGLinearTransportAssemble::InflowOutflow(const Vector &u) const {
  double inflow = 0.0;
  double outflow = 0.0;
  for (cell c = u.cells(); c != u.cells_end(); ++c) {
    BFParts bnd(u.GetMesh(), c);
    if (!bnd.onBnd()) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
201
    DGElement elem(u, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
202
203
204
205
206
207
    double U = 0;
    double area = 0;
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      U += w * elem.Value(q, u);
      area += w;
208
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
209
210
211
    U *= (1 / area);
    for (int f = 0; f < c.Faces(); ++f) {
      if (bnd[f] == -1) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
212
      DGFaceElement faceElem(u, c, f);
niklas.baumgarten's avatar
niklas.baumgarten committed
213
214
215
216
217
218
219
220
221
      for (int q = 0; q < faceElem.nQ(); ++q) {
        double w = faceElem.QWeight(q);
        VectorField Nq = faceElem.QNormal(q);
        Scalar BN = problem->FaceNormalFlux(c, f, Nq, elem.QPoint(q));
        if (BN > 0) outflow += w * BN * U;
        else inflow += w * BN * U;
      }
    }
  }
222
223
  inflow = PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit());
  outflow = PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit());
niklas.baumgarten's avatar
niklas.baumgarten committed
224
225
226
227
228
  if (inflow < 1e-7)
    inflow = 0.0;
  if (outflow < 1e-7)
    outflow = 0.0;
  return {inflow, outflow};
229
230
}

niklas.baumgarten's avatar
niklas.baumgarten committed
231
void DGLinearTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
232
233
234
  u_ex = 0;
  for (cell c = u_ex.cells(); c != u_ex.cells_end(); ++c) {
    row r = u_ex.find_row(c());
niklas.baumgarten's avatar
niklas.baumgarten committed
235
    DGElement Elem(u_ex, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
236
237
238
239
    for (int j = 0; j < Elem.NodalPoints(); ++j)
      u_ex(r, j) = problem->Solution(t, Elem.NodalPoint(j));
  }
  u_ex.Accumulate();
240
241
}

niklas.baumgarten's avatar
niklas.baumgarten committed
242
double DGLinearTransportAssemble::Mass(const Vector &u) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
243
244
  Scalar e = 0;
  for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
245
    DGElement elem(u, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
246
247
248
249
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      Scalar U = elem.Value(q, u);
      e += w * U;
250
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
251
  }
252
  return PPM->SumOnCommSplit(e, u.GetMesh().CommSplit());
253
254
}

niklas.baumgarten's avatar
niklas.baumgarten committed
255
void DGLinearTransportAssemble::SetInitialValue(Vector &u) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
256
257
258
  u = 0;
  for (cell c = u.cells(); c != u.cells_end(); ++c) {
    row r = u.find_row(c());
niklas.baumgarten's avatar
niklas.baumgarten committed
259
    DGElement Elem(u, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
260
261
262
263
264
    double lambda = 1e-9;
    for (int j = 0; j < Elem.NodalPoints(); ++j)
      u(r, j) = problem->Solution(0, lambda * c() + (1 - lambda) * Elem.NodalPoint(j));
  }
  u.Accumulate();
265
266
}

niklas.baumgarten's avatar
niklas.baumgarten committed
267
void DGLinearTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
268
  for (cell c = A.cells(); c != A.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
269
    DGElement elem(A, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
270
271
272
273
274
275
276
277
    DGRowEntries A_c(A, c, c);
    for (int i = 0; i < elem.NodalPoints(); ++i) {
      for (int j = 0; j < elem.NodalPoints(); ++j)
        mout << A_c(i, j) << " ";
      mout << endl;
    }
    mout << endl;
    if (diagonal) continue;
278

niklas.baumgarten's avatar
niklas.baumgarten committed
279
280
281
282
283
284
285
286
287
288
289
290
291
    for (int f = 0; f < c.Faces(); ++f) {
      cell cf = c;
      if (A.GetMesh().onBndDG(c, f))
        continue;
      else
        cf = A.GetMesh().find_neighbour_cell(c, f);
      DGRowEntries A_cf(A, c, cf);
      for (int i = 0; i < elem.NodalPoints(); ++i) {
        for (int j = 0; j < elem.NodalPoints(); ++j)
          mout << A_cf(i, j) << " ";
        mout << endl;
      }
      mout << endl;
292
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
293
294
    mout << "------------------------" << endl;
  }
295
296
}

niklas.baumgarten's avatar
niklas.baumgarten committed
297
298
299
void DGLinearTransportAssemble::VtkPlotting_cell(double t,
                                                 const Vector &u,
                                                 char *filename) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
300
301
  Vector u_tmp(u);
  for (cell c = u_tmp.cells(); c != u_tmp.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
302
    DGElement elem(u_tmp, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
303
304
305
306
307
308
    Scalar U = 0.0;
    double a = 0;
    for (int q = 0; q < elem.nQ(); ++q) {
      double w = elem.QWeight(q);
      U += w * elem.Value(q, u);
      a += w;
309
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
310
311
312
313
314
    U *= (1 / a);
    row r = u.find_row(c());
    u_tmp(r)[0] = U;
  }
  mpp::plot("U") << u_tmp << mpp::save_plot(NumberName("U", filename, step));
315
316
}

niklas.baumgarten's avatar
niklas.baumgarten committed
317
void DGLinearTransportAssemble::VtkPlotting(double t, const Vector &u) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
318
319
  char filename[128];
  VtkPlotting_cell(t, u, filename);
320
321
}

niklas.baumgarten's avatar
niklas.baumgarten committed
322
void DGLinearTransportAssemble::PrintInfo() const {
niklas.baumgarten's avatar
niklas.baumgarten committed
323
324
325
326
  mout.PrintInfo("Assemble", verbose,
                 PrintInfoEntry("Name", Name()),
                 PrintInfoEntry("Problem", problem->Name()),
                 PrintInfoEntry("Discretization", disc->DiscName()));
327
328
}

niklas.baumgarten's avatar
niklas.baumgarten committed
329
void DGLinearTransportAssemble::FinishTimeStep(double t, Vector &u) {
niklas.baumgarten's avatar
niklas.baumgarten committed
330
331
332
333
  step++;
  t = t;
  PrintInfo(u);
  VtkPlotting(t, u);
334
335
}

niklas.baumgarten's avatar
niklas.baumgarten committed
336
void DGLinearTransportAssemble::Initialize(Vector &u) {
niklas.baumgarten's avatar
niklas.baumgarten committed
337
338
339
340
  u = 0;
  SetInitialValue(u);
  PrintInfo(u);
  VtkPlotting(timeSeries.FirstTStep(), u);
341
342
}

niklas.baumgarten's avatar
niklas.baumgarten committed
343
double DGNonLinearTransportAssemble::Residual(const Vector &u, Vector &b) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
344
  b = 0;
niklas.baumgarten's avatar
niklas.baumgarten committed
345

niklas.baumgarten's avatar
niklas.baumgarten committed
346
347
  auto *massMatrix = new Matrix(u);
  MassMatrix(*massMatrix);
niklas.baumgarten's avatar
niklas.baumgarten committed
348

niklas.baumgarten's avatar
niklas.baumgarten committed
349
350
  auto *fluxMatrix = new Matrix(u);
  SystemMatrix(*fluxMatrix);
niklas.baumgarten's avatar
niklas.baumgarten committed
351

niklas.baumgarten's avatar
niklas.baumgarten committed
352
353
354
355
356
  Vector fluxMatrixU(u);
  fluxMatrixU = *fluxMatrix * u;
  fluxMatrixU *= -dt_;
  b = (*massMatrix * u + fluxMatrixU);
  b -= *massMatrix * U_old();
niklas.baumgarten's avatar
niklas.baumgarten committed
357

niklas.baumgarten's avatar
niklas.baumgarten committed
358
359
360
  Vector rhs(b);
  RHS(t_, rhs);
  b -= dt_ * rhs;
niklas.baumgarten's avatar
niklas.baumgarten committed
361

niklas.baumgarten's avatar
niklas.baumgarten committed
362
363
364
365
366
  b.ClearDirichletValues();
  b.Collect();
  delete fluxMatrix;
  delete massMatrix;
  return b.norm();
niklas.baumgarten's avatar
niklas.baumgarten committed
367
368
369
}

void DGNonLinearTransportAssemble::Jacobi(const Vector &u, Matrix &A) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
370
371
  Matrix *massMatrix = new Matrix(u);
  MassMatrix(*massMatrix);
niklas.baumgarten's avatar
niklas.baumgarten committed
372

niklas.baumgarten's avatar
niklas.baumgarten committed
373
374
375
376
377
378
  Matrix *fluxMatrix = new Matrix(u);
  SystemMatrix(*fluxMatrix);
  A = *massMatrix;
  A += -dt_ * (*fluxMatrix);
  delete fluxMatrix;
  delete massMatrix;
niklas.baumgarten's avatar
niklas.baumgarten committed
379
380
381
}