MonteCarloElliptic.C 2.94 KB
Newer Older
1
#include "MonteCarloElliptic.h"
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
49
50
51
52
53
54
55
56
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
82
83
84
85
86
87
88
89
90
91

using namespace std;

void MonteCarloElliptic::method() {
    Date Start;
    vout(1) << "  Start Monte Carlo method on level " << l;
    if (!baseLevel) vout(1) << " and level " << l - 1;
    vout(1) << " with " << dM << " samples" << endl;
    stochFields->setCurrentLevel(l);

    for (int m = M; m < M + dM; m++) {
        stochFields->generate();

        if (verbose >= 2) {
            vout(2) << "    Generated field:" << endl;
            pair<double, double> min_max_f, min_max_c;
            min_max_f = stochFields->getMinMaxValues(true);
            vout(2) << "      min_f=" << min_max_f.first
                    << " max_f=" << min_max_f.second << endl;
            if (!baseLevel) {
                min_max_c = stochFields->getMinMaxValues(false);
                vout(2) << "      min_c=" << min_max_c.first
                        << " max_c=" << min_max_c.second << endl;
            }
        }

        double Qf, Qc, costf, costc;;
        Qf = 0.0, Qc = 0.0, costf = 0.0, costc = 0.0;

        Vector uf(u_sum), uc_up(u_sum), uc((*graphs)[l - plevel - 1]);
        uf = 0.0, uc = 0.0, uc_up = 0.0;

        Date Start_solving;
        vout(2) << "    Solving PDE on level " << l;
        solvePDE(Qf, costf, uf);
        vout(2) << " took " << Start_solving - Date() << endl;

        if (!baseLevel) {
            Date Start_solving_c;
            vout(2) << "    Solving PDE on level " << l - 1;
            stochFields->setFineField(false);
            solvePDE(Qc, costc, uc);
            stochFields->setFineField(true);
            vout(2) << " took " << Start_solving_c - Date() << endl;
        } else {
            Qc = 0.0;
            costc = 0.0;
            uc = 0.0;
        }

        double sums[7] = {costf, Qf - Qc, pow((Qf - Qc), 2),
                          pow((Qf - Qc), 3), pow((Qf - Qc), 4),
                          Qf, pow(Qf, 2)};
        updateSum(sums);
        upscaleU(uc, uc_up);
        updateUSum(uc_up, uf);

        if (plotting >= 2)
            plotSolution(uf, m, "");
        if (plotting >= 3 && !baseLevel) {
            stochFields->setFineField(false);
            plotSolution(uc_up, m, "up");
            stochFields->setFineField(true);
        }
    }

    updateAvg();
    updateUAvg();
    updateStatistic();

    if (plotting >= 1) {
        plotSolution(u_avg_diff, -1, "avg_diff");
        plotSolution(u_avg, -1, "avg");
    }

    vout (1) << "  End after " << Date() - Start << ": |E(Y)|=" << avg_Y
             << " V(Y)=" << var_Y << endl;
}

void MonteCarloElliptic::solvePDE(double &Q, double &cost, Vector &u) {
    Preconditioner *preconditioner = GetPC(*graphs, *assemble, "LIB_PS");
    Solver solver(preconditioner, "GMRES");
    Newton newton(solver);
    newton(*assemble, u);
    cost = u.pSize();

    if (functional == "L2") Q = assemble->L2(u);
    if (functional == "Energy") Q = assemble->Energy(u);
    if (functional == "Outflow") Q = assemble->getInflowOutflow(u).second;
}