MonteCarloElliptic.C 2.69 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
        assemble->sampleDir = sampleDirName;
        assemble->plot = finePlot;
        assemble->problem->LoadNewSample(&finePerm);

35
        solvePDE(l, 0, true, fineQ, fineCost, fineSolution);
niklas.baumgarten's avatar
niklas.baumgarten committed
36
37

        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
            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);

51
            solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution);
niklas.baumgarten's avatar
niklas.baumgarten committed
52
53
54

            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

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

72
    newton(*assemble, solution);
niklas.baumgarten's avatar
niklas.baumgarten committed
73

74
    cost = solution.pSize();
75

76
77
78
79
80
    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
81
}