DGTransportAssemble.cpp 14.4 KB
Newer Older
1
2
3
4
5
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
#include "DGTransportAssemble.hpp"


double DGTransportAssemble::Residual(const Vector &u, Vector &b) const {
    b = 0;

    auto *massMatrix = new Matrix(u);
    MassMatrix(*massMatrix);

    auto *fluxMatrix = new Matrix(u);
    SystemMatrix(*fluxMatrix);

    Vector fluxMatrixU(u);
    fluxMatrixU = *fluxMatrix * u;
    fluxMatrixU *= -dt_;
    b = (*massMatrix * u + fluxMatrixU);
    b -= *massMatrix * U_old();

    Vector rhs(b);
    RHS(t_, rhs);
    b -= dt_ * rhs;

    b.ClearDirichletValues();
    Collect(b);
    delete fluxMatrix;
    delete massMatrix;
    return b.norm();
}

void DGTransportAssemble::Jacobi(const Vector &u, Matrix &A) const {
    Matrix *massMatrix = new Matrix(u);
    MassMatrix(*massMatrix);

    Matrix *fluxMatrix = new Matrix(u);
    SystemMatrix(*fluxMatrix);
    A = *massMatrix;
    A += -dt_ * (*fluxMatrix);
    delete fluxMatrix;
    delete massMatrix;
}

void DGTransportAssemble::MassMatrix(Matrix &massMatrix) const {
    massMatrix = 0;
    for (cell c = massMatrix.cells(); c != massMatrix.cells_end(); ++c) {
        DGElement elem(*disc, massMatrix, c);
        DGRowEntries M_c(massMatrix, c, c);
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
49
            for (int i = 0; i < elem.NodalPoints(); ++i) {
50
                Scalar Phi_i = elem.Value(q, i);
niklas.baumgarten's avatar
niklas.baumgarten committed
51
                for (int j = 0; j < elem.NodalPoints(); ++j) {
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
                    Scalar Phi_j = elem.Value(q, j);
                    M_c(i, j) += w * Phi_i * Phi_j;
                }
            }
        }
    }
}

void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
    systemMatrix = 0;
    for (cell c = systemMatrix.cells(); c != systemMatrix.cells_end(); ++c) {
        DGElement elem(*disc, systemMatrix, c);
        DGRowEntries A_c(systemMatrix, c, c);
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
67
            VectorField B = problem->CellFlux(c, elem.QPoint(q));
niklas.baumgarten's avatar
niklas.baumgarten committed
68
            for (int i = 0; i < elem.NodalPoints(); ++i) {
69
                Scalar Phi_i = elem.Value(q, i);
niklas.baumgarten's avatar
niklas.baumgarten committed
70
                for (int j = 0; j < elem.NodalPoints(); ++j) {
71
72
73
74
75
76
77
                    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) {
            DGFaceElement faceElem(*disc, systemMatrix, c, f);
78
            if (systemMatrix.GetMesh().onBndDG(c, f)) {
79
80
81
                for (int q = 0; q < faceElem.nQ(); ++q) {
                    const Point &Qf_c = faceElem.QPoint(q);
                    VectorField Nq = faceElem.QNormal(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
82
                    Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
83
84
                    if (BN > 0) continue;
                    double w = faceElem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
85
                    for (int i = 0; i < faceElem.NodalPoints(); ++i) {
86
                        Scalar phi_i = faceElem.Value(q, i);
niklas.baumgarten's avatar
niklas.baumgarten committed
87
                        for (int j = 0; j < faceElem.NodalPoints(); ++j) {
88
89
90
91
92
93
                            Scalar phi_j = faceElem.Value(q, j);
                            A_c(i, j) += w * BN * phi_j * phi_i;
                        }
                    }
                }
            } else {
94
95
                cell cf = systemMatrix.GetMesh().find_neighbour_cell(c, f);
                int f1 = systemMatrix.GetMesh().find_neighbour_face_id(c.Face(f), cf);
96
97
98
99
100
                DGFaceElement felem_1(*disc, systemMatrix, cf, f1);
                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);
niklas.baumgarten's avatar
niklas.baumgarten committed
101
                    Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
102
                    if (BN > 0) continue;
103
                    int q1 = felem_1.findQPointID(faceElem, Qf_c);
104
                    double w = faceElem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
105
                    for (int i = 0; i < faceElem.NodalPoints(); ++i) {
106
                        Scalar phi_i = faceElem.Value(q, i);
niklas.baumgarten's avatar
niklas.baumgarten committed
107
                        for (int j = 0; j < faceElem.NodalPoints(); ++j) {
108
109
110
                            Scalar phi_j = faceElem.Value(q, j);
                            A_c(i, j) += w * BN * phi_j * phi_i;
                        }
niklas.baumgarten's avatar
niklas.baumgarten committed
111
                        for (int j = 0; j < felem_1.NodalPoints(); ++j) {
112
113
114
115
116
117
118
119
120
121
122
123
124
                            Scalar phi_j = felem_1.Value(q1, j);
                            A_cf(i, j) -= w * BN * phi_j * phi_i;
                        }
                    }
                }
                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);
125
                    int q1 = felem_1.findQPointID(faceElem, Qf_c);
126
                    double s = 1;
niklas.baumgarten's avatar
niklas.baumgarten committed
127
                    for (int i = 0; i < faceElem.NodalPoints(); ++i) {
128
129
130
                        Scalar phi_i = faceElem.Value(q, i);
                        Scalar NDphi_i = diffusion *
                            (faceElem.Derivative(q, i) * N);
niklas.baumgarten's avatar
niklas.baumgarten committed
131
                        for (int j = 0; j < faceElem.NodalPoints(); ++j) {
132
133
134
135
136
137
138
                            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);
                        }
niklas.baumgarten's avatar
niklas.baumgarten committed
139
                        for (int j = 0; j < felem_1.NodalPoints(); ++j) {
140
141
142
143
144
145
146
147
148
149
150
                            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);
                        }
                    }
niklas.baumgarten's avatar
niklas.baumgarten committed
151
                    for (int i = 0; i < felem_1.NodalPoints(); ++i) {
152
153
154
                        Scalar phi_i = felem_1.Value(q1, i);
                        Scalar NDphi_i = diffusion
                            * (felem_1.Derivative(q1, i) * N);
niklas.baumgarten's avatar
niklas.baumgarten committed
155
                        for (int j = 0; j < felem_1.NodalPoints(); ++j) {
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
                            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);
                        }
                    }
                }
            }
        }
    }
}

void DGTransportAssemble::RHS(double t, Vector &rhs) const {
    rhs = 0;
niklas.baumgarten's avatar
niklas.baumgarten committed
172
    if (!problem->RHS()) return;
173
174
175
    for (cell c = rhs.cells(); c != rhs.cells_end(); ++c) {
        row r = rhs.find_row(c());
        for (int f = 0; f < c.Faces(); ++f) {
176
            if (!rhs.GetMesh().onBndDG(c, f)) continue;
177
178
179
180
181
            DGFaceElement felem(*disc, rhs, c, f);
            for (int q = 0; q < felem.nQ(); ++q) {
                const Point &z = felem.QPoint(q);
                double w = felem.QWeight(q);
                VectorField N = felem.QNormal(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
182
                Scalar BN = problem->FaceNormalFlux(c, f, N, z);
183
                if (BN > 0) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
184
                Scalar U = problem->Solution(t, z);
niklas.baumgarten's avatar
niklas.baumgarten committed
185
                for (int i = 0; i < felem.NodalPoints(); ++i) {
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
                    Scalar Phi_i = felem.Value(q, i);
                    rhs(r, i) -= w * U * BN * Phi_i;
                }
            }
        }
    }
}

double DGTransportAssemble::Energy(const Vector &u) const {
    double energy = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        DGElement elem(*disc, u, c);
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            Scalar U = elem.Value(q, u);
            energy += w * (U * U);
        }
    }
    return 0.5 * PPM->Sum(energy);
}

void DGTransportAssemble::Energy(const cell &c, const Vector &u, double &energy) const {
    DGElement elem(*disc, u, c);
    for (int q = 0; q < elem.nQ(); ++q) {
        double w = elem.QWeight(q);
        Scalar U = elem.Value(q, u);
        energy += w * (U * U);
    }
}

double DGTransportAssemble::Error(double t, const Vector &u) const {
    double err = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        DGElement elem(*disc, u, c);
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            Scalar U = elem.Value(q, u);
niklas.baumgarten's avatar
niklas.baumgarten committed
223
            Scalar Sol = problem->Solution(t, elem.QPoint(q));
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
            err += w * (U - Sol) * (U - Sol);
        }
    }
    return sqrt(PPM->Sum(err));
}

std::pair<double, double> DGTransportAssemble::InFlowOutFlowRate(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;
        DGElement elem(*disc, u, c);
        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;
        }
        U *= (1 / area);
        for (int f = 0; f < c.Faces(); ++f) {
            if (bnd[f] == -1) continue;
            DGFaceElement faceElem(*disc, u, c, f);
            for (int q = 0; q < faceElem.nQ(); ++q) {
                double w = faceElem.QWeight(q);
                VectorField Nq = faceElem.QNormal(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
251
                Scalar BN = problem->FaceNormalFlux(c, f, Nq, elem.QPoint(q));
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
                if (BN > 0) outflow += w * BN * U;
                else inflow += w * BN * U;
            }
        }
    }
    inflow = PPM->Sum(inflow);
    outflow = PPM->Sum(outflow);
    if (inflow < 1e-7)
        inflow = 0.0;
    if (outflow < 1e-7)
        outflow = 0.0;
    return {inflow, outflow};
}

void DGTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
    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
270
271
272
        DGElement Elem(*disc, u_ex, c);
        for (int j = 0; j < Elem.NodalPoints(); ++j)
            u_ex(r, j) = problem->Solution(t, Elem.NodalPoint(j));
273
274
275
276
277
278
279
280
281
282
283
284
    }
    Accumulate(u_ex);
}

void DGTransportAssemble::AssembleTransfer(TransferMatrix &TM) const {
    TM = 0;
    const matrixgraph &cg = TM.CoarseMatrixGraph();
    for (cell c = cg.cells(); c != cg.cells_end(); ++c)
        for (int k = 0; k < c.Children(); ++k)
            TM(c(), c.Child(k))[0] = 1;
}

285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
double DGTransportAssemble::Mass(const Vector &u) const {
    Scalar e = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        DGElement elem(*disc, u, c);
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            Scalar U = elem.Value(q, u);
            e += w * U;
        }
    }
    return PPM->Sum(e);
}

void DGTransportAssemble::SetInitialValue(Vector &u) const {
    u = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        row r = u.find_row(c());
        DGElement Elem(*disc, u, c);
        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));
    }
    Accumulate(u);
}

void DGTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
    for (cell c = A.cells(); c != A.cells_end(); ++c) {
        DGElement elem(*disc, A, c);
        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;

        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;
        }
        mout << "------------------------" << endl;
    }
}

void DGTransportAssemble::VtkPlotting_cell(double t,
                                           const Vector &u,
                                           char *filename) const {
    Vector u_tmp(u);
    for (cell c = u_tmp.cells(); c != u_tmp.cells_end(); ++c) {
        DGElement elem(*disc, u_tmp, c);
        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;
        }
        U *= (1 / a);
        row r = u.find_row(c());
        u_tmp(r)[0] = U;
    }
    plot->celldata(u_tmp, 1);
    NumberName("U", filename, step);
    plot->vtk_celldata(filename, 0);
    plot->gnu_celldata(filename, 0);
}

void DGTransportAssemble::VtkPlotting(double t, const Vector &u) const {
    // Todo, refactor plot class should be able to handle everything
    const Mesh &Mp = plot->GetMesh();
    const Mesh &Mu = u.GetMesh();
    if (Mp.CellCount() != Mu.CellCount()) return;
    vout (2) << "  => VtkPlotting result of step " << step << ".\n";
    char filename[128];
    VtkPlotting_cell(t, u, filename);
}

void DGTransportAssemble::PrintInfo() const {
    mout.PrintInfo("Assemble", verbose,
                   PrintInfoEntry("Name", Name()),
                   PrintInfoEntry("Problem", problem->Name()),
                   PrintInfoEntry("Discretization", disc->DiscName()));
}

void DGTransportAssemble::FinishTimeStep(double t, Vector &u, bool SpecialTimeStep) {
    step++;
    t = t;
    PrintInfo(u);
niklas.baumgarten's avatar
niklas.baumgarten committed
385
    VtkPlotting(t, u);
386
387
388
389
390
391
}

void DGTransportAssemble::Initialize(double t, Vector &u) {
    u = 0;
    SetInitialValue(u);
    PrintInfo(u);
niklas.baumgarten's avatar
niklas.baumgarten committed
392
    VtkPlotting(t, u);
393
394
}