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

3

4
5
using namespace std;

niklas.baumgarten's avatar
niklas.baumgarten committed
6
7
8
9
10
11
12
void MonteCarloElliptic::initialize() {
}

void MonteCarloElliptic::method() {
    Vector finePerm(cellMatrixGraphs[l - pLevel]);
    Vector coarsePerm(cellMatrixGraphs[l - pLevel - 1]);
    Vector coarsePermeabilityUp(finePerm);
13
14
15
16
    Vector fineSolution(solMatrixGraphs[l - pLevel]);
    Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
    Vector coarseSolutionUp(fineSolution);

niklas.baumgarten's avatar
niklas.baumgarten committed
17
    finePerm = 0.0, coarsePerm = 0.0, coarsePermeabilityUp = 0.0;
18
19
    fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;

20
    double fineQ, coarseQ, fineCost, coarseCost;
21
    fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
22
23

    for (int m = M; m < M + dM; m++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
24
25
26
27
28
29
        string sampleDirName = mkSampleDirName(m, true);
        mkSampleDir(m, sampleDirName);

        stochasticField->SetSampleDir(l, sampleDirName);
        stochasticField->SetPlot(l, finePlot);
        stochasticField->GetFineSample(l, finePerm);
30

niklas.baumgarten's avatar
niklas.baumgarten committed
31
32
33
34
35
36
37
        assemble->sampleDir = sampleDirName;
        assemble->plot = finePlot;
        assemble->problem->LoadNewSample(&finePerm);

        solvePDE(l, true, fineQ, fineCost, fineSolution);

        if (plotting) assemble->PlotResults(fineSolution);
38
39

        if (!onlyFineLevel) {
niklas.baumgarten's avatar
niklas.baumgarten committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
            sampleDirName = mkSampleDirName(m, false);
            mkSampleDir(m, sampleDirName);

            stochasticField->SetSampleDir(l, sampleDirName);
            stochasticField->SetPlot(l, coarsePlot);
            stochasticField->GetCoarseSample(l, finePerm, coarsePerm);

            assemble->sampleDir = sampleDirName;
            assemble->plot = coarsePlot;
            assemble->problem->LoadNewSample(&coarsePerm);

            solvePDE(l, false, coarseQ, coarseCost, coarseSolution);

            if (plotting) assemble->PlotResults(coarseSolution);

55
56
        } else {
            coarseQ = 0.0, coarseCost = 0.0, coarseSolutionUp = 0.0;
57
        }
58
59
60
        updateSums(fineCost, fineQ - coarseQ, fineQ);
    }
}
61

niklas.baumgarten's avatar
niklas.baumgarten committed
62
63
64
65
void MonteCarloElliptic::solvePDE(int l, bool fineLevel, double &valueQ, double &cost,
                                  Vector &solution) {
    Preconditioner *pc = GetPC("SuperLU");
    Solver solver(pc, "GMRES");
66
    Newton newton(solver);
niklas.baumgarten's avatar
niklas.baumgarten committed
67

68
    newton(*assemble, solution);
niklas.baumgarten's avatar
niklas.baumgarten committed
69

70
    cost = solution.pSize();
71

72
73
74
75
76
    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);
niklas.baumgarten's avatar
niklas.baumgarten committed
77
}