Skip to content

GitLab

Commit 9711af40 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

cleaning up transport problems and preparing them for tests

parent 4d714db3
Hide whitespace changes
Inline Side-by-side
mlmc/src/problems/StochasticTransportProblem.cpp
View file @ 9711af40
......@@ -8,22 +8,25 @@ CreateStochasticTransportProblem(std::string problemName, Meshes &meshes) {
if (problemName == "Pollution1D")
return new Pollution1D(meshes);
if (problemName == "CircleWave2D")
return new CircleWave2D(meshes);
if (problemName == "StochasticPollutionCosHat1D")
return new StochasticPollutionCosHat1D(meshes);
// if (problemName == "PollutionCosHat1D")
// return new PollutionCosHat1D(meshes);
if (problemName == "StochasticPollution2D")
return new StochasticPollution2D(meshes);
if (problemName == "StochasticPollution2D")
return new Pollution2D(meshes);
if (problemName == "StochasticPollutionCosHat2D")
return new StochasticPollutionCosHat2D(meshes);
if (problemName == "StochasticPollutionMollifiedBar2D")
return new StochasticPollutionMollifiedBar2D(meshes);
if (problemName == "PollutionCosHat2D")
return new PollutionCosHat2D(meshes);
if (problemName == "GaussHat2D")
return new GaussHat2D(meshes);
if (problemName == "StochasticGaussHat2D")
return new StochasticGaussHat2D(meshes);
if (problemName == "GaussHat2D")
return new GaussHat2D(meshes);
Exit(problemName + " not found")
}
mlmc/src/problems/StochasticTransportProblem.hpp
View file @ 9711af40
......@@ -14,9 +14,6 @@ public:
bool RHS() const { return rhs; }
// Todo rmv
virtual bool Dirichlet(const Point &x) const { return false; }
virtual Scalar Solution(double t, const Point &x) const = 0;
// Todo cell and point necessary?
......@@ -99,6 +96,11 @@ public:
string Name() const override { return "StochasticPollutionCosHat1D"; }
};
class PollutionCosHat1D : public IStochasticTransportProblem {
public:
PollutionCosHat1D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
};
class StochasticPollution2D : public IStochasticTransportProblem {
private:
double d = 1.0 / 16.0;
......@@ -125,45 +127,39 @@ public:
string Name() const override { return "StochasticPollution2D"; }
};
class CircleWave2D : public IStochasticTransportProblem {
double aa = 1.00;
double cc = 0.5;
class Pollution2D : public IStochasticTransportProblem {
private:
double d = 1.0 / 16.0;
public:
explicit CircleWave2D(Meshes &meshes) :
IStochasticTransportProblem(meshes) {
}
Pollution2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
double Solution(double t, const Point &x) const override {
Point midPoint = 5.0 * Point(cos(2. * Pi * t), sin(2. * Pi * t));
double r = dist(midPoint, x);
if (r < 1 / cc)
return aa * pow(cos(cc * Pi * r) + 1.0, 2.0);
Scalar Solution(double t, const Point &x) const override {
if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
return 1.0;
return 0.0;
}
bool Dirichlet(const Point &x) const override { return true; }
VectorField CircleVectorField(const Point &x) const {
return 2 * Pi * VectorField(-x[1], x[0], 0.0);
VectorField TransportFlux(const Point &x) const {
return VectorField(-1.0, 0.0, 0.0);
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return CircleVectorField(x);
};
return TransportFlux(x);
}
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override {
return CircleVectorField(x) * N;
};
return TransportFlux(x) * N;
}
string Name() const override { return "CircleWave2D"; }
string Name() const override { return "Pollution2D"; }
};
class StochasticPollutionCosHat2D : public IStochasticTransportProblem {
private:
double amplitude = 1.00;
double cc = 6.0;
public:
StochasticPollutionCosHat2D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator1D"}) {}
......@@ -188,77 +184,65 @@ public:
string Name() const override { return "StochasticPollutionCosHat2D"; }
};
class StochasticPollutionMollifiedBar2D : public IStochasticTransportProblem {
private:
double ee = 0.12, xl = 0.6, yl = 0.05;
Point midPoint = Point(0.5, 0.8);
double sh = 1 / (xl * yl + 2 * (xl + yl) * exp(1) * ee * 0.2219969080840397
+ exp(1) * 2 * Pi * ee * ee * 0.0742477533879610);
Point tl = Point(midPoint[0] - 0.5 * xl, midPoint[1] + 0.5 * yl);
Point tr = Point(midPoint[0] + 0.5 * xl, midPoint[1] + 0.5 * yl);
Point br = Point(midPoint[0] + 0.5 * xl, midPoint[1] - 0.5 * yl);
Point bl = Point(midPoint[0] - 0.5 * xl, midPoint[1] - 0.5 * yl);
double phit(const Point &x, const Point &y) const {
if (dist(x, y) < ee)
return sh * exp(1.0) * exp(-1 / (1 - pow(dist(x, y) / ee, 2)));
class PollutionCosHat2D : public IStochasticTransportProblem {
double amplitude = 1.00;
double cc = 0.5;
public:
explicit PollutionCosHat2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
double Solution(double t, const Point &x) const override {
Point midPoint = 5.0 * Point(cos(2. * Pi * t), sin(2. * Pi * t));
double r = dist(midPoint, x);
if (r < 1 / cc)
return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0;
}
public:
StochasticPollutionMollifiedBar2D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
if (abs(midPoint[0] - x[0]) <= 0.5 * xl &&
abs(midPoint[1] - x[1]) <= 0.5 * yl) {
return sh;
} else if (abs(midPoint[0] - x[0]) <= 0.5 * xl) {
Point hp1 = Point(x[0], midPoint[1] - 0.5 * yl);
Point hp2 = Point(x[0], midPoint[1] + 0.5 * yl);
return max(this->phit(x, hp1), this->phit(x, hp2));
} else if (abs(midPoint[1] - x[1]) <= 0.5 * yl) {
Point hp1 = Point(midPoint[0] - 0.5 * xl, x[1]);
Point hp2 = Point(midPoint[0] + 0.5 * xl, x[1]);
return max(this->phit(x, hp1), this->phit(x, hp2));
} else
return max(max(this->phit(x, tl), this->phit(x, tr)),
max(this->phit(x, br), this->phit(x, bl)));
VectorField CircleVectorField(const Point &x) const {
return 2 * Pi * VectorField(-x[1], x[0], 0.0);
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
return CircleVectorField(x);
};
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override {
return CircleVectorField(x) * N;
};
string Name() const override { return "StochasticPollutionMollifiedBar2D"; }
string Name() const override { return "PollutionCosHat2D"; }
};
class StochasticInitialConditionsTransport1D : public IStochasticTransportProblem {
class StochasticGaussHat2D : public IStochasticTransportProblem {
private:
double mu = 0.0;
double sigma = 1.0; // (sigma squared)
VectorField vectorField;
public:
StochasticInitialConditionsTransport1D(Meshes &meshes) :
IStochasticTransportProblem(meshes,
GeneratorNames{"NormalDistributionGenerator"}) {}
StochasticGaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0
// 1 / (sigma * sqrt(2 * PI)) exp ()
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0;
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return VectorField(1, 0, 0);
return vectorField;
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return VectorField(1, 0, 0) * N;
return vectorField * N;
}
string Name() const override { return "RiemannTransport1D"; }
string Name() const override { return "GaussHat2D"; }
};
class GaussHat2D : public IStochasticTransportProblem {
......@@ -300,75 +284,6 @@ public:
string Name() const override { return "GaussHat2D"; }
};
class StochasticGaussHat2D : public IStochasticTransportProblem {
private:
double mu = 0.0;
double sigma = 1.0; // (sigma squared)
VectorField vectorField;
public:
StochasticGaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0
// 1 / (sigma * sqrt(2 * PI)) exp ()
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0;
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return vectorField;
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return vectorField * N;
}
string Name() const override { return "GaussHat2D"; }
};
//class StochasticInflowTransport1D : public IStochasticTransportProblem {
//public:
// explicit StochasticInflowTransport1D() {}
//
// double ut(double t, const Point &x) const override {
// if (x[0] == 0.0 && x[1] >= 1.0) return 1.0;
// else return 0.0;
// }
//
// string Name() const override { return "InflowTest"; }
//};
//
//class StochasticCircleWave2D : public IStochasticTransportProblem {
// /*
// * Todo Idea: create stochastic flux
// */
// double aa = 1.00;
// double cc = 0.5;
//public:
// StochasticCircleWave2D() {}
//
// double ut(double t, const Point &x) const override {
// Point midPoint = 5.0 * Point(cos(2. * Pi * t), sin(2. * Pi * t));
// double r = dist(midPoint, x);
// if (r < 1 / cc)
// return aa * pow(cos(cc * Pi * r) + 1.0, 2.0);
// return 0.0;
// }
//
// bool Dirichlet(const Point &x) const override { return true; }
//
// VectorField B(const Point &x) const override {
// return 2 * Pi * VectorField(-x[1], x[0], 0.0);
// }
//
// string Name() const override { return "StochasticCircleWave2D"; }
//};
IStochasticTransportProblem *
CreateStochasticTransportProblem(std::string problemName, Meshes &meshes);
......
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

Imprint