Commit c13a9eb4 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

runnable elliptic, compilable transport

parent 9c771af7
......@@ -2,22 +2,23 @@ Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = DGFEM
Model = DGTransport
#Overlap = dG1
Overlap = dG1
#sign = 1
#penalty = 20
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
#Problem = StochasticLaplace2D
#Problem = StochasticPollution1D
#Problem = StochasticPollution2D
Problem = StochasticPollution2D
StochasticField = LogNormal
#StochasticField = Gaussian
#Experiment = ConvergenceTest
#Experiment = MLMCExperiment
Experiment = MLMCOverEpsilon
Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
initLevels = 3, 4, 5
initSampleAmount = 12, 6, 3
......@@ -30,7 +31,7 @@ mcOnly = false
epsilon_lst = 0.05, 0.01, 0.005
maxLevel = 8
pLevel = 2
plevel = 2
degree = 1
functional = L2
......
......@@ -11,6 +11,6 @@ set(MLMC_SRC
assemble/DGEllipticAssemble.C
stochastics/RNManager.C
stochastics/CirculantEmbedding.C
problem/StochasticEllipticProblem.h stochastics/SampleGenerator.h stochastics/HybridFluxCalculator.C stochastics/HybridFluxCalculator.h)
problem/StochasticEllipticProblem.h stochastics/SampleGenerator.h stochastics/HybridFluxGenerator.C stochastics/HybridFluxGenerator.h)
add_library(MLMC STATIC ${MLMC_SRC})
......@@ -6,18 +6,15 @@ using namespace std;
void MLMCMain::initialize() {
meshes = getMeshes();
disc = getDiscretization();
// fieldSample = getFieldSampleVector();
// stochasticField = getStochasticField();
stochasticField = new StochasticField(meshes);
sampleMatrixGraphs = getSampleMatrixGraphs();
stochasticField = getStochasticField();
problem = getStochasticProblem();
assemble = getAssemble();
graphs = getMatrixGraphs();
solutionMatrixGraphs = getSolutionMatrixGraphs();
mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes,
stochasticField, assemble, graphs);
stochasticField, assemble, solutionMatrixGraphs);
}
Meshes *MLMCMain::getMeshes() {
......@@ -54,26 +51,18 @@ Discretization *MLMCMain::getDiscretization() {
Exit("\nDiscretization not implemented yet\n")
}
Vector *MLMCMain::getFieldSampleVector() {
if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D") {
MatrixGraphs cellMatrixGraphs(*meshes, dof("cell", 1));
return new Vector(cellMatrixGraphs.fine());
}
if (problemName == "StochasticPollution1D" || "StochasticPollution2D") {
MatrixGraphs cellMatrixGraphs(*meshes, dof("cell", 3));
return new Vector(cellMatrixGraphs.fine());
}
MatrixGraphs *MLMCMain::getSampleMatrixGraphs() {
if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D")
return new MatrixGraphs(*meshes, dof("cell", 1));
if (problemName == "StochasticPollution1D" || "StochasticPollution2D")
return new MatrixGraphs(*meshes, dof("cell", 3));
}
StochasticField *MLMCMain::getStochasticField() {
if (problemName == "StochasticLaplace1D")
return new StochasticField(meshes);
if (problemName == "StochasticLaplace1D")
return new StochasticField(meshes);
if (problemName == "StochasticPollution1D")
return new StochasticVectorField(meshes);
if (problemName == "StochasticPollution2D")
return new StochasticField(meshes);
if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D")
return new StochasticField(*meshes, "CirculantEmbedding");
if (problemName == "StochasticPollution1D" || "StochasticPollution2D")
return new StochasticField(*meshes, "HybridFluxGenerator");
}
StochasticProblem *MLMCMain::getStochasticProblem() {
......@@ -107,12 +96,12 @@ Assemble *MLMCMain::getAssemble() {
Exit("\nAssembling for this problem not implemented yet\n")
}
MatrixGraphs *MLMCMain::getMatrixGraphs() {
MatrixGraphs *MLMCMain::getSolutionMatrixGraphs() {
if (modelName == "LagrangeFEM" || modelName == "MixedFEM")
return new MatrixGraphs(*meshes, *disc);
if (modelName == "HybridFEM")
return new MatrixGraphs(*meshes, dof("face", 1));
if (modelName == "DGFEM")
if (modelName == "DGFEM" || "DGTransport")
return new CellMatrixGraphs(*meshes, *disc);
Exit("\nMatrix graph not implemented yet\n")
}
......@@ -120,7 +109,6 @@ MatrixGraphs *MLMCMain::getMatrixGraphs() {
void MLMCMain::run() {
printParameters();
//Todo ConvergenceTest bool mcOnly
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
mlmc->Method();
......
......@@ -26,14 +26,15 @@ public:
StochasticField *stochasticField = nullptr;
StochasticProblem *problem = nullptr;
Assemble *assemble = nullptr;
MatrixGraphs *graphs = nullptr;
MatrixGraphs *solutionMatrixGraphs = nullptr;
MatrixGraphs *sampleMatrixGraphs = nullptr;
Plot *plot = nullptr;
MultilevelMonteCarlo *mlmc = nullptr;
MLMCMain() {
config.get("enablePython", enablePython);
config.get("epsilon", epsilon);
config.get("pLevel", pLevel);
config.get("plevel", pLevel);
config.get("maxLevel", maxLevel);
config.get("degree", degree);
config.get("initLevels", initLevels);
......@@ -54,7 +55,7 @@ public:
~MLMCMain() {
delete meshes, delete disc, delete stochasticField;
delete problem, delete assemble, delete graphs;
delete problem, delete assemble, delete solutionMatrixGraphs;
delete mlmc;
}
......@@ -64,7 +65,7 @@ public:
Discretization *getDiscretization();
Vector *getFieldSampleVector();
MatrixGraphs *getSampleMatrixGraphs();
StochasticField *getStochasticField();
......@@ -72,7 +73,7 @@ public:
Assemble *getAssemble();
MatrixGraphs *getMatrixGraphs();
MatrixGraphs *getSolutionMatrixGraphs();
void run();
......
......@@ -374,7 +374,6 @@ void EllipticAssemble::plotU(Plot &plot, const Vector &u, const char *name) {
}
void EllipticAssemble::plotFlux(Plot &plot, const Vector &flux, const char *name) {
plot.celldata(flux, 3);
plot.vtk_cellvector(name);
}
......
......@@ -21,11 +21,11 @@ public:
int M = 0, dM = 0;
double sumCost = 0.0, avgCost = 0.0;
double sumY = 0.0, avgY = 0.0;
double sumY1 = 0.0, avgY = 0.0;
double sumY2 = 0.0, avgY2 = 0.0, varY = 0.0;
double sumY3 = 0.0, avgY3 = 0.0;
double sumY4 = 0.0, avgY4 = 0.0;
double sumQ = 0.0, avgQ = 0.0;
double sumQ1 = 0.0, avgQ = 0.0;
double sumQ2 = 0.0, avgQ2 = 0.0, varQ = 0.0;
double kurtosis = 0.0;
......@@ -51,23 +51,23 @@ public:
protected:
void updateSums(double fineCost, double diffQ, double fineQ) {
sumCost += fineCost;
sumY += diffQ;
sumQ += fineQ;
sumY2 += sumY * diffQ;
sumY3 += sumY2 * diffQ;
sumY4 += sumY3 * diffQ;
sumQ2 += sumQ * fineQ;
sumQ1 += fineQ;
sumQ2 += fineQ * fineQ;
sumY1 += diffQ;
sumY2 += diffQ * diffQ;
sumY3 += diffQ * diffQ * diffQ;
sumY4 += diffQ * diffQ * diffQ * diffQ;
}
void updateAvg() {
M += dM;
dM = 0;
avgCost = sumCost / M;
avgY = abs(sumY / M);
avgY = abs(sumY1 / M);
avgY2 = sumY2 / M;
avgY3 = sumY3 / M;
avgY4 = sumY4 / M;
avgQ = abs(sumQ / M);
avgQ = abs(sumQ1 / M);
avgQ2 = sumQ2 / M;
}
......
......@@ -27,7 +27,7 @@ void MonteCarloElliptic::Method() {
fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;
fineFlux = 0.0, coarseFlux = 0.0, coarseFluxUp = 0.0;
double fineQ, coarseQ, fineCost, coarseCost;;
double fineQ, coarseQ, fineCost, coarseCost;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
for (int m = M; m < M + dM; m++) {
......
......@@ -3,25 +3,48 @@
using namespace std;
void MonteCarloTransport::Method() {
Date Start;
vout(1) << " Start Monte Carlo method on level " << l;
if (!onlyFineLevel) vout(1) << " and level " << l - 1;
vout(1) << " with " << dM << " samples" << endl;
logger->startMethodMSG();
if (!onlyFineLevel)
logger->logMSGV1("l: " + to_string(l)
+ " dM: " + to_string(dM));
else
logger->logMSGV1("l: " + to_string(l)
+ " dM: " + to_string(dM)
+ " onlyFineLevel: True");
MatrixGraphs fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3)));
Vector fineSolution(solMatrixGraphs[l - pLevel]);
Vector fineFlux(fluxMatrixGraphs[l - pLevel]);
Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
Vector coarseFlux(fluxMatrixGraphs[l - pLevel - 1]);
fineFlux = 0.0;
double fineQ, coarseQ, fineCost, coarseCost;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
for (int m = M; m < M + dM; m++) {
stochasticField->GetFineSample(l, fineFlux);
assemble->problem->LoadNewSample(&fineFlux);
solveTransport(l, fineQ, fineCost, fineSolution);
if (!onlyFineLevel) {
stochasticField->GetCoarseSample(l - 1, fineFlux, coarseFlux);
assemble->problem->LoadNewSample(&coarseFlux);
solveTransport(l, fineQ, fineCost, fineSolution);
}
}
}
void MonteCarloTransport::solveTransport(double &Q, double &cost, Vector &u) {
void MonteCarloTransport::solveTransport(int l,
double &valueQ,
double &cost,
Vector &u) {
Solver solver;
Preconditioner *BJ = GetPC("PointBlockJacobi");
Solver IM(BJ);
DGTimeIntegrator timeIntegrator(solver, solver, solver, IM);
TimeSeries timeSeries;
// transportAssemble.PrintInfo(u);
// assemble.PrintInfo(u);
timeIntegrator(*assemble, timeSeries, u);
}
\ No newline at end of file
......@@ -8,10 +8,11 @@
class MonteCarloTransport : public MonteCarlo {
private:
void solveTransport(double &Q, double &cost, Vector &u);
void solveTransport(int l, double &valueQ, double &cost, Vector &solution);
public:
DGTransportAssemble *assemble;
TimeSeries timeSeries;
CellMatrixGraphs solMatrixGraphs;
......@@ -23,8 +24,8 @@ public:
assemble(assemble),
solMatrixGraphs(*meshes, *assemble->disc) {
transfer = new MatrixTransfer(*assemble);
transfer->Construct((solMatrixGraphs)[l - pLevel], (solMatrixGraphs)[l - pLevel - 1]);
// transfer = new MatrixTransfer(*assemble);
// transfer->Construct((solMatrixGraphs)[l - pLevel], (solMatrixGraphs)[l - pLevel - 1]);
}
~MonteCarloTransport() {
......
......@@ -57,8 +57,11 @@ public:
}
VectorField cellB(const cell &c, const Point &x) const override {
VectorField localFlux = 0.0;
for (int d = 0; d < 3; d++)
localFlux[d] = (*fieldSample)(c(), d);
return localFlux;
// return hybrid->EvaluateCellFlux(*flux, c);
return Point();
}
double faceB(const cell &c, int f,
......
......@@ -16,8 +16,7 @@ public:
string fieldType = "LogNormal";
double mean = 0.0;
explicit CirculantEmbedding(int l)
: SampleGenerator(l) {
explicit CirculantEmbedding(int l, Meshes &meshes) : SampleGenerator(l, meshes) {
config.get("evtol", evtol);
config.get("StochasticFields", fieldType);
config.get("mean", mean);
......@@ -39,7 +38,8 @@ public:
double domain_start = 0;
double domain_end = 1;
explicit CirculantEmbedding1D(int l) : CirculantEmbedding(l) {
explicit CirculantEmbedding1D(int l, Meshes &meshes) :
CirculantEmbedding(l, meshes) {
numberOfCells = (int) pow(2, l);
sqrtEigenvalues = getSqrtEV();
numberOfCellsEmbedded = sqrtEigenvalues.size();
......@@ -69,7 +69,8 @@ private:
class CirculantEmbedding2D : public CirculantEmbedding, public CovarianceFunction2D {
public:
explicit CirculantEmbedding2D(int l) : CirculantEmbedding(l) {
explicit CirculantEmbedding2D(int l, Meshes &meshes) :
CirculantEmbedding(l, meshes) {
numberOfCells[0] = (int) pow(2, l);
numberOfCells[1] = (int) pow(2, l);
sqrtEigenvalues = getSqrtEV();
......
#include "HybridFluxCalculator.h"
void HybridPollutionFluxCalculator2D::GetFineSample(Vector &fineField) {
MatrixGraphs faceMatrixGraph(meshes, dof("face", 1));
MatrixGraphs cellMatrixGraph(meshes, dof("cell", 3));
Vector solution(faceMatrixGraph[l + meshes.pLevel()]);
Vector flux(cellMatrixGraph[l + meshes.pLevel()]);
Preconditioner *pc = GetPC(faceMatrixGraph, *assemble, "LIB_PS");
Solver solver(pc, "GMRES");
Newton newton(solver);
newton(*assemble, solution);
assemble->setFlux(solution, flux);
}
void HybridPollutionFluxCalculator2D::GetCoarseSample(const Vector &fineField, Vector &coarseField) {
}
#include "HybridFluxGenerator.h"
void HybridFluxGenerator1D::GetFineSample(Vector &fineField) {
}
void HybridFluxGenerator1D::GetCoarseSample(const Vector &fineField,
Vector &coarseField) {
}
void HybridFluxGenerator2D::GetFineSample(Vector &fineField) {
MatrixGraphs cellMatrixGraph(meshes, dof("cell", 1));
MatrixGraphs faceMatrixGraph(meshes, dof("face", 1));
Vector perm(cellMatrixGraph[l - meshes.pLevel()]);
Vector solution(faceMatrixGraph[l - meshes.pLevel()]);
permeabilityGenerator.GetFineSample(perm);
assemble.problem->LoadNewSample(&perm);
Preconditioner *pc = GetPC(faceMatrixGraph, assemble, "LIB_PS");
Solver solver(pc, "GMRES");
Newton newton(solver);
newton(assemble, solution);
assemble.setFlux(solution, fineField);
}
void HybridFluxGenerator2D::GetCoarseSample(const Vector &fineField,
Vector &coarseField) {
coarseField = 0;
mout << OUTX(fineField) << endl;
mout << OUTX(coarseField) << endl;
for (cell c = coarseField.cells(); c != coarseField.cells_end(); c++) {
row coarseRow = coarseField.find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = fineField.find_row(c.Child(k));
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) += fineField(fineRow, i);
}
for (int i = 0; i < coarseRow.n(); ++i)
coarseField(coarseRow, i) /= c.Children();
}
mout << OUTX(fineField) << endl;
mout << OUTX(coarseField) << endl;
}
\ No newline at end of file
#ifndef HYBRIDFLUXCALCULATOR_H
#define HYBRIDFLUXCALCULATOR_H
#ifndef HYBRIDFLUXGENERATOR_H
#define HYBRIDFLUXGENERATOR_H
#include "SampleGenerator.h"
#include "CirculantEmbedding.h"
#include "assemble/HybridEllipticAssemble.h"
class HybridPollutionFluxCalculator : public SampleGenerator {
class HybridFluxGenerator : public SampleGenerator {
public:
explicit HybridPollutionFluxCalculator(int l) :
SampleGenerator(l) {}
explicit HybridFluxGenerator(int l, Meshes &meshes) : SampleGenerator(l, meshes) {}
};
class HybridPollutionFluxCalculator1D : public HybridPollutionFluxCalculator {
public:
explicit HybridPollutionFluxCalculator1D(int l) :
HybridPollutionFluxCalculator(l) {
}
};
class HybridFluxGenerator1D : public HybridFluxGenerator {
private:
Meshes meshes;
Discretization *disc = nullptr;
StochasticEllipticProblem *problem = nullptr;
CirculantEmbedding2D *permeabilityGenerator = nullptr;
HybridEllipticAssemble *assemble = nullptr;
class HybridPollutionFluxCalculator2D : public HybridPollutionFluxCalculator {
public:
explicit HybridPollutionFluxCalculator2D(int l) : HybridPollutionFluxCalculator(l) {
explicit HybridFluxGenerator1D(int l, Meshes &meshes) :
HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, 2);
problem = new StochasticLaplace2D();
assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new CirculantEmbedding2D(l);
}
~HybridPollutionFluxCalculator2D() {
delete permeabilityGenerator;
delete assemble;
permeabilityGenerator = new CirculantEmbedding2D(l, meshes);
}
void GetFineSample(Vector &fineField) override;
void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
};
class HybridFluxGenerator2D : public HybridFluxGenerator {
private:
Meshes meshes;
Discretization *disc = nullptr;
StochasticEllipticProblem *problem = nullptr;
CirculantEmbedding2D *permeabilityGenerator = nullptr;
HybridEllipticAssemble *assemble = nullptr;
Discretization disc;
StochasticLaplace2D problem;
HybridEllipticAssemble assemble;
CirculantEmbedding2D permeabilityGenerator;
public:
explicit HybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes),
disc(Discretization("RT0_P0", 1, 2)),
problem(StochasticLaplace2D()),
assemble(HybridEllipticAssemble(&disc, &problem)),
permeabilityGenerator(CirculantEmbedding2D(l, meshes)) {}
void GetFineSample(Vector &fineField) override;
void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
};
#endif //HYBRIDFLUXCALCULATOR_H
#endif //HYBRIDFLUXGENERATOR_H
......@@ -7,8 +7,9 @@
class SampleGenerator {
public:
int l;
Meshes &meshes;
explicit SampleGenerator(int l) : l(l) {}
explicit SampleGenerator(int l, Meshes &meshes) : l(l), meshes(meshes) {}
virtual void GetFineSample(Vector &fineField) = 0;
......
#ifndef STOCHASTICFIELD_H
#define STOCHASTICFIELD_H
#include <utility>
#include "CirculantEmbedding.h"
#include "HybridFluxCalculator.h"
#include "HybridFluxGenerator.h"
class StochasticField {
private:
virtual SampleGenerator *getGenerator(int l) {
if (dim == 1)
return new CirculantEmbedding1D(l);
if (dim == 2)
return new CirculantEmbedding2D(l);
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")
}
public:
string generatorName = "";
map<int, SampleGenerator *> generators;
int pLevel = 0, maxLevel = 0, dim = 0;
Meshes &meshes;
explicit StochasticField(Meshes *meshes) {
pLevel = meshes->pLevel();
maxLevel = meshes->Level();
dim = meshes->dim();
for (int l = pLevel; l <= maxLevel; l++)
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);
}
......@@ -35,18 +45,4 @@ public:
}
};
class StochasticVectorField : public StochasticField {
private:
SampleGenerator *getGenerator(int l) {
if (dim == 1)
return nullptr;
// return new HybridPollutionFluxCalculator1D(l);
if (dim == 2)
return new HybridPollutionFluxCalculator2D(l);
}
public:
explicit StochasticVectorField(Meshes *meshes) : StochasticField(meshes) {}
};
#endif //STOCHASTICFIELD_H
Subproject commit 0d40b523a9d9bcba941b8208b2e11e0d010f388f
Subproject commit 5b33e85e111c8bdaf852e7f491526f821eeec01f