MonteCarloElliptic.C 3.64 KB
Newer Older
1
#include "MonteCarloElliptic.h"
2
3
4
5

using namespace std;

void MonteCarloElliptic::method() {
6
7
8
9
10
11
12
13
14
    logger->startMethodMSG();
    if (!baseLevel)
        logger->logMSGV1("l: " + to_string(l)
                         + " dM: " + to_string(dM));
    else
        logger->logMSGV1("l: " + to_string(l)
                         + " dM: " + to_string(dM)
                         + " baseLevel: True");

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
    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;
40
        logger->logMSGV2("solve PDE on level " + to_string(l));
41
42

        solveElliptic(Qf, costf, uf);
43
44
45
46
47
48
        vout(2) << " took " << Start_solving - Date() << endl;

        if (!baseLevel) {
            Date Start_solving_c;
            vout(2) << "    Solving PDE on level " << l - 1;
            stochFields->setFineField(false);
49
            solveElliptic(Qc, costc, uc);
50
51
52
53
54
55
56
57
58
59
60
            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)};
61

62
        updateSum(sums);
63
        upscale(uc, uc_up);
64
65
        updateUSum(uc_up, uf);

66
        if (plotting >= 2) {
67
68
            Vector fluxf(flux_sum), permf(flux_sum);;
            fluxf = 0.0, permf = 0.0;
69
70
71
72
            assemble->setFlux(uf, fluxf);
            assemble->setPerm(permf);
            plotSolution(uf, fluxf, permf, m, "");
        }
73
        if (plotting >= 3 && !baseLevel) {
74
75
            Vector fluxc_up(flux_sum), permc_up(flux_sum);
            fluxc_up = 0.0, permc_up = 0.0;
76
            stochFields->setFineField(false);
77
78
79
            assemble->setFlux(uc_up, fluxc_up);
            assemble->setPerm(permc_up);
            plotSolution(uc_up, fluxc_up, permc_up, m, "up");
80
81
82
83
84
85
86
87
88
            stochFields->setFineField(true);
        }
    }

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

    if (plotting >= 1) {
89
90
91
92
        plotSolution(u_avg_diff, flux_avg_diff, perm_avg_diff,
                     -1, "avg_diff");
        plotSolution(u_avg, flux_avg, perm_avg,
                     -1, "avg");
93
94
    }

95
96
97
    logger->logMSGV1("|E[Y]|: " + to_string(avg_Y)
                     + " V[Y]: " + to_string(var_Y));
    logger->endMethodMSG();
98
99
}

100
void MonteCarloElliptic::solveElliptic(double &Q, double &cost, Vector &u) {
101
102
103
104
105
106
107
108
109
    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;
niklas.baumgarten's avatar
niklas.baumgarten committed
110
    if (functional == "Goal") Q = assemble->getGoal(u);
111
}