MonteCarloTransport.cpp 2.97 KB
Newer Older
1
2
#include "MonteCarloTransport.h"

3

4
5
using namespace std;

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

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

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

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

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

28
        if (!onlyFineLevel) {
29

niklas.baumgarten's avatar
niklas.baumgarten committed
30
31
32
33
34
35
36
37
38
39
40
            sampleDirName = mkSampleDirName(m, false);
            mkSampleDir(m, sampleDirName);

            stochasticField->SetSampleDir(l, sampleDirName);
            stochasticField->SetPlot(l, coarsePlot);
            stochasticField->GetCoarseSample(l, fineNormalFlux, coarseNormalFlux);

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

niklas.baumgarten's avatar
niklas.baumgarten committed
41
            solvePDE(l - 1, m, false, coarseQ, coarseCost, coarseSolution);
niklas.baumgarten's avatar
niklas.baumgarten committed
42
        } else {
43

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

niklas.baumgarten's avatar
niklas.baumgarten committed
50
void MonteCarloTransport::solvePDE(int l,
51
                                   int m,
niklas.baumgarten's avatar
niklas.baumgarten committed
52
53
54
55
                                   bool fineLevel,
                                   double &valueQ,
                                   double &cost,
                                   Vector &solution) {
56
    logger->StartMethod("Start Solving PDE", verbose);
57
58
    logger->InnerMsg("l: " + to_string(l) +
        " m: " + to_string(m), verbose);
59

niklas.baumgarten's avatar
niklas.baumgarten committed
60
    Preconditioner *pc = GetPC("PointBlockGaussSeidel"); //TODO RMV after usage
61
    Solver solver(pc, "GMRES");
niklas.baumgarten's avatar
niklas.baumgarten committed
62
63
    Solver solver2(pc, "GMRES");  // TODO RMV
    DGTimeIntegrator timeIntegrator(solver, solver, solver, solver2);
64

65
66
    TimeSeries timeSeries = getTimeSeries(fineLevel);
    assemble->resetTAssemble();
67

niklas.baumgarten's avatar
niklas.baumgarten committed
68
    timeIntegrator(*assemble, timeSeries, solution);
69
70
71
72
73
74
    cost = solution.pSize() * assemble->getStep();

    if (functional == "Mass") valueQ = assemble->Mass(solution);
    if (functional == "Energy") valueQ = assemble->Energy(solution);
    if (functional == "Outflow")
        valueQ = assemble->InFlowOutFlowRate(solution).second;
niklas.baumgarten's avatar
niklas.baumgarten committed
75

76
    logger->EndMethod(verbose);
77
78
79
}

TimeSeries MonteCarloTransport::getTimeSeries(bool fineLevel) {
niklas.baumgarten's avatar
niklas.baumgarten committed
80
81
    int scaling = 1;
    config.get("scaling", scaling);
82
83
84
85
86
87
    if (fineLevel)
        return TimeSeries(startTime, endTime,
                          scaling * (int) (pow(2, l) * (endTime - startTime)));
    else
        return TimeSeries(startTime, endTime,
                          scaling * (int) (pow(2, l - 1) * (endTime - startTime)));
christian.wieners's avatar
.    
christian.wieners committed
88
}