MonteCarloElliptic.C 3.7 KB
Newer Older
1
#include "MonteCarloElliptic.h"
2

3

4
5
using namespace std;

6
void MonteCarloElliptic::Method() {
7
    logger->StartMethodMsg();
8
    if (!onlyFineLevel)
9
        logger->LogMsgv1("l: " + to_string(l)
10
                             + " dM: " + to_string(dM));
11
    else
12
        logger->LogMsgv1("l: " + to_string(l)
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
                             + " dM: " + to_string(dM)
                             + " onlyFineLevel: True");

    Vector finePermeability(permMatrixGraphs[l - pLevel]);
    Vector coarsePermeability(permMatrixGraphs[l - pLevel - 1]);
    Vector coarsePermeabilityUp(finePermeability);
    Vector fineSolution(solMatrixGraphs[l - pLevel]);
    Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
    Vector coarseSolutionUp(fineSolution);
    Vector fineFlux(fluxMatrixGraphs[l - pLevel]);
    Vector coarseFlux(fluxMatrixGraphs[l - pLevel - 1]);
    Vector coarseFluxUp(fineFlux);

    finePermeability = 0.0, coarsePermeability = 0.0, coarsePermeabilityUp = 0.0;
    fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;
    fineFlux = 0.0, coarseFlux = 0.0, coarseFluxUp = 0.0;

30
    double fineQ, coarseQ, fineCost, coarseCost;
31
    fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
32
33

    for (int m = M; m < M + dM; m++) {
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
        stochasticField->GetFineSample(l, finePermeability);
        assemble->problem->LoadNewSample(&finePermeability);
        solveElliptic(l, fineQ, fineCost, fineSolution);
        assemble->setFlux(fineSolution, fineFlux);

        if (plotting > 1)
            plotResults(l, m, "", fineSolution, fineFlux, finePermeability);

        if (!onlyFineLevel) {
            stochasticField->GetCoarseSample(l - 1, finePermeability, coarsePermeability);
            assemble->problem->LoadNewSample(&coarsePermeability);
            solveElliptic(l, coarseQ, coarseCost, coarseSolution);
            assemble->setFlux(coarseSolution, coarseFlux);

            UpscaleSolution(coarseSolution, coarseSolutionUp);
            UpscalePermeability(coarsePermeability, coarsePermeabilityUp);
            UpscaleFlux(coarseFlux, coarseFluxUp);

            if (plotting > 1)
                plotResults(l, m, "up", coarseSolutionUp,
                            coarseFluxUp, coarsePermeabilityUp);
        } else {
            coarseQ = 0.0, coarseCost = 0.0, coarseSolutionUp = 0.0;
57
        }
58
59
60
61
        updateSums(fineCost, fineQ - coarseQ, fineQ);
    }
    updateAvg();
    updateStatistic();
62

63
    logger->LogMsgv1("|E[Y]|: " + to_string(avgY)
64
                         + " V[Y]: " + to_string(varY));
65
    logger->EndMethodMsg();
66
}
67

68
69
70
void MonteCarloElliptic::solveElliptic(int l, double &valueQ, double &cost,
                                       Vector &solution) {
    Date startSolving;
71

72
73
74
75
    Solver solver(GetPC(solMatrixGraphs, *assemble, "LIB_PS"), "GMRES");
    Newton newton(solver);
    newton(*assemble, solution);
    cost = solution.pSize();
76

77
78
79
80
81
    if (functional == "L2") valueQ = assemble->L2(solution);
    if (functional == "Energy") valueQ = assemble->Energy(solution);
    if (functional == "Outflow")
        valueQ = assemble->getInflowOutflow(solution).second;
    if (functional == "Goal") valueQ = assemble->getGoal(solution);
82

83
    logger->LogMsgv2("Solving on level " + to_string(l) + " took " +
84
85
        to_string((startSolving - Date()).t) + "s");
}
86

87
88
89
void MonteCarloElliptic::plotResults(int l, int m, const string &suffix,
                                     const Vector &u, const Vector &flux,
                                     const Vector &perm) {
90

91
    string name = buildName(m, "u", suffix);
92
    assemble->plotU(plot, u, name);
93

94
    name = buildName(m, "flux", suffix);
95
    assemble->plotFlux(plot, flux, name);
96

97
    name = buildName(m, "perm", suffix);
98
    assemble->plotPerm(plot, perm, name);
99
}