Commit 018196fd authored by christian.wieners's avatar christian.wieners
Browse files

.

parent 9960fe59
......@@ -12,6 +12,7 @@ flux_alpha = 1 #Upwind: 1, Central: 0
#Problem = StochasticLaplace1D
#Problem = StochasticLaplace2D
#Problem = StochasticPollution1D
#Problem = Pollution2D
Problem = StochasticPollution2D
StochasticField = LogNormal
......@@ -21,8 +22,8 @@ StochasticField = LogNormal
Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
initLevels = 5, 6, 7
initSampleAmount = 12, 6, 3
initLevels = 4, 5,
initSampleAmount = 1, 12,
#l_init = 3, 4, 5
#M_init = 45, 15, 5
......@@ -31,7 +32,7 @@ epsilon = 0.01
mcOnly = false
epsilon_lst = 0.05, 0.01, 0.005
maxLevel = 8
maxLevel = 5
plevel = 2
degree = 2
......
......@@ -61,7 +61,7 @@ MatrixGraphs *MLMCMain::getSampleMatrixGraphs() {
StochasticField *MLMCMain::getStochasticField() {
if (problemName == "StochasticLaplace1D" || problemName == "StochasticLaplace2D")
return new StochasticField(*meshes, "CirculantEmbedding");
if (problemName == "StochasticPollution1D" || problemName == "StochasticPollution2D")
if (problemName == "StochasticPollution1D" || problemName == "StochasticPollution2D"|| problemName == "Pollution2D")
return new StochasticField(*meshes, "HybridFluxGenerator");
}
......@@ -74,6 +74,8 @@ StochasticProblem *MLMCMain::getStochasticProblem() {
return new StochasticPollution1D();
if (problemName == "StochasticPollution2D")
return new StochasticPollution2D();
if (problemName == "Pollution2D")
return new Pollution2D();
Exit("\nStochastic problem not implemented yet\n")
}
......
......@@ -216,7 +216,7 @@ public:
h = Scalar(0.0);
v[0] = u;
w[0] = M * u;
double beta = sqrt(::real(v[0] * w[0]));
double beta = sqrt(std::real(v[0] * w[0]));
v[0] *= 1 / beta;
w[0] *= 1 / beta;
s_old[0] = beta;
......@@ -232,7 +232,7 @@ public:
v[k + 1] -= h[j][k] * v[j];
}
w[k + 1] = M * v[k + 1];
h[k + 1][k] = sqrt(::real(v[k + 1] * w[k + 1]));
h[k + 1][k] = sqrt(std::real(v[k + 1] * w[k + 1]));
v[k + 1] *= (1.0 / h[k + 1][k]);
w[k + 1] *= (1.0 / h[k + 1][k]);
SmallMatrix S_k(k + 1, k + 1);
......@@ -281,7 +281,7 @@ public:
h = Scalar(0.0);
v[0] = u;
Mv = M * v[0];
double beta = sqrt(::real(v[0] * Mv));
double beta = sqrt(std::real(v[0] * Mv));
if (beta < Eps) {
du = 0;
return;
......@@ -297,7 +297,7 @@ public:
v[k + 1] -= h[j][k] * v[j];
}
Mv = M * v[k + 1];
h[k + 1][k] = sqrt(::real(v[k + 1] * Mv));
h[k + 1][k] = sqrt(std::real(v[k + 1] * Mv));
v[k + 1] *= (1.0 / h[k + 1][k]);
SmallMatrix H_k(k + 1, k + 1);
for (int j = 0; j <= k; ++j)
......@@ -338,7 +338,7 @@ public:
h = Scalar(0.0);
v[0] = u;
w[0] = M * u;
double beta = sqrt(::real(v[0] * w[0]));
double beta = sqrt(std::real(v[0] * w[0]));
v[0] *= 1 / beta;
w[0] *= 1 / beta;
s_old[0] = beta;
......@@ -354,7 +354,7 @@ public:
v[k + 1] -= h[j][k] * v[j];
}
w[k + 1] = M * v[k + 1];
h[k + 1][k] = sqrt(::real(v[k + 1] * w[k + 1]));
h[k + 1][k] = sqrt(std::real(v[k + 1] * w[k + 1]));
v[k + 1] *= (1.0 / h[k + 1][k]);
w[k + 1] *= (1.0 / h[k + 1][k]);
SmallMatrix S_k(k + 1, k + 1);
......@@ -492,7 +492,7 @@ public:
B2 = M;
B2 += Scalar(-dt / alpha) * A;
B3 = M;
B3 += Scalar(-dt / ::conj(alpha)) * A;
B3 += Scalar(-dt / std::conj(alpha)) * A;
mout << "Gauss (s=3)" << endl;
S(B, 1);
S2(B2, 1);
......
......@@ -65,9 +65,9 @@ void MonteCarloTransport::solveTransport(int l, double &valueQ, double &cost,
Preconditioner *BJ = GetPC("PointBlockJacobi");
Solver IM(BJ);
DGTimeIntegrator timeIntegrator(solver, solver, solver, IM);
// assemble.PrintInfo(solution);
assemble->PrintInfo(solution);
timeIntegrator(*assemble, timeSeries, solution);
logger->logMSGV2("Solving on level " + to_string(l) + " took " +
to_string((startSolving - Date()).t) + " seconds");
}
\ No newline at end of file
}
......@@ -57,7 +57,7 @@ public:
}
VectorField cellB(const cell &c, const Point &x) const override {
VectorField localFlux = 0.0;
VectorField localFlux;
for (int d = 0; d < 3; d++)
localFlux[d] = (*fieldSample)(c(), d);
return localFlux;
......@@ -67,12 +67,35 @@ public:
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
// return hybrid->EvaluateNormalFlux(*flux, c, f);
return 0.0;
return N * VectorField(0.0,-1.0);
}
string Name() const override { return "Stochastic Pollution 2D"; }
};
class Pollution2D : public StochasticTransportProblem {
public:
Pollution2D() = default;
double ut(double t, const Point &x) const override {
double d = 1.0 / 16.0;
if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
return 1.0;
else return 0.0;
}
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"; }
};
//HybridEllipticAssemble *hybrid;
//Vector *flux;
......
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