DGTransportAssemble.h 15.3 KB
Newer Older
1
2
3
#ifndef TUTORIAL_DGTRANSPORT_HPP
#define TUTORIAL_DGTRANSPORT_HPP

4
#include "DGTAssemble.h"
5
#include "utils/IndentedLogger.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
6
#include <memory>
7

8

9
10
class DGTransportAssemble : public DGTAssemble {
public:
11
12
//    std::shared_ptr<IDiscretizationT<>> disc;
    std::shared_ptr<StochasticTransportProblem> problem;
13

14
    IndentedLogger *logger;
15

16
17
18
    double flux_alpha = 1.0;
    double diffusion = 0.0;

19
20
21
22
23
    // Todo Super Hacky!
    DGTransportAssemble(std::shared_ptr<IDiscretizationT<>> disc,
                        std::shared_ptr<StochasticTransportProblem> problem)
        : DGTAssemble(dynamic_cast<DGDiscretization *>(disc.get())),
          problem(move(problem)) {
24
        config.get("flux_alpha", flux_alpha);
25
        logger = IndentedLogger::GetInstance();
26
27
28
29
    }

    const char *Name() const override { return "DGTransportAssemble"; }

30
    double Residual(const Vector &u, Vector &b) const override {
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
57
58
59
60
61
62
63
64
65
66
67
        b = 0;

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

        Matrix *fluxMatrix = new Matrix(u);
        FluxMatrix(*fluxMatrix);

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

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

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

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

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

68
69
70
    void MassMatrix(Matrix &M) const override {
        M = 0;
        for (cell c = M.cells(); c != M.cells_end(); ++c) {
71
            DGElement elem(*disc, M, c);
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
            DGRowEntries M_c(M, c, c);
            for (int q = 0; q < elem.nQ(); ++q) {
                double w = elem.QWeight(q);
                for (int i = 0; i < elem.dimension(); ++i) {
                    Scalar Phi_i = elem.Value(q, i);
                    for (int j = 0; j < elem.dimension(); ++j) {
                        Scalar Phi_j = elem.Value(q, j);
                        M_c(i, j) += w * Phi_i * Phi_j;
                    }
                }
            }
        }
    }

    void FluxMatrix(Matrix &A) const override {
        A = 0;
        for (cell c = A.cells(); c != A.cells_end(); ++c) {
89
            DGElement elem(*disc, A, c);
90
91
92
            DGRowEntries A_c(A, c, c);
            for (int q = 0; q < elem.nQ(); ++q) {
                double w = elem.QWeight(q);
93
                VectorField B = problem->cellB(c, elem.QPoint(q));
94
95
96
97
98
99
100
101
102
                for (int i = 0; i < elem.dimension(); ++i) {
                    Scalar Phi_i = elem.Value(q, i);
                    for (int j = 0; j < elem.dimension(); ++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) {
103
                DGFaceElement faceElem(*disc, A, c, f);
104
105
                if (checkIfBND(A, c, f)) {
                    for (int q = 0; q < faceElem.nQ(); ++q) {
106
                        const Point &Qf_c = faceElem.QPoint(q);
107
                        VectorField Nq = faceElem.QNormal(q);
108
                        Scalar BN = problem->faceB(c, f, Nq, Qf_c);
109
110
111
112
113
114
115
116
117
118
119
120
121
                        if (BN > 0) continue;
                        double w = faceElem.QWeight(q);
                        for (int i = 0; i < faceElem.dimension(); ++i) {
                            Scalar phi_i = faceElem.Value(q, i);
                            for (int j = 0; j < faceElem.dimension(); ++j) {
                                Scalar phi_j = faceElem.Value(q, j);
                                A_c(i, j) += w * BN * phi_j * phi_i;
                            }
                        }
                    }
                } else {
                    cell cf = findNBCell(A, c, f);
                    int f1 = findNBFaceID(A, c.Face(f), cf);
122
                    DGFaceElement felem_1(*disc, A, cf, f1);
123
124
125
126
                    DGRowEntries A_cf(A, c, cf);
                    for (int q = 0; q < faceElem.nQ(); ++q) {
                        const Point &Qf_c = faceElem.QPoint(q);
                        VectorField Nq = faceElem.QNormal(q);
127
                        Scalar BN = problem->faceB(c, f, Nq, Qf_c);
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
                        if (BN > 0) continue;
                        int q1 = findQPointID(faceElem, felem_1, Qf_c);
                        double w = faceElem.QWeight(q);
                        for (int i = 0; i < faceElem.dimension(); ++i) {
                            Scalar phi_i = faceElem.Value(q, i);
                            for (int j = 0; j < faceElem.dimension(); ++j) {
                                Scalar phi_j = faceElem.Value(q, j);
                                A_c(i, j) += w * BN * phi_j * phi_i;
                            }
                            for (int j = 0; j < felem_1.dimension(); ++j) {
                                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(A, cf, c);
                    DGRowEntries A_ff(A, 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 = findQPointID(faceElem, felem_1, Qf_c);
                        double s = 1;
                        for (int i = 0; i < faceElem.dimension(); ++i) {
                            Scalar phi_i = faceElem.Value(q, i);
                            Scalar NDphi_i = diffusion *
156
                                (faceElem.Derivative(q, i) * N);
157
158
159
                            for (int j = 0; j < faceElem.dimension(); ++j) {
                                Scalar phi_j = faceElem.Value(q, j);
                                Scalar NDphi_j =
160
                                    diffusion * (faceElem.Derivative(q, j) * N);
161
                                A_c(i, j) -= w * (-0.5 * NDphi_i * phi_j
162
163
                                    - 0.5 * phi_i * NDphi_j
                                    + s * phi_i * phi_j);
164
165
166
167
                            }
                            for (int j = 0; j < felem_1.dimension(); ++j) {
                                Scalar phi_j = felem_1.Value(q1, j);
                                Scalar NDphi_j =
168
                                    diffusion * (felem_1.Derivative(q1, j) * N);
169
                                A_cf(i, j) -= w * (0.5 * NDphi_i * phi_j
170
171
                                    - phi_i * 0.5 * NDphi_j
                                    - s * phi_i * phi_j);
172
                                A_fc(j, i) -= w * (0.5 * NDphi_i * phi_j
173
174
                                    - phi_i * 0.5 * NDphi_j
                                    - s * phi_i * phi_j);
175
176
177
178
179
                            }
                        }
                        for (int i = 0; i < felem_1.dimension(); ++i) {
                            Scalar phi_i = felem_1.Value(q1, i);
                            Scalar NDphi_i = diffusion
180
                                * (felem_1.Derivative(q1, i) * N);
181
182
183
                            for (int j = 0; j < felem_1.dimension(); ++j) {
                                Scalar phi_j = felem_1.Value(q1, j);
                                Scalar NDphi_j = diffusion
184
                                    * (felem_1.Derivative(q1, j) * N);
185
                                A_ff(i, j) -= w * (0.5 * NDphi_i * phi_j
186
187
                                    + phi_i * 0.5 * NDphi_j
                                    + s * phi_i * phi_j);
188
189
190
191
192
193
194
195
196
197
                            }
                        }
                    }
                }
            }
        }
    }

    void RHS(Vector &b, double t) const override {
        b = 0;
198
        if (!problem->RHS()) return;
199
200
201
202
        for (cell c = b.cells(); c != b.cells_end(); ++c) {
            row r = b.find_row(c());
            for (int f = 0; f < c.Faces(); ++f) {
                if (!checkIfBND(b, c, f)) continue;
203
                DGFaceElement felem(*disc, b, c, f);
204
                for (int q = 0; q < felem.nQ(); ++q) {
205
                    const Point &z = felem.QPoint(q);
206
207
                    double w = felem.QWeight(q);
                    VectorField N = felem.QNormal(q);
208
                    Scalar BN = problem->faceB(c, f, N, z);
209
                    if (BN > 0) continue;
210
                    Scalar U = problem->ut(t, z);
211
212
213
214
215
216
217
218
219
220
221
222
223
224
                    for (int i = 0; i < felem.dimension(); ++i) {
                        Scalar Phi_i = felem.Value(q, i);
                        b(r, i) -= w * U * BN * Phi_i;
                    }
                }
            }
        }
    }

    void SetInitialValue(Vector &u) const override {
        u = 0;
        Point z[MaxNodalPoints];
        for (cell c = u.cells(); c != u.cells_end(); ++c) {
            row r = u.find_row(c());
225
            DGElement Elem(*disc, u, c);
226
227
228
            Elem.NodalPoints(z);
            double lambda = 1e-9;
            for (int j = 0; j < Elem.dimension(); ++j)
229
                u(r, j) = problem->ut(0, lambda * c() + (1 - lambda) * z[j]);
230
231
232
233
234
235
        }
        Accumulate(u);
    }

    void PrintStep(double t, const Vector &u) const override {
        pair<double, double> rate = InFlowOutFlowRate(u);
niklas.baumgarten's avatar
niklas.baumgarten committed
236
237
238
239
240
241

        string msg = "Step(";
        msg = msg.append(to_string(step)).append("), t(");
        msg = msg.append(to_string(t)).append("): Energy(u) = ");
        msg = msg.append(to_string(Energy(u))).append(" Mass(u) = ");
        msg = msg.append(to_string(Mass(u)));
242
        if (abs(rate.first) > 1e-7)
niklas.baumgarten's avatar
niklas.baumgarten committed
243
            msg = msg.append(" InFlowRate(u) = ").append(to_string(rate.first));
244
        if (abs(rate.second) > 1e-7)
niklas.baumgarten's avatar
niklas.baumgarten committed
245
246
            msg = msg.append(" OutFlowRate(u) = ").append(to_string(rate.second));

247
        logger->LogMsgFlush(msg, verbose);
248
249
250
    }

    void VtkPlotting(double t, const Vector &u) const override {
251
252
253
254
255
256
257
        // Todo
//        const Mesh &Mp = plot->GetMesh();
//        const Mesh &Mu = u.GetMesh();
//        if (Mp.Cells::size() != Mu.Cells::size()) return;
//        vout (2) << "  => VtkPlotting result of step " << step << ".\n";
//        char filename[128];
//        VtkPlotting_cell(t, u, filename);
258
259
    }

260
    void VtkPlotting_cell(double t, const Vector &u, char *filename) const {
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
        // Todo
//        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);
//        string name = "";
//        name = name.append(sampleDir).append("/U");
//        NumberName(const_cast<char *>(name.c_str()), filename, step);
//        plot->vtk_celldata(filename, 0);
//        plot->gnu_celldata(filename, 0);
282
283
    }

284
285
    void PrintMatrixInfo(Matrix &A, int diagonal) const override {
        for (cell c = A.cells(); c != A.cells_end(); ++c) {
286
            DGElement elem(*disc, A, c);
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
            DGRowEntries A_c(A, c, c);
            for (int i = 0; i < elem.dimension(); ++i) {
                for (int j = 0; j < elem.dimension(); ++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 (checkIfBND(A, c, f))
                    continue;
                else
                    cf = findNBCell(A, c, f);
                DGRowEntries A_cf(A, c, cf);
                for (int i = 0; i < elem.dimension(); ++i) {
                    for (int j = 0; j < elem.dimension(); ++j)
                        mout << A_cf(i, j) << " ";
                    mout << endl;
                }
                mout << endl;
            }
            mout << "------------------------" << endl;
        }
    }

    void AssembleTransfer(TransferMatrix &TM) const override {
        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;
    }

322
323
324
325
326
327
328
    void PrintInfo(const Vector &u) const override {
        Vout (1) << endl
                 << "Assemble Info:" << endl
                 << "  Name:            " << Name() << endl
                 << "  Problem:         " << problem->Name() << endl
                 << "  Problem size:    " << u.pSize() << endl
                 << "  Discretization:  " << disc->DiscName() << endl
329
330
331
332
333
334
                 << endl;
    }

    double Energy(const Vector &u) const override {
        double energy = 0;
        for (cell c = u.cells(); c != u.cells_end(); ++c) {
335
            DGElement elem(*disc, u, c);
336
337
338
339
340
341
            for (int q = 0; q < elem.nQ(); ++q) {
                double w = elem.QWeight(q);
                Scalar U = elem.Value(q, u);
                energy += w * (U * U);
            }
        }
342
        return 0.5 * PPM->Sum(energy);
343
344
345
346
347
    }

    double Error(double t, const Vector &u) const {
        Scalar e = 0;
        for (cell c = u.cells(); c != u.cells_end(); ++c) {
348
            DGElement elem(*disc, u, c);
349
350
            for (int q = 0; q < elem.nQ(); ++q) {
                double w = elem.QWeight(q);
351
                Scalar U = elem.Value(q, u) - problem->ut(t, elem.QPoint(q));
352
353
354
355
356
357
358
359
360
361
362
363
                e += w * U * U;
            }
        }
        return sqrt(PPM->Sum(e));
    }

    pair<double, double> InFlowOutFlowRate(const Vector &u) const {
        double inflow = 0;
        double outflow = 0;
        for (cell c = u.cells(); c != u.cells_end(); ++c) {
            BFParts bnd(u.GetMesh(), c);
            if (!bnd.onBnd()) continue;
364
            DGElement elem(*disc, u, c);
365
366
367
368
369
370
371
372
373
374
            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;
375
                DGFaceElement faceElem(*disc, u, c, f);
376
377
378
                for (int q = 0; q < faceElem.nQ(); ++q) {
                    double w = faceElem.QWeight(q);
                    VectorField Nq = faceElem.QNormal(q);
379
                    Scalar BN = problem->faceB(c, f, Nq, elem.QPoint(q));
380
381
382
383
384
385
386
387
388
389
                    if (BN > 0) outflow += w * BN * U;
                    else inflow += w * BN * U;
                }
            }
        }
        return {PPM->Sum(inflow), PPM->Sum(outflow)};
    }
};

#endif //TUTORIAL_DGTRANSPORT_HPP