Commit 36f64c32 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

added deterministic problems by deterministic generators

parent c485c506
...@@ -10,12 +10,14 @@ flux_alpha = 1 #Upwind: 1, Central: 0 ...@@ -10,12 +10,14 @@ flux_alpha = 1 #Upwind: 1, Central: 0
#penalty = 20 #penalty = 20
#Problem = StochasticLaplace1D #Problem = StochasticLaplace1D
#Problem = StochasticLaplace2D
#Problem = StochasticPollution1D #Problem = StochasticPollution1D
#Problem = DeterministicPollution1D
#Problem = StochasticLaplace2D
#Problem = StochasticPollution2D #Problem = StochasticPollution2D
Problem = StochasticPollution2DCosHat
#Problem = DeterministicPollution2D #Problem = DeterministicPollution2D
#Problem = DeterministicPollution2DCosHat #Problem = StochasticPollutionCosHat2D
#Problem = DeterministicPollutionCosHat2D
Problem = DeterministicPollutionCosHatSquare500
StochasticField = LogNormal StochasticField = LogNormal
#StochasticField = Gaussian #StochasticField = Gaussian
...@@ -24,14 +26,14 @@ StochasticField = LogNormal ...@@ -24,14 +26,14 @@ StochasticField = LogNormal
Experiment = MLMCExperiment Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon #Experiment = MLMCOverEpsilon
initLevels = 4, 5, 6 initLevels = 3, 4
initSampleAmount = 0, 0, 4 initSampleAmount = 1, 0
epsilon = 0.01 epsilon = 0.01
mcOnly = true mcOnly = true
epsilon_lst = 0.05, 0.01, 0.005 epsilon_lst = 0.05, 0.01, 0.005
maxLevel = 8 maxLevel = 4
plevel = 2 plevel = 2
degree = 2 degree = 2
...@@ -74,6 +76,7 @@ LinearSteps = 3000 ...@@ -74,6 +76,7 @@ LinearSteps = 3000
NewtonSteps = 1 NewtonSteps = 1
NewtonLineSearchSteps = 0 NewtonLineSearchSteps = 0
ConfigVerbose = 0
TimeLevel = -1 TimeLevel = -1
LinearVerbose = -1 LinearVerbose = -1
BaseSolverVerbose = -1 BaseSolverVerbose = -1
...@@ -99,4 +102,4 @@ rkorder = -2 #Implicit MP-rule ...@@ -99,4 +102,4 @@ rkorder = -2 #Implicit MP-rule
gamma = 0.01 gamma = 0.01
plotPermForTransport = false plotPermForTransport = true
\ No newline at end of file \ No newline at end of file
#include "MLMCMain.h" #include "MLMCMain.h"
using namespace std; using namespace std;
void MLMCMain::initialize() { void MLMCMain::initialize() {
...@@ -21,6 +22,8 @@ Meshes *MLMCMain::getMeshes() { ...@@ -21,6 +22,8 @@ Meshes *MLMCMain::getMeshes() {
return new Meshes("Line", pLevel, maxLevel); return new Meshes("Line", pLevel, maxLevel);
if (problemName.find("2D") != string::npos) if (problemName.find("2D") != string::npos)
return new Meshes("UnitSquare", pLevel, maxLevel); return new Meshes("UnitSquare", pLevel, maxLevel);
if (problemName.find("Square500") != string::npos)
return new Meshes("Square500", pLevel, maxLevel);
Exit("\nMesh name not found\n") Exit("\nMesh name not found\n")
} }
...@@ -51,20 +54,25 @@ Discretization *MLMCMain::getDiscretization() { ...@@ -51,20 +54,25 @@ Discretization *MLMCMain::getDiscretization() {
} }
MatrixGraphs *MLMCMain::getSampleMatrixGraphs() { MatrixGraphs *MLMCMain::getSampleMatrixGraphs() {
if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D") if (problemName.find("Laplace") != string::npos)
return new MatrixGraphs(*meshes, dof("cell", 1)); return new MatrixGraphs(*meshes, dof("cell", 1));
if (problemName == "StochasticPollution1D" || "StochasticPollution2D") if (problemName.find("Pollution") != string::npos)
return new MatrixGraphs(*meshes, dof("cell", 3)); return new MatrixGraphs(*meshes, dof("cell", 3));
} }
StochasticField *MLMCMain::getStochasticField() { StochasticField *MLMCMain::getStochasticField() {
if (problemName == "StochasticLaplace1D" || problemName == "StochasticLaplace2D") if (problemName.find("Stochastic") != string::npos) {
return new StochasticField(*meshes, "CirculantEmbedding"); if (problemName.find("Laplace") != string::npos)
if (problemName == "StochasticPollution1D" || problemName == "StochasticPollution2D" return new StochasticField(*meshes, "CirculantEmbedding");
|| problemName == "StochasticPollution2DCosHat" if (problemName.find("Pollution") != string::npos)
|| problemName == "DeterministicPollution2D" return new StochasticField(*meshes, "HybridFluxGenerator");
|| problemName == "DeterministicPollution2DCosHat") }
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");
}
} }
StochasticProblem *MLMCMain::getStochasticProblem() { StochasticProblem *MLMCMain::getStochasticProblem() {
...@@ -72,16 +80,16 @@ StochasticProblem *MLMCMain::getStochasticProblem() { ...@@ -72,16 +80,16 @@ StochasticProblem *MLMCMain::getStochasticProblem() {
return new StochasticLaplace1D(); return new StochasticLaplace1D();
if (problemName == "StochasticLaplace2D") if (problemName == "StochasticLaplace2D")
return new StochasticLaplace2D(); return new StochasticLaplace2D();
if (problemName == "StochasticPollution1D") if (problemName == "StochasticPollution1D"
|| problemName == "DeterministicPollution1D")
return new StochasticPollution1D(); return new StochasticPollution1D();
if (problemName == "StochasticPollution2D") if (problemName == "StochasticPollution2D"
|| problemName == "DeterministicPollution2D")
return new StochasticPollution2D(); return new StochasticPollution2D();
if (problemName == "StochasticPollution2DCosHat") if (problemName == "StochasticPollutionCosHat2D"
return new StochasticPollution2DCosHat(); || problemName == "DeterministicPollutionCosHat2D"
if (problemName == "DeterministicPollution2D") || problemName == "DeterministicPollutionCosHatSquare500")
return new DeterministicPollution2D(); return new StochasticPollutionCosHat2D();
if (problemName == "DeterministicPollution2DCosHat")
return new DeterministicPollution2DCosHat();
Exit("\nStochastic problem not implemented yet\n") Exit("\nStochastic problem not implemented yet\n")
} }
...@@ -117,6 +125,9 @@ MatrixGraphs *MLMCMain::getSolutionMatrixGraphs() { ...@@ -117,6 +125,9 @@ MatrixGraphs *MLMCMain::getSolutionMatrixGraphs() {
void MLMCMain::run() { void MLMCMain::run() {
printParameters(); printParameters();
assemble->PrintInfo(solutionMatrixGraphs->fine());
meshes->PrintInfo();
if (experimentName == "ConvergenceTest") { if (experimentName == "ConvergenceTest") {
headConvergenceTest(); headConvergenceTest();
mlmc->Method(); mlmc->Method();
......
...@@ -45,9 +45,9 @@ public: ...@@ -45,9 +45,9 @@ public:
string Name() const override { return "Stochastic Pollution 1D"; } string Name() const override { return "Stochastic Pollution 1D"; }
}; };
class Pollution2D : public StochasticTransportProblem { class StochasticPollution2D : public StochasticTransportProblem {
public: public:
Pollution2D() = default; StochasticPollution2D() = default;
double ut(double t, const Point &x) const override { double ut(double t, const Point &x) const override {
double d = 1.0 / 16.0; double d = 1.0 / 16.0;
...@@ -55,27 +55,6 @@ public: ...@@ -55,27 +55,6 @@ public:
return 1.0; return 1.0;
else return 0.0; else return 0.0;
} }
};
class Pollution2DCosHat : public StochasticTransportProblem {
private:
double aa = 1.00;
double cc = 5.0;
public:
Pollution2DCosHat() = default;
double ut(double t, const Point &x) const override {
Point midPoint = Point(0.5, 0.8);
double r = dist(midPoint, x);
if (r < 1 / cc)
return aa * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0;
}
};
class StochasticPollution2D : public Pollution2D {
public:
StochasticPollution2D() = default;
VectorField cellB(const cell &c, const Point &x) const override { VectorField cellB(const cell &c, const Point &x) const override {
VectorField localFlux; VectorField localFlux;
...@@ -92,9 +71,21 @@ public: ...@@ -92,9 +71,21 @@ public:
string Name() const override { return "Stochastic Pollution 2D"; } string Name() const override { return "Stochastic Pollution 2D"; }
}; };
class StochasticPollution2DCosHat : public Pollution2DCosHat { class StochasticPollutionCosHat2D : public StochasticTransportProblem {
private:
double aa = 1.00;
double cc = 5.0;
public: public:
StochasticPollution2DCosHat() = default; StochasticPollutionCosHat2D() = default;
double ut(double t, const Point &x) const override {
Point midPoint = Point(0.5, 0.8);
double r = dist(midPoint, x);
if (r < 1 / cc)
return aa * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0;
}
VectorField cellB(const cell &c, const Point &x) const override { VectorField cellB(const cell &c, const Point &x) const override {
VectorField localFlux; VectorField localFlux;
...@@ -111,36 +102,4 @@ public: ...@@ -111,36 +102,4 @@ public:
string Name() const override { return "Stochastic Pollution 2D"; } string Name() const override { return "Stochastic Pollution 2D"; }
}; };
class DeterministicPollution2D : public Pollution2D {
public:
DeterministicPollution2D() = default;
VectorField cellB(const cell &c, const Point &x) const override {
return VectorField(0.0, -1.0);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return N * VectorField(0.0, -1.0);
}
string Name() const override { return "Pollution 2D"; }
};
class DeterministicPollution2DCosHat : public Pollution2DCosHat {
public:
DeterministicPollution2DCosHat() = default;
VectorField cellB(const cell &c, const Point &x) const override {
return VectorField(0.0, -1.0);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return N * VectorField(0.0, -1.0);
}
string Name() const override { return "Pollution 2D"; }
};
#endif //MLMC_STOCHASTICTRANSPORTPROBLEM_HPP #endif //MLMC_STOCHASTICTRANSPORTPROBLEM_HPP
...@@ -25,7 +25,21 @@ public: ...@@ -25,7 +25,21 @@ public:
rnmg = new RNManager(streamnum, nstreams, gtype); rnmg = new RNManager(streamnum, nstreams, gtype);
} }
~CirculantEmbedding() { delete rnmg; }; virtual ~CirculantEmbedding() { delete rnmg; };
};
class DeterministicPermeabilityGenerator : public CirculantEmbedding {
public:
DeterministicPermeabilityGenerator(int l, Meshes &meshes)
: CirculantEmbedding(l, meshes) {}
void GetFineSample(Vector &fineField) override {
fineField = 1.0;
};
void GetCoarseSample(const Vector &fineField, Vector &coarseField) override {
coarseField = 1.0;
}
}; };
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D { class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
......
#include "HybridFluxGenerator.h" #include "HybridFluxGenerator.h"
void HybridFluxGenerator1D::GetFineSample(Vector &fineField) { void HybridFluxGenerator::GetFineSample(Vector &fineField) {
}
void HybridFluxGenerator1D::GetCoarseSample(const Vector &fineField,
Vector &coarseField) {
}
void HybridFluxGenerator2D::GetFineSample(Vector &fineField) {
MatrixGraphs cellMatrixGraph(meshes, dof("cell", 1)); MatrixGraphs cellMatrixGraph(meshes, dof("cell", 1));
MatrixGraphs faceMatrixGraph(meshes, dof("face", 1)); MatrixGraphs faceMatrixGraph(meshes, dof("face", 1));
Vector perm(cellMatrixGraph[l - meshes.pLevel()]); Vector perm(cellMatrixGraph[l - meshes.pLevel()]);
Vector solution(faceMatrixGraph[l - meshes.pLevel()]); Vector solution(faceMatrixGraph[l - meshes.pLevel()]);
permeabilityGenerator.GetFineSample(perm); permeabilityGenerator->GetFineSample(perm);
bool plotPermForTransport = false; bool plotPermForTransport = false;
config.get("plotPermForTransport", plotPermForTransport); config.get("plotPermForTransport", plotPermForTransport);
...@@ -23,20 +16,20 @@ void HybridFluxGenerator2D::GetFineSample(Vector &fineField) { ...@@ -23,20 +16,20 @@ void HybridFluxGenerator2D::GetFineSample(Vector &fineField) {
plot.vtk_celldata("permeability"); plot.vtk_celldata("permeability");
} }
assemble.problem->LoadNewSample(&perm); assemble->problem->LoadNewSample(&perm);
Preconditioner *pc = GetPC(faceMatrixGraph, assemble, "GaussSeidel"); Preconditioner *pc = GetPC(faceMatrixGraph, *assemble, "SuperLU");
Solver solver(pc, "GMRES"); Solver solver(pc, "GMRES");
Newton newton(solver); Newton newton(solver);
newton(assemble, solution); newton(*assemble, solution);
assemble.setFlux(solution, fineField); assemble->setFlux(solution, fineField);
} }
void HybridFluxGenerator2D::GetCoarseSample(const Vector &fineField, void HybridFluxGenerator::GetCoarseSample(const Vector &fineField,
Vector &coarseField) { Vector &coarseField) {
coarseField = 0; coarseField = 0;
for (cell c = coarseField.cells(); c != coarseField.cells_end(); c++) { for (cell c = coarseField.cells(); c != coarseField.cells_end(); c++) {
row coarseRow = coarseField.find_row(c()); row coarseRow = coarseField.find_row(c());
......
...@@ -7,50 +7,70 @@ ...@@ -7,50 +7,70 @@
class HybridFluxGenerator : public SampleGenerator { class HybridFluxGenerator : public SampleGenerator {
public: protected:
explicit HybridFluxGenerator(int l, Meshes &meshes) : SampleGenerator(l, meshes) {}
};
class HybridFluxGenerator1D : public HybridFluxGenerator {
private:
Meshes meshes;
Discretization *disc = nullptr; Discretization *disc = nullptr;
StochasticEllipticProblem *problem = nullptr; StochasticEllipticProblem *problem = nullptr;
CirculantEmbedding2D *permeabilityGenerator = nullptr;
HybridEllipticAssemble *assemble = nullptr; HybridEllipticAssemble *assemble = nullptr;
CirculantEmbedding *permeabilityGenerator = nullptr;
public: public:
explicit HybridFluxGenerator1D(int l, Meshes &meshes) : explicit HybridFluxGenerator(int l, Meshes &meshes)
HybridFluxGenerator(l, meshes) { : SampleGenerator(l, meshes) {}
disc = new Discretization("RT0_P0", 1, 2);
problem = new StochasticLaplace2D(); ~HybridFluxGenerator() {
assemble = new HybridEllipticAssemble(disc, problem); delete disc;
permeabilityGenerator = new CirculantEmbedding2D(l, meshes); delete problem;
} delete assemble;
delete permeabilityGenerator;
};
void GetFineSample(Vector &fineField) override; void GetFineSample(Vector &fineField) override;
void GetCoarseSample(const Vector &fineField, Vector &coarseField) override; void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
}; };
class HybridFluxGenerator2D : public HybridFluxGenerator { class HybridFluxGenerator1D : public HybridFluxGenerator {
private: public:
Discretization disc; explicit HybridFluxGenerator1D(int l, Meshes &meshes) :
StochasticLaplace2D problem; HybridFluxGenerator(l, meshes) {
HybridEllipticAssemble assemble; disc = new Discretization("RT0_P0", 1, meshes.dim());
CirculantEmbedding2D permeabilityGenerator; problem = new StochasticLaplace1D();
assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new CirculantEmbedding1D(l, meshes);
}
};
class HybridFluxGenerator2D : public HybridFluxGenerator {
public: public:
explicit HybridFluxGenerator2D(int l, Meshes &meshes) explicit HybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes), : HybridFluxGenerator(l, meshes) {
disc(Discretization("RT0_P0", 1, 2)), disc = new Discretization("RT0_P0", 1, meshes.dim());
problem(StochasticLaplace2D()), problem = new StochasticLaplace2D();
assemble(HybridEllipticAssemble(&disc, &problem)), assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator(CirculantEmbedding2D(l, meshes)) {} permeabilityGenerator = new CirculantEmbedding2D(l, meshes);
}
};
void GetFineSample(Vector &fineField) override; class DeterministicHybridFluxGenerator1D : public HybridFluxGenerator {
public:
explicit DeterministicHybridFluxGenerator1D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, meshes.dim());
problem = new StochasticLaplace1D();
assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new DeterministicPermeabilityGenerator(l, meshes);
}
};
void GetCoarseSample(const Vector &fineField, Vector &coarseField) override; class DeterministicHybridFluxGenerator2D : public HybridFluxGenerator {
public:
explicit DeterministicHybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, meshes.dim());
problem = new StochasticLaplace2D();
assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new DeterministicPermeabilityGenerator(l, meshes);
}
}; };
#endif //HYBRIDFLUXGENERATOR_H #endif //HYBRIDFLUXGENERATOR_H
...@@ -22,6 +22,12 @@ private: ...@@ -22,6 +22,12 @@ private:
if (meshes.dim() == 2) if (meshes.dim() == 2)
return new HybridFluxGenerator2D(l, meshes); return new HybridFluxGenerator2D(l, meshes);
} }
if (generatorName == "DeterministicHybridFluxGenerator") {
if (meshes.dim() == 1)
return new DeterministicHybridFluxGenerator1D(l, meshes);
if (meshes.dim() == 2)
return new DeterministicHybridFluxGenerator2D(l, meshes);
}
Exit("Generator does not exist") Exit("Generator does not exist")
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment