Commit 2ce4e6b6 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

added GetStep size and fixed GaussHat

parent eb5293d4
...@@ -6,306 +6,308 @@ ...@@ -6,306 +6,308 @@
class IStochasticTransportProblem : public IStochasticProblem { class IStochasticTransportProblem : public IStochasticProblem {
private: private:
bool rhs = true; bool rhs = true;
// tau = CFL * h_max // tau = CFL * h_max
double CFL = 0.25; double CFL = 0.25;
double t0 = 0.0; double t0 = 0.0;
double T = 1.0; double T = 1.0;
public: public:
IStochasticTransportProblem(Meshes &meshes, GeneratorNames genNames = {}) IStochasticTransportProblem(Meshes &meshes, GeneratorNames genNames = {})
: IStochasticProblem(meshes, genNames) { : IStochasticProblem(meshes, genNames) {
config.get("CFL", CFL); config.get("CFL", CFL);
config.get("t0", t0); config.get("t0", t0);
config.get("T", T); config.get("T", T);
} }
double GetStartTime() { return t0; } double GetStepSize() { return CFL * meshes.fineMesh().MaxMeshWidth(); }
double GetEndTime() { return T; } double GetStartTime() { return t0; }
double GetCFL() { return CFL; } double GetEndTime() { return T; }
bool RHS() const { return rhs; } double GetCFL() { return CFL; }
virtual Scalar Solution(double t, const Point &x) const = 0; bool RHS() const { return rhs; }
// Todo cell and point necessary? virtual Scalar Solution(double t, const Point &x) const = 0;
virtual VectorField CellFlux(const cell &c, const Point &x) const = 0;
// Todo cell and point necessary? // Todo cell and point necessary?
virtual Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N, virtual VectorField CellFlux(const cell &c, const Point &x) const = 0;
const Point &x) const = 0;
// Todo cell and point necessary?
virtual Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const = 0;
}; };
class StochasticPollution1D : public IStochasticTransportProblem { class StochasticPollution1D : public IStochasticTransportProblem {
public: public:
StochasticPollution1D(Meshes &meshes) : StochasticPollution1D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {} IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0; if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0; return 1.0;
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return this->genContainer.vectorFieldGenerator->EvalSample(c); return this->genContainer.vectorFieldGenerator->EvalSample(c);
} }
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N, Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override { const Point &x) const override {
return this->genContainer.scalarGenerator->EvalSample(face, c); return this->genContainer.scalarGenerator->EvalSample(face, c);
} }
string Name() const override { return "StochasticPollution1D"; } string Name() const override { return "StochasticPollution1D"; }
}; };
class Pollution1D : public IStochasticTransportProblem { class Pollution1D : public IStochasticTransportProblem {
public: public:
Pollution1D(Meshes &meshes) : IStochasticTransportProblem(meshes) {} Pollution1D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0; if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0; return 1.0;
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return VectorField(-1.0, 0.0); return VectorField(-1.0, 0.0);
} }
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N, Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override { const Point &x) const override {
return 1.0; return 1.0;
} }
string Name() const override { return "Pollution1D"; } string Name() const override { return "Pollution1D"; }
}; };
class StochasticPollutionCosHat1D : public IStochasticTransportProblem { class StochasticPollutionCosHat1D : public IStochasticTransportProblem {
private: private:
double amplitude = 1.00; double amplitude = 1.00;
double cc = 6.0; double cc = 6.0;
public: public:
StochasticPollutionCosHat1D(Meshes &meshes) : StochasticPollutionCosHat1D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {} IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
Point midPoint = Point(0.2, 0.0); Point midPoint = Point(0.2, 0.0);
double r = dist(midPoint, x); double r = dist(midPoint, x);
if (r < 1 / cc) if (r < 1 / cc)
return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0); return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0; return 0.0;
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return this->genContainer.vectorFieldGenerator->EvalSample(c); return this->genContainer.vectorFieldGenerator->EvalSample(c);
} }
Scalar FaceNormalFlux(const cell &c, int face, Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override { const VectorField &N, const Point &x) const override {
return this->genContainer.scalarGenerator->EvalSample(face, c); return this->genContainer.scalarGenerator->EvalSample(face, c);
} }
string Name() const override { return "StochasticPollutionCosHat1D"; } string Name() const override { return "StochasticPollutionCosHat1D"; }
}; };
class PollutionCosHat1D : public IStochasticTransportProblem { class PollutionCosHat1D : public IStochasticTransportProblem {
public: public:
PollutionCosHat1D(Meshes &meshes) : IStochasticTransportProblem(meshes) {} PollutionCosHat1D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
}; };
class StochasticPollution2D : public IStochasticTransportProblem { class StochasticPollution2D : public IStochasticTransportProblem {
private: private:
double d = 1.0 / 16.0; double d = 1.0 / 16.0;
public: public:
StochasticPollution2D(Meshes &meshes) : StochasticPollution2D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {} IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override { 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)) if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
return 1.0; return 1.0;
return 0.0; return 0.0;
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return this->genContainer.vectorFieldGenerator->EvalSample(c); return this->genContainer.vectorFieldGenerator->EvalSample(c);
} }
Scalar FaceNormalFlux(const cell &c, int face, Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override { const VectorField &N, const Point &x) const override {
return this->genContainer.scalarGenerator->EvalSample(face, c); return this->genContainer.scalarGenerator->EvalSample(face, c);
} }
string Name() const override { return "StochasticPollution2D"; } string Name() const override { return "StochasticPollution2D"; }
}; };
class Pollution2D : public IStochasticTransportProblem { class Pollution2D : public IStochasticTransportProblem {
private: private:
double d = 1.0 / 16.0; double d = 1.0 / 16.0;
public: public:
Pollution2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {} Pollution2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 4 * d)) if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 4 * d))
return 32.0; return 32.0;
return 0.0; return 0.0;
} }
VectorField TransportFlux(const Point &x) const { VectorField TransportFlux(const Point &x) const {
return VectorField(0.0, -1.0); return VectorField(0.0, -1.0);
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return TransportFlux(x); return TransportFlux(x);
} }
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N, Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override { const Point &x) const override {
return TransportFlux(x) * N; return TransportFlux(x) * N;
} }
string Name() const override { return "Pollution2D"; } string Name() const override { return "Pollution2D"; }
}; };
class StochasticPollutionCosHat2D : public IStochasticTransportProblem { class StochasticPollutionCosHat2D : public IStochasticTransportProblem {
private: private:
double amplitude = 1.00; double amplitude = 1.00;
double cc = 6.0; double cc = 6.0;
public: public:
StochasticPollutionCosHat2D(Meshes &meshes) : StochasticPollutionCosHat2D(Meshes &meshes) :
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator1D"}) {} IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator1D"}) {}
double Solution(double t, const Point &x) const override { double Solution(double t, const Point &x) const override {
Point midPoint = Point(0.5, 0.8); Point midPoint = Point(0.5, 0.8);
double r = dist(midPoint, x); double r = dist(midPoint, x);
if (r < 1 / cc) if (r < 1 / cc)
return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0); return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0; return 0.0;
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return this->genContainer.vectorFieldGenerator->EvalSample(c); return this->genContainer.vectorFieldGenerator->EvalSample(c);
} }
Scalar FaceNormalFlux(const cell &c, int face, Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override { const VectorField &N, const Point &x) const override {
return this->genContainer.scalarGenerator->EvalSample(face, c); return this->genContainer.scalarGenerator->EvalSample(face, c);
} }
string Name() const override { return "StochasticPollutionCosHat2D"; } string Name() const override { return "StochasticPollutionCosHat2D"; }
}; };
class PollutionCosHat2D : public IStochasticTransportProblem { class PollutionCosHat2D : public IStochasticTransportProblem {
double amplitude = 2.0 / Pi; double amplitude = 2.0 / Pi;
public: public:
explicit PollutionCosHat2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {} explicit PollutionCosHat2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
double Solution(double t, const Point &x) const override { double Solution(double t, const Point &x) const override {
Point midPoint = 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t)) Point midPoint = 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t))
+ Point(1.0, 1.0)); + Point(1.0, 1.0));
double r = dist(midPoint, x); double r = dist(midPoint, x);
if (r < 1.0 / (2.0 * Pi)) return amplitude * pow(cos(r / 2.0), 2.0); if (r < 1.0 / (2.0 * Pi)) return amplitude * pow(cos(r / 2.0), 2.0);
return 0.0; return 0.0;
} }
VectorField CircleVectorField(const Point &x) const { VectorField CircleVectorField(const Point &x) const {
return 2 * Pi * VectorField(-x[1], x[0], 0.0); return 2 * Pi * VectorField(-x[1], x[0], 0.0);
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return CircleVectorField(x); return CircleVectorField(x);
}; };
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N, Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override { const Point &x) const override {
return CircleVectorField(x) * N; return CircleVectorField(x) * N;
}; };
string Name() const override { return "PollutionCosHat2D"; } string Name() const override { return "PollutionCosHat2D"; }
}; };
class StochasticGaussHat2D : public IStochasticTransportProblem { class StochasticGaussHat2D : public IStochasticTransportProblem {
private: private:
double mu = 0.0; double mu = 0.0;
double sigma = 1.0; // (sigma squared) double sigma = 1.0; // (sigma squared)
VectorField vectorField; VectorField vectorField;
public: public:
StochasticGaussHat2D(Meshes &meshes) : StochasticGaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)), vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(meshes) {} IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0 // Todo 2D Gaussfunction for t = 0
// 1 / (sigma * sqrt(2 * PI)) exp () // 1 / (sigma * sqrt(2 * PI)) exp ()
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0; if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0; return 1.0;
} }
VectorField CircleVectorField(const Point &x) const { VectorField CircleVectorField(const Point &x) const {
return 2 * Pi * VectorField(-x[1], x[0], 0.0); return 2 * Pi * VectorField(-x[1], x[0], 0.0);
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return CircleVectorField(x); return CircleVectorField(x);
} }
Scalar FaceNormalFlux(const cell &c, int face, Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override { const VectorField &N, const Point &x) const override {
return CircleVectorField(x) * N; return CircleVectorField(x) * N;
} }
string Name() const override { return "GaussHat2D"; } string Name() const override { return "GaussHat2D"; }
}; };
class GaussHat2D : public IStochasticTransportProblem { class GaussHat2D : public IStochasticTransportProblem {
private: private:
Tensor sigma; Tensor sigma;
public: public:
GaussHat2D(Meshes &meshes) : GaussHat2D(Meshes &meshes) :
sigma(Tensor()), sigma(Tensor()),
IStochasticTransportProblem(meshes) { IStochasticTransportProblem(meshes) {
sigma[0][0] = 0.005; sigma[0][0] = 0.005;
sigma[0][1] = 0.0; sigma[0][1] = 0.0;
sigma[1][0] = 0.0; sigma[1][0] = 0.0;
sigma[1][1] = 0.005; sigma[1][1] = 0.005;
// sigma[2][2] = 1.0; // Todo might not be working for 2D Point sigma[2][2] = 1.0; // Todo this is needed even in 2D case which is not quite clean
} }
Point mu(double t) const { Point mu(double t) const {
return 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t)) + Point(1.0, 1.0)); return 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t)) + Point(1.0, 1.0));
} }
Scalar Solution(double t, const Point &x) const override { Scalar Solution(double t, const Point &x) const override {
double factor = 1 / sqrt(pow(2 * Pi, 2) * sigma.det()); double factor = 1 / sqrt(pow(2 * Pi, 2) * sigma.det());
VectorField diff = (x - mu(t)); VectorField diff = (x - mu(t));
VectorField temp = Invert(sigma) * diff; VectorField temp = Invert(sigma) * diff;
Scalar exponent = -0.5 * diff * temp; Scalar exponent = -0.5 * diff * temp;
return factor * exp(exponent); return factor * exp(exponent);
} }
VectorField CircleVectorField(const Point &x) const { VectorField CircleVectorField(const Point &x) const {
return 2 * Pi * VectorField(-x[1], x[0], 0.0); return 2 * Pi * VectorField(-x[1], x[0], 0.0);
} }
VectorField CellFlux(const cell &c, const Point &x) const override { VectorField CellFlux(const cell &c, const Point &x) const override {
return CircleVectorField(x); return CircleVectorField(x);
} }
Scalar FaceNormalFlux(const cell &c, int face, Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override { const VectorField &N, const Point &x) const override {
return CircleVectorField(x) * N; return CircleVectorField(x) * N;
} }
string Name() const override { return "GaussHat2D"; } string Name() const override { return "GaussHat2D"; }
}; };
IStochasticTransportProblem * IStochasticTransportProblem *
......
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