MLMCMain.C 6.89 KB
Newer Older
1
2
#include "MLMCMain.h"

3

4
using namespace std;
5

6
7
8
void MLMCMain::initialize() {
    meshes = getMeshes();
    disc = getDiscretization();
9
10
    sampleMatrixGraphs = getSampleMatrixGraphs();
    stochasticField = getStochasticField();
11

12
13
    problem = getStochasticProblem();
    assemble = getAssemble();
14
    solutionMatrixGraphs = getSolutionMatrixGraphs();
15

16
    mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes,
17
                                    stochasticField, assemble, solutionMatrixGraphs);
18
}
19

20
21
Meshes *MLMCMain::getMeshes() {
    if (problemName.find("1D") != string::npos)
22
        return new Meshes("Line", pLevel, maxLevel);
23
    if (problemName.find("2D") != string::npos)
24
        return new Meshes("UnitSquare", pLevel, maxLevel);
25
26
    if (problemName.find("Square500") != string::npos)
        return new Meshes("Square500", pLevel, maxLevel);
niklas.baumgarten's avatar
niklas.baumgarten committed
27
    Exit("\nMesh name not found\n")
28
}
29

30
31
32
33
34
35
36
37
Discretization *MLMCMain::getDiscretization() {
    if (modelName == "LagrangeFEM") {
        if (degree == 1)
            return new Discretization("linear", 1, meshes->dim());
        if (degree == 2)
            return new Discretization("serendipity", 1, meshes->dim());
    }
    if (modelName == "MixedFEM" || modelName == "HybridFEM")
38
        return new Discretization("RT0_P0", 1, meshes->dim());
39
40
41
42
43
44
45
46
47
48
49
50
51
52
    if (modelName == "DGFEM")
        return new DGDiscretization(new DGDoF(degree));
    if (modelName == "EGFEM")
        return new EGDiscretization(new EGDoF(degree));
    if (modelName == "LinearWC")
        return nullptr;
    if (modelName == "QuadraticWC")
        return nullptr;
    if (modelName == "LinearDPG")
        return nullptr;
    if (modelName == "QuadraticDPG")
        return nullptr;
    if (modelName == "DGTransport")
        return new DGDiscretization(new DGDoF(degree));
niklas.baumgarten's avatar
niklas.baumgarten committed
53
    Exit("\nDiscretization not implemented yet\n")
niklas.baumgarten's avatar
niklas.baumgarten committed
54
}
55

56
MatrixGraphs *MLMCMain::getSampleMatrixGraphs() {
57
    if (problemName.find("Laplace") != string::npos)
58
        return new MatrixGraphs(*meshes, dof("cell", 1));
59
    if (problemName.find("Pollution") != string::npos)
60
        return new MatrixGraphs(*meshes, dof("cell", 3));
61
}
62

63
StochasticField *MLMCMain::getStochasticField() {
64
65
66
67
68
69
70
71
72
73
74
75
    if (problemName.find("Stochastic") != string::npos) {
        if (problemName.find("Laplace") != string::npos)
            return new StochasticField(*meshes, "CirculantEmbedding");
        if (problemName.find("Pollution") != string::npos)
            return new StochasticField(*meshes, "HybridFluxGenerator");
    }
    if (problemName.find("Deterministic") != string::npos) {
        if (problemName.find("Laplace") != string::npos)
            return nullptr;
        if (problemName.find("Pollution") != string::npos)
            return new StochasticField(*meshes, "DeterministicHybridFluxGenerator");
    }
76
}
77

78
79
80
81
82
StochasticProblem *MLMCMain::getStochasticProblem() {
    if (problemName == "StochasticLaplace1D")
        return new StochasticLaplace1D();
    if (problemName == "StochasticLaplace2D")
        return new StochasticLaplace2D();
83
84
    if (problemName == "StochasticPollution1D"
        || problemName == "DeterministicPollution1D")
85
        return new StochasticPollution1D();
86
87
    if (problemName == "StochasticPollution2D"
        || problemName == "DeterministicPollution2D")
88
        return new StochasticPollution2D();
89
90
91
92
    if (problemName == "StochasticPollutionCosHat2D"
        || problemName == "DeterministicPollutionCosHat2D"
        || problemName == "DeterministicPollutionCosHatSquare500")
        return new StochasticPollutionCosHat2D();
niklas.baumgarten's avatar
niklas.baumgarten committed
93
    Exit("\nStochastic problem not implemented yet\n")
94
95
}

96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
Assemble *MLMCMain::getAssemble() {
    if (problemName.find("Laplace") != string::npos) {
        auto ellipticProblem = dynamic_cast<StochasticEllipticProblem *>(problem);
        if (modelName == "LagrangeFEM")
            return new EllipticAssemble(disc, ellipticProblem);
        if (modelName == "MixedFEM")
            return new MixedEllipticAssemble(disc, ellipticProblem);
        if (modelName == "HybridFEM")
            return new HybridEllipticAssemble(disc, ellipticProblem);
        if (modelName == "DGFEM")
            return new DGEllipticAssemble(disc, ellipticProblem);
    } else if (problemName.find("Pollution") != string::npos) {
        auto transportProblem = dynamic_cast<StochasticTransportProblem *>(problem);
        Plot plot(meshes->fine(), meshes->dim(), 2);
        return new DGTransportAssemble(disc, transportProblem, &plot);
niklas.baumgarten's avatar
niklas.baumgarten committed
111
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
112
    Exit("\nAssembling for this problem not implemented yet\n")
niklas.baumgarten's avatar
niklas.baumgarten committed
113
}
114

115
MatrixGraphs *MLMCMain::getSolutionMatrixGraphs() {
116
117
118
119
    if (modelName == "LagrangeFEM" || modelName == "MixedFEM")
        return new MatrixGraphs(*meshes, *disc);
    if (modelName == "HybridFEM")
        return new MatrixGraphs(*meshes, dof("face", 1));
120
    if (modelName == "DGFEM" || "DGTransport")
121
122
        return new CellMatrixGraphs(*meshes, *disc);
    Exit("\nMatrix graph not implemented yet\n")
123
124
}

125
void MLMCMain::run() {
niklas.baumgarten's avatar
niklas.baumgarten committed
126
127
    printParameters();

128
129
130
    assemble->PrintInfo(solutionMatrixGraphs->fine());
    meshes->PrintInfo();

131
132
    if (experimentName == "ConvergenceTest") {
        headConvergenceTest();
133
134
135
136
        mlmc->Method();
        mlmc->ShowMCTable();
        mlmc->ShowKurtosisWarning();
        mlmc->ShowExponentResults();
137
138
139
        if (enablePython == 1) {
            system("python3 ../tools/plot_statistics.py ConvergenceTest");
        }
140
141
    } else if (experimentName == "MLMCExperiment") {
        headMLMCExperiment();
142
143
144
145
146
        mlmc->Method(epsilon);
        mlmc->ShowMCTable();
        mlmc->ShowResultsMLMCRun(epsilon);
        mlmc->ShowKurtosisWarning();
        mlmc->ShowExponentResults();
147
148
149
150
        if (enablePython == 1 && PPM->master()) {
            system("python3 ../tools/plot_statistics.py MLMCExperiment");
            system("python3 ../tools/plot_mlmc.py MLMCExperiment");
        }
151
152
    } else if (experimentName == "MLMCOverEpsilon") {
        headMLMCOverEpsilon();
153
154
155
        mlmc->Method(epsilonList);
        mlmc->ShowMCTable();
        mlmc->ShowResultsOverEpsilon(epsilonList);
156
157
158
159
        if (enablePython == 1 && PPM->master()) {
            system("python3 ../tools/plot_statistics.py MLMCOverEpsilon");
            system("python3 ../tools/plot_mlmc.py MLMCOverEpsilon");
        }
160
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
161
162
163
}

void MLMCMain::printParameters() {
niklas.baumgarten's avatar
niklas.baumgarten committed
164
    mout << endl;
165
166
167
168
    mout << "Problem   \t=\t" << problemName << endl;
    mout << "Model     \t=\t" << modelName << endl;
    mout << "Experiment\t=\t" << experimentName << endl;
    mout << "StochField\t=\t" << fieldName << endl;
169
170
    mout << "initLevels\t=\t" << vec2str(initLevels) << endl;
    mout << "initSample\t=\t" << vec2str(initSampleAmount) << endl;
niklas.baumgarten's avatar
niklas.baumgarten committed
171
    if (experimentName == "MLMCExperiment")
172
        mout << "eps       \t=\t" << epsilon << endl;
niklas.baumgarten's avatar
niklas.baumgarten committed
173
    if (experimentName == "MLMCOverEpsilon")
174
175
176
        mout << "eps_lst   \t=\t" << vec2str(epsilonList) << endl;
    mout << "maxLevel  \t=\t" << maxLevel << endl;
    mout << "pLevel    \t=\t" << pLevel << endl;
177
    mout << "enablePy  \t=\t" << enablePython << endl;
niklas.baumgarten's avatar
niklas.baumgarten committed
178
}