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

refactored problems

parent 4ff5b5aa
......@@ -2,7 +2,7 @@
#define STOCHASTICREACTIONPROBLEM_HPP
#include "StochasticProblem.hpp"
#include "HybridFluxTransport.hpp"
#include "StochasticHybridFlux.hpp"
class StochasticReactionProblem : public StochasticProblem {
......@@ -145,7 +145,7 @@ public:
};
class HybridExponentialReaction2D : public ExponentialReaction2D,
public HybridFluxTransport {
public StochasticHybridFlux {
VectorField CellConvection(int sd, const cell &c, const Point &x) const override {
return convection * hybrid->EvaluateCellFlux(*flux, c);
......@@ -160,7 +160,7 @@ class HybridExponentialReaction2D : public ExponentialReaction2D,
};
class HybridLogisticReaction2D : public LogisticReaction2D,
public HybridFluxTransport {
public StochasticHybridFlux {
VectorField CellConvection(int sd, const cell &c,
const Point &x) const override {
......
......@@ -2,19 +2,33 @@
#define MLMC_STOCHASTICTRANSPORTPROBLEM_HPP
#include "StochasticProblem.hpp"
#include "assemble/HybridEllipticAssemble.hpp"
#include "StochasticHybridFlux.hpp"
class StochasticTransportProblem : public StochasticProblem {
/*
* Stochastic Transport Problem can handle:
* - Random right hand sides ->
* - Random initial conditions ->
* - Random inflow conditions ->
* - Random coefficients -> StochasticHybridFluxTransport
* -> warum hat der hier nicht einfach einen generator
*
* Design your problem by deriving from
*
*/
protected:
std::shared_ptr<Discretization> disc;
std::shared_ptr<HybridEllipticAssemble> assemble;
VectorField fluxField = Point(0.0, -1.0, 0.0);
bool exactAverage = false; // Todo maybe
public:
bool rhs = false;
VectorField b = Point(0.0, -1.0, 0.0);
explicit StochasticTransportProblem() {
// Deal with b
}
bool rhs = false;
virtual ~StochasticTransportProblem() = default;
bool RHS() const { return rhs; }
......@@ -24,81 +38,77 @@ public:
virtual bool Dirichlet(const Point &x) const { return false; }
virtual VectorField B(const Point &x) const { return b; }
virtual VectorField B(const Point &x) const {
return fluxField;
}
virtual VectorField cellB(const cell &c, const Point &x) const {
return assemble->EvaluateCellFlux(*fieldSample, c);
return B(x);
}
virtual double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const {
return assemble->EvaluateNormalFlux(*fieldSample, c, f);
virtual double faceB(const cell &c, int f, const VectorField &N,
const Point &x) const {
return B(x) * N;
}
};
class StochasticPollution1D : public StochasticTransportProblem {
private:
double c = 1.0;
class StochasticPollution1D : public StochasticTransportProblem,
public StochasticHybridFlux {
public:
StochasticPollution1D() {
disc = std::make_shared<Discretization>("RT0_P0", 1, 1);
assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), nullptr);
}
StochasticPollution1D() {}
double ut(double t, const Point &x) const override {
if (abs(c * t - x[0] + 0.5) > 0.06251)
if (abs(1.0 * t - x[0] + 0.5) > 0.06251)
return 0.0;
else
return 1.0;
}
// VectorField cellB(const cell &c, const Point &x) const override {
// return VectorField(1.0, 0.0);
// }
//
// double faceB(const cell &c, int f,
// const VectorField &N, const Point &x) const override {
// return VectorField(1.0, 0.0) * N;
// }
VectorField cellB(const cell &c, const Point &x) const override {
return hybrid->EvaluateCellFlux(*flux, c);
}
string Name() const override { return "Stochastic Pollution 1D"; }
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return hybrid->EvaluateNormalFlux(*flux, c, f);
}
string Name() const override { return "StochasticPollution1D"; }
};
class StochasticPollutionCosHat1D : public StochasticTransportProblem {
class StochasticPollutionCosHat1D : public StochasticTransportProblem,
public StochasticHybridFlux {
private:
double amplitude = 1.00;
double cc = 6.0;
public:
StochasticPollutionCosHat1D() {
disc = std::make_shared<Discretization>("RT0_P0", 1, 1);
assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), nullptr);
}
StochasticPollutionCosHat1D() {}
double ut(double t, const Point &x) const override {
Point midPoint = Point(0.2, 0.0);
double r = dist(midPoint, x);
if (r < 1 / cc)
return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0;
}
if (r < 0.1)
return 1.0;
else
return 0.0;
// if (r < 1 / cc)
// return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
// return 0.0;
VectorField cellB(const cell &c, const Point &x) const override {
return hybrid->EvaluateCellFlux(*flux, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return hybrid->EvaluateNormalFlux(*flux, c, f);
}
string Name() const override { return "Stochastic Pollution CosHat 2D"; }
string Name() const override { return "StochasticPollutionCosHat1D"; }
};
class StochasticPollution2D : public StochasticTransportProblem {
class StochasticPollution2D : public StochasticTransportProblem,
public StochasticHybridFlux {
public:
StochasticPollution2D() {
disc = std::make_shared<Discretization>("RT0_P0", 1, 2);
assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), nullptr);
}
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))
......@@ -106,19 +116,26 @@ public:
else return 0.0;
}
VectorField cellB(const cell &c, const Point &x) const override {
return hybrid->EvaluateCellFlux(*flux, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return hybrid->EvaluateNormalFlux(*flux, c, f);
}
string Name() const override { return "Stochastic Pollution 2D"; }
};
class StochasticPollutionCosHat2D : public StochasticTransportProblem {
class StochasticPollutionCosHat2D : public StochasticTransportProblem,
public StochasticHybridFlux {
private:
double amplitude = 1.00;
double cc = 6.0;
public:
StochasticPollutionCosHat2D() {
disc = std::make_shared<Discretization>("RT0_P0", 1, 2);
assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), nullptr);
}
StochasticPollutionCosHat2D() {}
double ut(double t, const Point &x) const override {
Point midPoint = Point(0.5, 0.8);
......@@ -128,10 +145,20 @@ public:
return 0.0;
}
string Name() const override { return "Stochastic Pollution 2D"; }
VectorField cellB(const cell &c, const Point &x) const override {
return hybrid->EvaluateCellFlux(*flux, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return hybrid->EvaluateNormalFlux(*flux, c, f);
}
string Name() const override { return "StochasticPollutionCosHat2D"; }
};
class StochasticPollutionMollifiedBar2D : public StochasticTransportProblem {
class StochasticPollutionMollifiedBar2D : public StochasticTransportProblem,
public StochasticHybridFlux {
private:
double ee = 0.12;
double xl = 0.6;
......@@ -151,10 +178,7 @@ private:
}
public:
StochasticPollutionMollifiedBar2D() {
disc = std::make_shared<Discretization>("RT0_P0", 1, 2);
assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), nullptr);
}
StochasticPollutionMollifiedBar2D() {}
double ut(double t, const Point &x) const override {
if (abs(midPoint[0] - x[0]) <= 0.5 * xl && abs(midPoint[1] - x[1]) <= 0.5 * yl) {
......@@ -172,7 +196,77 @@ public:
max(this->phit(x, br), this->phit(x, bl)));
}
VectorField cellB(const cell &c, const Point &x) const override {
return hybrid->EvaluateCellFlux(*flux, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return hybrid->EvaluateNormalFlux(*flux, c, f);
}
string Name() const override { return "Stochastic Pollution 2D"; }
};
#endif //MLMC_STOCHASTICTRANSPORTPROBLEM_HPP
class StochasticInitialConditionsTransport1D : public StochasticTransportProblem {
double c = 1.0;
public:
explicit StochasticInitialConditionsTransport1D() {
rhs = false;
}
double Amplitude(double r) const {
if (abs(r) > 0.06251) return 0.0;
return 1.0;
}
double ut(double t, const Point &x) const override {
return Amplitude(c * t - x[0] + 0.5);
}
VectorField B(const Point &) const override {
return VectorField(1.0, 0.0);
}
string Name() const override { return "RiemannTransport1D"; }
};
class StochasticInflowTransport1D : public StochasticTransportProblem {
public:
explicit StochasticInflowTransport1D() : StochasticTransportProblem() {}
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 StochasticTransportProblem {
/*
* 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"; }
};
#endif
\ No newline at end of file
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