MonteCarloElliptic.cpp 2.32 KB
Newer Older
1
#include "MonteCarloElliptic.h"
2

3

4
5
using namespace std;

niklas.baumgarten's avatar
niklas.baumgarten committed
6
void MonteCarloElliptic::initialize() {
niklas.baumgarten's avatar
niklas.baumgarten committed
7
8
9
    finePermeability = 0.0, coarsePermeability = 0.0;
    fineSolution = 0.0, coarseSolution = 0.0;
    fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
10
11
12
}

void MonteCarloElliptic::method() {
13
    for (int m = M; m < M + dM; m++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
14
15
16
17
18
        string sampleDirName = mkSampleDirName(m, true);
        mkSampleDir(m, sampleDirName);

        stochasticField->SetSampleDir(l, sampleDirName);
        stochasticField->SetPlot(l, finePlot);
niklas.baumgarten's avatar
niklas.baumgarten committed
19
        stochasticField->GetFineSample(l, finePermeability);
20

niklas.baumgarten's avatar
niklas.baumgarten committed
21
22
        assemble->sampleDir = sampleDirName;
        assemble->plot = finePlot;
niklas.baumgarten's avatar
niklas.baumgarten committed
23
24
25
        assemble->problem->LoadNewSample(&finePermeability);

        solvePDE(l, m, true, fineQ, fineCost, fineSolution);
niklas.baumgarten's avatar
niklas.baumgarten committed
26
27

        if (plotting) assemble->PlotResults(fineSolution);
28
29

        if (!onlyFineLevel) {
niklas.baumgarten's avatar
niklas.baumgarten committed
30
31
32
33
34
            sampleDirName = mkSampleDirName(m, false);
            mkSampleDir(m, sampleDirName);

            stochasticField->SetSampleDir(l, sampleDirName);
            stochasticField->SetPlot(l, coarsePlot);
niklas.baumgarten's avatar
niklas.baumgarten committed
35
            stochasticField->GetCoarseSample(l, finePermeability, coarsePermeability);
niklas.baumgarten's avatar
niklas.baumgarten committed
36
37
38

            assemble->sampleDir = sampleDirName;
            assemble->plot = coarsePlot;
niklas.baumgarten's avatar
niklas.baumgarten committed
39
            assemble->problem->LoadNewSample(&coarsePermeability);
niklas.baumgarten's avatar
niklas.baumgarten committed
40

41
            solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution);
niklas.baumgarten's avatar
niklas.baumgarten committed
42
43
44

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

45
        } else {
niklas.baumgarten's avatar
niklas.baumgarten committed
46
            coarseQ = 0.0, coarseCost = 0.0;
47
        }
48
49
50
        updateSums(fineCost, fineQ - coarseQ, fineQ);
    }
}
51

52
53
54
55
56
void MonteCarloElliptic::solvePDE(int l,
                                  int m,
                                  bool fineLevel,
                                  double &valueQ,
                                  double &cost,
niklas.baumgarten's avatar
niklas.baumgarten committed
57
58
59
                                  Vector &solution) {
    Preconditioner *pc = GetPC("SuperLU");
    Solver solver(pc, "GMRES");
60
    Newton newton(solver);
niklas.baumgarten's avatar
niklas.baumgarten committed
61

62
    newton(*assemble, solution);
niklas.baumgarten's avatar
niklas.baumgarten committed
63

64
    cost = solution.pSize();
65

66
67
68
69
70
    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
71
}