Commit 2027608e authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored plotter

parent d6d9f6dd
...@@ -14,10 +14,10 @@ degree = 1 ...@@ -14,10 +14,10 @@ degree = 1
plevel = 2 plevel = 2
Functional = L2 Functional = L2
#functional = Energy #Functional = Energy
#funcitonal = H1 #Funcitonal = H1
#functional = Outflow #Functional = Outflow
#functional = Goal #Functional = Goal
# ----- Multilevel Monte Carlo ----- # ----- Multilevel Monte Carlo -----
maxLevel = 9 maxLevel = 9
......
...@@ -10,6 +10,7 @@ void MainProgram::Initialize() { ...@@ -10,6 +10,7 @@ void MainProgram::Initialize() {
problem = createStochasticProblem(); problem = createStochasticProblem();
assemble = createAssemble(); assemble = createAssemble();
plotter = make_shared<MultilevelPlotter>(*meshes);
mlmc = make_unique<MultilevelMonteCarlo>(initLevels, initSampleAmount, mlmc = make_unique<MultilevelMonteCarlo>(initLevels, initSampleAmount,
meshes.get(), stochField.get(), meshes.get(), stochField.get(),
assemble.get()); assemble.get());
......
...@@ -58,7 +58,7 @@ void MonteCarlo::mkSampleDir(int m, const string &sampleDirName) { ...@@ -58,7 +58,7 @@ void MonteCarlo::mkSampleDir(int m, const string &sampleDirName) {
string pathToData = "data/vtk/"; string pathToData = "data/vtk/";
char *sampleDirNameCharArr = char *sampleDirNameCharArr =
const_cast<char *> ((pathToData.append(sampleDirName)).c_str()); const_cast<char *> ((pathToData.append(sampleDirName)).c_str());
int status = mkdir(sampleDirNameCharArr, 0777); mkdir(sampleDirNameCharArr, 0777);
} }
string MonteCarlo::buildName(int m, const string &name, const string &suffix) { string MonteCarlo::buildName(int m, const string &name, const string &suffix) {
......
...@@ -4,7 +4,6 @@ ...@@ -4,7 +4,6 @@
#include "m++.h" #include "m++.h"
#include "stochastics/StochasticField.h" #include "stochastics/StochasticField.h"
#include "utils/IndentedLogger.hpp" #include "utils/IndentedLogger.hpp"
#include <sys/stat.h>
#include <utility> #include <utility>
#include <utils/MultilevelPlotter.hpp> #include <utils/MultilevelPlotter.hpp>
......
...@@ -4,69 +4,57 @@ ...@@ -4,69 +4,57 @@
using namespace std; using namespace std;
void MonteCarloElliptic::Initialize() { void MonteCarloElliptic::Initialize() {
finePermeability = 0.0, coarsePermeability = 0.0; fineInput = 0.0, coarseInput = 0.0;
fineSolution = 0.0, coarseSolution = 0.0; fineSolution = 0.0, coarseSolution = 0.0;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0; fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
} }
void MonteCarloElliptic::method() { void MonteCarloElliptic::method() {
Initialize();
for (int m = M; m < M + dM; m++) { for (int m = M; m < M + dM; m++) {
// Refactor dir stuff
string sampleDirName = mkSampleDirName(m, true); string sampleDirName = mkSampleDirName(m, true);
mkSampleDir(m, sampleDirName); plotter->SetDirectory(sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->GetFineSample(l, finePermeability);
// assemble->sampleDir = sampleDirName;
// assemble->plot = finePlot;
assemble->problem->LoadNewSample(&finePermeability);
stochasticField->GetFineSample(l, fineInput);
assemble->problem->SetSample(&fineInput);
solvePDE(l, m, true, fineQ, fineCost, fineSolution); solvePDE(l, m, true, fineQ, fineCost, fineSolution);
// if (plotting) assemble->PlotResults(fineSolution); if (plotting)
plotter->PlotVector("u", fineSolution,
if (!onlyFineLevel) { 1, l, "VertexData");
sampleDirName = mkSampleDirName(m, false);
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
// stochasticField->SetPlot(l, coarsePlot);
stochasticField->GetCoarseSample(l, finePermeability, coarsePermeability);
// assemble->sampleDir = sampleDirName; if (onlyFineLevel)
// assemble->plot = coarsePlot; coarseQ = 0.0, coarseCost = 0.0;
assemble->problem->LoadNewSample(&coarsePermeability);
solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution); else {
sampleDirName = mkSampleDirName(m, false);
plotter->SetDirectory(sampleDirName);
// if (plotting) assemble->PlotResults(coarseSolution); stochasticField->GetCoarseSample(l, fineInput, coarseInput);
assemble->problem->SetSample(&coarseInput);
solvePDE(l, 0, false,
coarseQ, coarseCost, coarseSolution);
} else { if (plotting)
coarseQ = 0.0, coarseCost = 0.0; plotter->PlotVector("u", coarseSolution,
1, l - 1, "VertexData");
} }
updateSums(fineCost, fineQ - coarseQ, fineQ); updateSums(fineCost, fineQ - coarseQ, fineQ);
} }
} }
void MonteCarloElliptic::solvePDE(int l, void MonteCarloElliptic::solvePDE(int l, int m, bool fineLevel,
int m, double &valueQ, double &cost, Vector &solution) {
bool fineLevel, // Todo refactor
double &valueQ,
double &cost,
Vector &solution) {
Preconditioner *pc = GetPC("SuperLU"); Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES"); Solver solver(pc, "GMRES");
Newton newton(solver); Newton newton(solver);
newton(*assemble, solution); newton(*assemble, solution);
cost = solution.pSize(); cost = solution.pSize();
if (functional == "L2") valueQ = assemble->L2(solution); if (functional == "L2") valueQ = assemble->L2(solution);
if (functional == "Energy") valueQ = assemble->Energy(solution); if (functional == "Energy") valueQ = assemble->Energy(solution);
if (functional == "Outflow") if (functional == "Outflow") valueQ = assemble->getInflowOutflow(solution).second;
valueQ = assemble->getInflowOutflow(solution).second;
if (functional == "Goal") valueQ = assemble->getGoal(solution); if (functional == "Goal") valueQ = assemble->getGoal(solution);
} }
\ No newline at end of file
...@@ -24,8 +24,8 @@ public: ...@@ -24,8 +24,8 @@ public:
MatrixGraphs cellMatrixGraphs; MatrixGraphs cellMatrixGraphs;
MatrixGraphs solMatrixGraphs; MatrixGraphs solMatrixGraphs;
Vector finePermeability; Vector fineInput;
Vector coarsePermeability; Vector coarseInput;
Vector fineSolution; Vector fineSolution;
Vector coarseSolution; Vector coarseSolution;
...@@ -39,8 +39,8 @@ public: ...@@ -39,8 +39,8 @@ public:
assemble(assemble), assemble(assemble),
cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))), cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)), solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)),
finePermeability(Vector(cellMatrixGraphs[l - pLevel])), fineInput(Vector(cellMatrixGraphs[l - pLevel])),
coarsePermeability(Vector(cellMatrixGraphs[l - pLevel - 1])), coarseInput(Vector(cellMatrixGraphs[l - pLevel - 1])),
fineSolution(Vector(solMatrixGraphs[l - pLevel])), fineSolution(Vector(solMatrixGraphs[l - pLevel])),
coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) { coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) {
......
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
using namespace std; using namespace std;
void MonteCarloTransport::Initialize() { void MonteCarloTransport::Initialize() {
fineNormalFlux = 0.0, coarseNormalFlux = 0.0; fineInput = 0.0, coarseInput = 0.0;
fineSolution = 0.0, coarseSolution = 0.0; fineSolution = 0.0, coarseSolution = 0.0;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0; fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
} }
...@@ -13,40 +13,38 @@ void MonteCarloTransport::method() { ...@@ -13,40 +13,38 @@ void MonteCarloTransport::method() {
Initialize(); Initialize();
for (int m = M; m < M + dM; m++) { for (int m = M; m < M + dM; m++) {
string sampleDirName = mkSampleDirName(m, true); string sampleDirName = mkSampleDirName(m, true);
mkSampleDir(m, sampleDirName); plotter->SetDirectory(sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->GetFineSample(l, fineNormalFlux);
assemble->problem->LoadNewSample(&fineNormalFlux);
stochasticField->GetFineSample(l, fineInput);
assemble->problem->SetSample(&fineInput);
solvePDE(l, m, true, fineQ, fineCost, fineSolution); solvePDE(l, m, true, fineQ, fineCost, fineSolution);
if (!onlyFineLevel) { // if (plotting)
// plotter->PlotVector("u", fineSolution,
sampleDirName = mkSampleDirName(m, false); // 1, l, "VertexData");
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName); if (onlyFineLevel)
stochasticField->GetCoarseSample(l, fineNormalFlux, coarseNormalFlux); coarseQ = 0.0, coarseCost = 0.0;
assemble->problem->LoadNewSample(&coarseNormalFlux); else {
sampleDirName = mkSampleDirName(m, false);
plotter->SetDirectory(sampleDirName);
solvePDE(l - 1, m, false, coarseQ, coarseCost, coarseSolution); stochasticField->GetCoarseSample(l, fineInput, coarseInput);
} else { assemble->problem->SetSample(&coarseInput);
solvePDE(l - 1, m, false,
coarseQ, coarseCost, coarseSolution);
coarseQ = 0.0, coarseCost = 0.0; // if (plotting)
// plotter->PlotVector("u", coarseSolution,
// 1, l - 1, "VertexData");
} }
updateSums(fineCost, fineQ - coarseQ, fineQ); updateSums(fineCost, fineQ - coarseQ, fineQ);
} }
} }
void MonteCarloTransport::solvePDE(int l, void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel,
int m, double &valueQ, double &cost, Vector &solution) {
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) {
logger->StartMethod("Start Solving PDE", verbose); logger->StartMethod("Start Solving PDE", verbose);
logger->InnerMsg("l: " + to_string(l) + logger->InnerMsg("l: " + to_string(l) +
" m: " + to_string(m), verbose); " m: " + to_string(m), verbose);
......
...@@ -33,8 +33,8 @@ public: ...@@ -33,8 +33,8 @@ public:
MatrixGraphs faceMatrixGraphs; MatrixGraphs faceMatrixGraphs;
CellMatrixGraphs solMatrixGraphs; CellMatrixGraphs solMatrixGraphs;
Vector fineNormalFlux; Vector fineInput;
Vector coarseNormalFlux; Vector coarseInput;
Vector fineSolution; Vector fineSolution;
Vector coarseSolution; Vector coarseSolution;
...@@ -49,8 +49,8 @@ public: ...@@ -49,8 +49,8 @@ public:
cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))), cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
faceMatrixGraphs(MatrixGraphs(*meshes, dof("face", 3))), faceMatrixGraphs(MatrixGraphs(*meshes, dof("face", 3))),
solMatrixGraphs(*meshes, *assemble->disc), solMatrixGraphs(*meshes, *assemble->disc),
fineNormalFlux(Vector(faceMatrixGraphs[l - pLevel])), fineInput(Vector(faceMatrixGraphs[l - pLevel])),
coarseNormalFlux(Vector(faceMatrixGraphs[l - pLevel - 1])), coarseInput(Vector(faceMatrixGraphs[l - pLevel - 1])),
fineSolution(Vector(solMatrixGraphs[l - pLevel])), fineSolution(Vector(solMatrixGraphs[l - pLevel])),
coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) { coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) {
......
...@@ -12,7 +12,7 @@ public: ...@@ -12,7 +12,7 @@ public:
virtual ~StochasticProblem() = default;; virtual ~StochasticProblem() = default;;
void LoadNewSample(Vector *sample) { void SetSample(Vector *sample) {
fieldSample = sample; fieldSample = sample;
} }
......
...@@ -22,7 +22,9 @@ void CirculantEmbedding1D::generateFineSample(Vector &fineField) { ...@@ -22,7 +22,9 @@ void CirculantEmbedding1D::generateFineSample(Vector &fineField) {
internalCounter += 1; internalCounter += 1;
internalCounter %= 2; internalCounter %= 2;
if (plotting) generatePlot(fineField); if (plotting)
plotter->PlotVector("permeability", fineField,
1, l, "CellData");
} }
void CirculantEmbedding1D::generateCoarseSample(const Vector &fineField, void CirculantEmbedding1D::generateCoarseSample(const Vector &fineField,
...@@ -39,7 +41,9 @@ void CirculantEmbedding1D::generateCoarseSample(const Vector &fineField, ...@@ -39,7 +41,9 @@ void CirculantEmbedding1D::generateCoarseSample(const Vector &fineField,
coarseField(coarseRow, i) /= c.Children(); coarseField(coarseRow, i) /= c.Children();
} }
if (plotting) generatePlot(coarseField); if (plotting)
plotter->PlotVector("permeability", coarseField, 1,
l - 1, "CellData");
} }
vector<complex<double>> CirculantEmbedding1D::generateField() { vector<complex<double>> CirculantEmbedding1D::generateField() {
...@@ -167,7 +171,9 @@ void CirculantEmbedding2D::generateFineSample(Vector &fineField) { ...@@ -167,7 +171,9 @@ void CirculantEmbedding2D::generateFineSample(Vector &fineField) {
internalCounter += 1; internalCounter += 1;
internalCounter %= 2; internalCounter %= 2;
if (plotting) generatePlot(fineField); if (plotting)
plotter->PlotVector("permeability", fineField,
1, l, "CellData");
} }
void CirculantEmbedding2D::generateCoarseSample(const Vector &fineField, void CirculantEmbedding2D::generateCoarseSample(const Vector &fineField,
...@@ -184,7 +190,9 @@ void CirculantEmbedding2D::generateCoarseSample(const Vector &fineField, ...@@ -184,7 +190,9 @@ void CirculantEmbedding2D::generateCoarseSample(const Vector &fineField,
coarseField(coarseRow, i) /= c.Children(); coarseField(coarseRow, i) /= c.Children();
} }
if (plotting) generatePlot(coarseField); if (plotting)
plotter->PlotVector("permeability", coarseField,
1, l - 1, "CellData");
} }
vector<vector<complex<double>>> CirculantEmbedding2D::generateField() { vector<vector<complex<double>>> CirculantEmbedding2D::generateField() {
......
...@@ -12,15 +12,6 @@ protected: ...@@ -12,15 +12,6 @@ protected:
int internalCounter = 0; int internalCounter = 0;
RNManager *rnmg; RNManager *rnmg;
void generatePlot(const Vector &sample) override {
string namePerm = sampleDir + "/permeability";
// mout << OUT(sample) << endl;
plotter->PlotVector(namePerm, sample, 1, l, "CellData");
}
public: public:
double evtol = 1e-10; double evtol = 1e-10;
string fieldType = "LogNormal"; string fieldType = "LogNormal";
......
...@@ -5,7 +5,7 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) { ...@@ -5,7 +5,7 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) {
initialize(true); initialize(true);
permeabilityGenerator->SetSampleDir(sampleDir); permeabilityGenerator->SetSampleDir(sampleDir);
permeabilityGenerator->GetFineSample(*finePerm); permeabilityGenerator->GetFineSample(*finePerm);
assemble->problem->LoadNewSample(finePerm); assemble->problem->SetSample(finePerm);
Preconditioner *pc = GetPC("SuperLU"); Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES"); Solver solver(pc, "GMRES");
...@@ -15,14 +15,16 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) { ...@@ -15,14 +15,16 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) {
assemble->setFlux(*u, *flux); assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, fineField); assemble->SetNormalFlux(*u, fineField);
if (plotting) generatePlot(*flux); if (plotting)
plotter->PlotVector("flux", *flux,
3, l, "CellData");
} }
void HybridFluxGenerator::generateCoarseSample(const Vector &fineField, void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
Vector &coarseField) { Vector &coarseField) {
initialize(false); initialize(false);
permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm); permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm);
assemble->problem->LoadNewSample(coarsePerm); assemble->problem->SetSample(coarsePerm);
Preconditioner *pc = GetPC("SuperLU"); Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES"); Solver solver(pc, "GMRES");
...@@ -32,7 +34,9 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField, ...@@ -32,7 +34,9 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
assemble->setFlux(*u, *flux); assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, coarseField); assemble->SetNormalFlux(*u, coarseField);
if (plotting) generatePlot(*flux); if (plotting)
plotter->PlotVector("flux", *flux,
3, l - 1, "CellData");
} }
void HybridFluxGenerator::initialize(bool fineField) { void HybridFluxGenerator::initialize(bool fineField) {
......
...@@ -29,11 +29,6 @@ protected: ...@@ -29,11 +29,6 @@ protected:
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override; void generateCoarseSample(const Vector &fineField, Vector &coarseField) override;
void generatePlot(const Vector &sample) override {
string nameFlux = sampleDir + "/flux";
plotter->PlotVector(nameFlux, sample, 3, l, "CellData3");
}
public: public:
explicit HybridFluxGenerator(int l, Meshes &meshes) explicit HybridFluxGenerator(int l, Meshes &meshes)
: SampleGenerator(l, meshes), : SampleGenerator(l, meshes),
......
...@@ -11,7 +11,6 @@ ...@@ -11,7 +11,6 @@
class SampleGenerator { class SampleGenerator {
protected: protected:
int plotting = 0; int plotting = 0;
MultilevelPlotter *plotter = nullptr;
int verbose = 0; int verbose = 0;
IndentedLogger *logger = nullptr; IndentedLogger *logger = nullptr;
...@@ -22,18 +21,14 @@ protected: ...@@ -22,18 +21,14 @@ protected:
virtual void generateCoarseSample(const Vector &fineField, Vector &coarseField) = 0; virtual void generateCoarseSample(const Vector &fineField, Vector &coarseField) = 0;
virtual void generatePlot(const Vector &sample) = 0;
public: public:
int l; int l;
Meshes &meshes; Meshes &meshes;
explicit SampleGenerator(int l, Meshes &meshes) : l(l), meshes(meshes) { explicit SampleGenerator(int l, Meshes &meshes) : l(l), meshes(meshes) {
config.get("GeneratorPlotting", plotting); config.get("GeneratorPlotting", plotting);
config.get("GeneratorVerbose", verbose); config.get("GeneratorVerbose", verbose);
logger = IndentedLogger::GetInstance(); logger = IndentedLogger::GetInstance();
// plotter = MultilevelPlotter::GetInstance(meshes);
} }
void GetFineSample(Vector &fineField) { void GetFineSample(Vector &fineField) {
......
...@@ -27,10 +27,6 @@ public: ...@@ -27,10 +27,6 @@ public:
generators[l]->GetCoarseSample(fineField, coarseField); generators[l]->GetCoarseSample(fineField, coarseField);
} }
void SetSampleDir(int l, const string &sampleDirName) {
generators[l]->SetSampleDir(sampleDirName);
}
private: private:
SampleGenerator *getGenerator(int l) { SampleGenerator *getGenerator(int l) {
if (generatorName == "CirculantEmbedding") { if (generatorName == "CirculantEmbedding") {
......
...@@ -3,11 +3,16 @@ ...@@ -3,11 +3,16 @@
#include "Plot.h" #include "Plot.h"
#include <memory> #include <memory>
#include <sys/stat.h>
//#include "plot/VtuPlot.h" //#include "plot/VtuPlot.h"
class MultilevelPlotter { class MultilevelPlotter {
private: private:
const string defaultDir = "data/vtk/";
string currentDir = defaultDir;
string currentDirShort = "";
void initializeMap(const vector<int> &levels, const Meshes &meshes) { void initializeMap(const vector<int> &levels, const Meshes &meshes) {
for (auto &l: levels) { for (auto &l: levels) {
PlotMap[l] = new Plot(meshes[l], meshes.dim()); PlotMap[l] = new Plot(meshes[l], meshes.dim());
...@@ -23,34 +28,35 @@ public: ...@@ -23,34 +28,35 @@ public:
initializeMap(levels, meshes); initializeMap(levels, meshes);
} }