Commit 9fdb199a authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

template based sample generator adapted to stochastic problems and removed stocahstic field

parent 5f7070c3
......@@ -2,8 +2,8 @@ set(assemble
assemble/LagrangeEllipticAssemble.cpp
assemble/MixedEllipticAssemble.cpp
assemble/HybridEllipticAssemble.cpp
assemble/DGEllipticAssemble.cpp
assemble/DGTransportAssemble.cpp
assemble/PGReactionAssemble.cpp
assemble/DGReactionAssemble.cpp
# assemble/DGEllipticAssemble.cpp
# assemble/DGTransportAssemble.cpp
# assemble/PGReactionAssemble.cpp
# assemble/DGReactionAssemble.cpp
)
\ No newline at end of file
#include "CirculantEmbedding.hpp"
#include "main/MultilevelPlotter.hpp"
#include "dof/BasicDoFs.hpp"
void CirculantEmbedding1D::generateFineSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
MatrixGraphs cellMatrixGraphs(this->meshes, dof(new CellDoF(1)));
fineSample = new Vector(cellMatrixGraphs[id.level - this->meshes.pLevel()]);
coarseSample = nullptr;
void CirculantEmbedding1D::generateFineSample(Vector &fineField) {
if (internalCounter == 0)
fineComplexField = generateField();
for (cell c = fineField.cells(); c != fineField.cells_end(); c++) {
int i = floor(c()[0] * sqrt(fineField.size()));
for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
int i = floor(c()[0] * sqrt(fineSample->size()));
if (fieldType == "Gaussian") {
if (internalCounter == 0)
fineField(c(), 0) = fineComplexField[i].real() + mean;
fineSample->operator()(c(), 0) = fineComplexField[i].real() + mean;
else
fineField(c(), 0) = fineComplexField[i].imag() + mean;
fineSample->operator()(c(), 0) = fineComplexField[i].imag() + mean;
} else if (fieldType == "LogNormal") {
if (internalCounter == 0)
fineField(c(), 0) = exp(fineComplexField[i].real() + mean);
fineSample->operator()(c(), 0) = exp(fineComplexField[i].real() + mean);
else
fineField(c(), 0) = exp(fineComplexField[i].imag() + mean);
fineSample->operator()(c(), 0) = exp(fineComplexField[i].imag() + mean);
} else Exit("Has to be Gaussian or LogNormal")
}
internalCounter += 1;
internalCounter %= 2;
if (plotting)
plotter->PlotVector("kappa", fineField,
plotter->PlotVector("kappa", *fineSample,
1, l, "CellData");
}
void CirculantEmbedding1D::generateCoarseSample(const Vector &fineField,
Vector &coarseField) {
coarseField = 0;
for (cell c = coarseField.cells(); c != coarseField.cells_end(); c++) {
row coarseRow = coarseField.find_row(c());
void CirculantEmbedding1D::generateCoarseSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
MatrixGraphs cellMatrixGraphs(this->meshes, dof(new CellDoF(1)));
coarseSample = new Vector(cellMatrixGraphs[l - this->meshes.pLevel() - 1]);
(*coarseSample) = 0;
for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); c++) {
row coarseRow = coarseSample->find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = fineField.find_row(c.Child(k));
row fineRow = fineSample->find_row(c.Child(k));
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) += fineField(fineRow, i);
coarseSample->operator()(coarseRow, i) +=
fineSample->operator()(fineRow, i);
}
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) /= c.Children();
coarseSample->operator()(coarseRow, i) /= c.Children();
}
if (plotting)
plotter->PlotVector("kappa", coarseField, 1,
plotter->PlotVector("kappa", *coarseSample, 1,
l - 1, "CellData");
}
......@@ -149,49 +161,53 @@ vector<double> CirculantEmbedding1D::getSqrtEV() {
return sqrt_ev_return;
}
void CirculantEmbedding2D::generateFineSample(Vector &fineField) {
void CirculantEmbedding2D::generateFineSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
if (internalCounter == 0)
fineComplexField = generateField();
for (cell c = fineField.cells(); c != fineField.cells_end(); c++) {
int i = floor(c()[0] * sqrt(fineField.size()));
int j = floor(c()[1] * sqrt(fineField.size()));
for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
int i = floor(c()[0] * sqrt(fineSample->size()));
int j = floor(c()[1] * sqrt(fineSample->size()));
if (fieldType == "Gaussian") {
if (internalCounter == 0)
fineField(c(), 0) = fineComplexField[i][j].real() + mean;
fineSample->operator()(c(), 0) = fineComplexField[i][j].real() + mean;
else
fineField(c(), 0) = fineComplexField[i][j].imag() + mean;
fineSample->operator()(c(), 0) = fineComplexField[i][j].imag() + mean;
} else if (fieldType == "LogNormal") {
if (internalCounter == 0)
fineField(c(), 0) = exp(fineComplexField[i][j].real() + mean);
fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].real() + mean);
else
fineField(c(), 0) = exp(fineComplexField[i][j].imag() + mean);
fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].imag() + mean);
} else Exit("Has to be Gaussian or LogNormal")
}
internalCounter += 1;
internalCounter %= 2;
if (plotting)
plotter->PlotVector("kappa", fineField,
plotter->PlotVector("kappa", *fineSample,
1, l, "CellData");
}
void CirculantEmbedding2D::generateCoarseSample(const Vector &fineField,
Vector &coarseField) {
coarseField = 0;
for (cell c = coarseField.cells(); c != coarseField.cells_end(); ++c) {
row coarseRow = coarseField.find_row(c());
void CirculantEmbedding2D::generateCoarseSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
*coarseSample = 0;
for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
row coarseRow = coarseSample->find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = fineField.find_row(c.Child(k));
row fineRow = fineSample->find_row(c.Child(k));
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) += fineField(fineRow, i);
coarseSample->operator()(coarseRow, i) +=
fineSample->operator()(fineRow, i);
}
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) /= c.Children();
coarseSample->operator()(coarseRow, i) /= c.Children();
}
if (plotting)
plotter->PlotVector("kappa", coarseField,
plotter->PlotVector("kappa", *coarseSample,
1, l - 1, "CellData");
}
......
......@@ -8,7 +8,7 @@
#include "SampleGenerator.hpp"
class CirculantEmbedding : public SampleGenerator {
class CirculantEmbedding : public SampleGenerator<Vector> {
protected:
int internalCounter = 0;
RNManager *rnmg;
......@@ -30,20 +30,6 @@ public:
virtual ~CirculantEmbedding() { delete rnmg; };
};
class DeterministicPermeabilityGenerator : public CirculantEmbedding {
public:
DeterministicPermeabilityGenerator(int l, Meshes &meshes)
: CirculantEmbedding(l, meshes) {}
void generateFineSample(Vector &fineField) override {
fineField = 1.0;
};
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override {
coarseField = 1.0;
}
};
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public:
vector<complex<double>> fineComplexField;
......@@ -61,9 +47,13 @@ public:
numberOfCellsEmbedded = sqrtEigenvalues.size();
}
void generateFineSample(Vector &fineField) override;
void generateFineSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) override;
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override;
void generateCoarseSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) override;
private:
vector<complex<double>> generateField();
......@@ -92,9 +82,13 @@ public:
numberOfCellsEmbedded[1] = sqrtEigenvalues[0].size();
}
void generateFineSample(Vector &fineField) override;
void generateFineSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) override;
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override;
void generateCoarseSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) override;
private:
vector<vector<complex<double>>> fineComplexField;
......
#include "HybridFluxGenerator.hpp"
void HybridFluxGenerator::generateFineSample(Vector &fineField) {
template<typename Sample>
void HybridFluxGenerator<Sample>::generateFineSample(SampleID id,
Vector *fineSample,
Vector *coarseSample) {
initialize(true);
permeabilityGenerator->SetSampleDir(sampleDir);
permeabilityGenerator->GetFineSample(*finePerm);
assemble->problem->SetSample(finePerm.get());
newton(*assemble, *u);
assemble->problem->DrawSample(id);
newton->operator()(*assemble, *u);
assemble->SetFlux(*u, *flux);
assemble->SetNormalFlux(*u, fineField);
assemble->SetNormalFlux(*u, *fineSample);
// if (plotting)
// plotter->PlotVector("flux", *flux,
// 3, l, "CellData");
}
void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
Vector &coarseField) {
initialize(false);
permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm);
assemble->problem->SetSample(coarsePerm.get());
newton(*assemble, *u);
template<typename Sample>
void HybridFluxGenerator<Sample>::generateCoarseSample(SampleID id,
Vector *fineSample,
Vector *coarseSample) {
initialize(true);
assemble->problem->DrawSample(id);
newton->operator()(*assemble, *u);
assemble->SetFlux(*u, *flux);
assemble->SetNormalFlux(*u, coarseField);
assemble->SetNormalFlux(*u, *coarseSample);
// if (plotting)
// plotter->PlotVector("flux", *flux,
// 3, l - 1, "CellData");
}
void HybridFluxGenerator::initialize(bool fineField) {
template<typename Sample>
void HybridFluxGenerator<Sample>::initialize(bool fineField) {
if (fineField) {
flux = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel()]);
u = make_unique<Vector>(faceMatrixGraphs[l - meshes.pLevel()]);
flux = new Vector(cellMatrixGraphs[this->l - this->meshes.pLevel()]);
u = new Vector(faceMatrixGraphs[this->l - this->meshes.pLevel()]);
} else {
flux = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel() - 1]);
u = make_unique<Vector>(faceMatrixGraphs[l - meshes.pLevel() - 1]);
flux = new Vector(cellMatrixGraphs[this->l - this->meshes.pLevel() - 1]);
u = new Vector(faceMatrixGraphs[this->l - this->meshes.pLevel() - 1]);
}
*flux = 0.0, *u = 0.0;
}
......
......@@ -12,113 +12,75 @@
#include <memory>
class HybridFluxGenerator : public SampleGenerator {
template<typename Sample>
class HybridFluxGenerator : public SampleGenerator<Sample> {
private:
MatrixGraphs cellMatrixGraphs;
MatrixGraphs faceMatrixGraphs;
unique_ptr<Vector> finePerm;
unique_ptr<Vector> coarsePerm;
unique_ptr<Vector> flux;
unique_ptr<Vector> u;
Vector *finePerm;
Vector *coarsePerm;
Vector *flux;
Vector *u;
void initialize(bool fineField);
protected:
shared_ptr<RTLagrangeDiscretization> disc;
shared_ptr<StochasticEllipticProblem> problem;
unique_ptr<HybridEllipticAssemble> assemble;
unique_ptr<CirculantEmbedding> permeabilityGenerator;
unique_ptr<Preconditioner> pc;
IStochasticEllipticProblemT<Sample> *problem;
HybridEllipticAssemble *assemble;
RTLagrangeDiscretization *disc;
Preconditioner *pc;
Solver *solver;
Newton *newton;
Solver solver;
Newton newton;
void generateFineSample(SampleID id,
Vector *fineSample,
Vector *coarseSample) override;
void generateFineSample(Vector &fineField) override;
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override;
void generateCoarseSample(SampleID id,
Vector *fineSample,
Vector *coarseSample) override;
public:
explicit HybridFluxGenerator(int l, Meshes &meshes) :
SampleGenerator(l, meshes),
SampleGenerator<Sample>(l, meshes),
cellMatrixGraphs(MatrixGraphs(meshes, dof("cell", 3))),
faceMatrixGraphs(MatrixGraphs(meshes, dof("face", 1))),
pc(unique_ptr<Preconditioner>(GetPC("SuperLU"))),
solver(Solver(pc.get(), "GMRES")),
newton(Newton(solver)) {
disc(new RTLagrangeDiscretization(meshes, 0, 0)),
pc(GetPC("SuperLU")),
solver(new Solver(pc, "GMRES")),
newton(new Newton(*solver)) {
finePerm = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel()]);
finePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel()]);
if (l != meshes.pLevel())
coarsePerm = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel() - 1]);
coarsePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]);
}
};
class HybridFluxGenerator1D : public HybridFluxGenerator {
public:
explicit HybridFluxGenerator1D(int l, Meshes &meshes) :
HybridFluxGenerator(l, meshes) {
virtual ~HybridFluxGenerator() {
disc = make_shared<RTLagrangeDiscretization>(meshes, 0, 0);
problem = make_shared<StochasticLaplace1D>();
assemble = make_unique<HybridEllipticAssemble>(disc.get(), problem.get());
permeabilityGenerator = make_unique<CirculantEmbedding1D>(l, meshes);
}
void PrintInfo() const;
};
class HybridFluxGenerator2D : public HybridFluxGenerator {
template<typename Sample=Vector>
class HybridFluxGenerator1D : public HybridFluxGenerator<Sample> {
public:
explicit HybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes) {
disc = make_shared<RTLagrangeDiscretization>(meshes, 0,0);
problem = make_shared<StochasticLaplace2D>();
assemble = make_unique<HybridEllipticAssemble>(disc.get(), problem.get());
permeabilityGenerator = make_unique<CirculantEmbedding2D>(l, meshes);
explicit HybridFluxGenerator1D(int l, Meshes &meshes)
: HybridFluxGenerator<Sample>(l, meshes) {
this->problem = new StochasticLaplace1D(meshes);
this->assemble = new HybridEllipticAssemble(this->disc, this->problem);
}
};
/*
* This is class just serves as a helper to equip
* advection problems with the results of an hybrid
* flux computation.
*/
class HybridMainForPollutionProblem {
private:
std::shared_ptr<RTLagrangeDiscretization> disc;
std::shared_ptr<StochasticLaplace2D> problem;
std::shared_ptr<MatrixGraphs> mGraphs;
std::shared_ptr<Preconditioner> pc;
std::shared_ptr<Solver> solver;
std::shared_ptr<Newton> newton;
template<typename Sample=Vector>
class HybridFluxGenerator2D : public HybridFluxGenerator<Sample> {
public:
HybridMainForPollutionProblem(Meshes &meshes);
void Initialize(Meshes &meshes);
void PrintInfo() const;
void ResetVectors() const;
void Solve() const;
void EquipHybridFluxTransport(StochasticHybridFlux *hybridFluxTransport);
std::shared_ptr<Vector> u;
std::shared_ptr<Vector> normalFlux;
std::shared_ptr<HybridEllipticAssemble> assemble;
explicit HybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator<Sample>(l, meshes) {
this->problem = new StochasticLaplace2D(meshes);
this->assemble = new HybridEllipticAssemble(this->disc, this->problem);
}
};
#endif //HYBRIDFLUXGENERATOR_HPP
......@@ -4,19 +4,23 @@
#include <utility>
#include "utility/Config.hpp"
#include "Algebra.h"
#include "montecarlo/SampleID.hpp"
template<typename Sample>
class SampleGenerator {
protected:
int plotting = 0;
int verbose = 0;
string sampleDir = "";
virtual void generateFineSample(SampleID id,
Sample *&fineSample,
Sample *&coarseSample) = 0;
virtual void generateFineSample(Vector &fineField) = 0;
virtual void generateCoarseSample(const Vector &fineField, Vector &coarseField) = 0;
virtual void generateCoarseSample(SampleID id,
Sample *&fineSample,
Sample *&coarseSample) = 0;
public:
int l;
......@@ -27,23 +31,15 @@ public:
config.get("GeneratorVerbose", verbose);
}
void GetFineSample(Vector &fineField) {
void DrawSample(SampleID id, Sample *&fineSample, Sample *&coarseSample) {
mout.StartBlock("Sample Generator");
vout(1) << "fine=True" << endl;
generateFineSample(fineField);
vout(1) << id.Str() << endl;
if (!id.coarse)
generateFineSample(id, fineSample, coarseSample);
else
generateCoarseSample(id, fineSample, coarseSample);
mout.EndBlock(verbose == 0);
}
void GetCoarseSample(const Vector &fineField, Vector &coarseField) {
mout.StartBlock("Sample Generator");
vout(1) << "fine=False" << endl;
generateCoarseSample(fineField, coarseField);
mout.EndBlock(verbose == 0);
}
void SetSampleDir(const string &sampleDirName) {
sampleDir = sampleDirName;
}
};
#endif //SAMPLEGENERATOR_HPP
#ifndef STOCHASTICFIELD_HPP
#define STOCHASTICFIELD_HPP
#include <utility>
#include "CirculantEmbedding.hpp"
#include "HybridFluxGenerator.hpp"
class StochasticField {
public:
string generatorName = "";
map<int, SampleGenerator *> generators;
Meshes &meshes;
explicit StochasticField(Meshes &meshes, string generatorName) :
meshes(meshes), generatorName(std::move(generatorName)) {
for (int l = meshes.pLevel(); l <= meshes.Level(); l++)
generators[l] = getGenerator(l);
}
void PrintInfo() {
// Todo refactor
string stochField = "";
config.get("StochasticField", stochField);
string mean = "";
config.get("mean", mean);
string sigma = "";
config.get("sigma", sigma);
vector<double> lambda;
config.get("lambda", lambda);
string smoothing = "";
config.get("smoothing", smoothing);
mout.PrintInfo("Stochastic Field", 1,
PrintInfoEntry("Stochastic Field", stochField),
PrintInfoEntry("Generator Name", generatorName),
PrintInfoEntry("mean", mean),
PrintInfoEntry("sigma", sigma),
PrintInfoEntry("lambda", vec2str(lambda)),
PrintInfoEntry("smoothing", smoothing));
}
void GetFineSample(int l, Vector &fineField) {
generators[l]->GetFineSample(fineField);
}
void GetCoarseSample(int l, const Vector &fineField, Vector &coarseField) {
generators[l]->GetCoarseSample(fineField, coarseField);
}
private:
SampleGenerator *getGenerator(int l) {
if (generatorName == "CirculantEmbedding") {
if (meshes.dim() == 1)
return new CirculantEmbedding1D(l, meshes);
if (meshes.dim() == 2)
return new CirculantEmbedding2D(l, meshes);
}
if (generatorName == "HybridFluxGenerator") {
if (meshes.dim() == 1)
return new HybridFluxGenerator1D(l, meshes);
if (meshes.dim() == 2)
return new HybridFluxGenerator2D(l, meshes);
}
Exit("Generator does not exist")
}
};
#endif //STOCHASTICFIELD_HPP
Markdown is supported
0% or .