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

added variables for flux and permeability

parent ba749830
Todo Circulant Embedding & Stochastic fields:
- Use Information from Mesh to be able to initialize the right domain for
circulant embedding -> CoarseGeometry knows coordinates
- Remove hack in Stochastic fields for problem size
- Add permeability average
Todo Plotting routines:
- Averaged Flux & Averaged Permeability -> Upscaling has to work for vector
valued data
- Averaged u is not implemented due to refactoring
\ No newline at end of file
......@@ -16,9 +16,9 @@ Experiment = ConvergenceTest
#Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
size_init = 7
l_init = 3 4 5 6 7 8 9
M_init = 100 100 100 100 100 100 100
size_init = 3
l_init = 3 4 5
M_init = 5 5 5
Lmax = 9
epsilon = 0.01
......@@ -33,7 +33,7 @@ functional = L2
MCVerbose = 1
MLMCVerbose = 2
MCPlotting = 0
MCPlotting = 3
MLMCPlotting = 0
enablePython = 1
\ No newline at end of file
......@@ -84,7 +84,7 @@ void MultilevelMonteCarlo::method(const double eps) {
u = 0.0;
wrapUpResults(u, value, cost, levels, numsamples);
if (plotting >= 1) {
map_mc.rbegin()->second->plotSolution(u, -1, "");
}
if (plotting >= 2) {
......@@ -178,7 +178,7 @@ void MultilevelMonteCarlo::wrapUpResults(Vector &u, double &val, double &mlmc_co
for (int l = mc.first; l < L; l++) {
if (l == mc.first)
*map_u_temp[l] = mc.second->u_avg_diff;
map_mc[l + 1]->upscaleU(*map_u_temp[l], *map_u_temp[l + 1]);
map_mc[l + 1]->upscale(*map_u_temp[l], *map_u_temp[l + 1]);
if (l == L - 1)
u += *map_u_temp[l + 1];
}
......
......@@ -220,7 +220,7 @@ void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
}
}
void DGEllipticAssemble::setExactSolution(Vector &u) const {
void DGEllipticAssemble::setExactU(Vector &u) const {
u = 0;
Point z[MaxNodalPoints];
for (cell c = u.cells(); c != u.cells_end(); ++c) {
......@@ -243,8 +243,7 @@ void DGEllipticAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM(c(), c.Child(k))[0] = 1;
}
void DGEllipticAssemble::plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
void DGEllipticAssemble::plotU(Plot &plot, const Vector &u, const char *name) {
Vector plot_u(u);
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
......
......@@ -51,12 +51,11 @@ public:
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
void setExactSolution(Vector &u) const override;
void setExactU(Vector &u) const override;
void AssembleTransfer(TransferMatrix &TM) const override;
void plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) override;
void plotU(Plot &plot, const Vector &u, const char *name) override;
};
#endif
......@@ -343,12 +343,23 @@ double EllipticAssemble::getGoal(const Vector &u) const {
return PPM->Sum(goal);
}
void EllipticAssemble::setExactSolution(Vector &u) const {
void EllipticAssemble::setPerm(Vector &perm) const {
for (cell c = perm.cells(); c != perm.cells_end(); ++c) {
Tensor K = problem->Permeability(c.Subdomain(), c());
perm(c(), 0) = K[0][0];
}
}
void EllipticAssemble::setU(Vector &u) const {
// Nothing to do here
}
void EllipticAssemble::setExactU(Vector &u) const {
for (row r = u.rows(); r != u.rows_end(); ++r)
u(r, 0) = problem->Solution(r());
}
void EllipticAssemble::setFlux(const Vector &u, Vector &flux) {
void EllipticAssemble::setFlux(const Vector &u, Vector &flux) const {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
auto *elem = getElement(c, u);
VectorField F = zero;
......@@ -366,32 +377,18 @@ void EllipticAssemble::setFlux(const Vector &u, Vector &flux) {
}
}
void EllipticAssemble::plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
void EllipticAssemble::plotU(Plot &plot, const Vector &u, const char *name) {
plot.vertexdata(u);
plot.vtk_vertexdata(name);
}
void EllipticAssemble::plotFlux(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
MatrixGraphs graphs(meshes, dof("cell", 3));
if (l == -1) l = meshes.Level();
Vector flux(graphs[l - meshes.pLevel()]);
setFlux(u, flux);
void EllipticAssemble::plotFlux(Plot &plot, const Vector &flux, const char *name) {
plot.celldata(flux, 3);
plot.vtk_celldata(name);
plot.vtk_cellvector(name);
}
void EllipticAssemble::plotPermeability(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
MatrixGraphs graphs(meshes, dof("cell", 3));
if (l == -1) l = meshes.Level();
Vector permeability(graphs[l - meshes.pLevel()]);
for (cell c = permeability.cells(); c != permeability.cells_end(); ++c) {
Tensor K = problem->Permeability(c.Subdomain(), c());
permeability(c(), 0) = K[0][0];
}
plot.celldata(permeability);
void EllipticAssemble::plotPerm(Plot &plot, const Vector &perm, const char *name) {
plot.celldata(perm);
plot.vtk_celldata(name);
}
......@@ -435,21 +432,22 @@ void EllipticAssemble::printSolutionInfo(const Vector &u) {
}
}
void EllipticAssemble::plotSolution(Plot &plot, const Vector &u, Meshes &meshes,
int l = -1, int m = -1, const string &suffix = "") {
void EllipticAssemble::plotSolution(Plot &plot, int l, int m, const string &suffix,
const Vector &u, const Vector &flux,
const Vector &perm) {
const char *name = buildString("u", l, m, suffix);
plotU(plot, u, meshes, name, l);
plotU(plot, u, name);
name = buildString("flux", l, m, suffix);
plotFlux(plot, u, meshes, name, l);
plotFlux(plot, flux, name);
name = buildString("perm", l, m, suffix);
plotPermeability(plot, u, meshes, name, l);
plotPerm(plot, perm, name);
}
const char *
EllipticAssemble::buildString(string name, int l, int m, const string &suffix) {
const char *EllipticAssemble::buildString(string name, int l, int m,
const string &suffix) {
if (!suffix.empty())
name += "_" + suffix;
if (l != -1)
......
......@@ -65,18 +65,19 @@ public:
virtual double getGoal(const Vector &u) const;
virtual void setExactSolution(Vector &u) const;
virtual void setPerm(Vector &perm) const;
virtual void setFlux(const Vector &u, Vector &flux);
virtual void setU(Vector &u) const;
virtual void plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l);
virtual void setExactU(Vector &u) const;
virtual void plotFlux(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l);
virtual void setFlux(const Vector &u, Vector &flux) const;
virtual void plotPermeability(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l);
virtual void plotU(Plot &plot, const Vector &u, const char *name);
virtual void plotFlux(Plot &plot, const Vector &flux, const char *name);
virtual void plotPerm(Plot &plot, const Vector &perm, const char *name);
void AssembleTransfer(TransferMatrix &TM) const override;
......@@ -84,8 +85,9 @@ public:
static const char *buildString(string name, int l, int m, const string &suffix);
virtual void plotSolution(Plot &plot, const Vector &u, Meshes &meshes,
int l, int m, const string &suffix);
virtual void plotSolution(Plot &plot, int l, int m, const string &suffix,
const Vector &u, const Vector &flux,
const Vector &perm);
};
#endif
......@@ -263,7 +263,29 @@ pair<double, double> HybridEllipticAssemble::getOutflowLeftRight(const Vector &p
return {PPM->Sum(outflowLeft), PPM->Sum(outflowRight)};
}
void HybridEllipticAssemble::setFlux(const Vector &u, Vector &flux) {
void HybridEllipticAssemble::setPerm(Vector &perm) const {
}
void HybridEllipticAssemble::setU(Vector &u) const {
// dof dofCell("cell", 3);
// MatrixGraphs GDcell(meshes, dofCell);
// Vector p(GDcell[l - meshes.pLevel()]);
// for (cell c = u.cells(); c != u.cells_end(); ++c) {
// double U = 0;
// for (int k = 0; k < c.Faces(); ++k)
// U += u(c.Face(k), 0);
// p(c(), 0) = U / c.Faces();
// }
}
void HybridEllipticAssemble::setExactU(Vector &u) const {
const Mesh &M = u.GetMesh();
for (face f = M.faces(); f != M.faces_end(); ++f)
u(f(), 0) = problem->Solution(f());
}
void HybridEllipticAssemble::setFlux(const Vector &u, Vector &flux) const {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
int N = c.Faces();
SmallMatrix A(N);
......@@ -301,26 +323,9 @@ void HybridEllipticAssemble::setFlux(const Vector &u, Vector &flux) {
}
}
void HybridEllipticAssemble::setExactSolution(Vector &u) const {
const Mesh &M = u.GetMesh();
for (face f = M.faces(); f != M.faces_end(); ++f)
u(f(), 0) = problem->Solution(f());
}
void HybridEllipticAssemble::plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
dof dofCell("cell", 3);
MatrixGraphs GDcell(meshes, dofCell);
Vector p(GDcell[l - meshes.pLevel()]);
for (cell c = u.cells(); c != u.cells_end(); ++c) {
double U = 0;
for (int k = 0; k < c.Faces(); ++k)
U += u(c.Face(k), 0);
p(c(), 0) = U / c.Faces();
}
plot.celldata(p);
void HybridEllipticAssemble::plotU(Plot &plot, const Vector &u, const char *name) {
plot.celldata(u);
plot.vtk_celldata("u");
plot.gnu_celldata("u");
}
void HybridEllipticAssemble::SetNormalFlux(const Vector &u, Vector &flux) {
......
......@@ -35,12 +35,15 @@ public:
pair<double, double> getOutflowLeftRight(const Vector &p) const override;
void setFlux(const Vector &u, Vector &flux) override;
void setPerm(Vector &perm) const override;
void setExactSolution(Vector &u) const override;
void setU(Vector &u) const override;
void plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) override;
void setExactU(Vector &u) const override;
void setFlux(const Vector &u, Vector &flux) const override;
void plotU(Plot &plot, const Vector &u, const char *name) override;
void SetNormalFlux(const Vector &u, Vector &flux);
......
......@@ -145,7 +145,7 @@ double MixedEllipticAssemble::DualPrimalError(const Discretization &discretizati
return sqrt(PPM->Sum(err));
}
void MixedEllipticAssemble::setExactSolution(Vector &u) const {
void MixedEllipticAssemble::setExactU(Vector &u) const {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
u(c(), 0) = problem->Solution(c());
for (int i = 0; i < c.Faces(); ++i) {
......@@ -158,8 +158,7 @@ void MixedEllipticAssemble::setExactSolution(Vector &u) const {
}
}
void MixedEllipticAssemble::plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) {
void MixedEllipticAssemble::plotU(Plot &plot, const Vector &u, const char *name) {
plot.celldata(u);
plot.vtk_celldata("u");
}
......
......@@ -38,10 +38,9 @@ public:
virtual double DualPrimalError(const Discretization &discretization,
const Vector &u_ref, const Vector &u) const;
void setExactSolution(Vector &u) const override;
void setExactU(Vector &u) const override;
void plotU(Plot &plot, const Vector &u, Meshes &meshes,
const char *name, int l) override;
void plotU(Plot &plot, const Vector &u, const char *name) override;
void AssembleTransfer(TransferMatrix &TM) const override;
};
......
......@@ -80,17 +80,10 @@ public:
3.0 * avg_Y * avg_Y * avg_Y * avg_Y) / var_Y / var_Y;
}
virtual void upscaleU(const Vector &uc, Vector &uc_up) = 0;
virtual void updateUSum(const Vector &uc_up, const Vector &uf) = 0;
virtual void updateUAvg() = 0;
virtual void plotSolution(Vector &u, int m, string suffix) = 0;
virtual void upscale(const Vector &uc, Vector &uc_up) = 0;
virtual void method() = 0;
virtual void solvePDE(double &Q, double &cost, Vector &u) {};
};
#endif //MLMC_MC_HPP
\ No newline at end of file
......@@ -31,16 +31,23 @@ void MonteCarloElliptic::method() {
Vector uf(u_sum), uc_up(u_sum), uc((*graphs)[l - plevel - 1]);
uf = 0.0, uc = 0.0, uc_up = 0.0;
Vector fluxf(flux_sum), fluxc_up(flux_sum);
fluxf = 0.0, fluxc_up = 0.0;
Vector permf(flux_sum), permc_up(flux_sum);
permf = 0.0, permc_up = 0.0;
Date Start_solving;
vout(2) << " Solving PDE on level " << l;
solvePDE(Qf, costf, uf);
solveElliptic(Qf, costf, uf);
vout(2) << " took " << Start_solving - Date() << endl;
if (!baseLevel) {
Date Start_solving_c;
vout(2) << " Solving PDE on level " << l - 1;
stochFields->setFineField(false);
solvePDE(Qc, costc, uc);
solveElliptic(Qc, costc, uc);
stochFields->setFineField(true);
vout(2) << " took " << Start_solving_c - Date() << endl;
} else {
......@@ -52,15 +59,21 @@ void MonteCarloElliptic::method() {
double sums[7] = {costf, Qf - Qc, pow((Qf - Qc), 2),
pow((Qf - Qc), 3), pow((Qf - Qc), 4),
Qf, pow(Qf, 2)};
updateSum(sums);
upscaleU(uc, uc_up);
upscale(uc, uc_up);
updateUSum(uc_up, uf);
if (plotting >= 2)
plotSolution(uf, m, "");
if (plotting >= 2) {
assemble->setFlux(uf, fluxf);
assemble->setPerm(permf);
plotSolution(uf, fluxf, permf, m, "");
}
if (plotting >= 3 && !baseLevel) {
stochFields->setFineField(false);
plotSolution(uc_up, m, "up");
assemble->setFlux(uc_up, fluxc_up);
assemble->setPerm(permc_up);
plotSolution(uc_up, fluxc_up, permc_up, m, "up");
stochFields->setFineField(true);
}
}
......@@ -70,15 +83,17 @@ void MonteCarloElliptic::method() {
updateStatistic();
if (plotting >= 1) {
plotSolution(u_avg_diff, -1, "avg_diff");
plotSolution(u_avg, -1, "avg");
plotSolution(u_avg_diff, flux_avg_diff, perm_avg_diff,
-1, "avg_diff");
plotSolution(u_avg, flux_avg, perm_avg,
-1, "avg");
}
vout (1) << " End after " << Date() - Start << ": |E(Y)|=" << avg_Y
<< " V(Y)=" << var_Y << endl;
}
void MonteCarloElliptic::solvePDE(double &Q, double &cost, Vector &u) {
void MonteCarloElliptic::solveElliptic(double &Q, double &cost, Vector &u) {
Preconditioner *preconditioner = GetPC(*graphs, *assemble, "LIB_PS");
Solver solver(preconditioner, "GMRES");
Newton newton(solver);
......@@ -88,4 +103,6 @@ void MonteCarloElliptic::solvePDE(double &Q, double &cost, Vector &u) {
if (functional == "L2") Q = assemble->L2(u);
if (functional == "Energy") Q = assemble->Energy(u);
if (functional == "Outflow") Q = assemble->getInflowOutflow(u).second;
}
......@@ -9,15 +9,27 @@ public:
EllipticAssemble *assemble;
Transfer *transfer;
MatrixGraphs cellGraphs;
Vector flux_sum, flux_avg, flux_sum_diff, flux_avg_diff;
Vector perm_sum, perm_avg, perm_sum_diff, perm_avg_diff;
MonteCarloElliptic(int l, int dM, bool baseLevel, Meshes *meshes,
StochasticFields *stochFields, MatrixGraphs *graphs,
EllipticAssemble *assemble) :
assemble(assemble),
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs) {
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs),
cellGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
flux_sum(cellGraphs[l - plevel]), flux_avg(flux_sum),
flux_sum_diff(flux_sum), flux_avg_diff(flux_sum),
perm_sum(flux_sum), perm_avg(flux_sum),
perm_sum_diff(flux_sum), perm_avg_diff(flux_sum),
assemble(assemble) {
transfer = new MatrixTransfer(*assemble);
transfer->Construct((*graphs)[l - plevel], (*graphs)[l - plevel - 1]);
flux_sum = 0.0, flux_avg = 0.0, flux_sum_diff = 0.0, flux_avg_diff = 0.0;
perm_sum = 0.0, perm_avg = 0.0, perm_sum_diff = 0.0, perm_avg_diff = 0.0;
}
~MonteCarloElliptic() override {
......@@ -25,30 +37,30 @@ public:
delete transfer;
}
void upscaleU(const Vector &uc, Vector &uc_up) override {
void upscale(const Vector &uc, Vector &uc_up) override {
uc_up = (*transfer) * uc;
}
void updateUSum(const Vector &uc_up, const Vector &uf) override {
void updateUSum(const Vector &uc_up, const Vector &uf) {
u_sum += uf;
u_sum_diff += uf;
u_sum_diff -= uc_up;
}
void updateUAvg() override {
void updateUAvg() {
Scalar M_inverted = 1.0 / M;
u_avg = M_inverted * u_sum;
u_avg_diff = M_inverted * u_sum_diff;
}
void plotSolution(Vector &u, int m, string suffix) override {
assemble->plotSolution(plot, u, *meshes, l, m, suffix);
void plotSolution(Vector &u, Vector &flux, Vector &perm,
int m, const string &suffix) {
assemble->plotSolution(plot, l, m, suffix, u, flux, perm);
}
void method() override;
void solvePDE(double &Q, double &cost, Vector &u) override;
void solveElliptic(double &Q, double &cost, Vector &u);
};
#endif //MLMC_MONTECARLOELLIPTIC_H
......@@ -6,6 +6,5 @@ void MonteCarloTransport::method() {
}
void MonteCarloTransport::solvePDE(double &Q, double &cost, Vector &u) {
MonteCarlo::solvePDE(Q, cost, u);
void MonteCarloTransport::solveTransport(double &Q, double &cost, Vector &u) {
}
......@@ -25,21 +25,18 @@ public:
delete transfer;
}
void upscaleU(const Vector &uc, Vector &uc_up) override {
void upscale(const Vector &uc, Vector &uc_up) override {
}
void updateUSum(const Vector &uc_up, const Vector &uf) override {
void updateUSum(const Vector &uc_up, const Vector &uf) {
}
void updateUAvg() override {
}
void plotSolution(Vector &u, int m, string suffix) override {
void updateUAvg() {
}
void method() override;
void solvePDE(double &Q, double &cost, Vector &u) override;
void solveTransport(double &Q, double &cost, Vector &u);
};
......
......@@ -5,18 +5,6 @@
#include "CovarianceFunction.h"
#include "Utils.h"
using namespace std;
/*
* Todo:
* - Use Information from Mesh to be able to initialize the right domain for
* circulant embedding -> CoarseGeometry knows coordinates
* - Remove hack in Stochastic fields for problem size
* - Add permeability average
*
*
*/
class CirculantEmbedding {
/* Circulant Embedding
......
#include "StochasticFields.h"
using namespace std;
class StochasticFieldPair1D : public StochasticFieldPair {
int N;
CirculantEmbedding1D *circEmbedding1D;
......@@ -70,8 +72,8 @@ class StochasticFieldPair2D : public StochasticFieldPair {
int N1;
int N2;
CirculantEmbedding2D *circEmbedding2D;
std::vector<std::vector<double>> fieldf;
std::vector<std::vector<double>> fieldc;
vector<vector<double>> fieldf;
vector<vector<double>> fieldc;
public:
explicit StochasticFieldPair2D(int l) : StochasticFieldPair() {
N1 = (int) pow(2, l); // Hack
......
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