Commit 636438b2 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored MC classes

parent a876ba72
......@@ -12,14 +12,17 @@
class MonteCarlo {
public:
int plotting = 0;
MonteCarloLogger *logger;
Plot *coarsePlot = nullptr;
Plot *finePlot = nullptr;
MonteCarloLogger *logger = nullptr;
int l = 0, pLevel = 0;
int l = 0;
int pLevel = 0;
bool onlyFineLevel = true;
string functional = "L2";
Meshes *meshes;
StochasticField *stochasticField;
Meshes *meshes = nullptr;
StochasticField *stochasticField = nullptr;
int M = 0, dM = 0;
......@@ -34,7 +37,6 @@ public:
MonteCarlo(int l, int dM, bool onlyFineLevel,
Meshes *meshes, StochasticField *stochFields) :
l(l), dM(dM), onlyFineLevel(onlyFineLevel),
meshes(meshes), pLevel(meshes->pLevel()),
stochasticField(stochFields) {
......@@ -42,16 +44,45 @@ public:
config.get("MCPlotting", plotting);
config.get("functional", functional);
coarsePlot = new Plot((*meshes)[l - 1], meshes->dim());
finePlot = new Plot((*meshes)[l], meshes->dim());
logger = new MonteCarloLogger();
}
~MonteCarlo() {
delete coarsePlot;
delete finePlot;
delete logger;
}
virtual void Method() = 0;
void Method() {
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");
method();
updateAvg();
updateStatistic();
logger->LogMsgv1("|E[Y]|: " + to_string(avgY)
+ " V[Y]: " + to_string(varY));
logger->EndMethodMsg();
}
protected:
virtual void initialize() = 0;
virtual void method() = 0;
virtual void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
Vector &solution) = 0;
void updateSums(double fineCost, double diffQ, double fineQ) {
sumCost += fineCost;
sumQ1 += fineQ;
......@@ -81,18 +112,22 @@ protected:
3.0 * avgY * avgY * avgY * avgY) / varY / varY;
}
string mkSampleDirName(int m) {
return buildName(m, "sample");
string mkSampleDirName(int m, bool fineLevel) {
if (fineLevel)
return buildName(m, "sample");
else
return buildName(m, "sample", "coarse");
}
void mkSampleDir(int m) {
static void mkSampleDir(int m, const string &sampleDirName) {
string pathToData = "data/vtk/";
string sampleDirName = mkSampleDirName(m);
int status = mkdir(const_cast<char *> ((pathToData.append(sampleDirName)).c_str()),0777);
char *sampleDirNameCharArr =
const_cast<char *> ((pathToData.append(sampleDirName)).c_str());
int status = mkdir(sampleDirNameCharArr, 0777);
}
string buildName(int m, string name, const string &suffix = "") {
string returnName = std::move(name);
string buildName(int m, const string &name, const string &suffix = "") {
string returnName = name;
if (!suffix.empty())
returnName = returnName.append("_").append(suffix);
if (l != -1)
......
......@@ -3,75 +3,70 @@
using namespace std;
void MonteCarloElliptic::Method() {
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");
Vector finePermeability(permMatrixGraphs[l - pLevel]);
Vector coarsePermeability(permMatrixGraphs[l - pLevel - 1]);
Vector coarsePermeabilityUp(finePermeability);
void MonteCarloElliptic::initialize() {
}
void MonteCarloElliptic::method() {
Vector finePerm(cellMatrixGraphs[l - pLevel]);
Vector coarsePerm(cellMatrixGraphs[l - pLevel - 1]);
Vector coarsePermeabilityUp(finePerm);
Vector fineSolution(solMatrixGraphs[l - pLevel]);
Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
Vector coarseSolutionUp(fineSolution);
Vector fineFlux(fluxMatrixGraphs[l - pLevel]);
Vector coarseFlux(fluxMatrixGraphs[l - pLevel - 1]);
Vector coarseFluxUp(fineFlux);
finePermeability = 0.0, coarsePermeability = 0.0, coarsePermeabilityUp = 0.0;
finePerm = 0.0, coarsePerm = 0.0, coarsePermeabilityUp = 0.0;
fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;
fineFlux = 0.0, coarseFlux = 0.0, coarseFluxUp = 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, finePermeability);
assemble->problem->LoadNewSample(&finePermeability);
solveElliptic(l, fineQ, fineCost, fineSolution);
assemble->setFlux(fineSolution, fineFlux);
string sampleDirName = mkSampleDirName(m, true);
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->SetPlot(l, finePlot);
stochasticField->GetFineSample(l, finePerm);
if (plotting > 1)
plotResults(l, m, "", fineSolution, fineFlux, finePermeability);
assemble->sampleDir = sampleDirName;
assemble->plot = finePlot;
assemble->problem->LoadNewSample(&finePerm);
solvePDE(l, true, fineQ, fineCost, fineSolution);
if (plotting) assemble->PlotResults(fineSolution);
if (!onlyFineLevel) {
stochasticField->GetCoarseSample(l - 1, finePermeability, coarsePermeability);
assemble->problem->LoadNewSample(&coarsePermeability);
solveElliptic(l, coarseQ, coarseCost, coarseSolution);
assemble->setFlux(coarseSolution, coarseFlux);
UpscaleSolution(coarseSolution, coarseSolutionUp);
UpscalePermeability(coarsePermeability, coarsePermeabilityUp);
UpscaleFlux(coarseFlux, coarseFluxUp);
if (plotting > 1)
plotResults(l, m, "up", coarseSolutionUp,
coarseFluxUp, coarsePermeabilityUp);
sampleDirName = mkSampleDirName(m, false);
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->SetPlot(l, coarsePlot);
stochasticField->GetCoarseSample(l, finePerm, coarsePerm);
assemble->sampleDir = sampleDirName;
assemble->plot = coarsePlot;
assemble->problem->LoadNewSample(&coarsePerm);
solvePDE(l, false, coarseQ, coarseCost, coarseSolution);
if (plotting) assemble->PlotResults(coarseSolution);
} else {
coarseQ = 0.0, coarseCost = 0.0, coarseSolutionUp = 0.0;
}
updateSums(fineCost, fineQ - coarseQ, fineQ);
}
updateAvg();
updateStatistic();
logger->LogMsgv1("|E[Y]|: " + to_string(avgY)
+ " V[Y]: " + to_string(varY));
logger->EndMethodMsg();
}
void MonteCarloElliptic::solveElliptic(int l, double &valueQ, double &cost,
Vector &solution) {
Date startSolving;
Solver solver(GetPC(solMatrixGraphs, *assemble, "LIB_PS"), "GMRES");
void MonteCarloElliptic::solvePDE(int l, bool fineLevel, double &valueQ, double &cost,
Vector &solution) {
Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES");
Newton newton(solver);
newton(*assemble, solution);
cost = solution.pSize();
if (functional == "L2") valueQ = assemble->L2(solution);
......@@ -79,21 +74,4 @@ void MonteCarloElliptic::solveElliptic(int l, double &valueQ, double &cost,
if (functional == "Outflow")
valueQ = assemble->getInflowOutflow(solution).second;
if (functional == "Goal") valueQ = assemble->getGoal(solution);
logger->LogMsgv2("Solving on level " + to_string(l) + " took " +
to_string((startSolving - Date()).t) + "s");
}
void MonteCarloElliptic::plotResults(int l, int m, const string &suffix,
const Vector &u, const Vector &flux,
const Vector &perm) {
string name = buildName(m, "u", suffix);
assemble->plotU(plot, u, name);
name = buildName(m, "flux", suffix);
assemble->plotFlux(plot, flux, name);
name = buildName(m, "perm", suffix);
assemble->plotPerm(plot, perm, name);
}
}
\ No newline at end of file
......@@ -6,84 +6,32 @@
class MonteCarloElliptic : public MonteCarlo {
private:
Plot plot;
protected:
void initialize() override;
void solveElliptic(int l, double &valueQ, double &cost,
Vector &solution);
void method() override;
virtual void plotResults(int l, int m, const string &suffix,
const Vector &u, const Vector &flux,
const Vector &perm);
void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
Vector &solution) override;
public:
EllipticAssemble *assemble;
MatrixGraphs permMatrixGraphs;
MatrixGraphs fluxMatrixGraphs;
MatrixGraphs cellMatrixGraphs;
MatrixGraphs solMatrixGraphs;
Transfer *transfer;
MonteCarloElliptic(int l, int dM, bool onlyFineLevel, Meshes *meshes,
StochasticField *stochFields, EllipticAssemble *assemble) :
MonteCarlo(l, dM, onlyFineLevel, meshes, stochFields),
assemble(assemble),
permMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 1))),
fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)),
plot(Plot((*meshes)[l], meshes->dim(), 2)) {
cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)) {
transfer = new MatrixTransfer(*assemble);
transfer->Construct(solMatrixGraphs[l - pLevel], solMatrixGraphs[l - pLevel - 1]);
}
~MonteCarloElliptic() {
transfer->Destruct();
delete transfer;
}
void UpscaleSolution(const Vector &coarseSol, Vector &coarseSolUp) {
coarseSolUp = (*transfer) * coarseSol;
}
static void UpscalePermeability(const Vector &coarsePerm, Vector &coarsePermUp) {
coarsePermUp = 0;
for (cell c = coarsePerm.cells(); c != coarsePerm.cells_end(); c++) {
row coarseRow = coarsePerm.find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = coarsePermUp.find_row(c.Child(k));
for (int i = 0; i < fineRow.n(); ++i)
coarsePermUp(fineRow, i) += coarsePerm(coarseRow, i);
}
}
}
static void UpscaleFlux(const Vector &coarseFlux, Vector &coarseFluxUp) {
coarseFluxUp = 0;
for (cell c = coarseFlux.cells(); c != coarseFlux.cells_end(); c++) {
row coarseRow = coarseFlux.find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = coarseFluxUp.find_row(c.Child(k));
for (int i = 0; i < fineRow.n(); ++i)
coarseFluxUp(fineRow, i) += coarseFlux(coarseRow, i);
}
}
}
// void updateUSum(const Vector &uc_up, const Vector &uf) {
// u_sum += uf;
// u_sum_diff += uf;
// u_sum_diff -= uc_up;
// }
//
// void updateUAvg() {
// Scalar M_inverted = 1.0 / M;
// u_avg = M_inverted * u_sum;
// u_avg_diff = M_inverted * u_sum_diff;
// }
void Method() override;
};
#endif //MLMC_MONTECARLOELLIPTIC_H
......@@ -2,86 +2,67 @@
using namespace std;
void MonteCarloTransport::Method() {
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");
void MonteCarloTransport::initialize() {
}
Vector fineFlux(faceMatrixGraphs[l - pLevel]);
Vector coarseFlux(faceMatrixGraphs[l - pLevel - 1]);
Vector coarseFluxUp(fineFlux);
void MonteCarloTransport::method() {
Vector fineNormalFlux(faceMatrixGraphs[l - pLevel]);
Vector coarseNormalFlux(faceMatrixGraphs[l - pLevel - 1]);
Vector fineSolution(solMatrixGraphs[l - pLevel]);
Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
Vector coarseSolutionUp(fineSolution);
fineFlux = 0.0, coarseFlux = 0.0, coarseFluxUp = 0.0;
fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;
fineNormalFlux = 0.0, coarseNormalFlux = 0.0;
fineSolution = 0.0, coarseSolution = 0.0;
double fineQ, coarseQ, fineCost, coarseCost;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
assemble->plot = &plot;
for (int m = M; m < M + dM; m++) {
string sampleDirName = mkSampleDirName(m);
mkSampleDir(m);
string sampleDirName = mkSampleDirName(m, true);
mkSampleDir(m, sampleDirName);
stochasticField->generators[l]->SetSampleDir(sampleDirName);
stochasticField->GetFineSample(l, fineFlux);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->SetPlot(l, finePlot);
stochasticField->GetFineSample(l, fineNormalFlux);
assemble->sampleDir = sampleDirName;
assemble->problem->LoadNewSample(&fineFlux);
assemble->plot = finePlot;
assemble->problem->LoadNewSample(&fineNormalFlux);
// if (plotting > 1) {
// string name = buildName(m, mkSampleDirName(m).append("/flux"));
// plot.celldata(fineFlux, 3);
// plot.vtk_cellvector(name);
// }
solveTransport(l, true, fineQ, fineCost, fineSolution);
solvePDE(l, true, fineQ, fineCost, fineSolution);
if (!onlyFineLevel) {
stochasticField->GetCoarseSample(l - 1, fineFlux, coarseFlux);
assemble->problem->LoadNewSample(&coarseFlux);
solveTransport(l, false, coarseQ, coarseCost, coarseSolution);
UpscaleSolution(coarseSolution, coarseSolutionUp);
UpscaleFlux(coarseFlux, coarseFluxUp);
if (plotting > 1) {
string name = buildName(m, "flux_up");
plot.celldata(coarseFluxUp, 3);
plot.vtk_cellvector(name);
}
sampleDirName = mkSampleDirName(m, false);
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->SetPlot(l, coarsePlot);
stochasticField->GetCoarseSample(l, fineNormalFlux, coarseNormalFlux);
assemble->sampleDir = sampleDirName;
assemble->plot = coarsePlot;
assemble->problem->LoadNewSample(&coarseNormalFlux);
solvePDE(l, false, coarseQ, coarseCost, coarseSolution);
} else {
coarseQ = 0.0, coarseCost = 0.0, coarseSolutionUp = 0.0;
coarseQ = 0.0, coarseCost = 0.0;
}
updateSums(fineCost, fineQ - coarseQ, fineQ);
}
updateAvg();
updateStatistic();
logger->LogMsgv1("|E[Y]|: " + to_string(avgY)
+ " V[Y]: " + to_string(varY));
logger->EndMethodMsg();
}
void MonteCarloTransport::solveTransport(int l,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) {
void MonteCarloTransport::solvePDE(int l,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) {
assemble->logger->StartMethodMsg();
Solver s; // TODO RMV
Preconditioner *pc = GetPC("PointBlockGaussSeidel"); // TODO ADD MG AND ASSEMBLE
Preconditioner *pc = GetPC("PointBlockGaussSeidel"); //TODO RMV after usage
Solver solver(pc, "GMRES");
DGTimeIntegrator timeIntegrator(s, s, s, solver);
Solver solver2(pc, "GMRES"); // TODO RMV
DGTimeIntegrator timeIntegrator(solver, solver, solver, solver2);
TimeSeries timeSeries = getTimeSeries(fineLevel);
assemble->resetTAssemble();
......
......@@ -8,67 +8,42 @@
class MonteCarloTransport : public MonteCarlo {
private:
Plot plot;
protected:
void initialize() override;
void method() override;
void solveTransport(int l,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution);
void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
Vector &solution) override;
private:
double startTime = 0.0;
double endTime = 0.0;
TimeSeries getTimeSeries(bool fineLevel);
public:
DGTransportAssemble *assemble;
double startTime = 0.0;
double endTime = 0.0;
MatrixGraphs fluxMatrixGraphs;
MatrixGraphs cellMatrixGraphs;
MatrixGraphs faceMatrixGraphs;
CellMatrixGraphs solMatrixGraphs;
Transfer *transfer;
MonteCarloTransport(int l, int dM, bool baseLevel, Meshes *meshes,
StochasticField *stochFields, DGTransportAssemble *assemble) :
MonteCarlo(l, dM, baseLevel, meshes, stochFields),
assemble(assemble),
fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
faceMatrixGraphs(MatrixGraphs(*meshes, dof("face", 3))),
solMatrixGraphs(*meshes, *assemble->disc),
plot(Plot((*meshes)[l], meshes->dim())) {
solMatrixGraphs(*meshes, *assemble->disc) {
config.get("startTime", startTime);
config.get("endTime", endTime);
transfer = new MatrixTransfer(*assemble);
if (meshes->dim() != 1)
transfer->Construct((solMatrixGraphs)[l - pLevel], (solMatrixGraphs)[l - pLevel - 1]);
}
~MonteCarloTransport() {
transfer->Destruct();
delete transfer;
}
void UpscaleSolution(const Vector &coarseSol, Vector &coarseSolUp) {
coarseSolUp = (*transfer) * coarseSol;
}
static void UpscaleFlux(const Vector &coarseFlux, Vector &coarseFluxUp) {
coarseFluxUp = 0;
for (cell c = coarseFlux.cells(); c != coarseFlux.cells_end(); c++) {
row coarseRow = coarseFlux.find_row(c());
for (int k = 0; k < c.Children(); ++k) {
row fineRow = coarseFluxUp.find_row(c.Child(k));
for (int i = 0; i < fineRow.n(); ++i)
coarseFluxUp(fineRow, i) += coarseFlux(coarseRow, i);
}
}
}
void Method() override;
};
#endif //MLMC_MONTECARLOTRANSPORT_HPP
Supports Markdown
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