StochasticReactionProblem.hpp 5.42 KB
Newer Older
niklas.baumgarten's avatar
niklas.baumgarten committed
1
2
#ifndef STOCHASTICREACTIONPROBLEM_HPP
#define STOCHASTICREACTIONPROBLEM_HPP
3

4
#include "IStochasticProblem.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30


// Todo rmv
#include "pdesolver/assembling/HybridEllipticAssemble.hpp"


class StochasticHybridFlux {
public:
    HybridEllipticAssemble *hybrid = nullptr;
    Vector *flux = nullptr;

    explicit StochasticHybridFlux() = default;;

    StochasticHybridFlux(HybridEllipticAssemble *hybrid, Vector *flux) :
        hybrid(hybrid), flux(flux) {}

    void SetHybridEllipticAssemble(HybridEllipticAssemble *hybridEllipticAssemble) {
        this->hybrid = hybridEllipticAssemble;
    }

    void SetFlux(Vector *normalFlux) {
        this->flux = normalFlux;
    }
};


31
32


niklas.baumgarten's avatar
niklas.baumgarten committed
33
class IStochasticReactionProblem : public IStochasticProblem {
34
35
36
37
protected:
    double convection = 1.0;
    double diffusion = 0.001;
    double reaction = 1.0;
38

39
public:
40
41
    IStochasticReactionProblem(Meshes &meshes, GeneratorNames genNames={}) :
        IStochasticProblem(meshes, genNames) {
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
        config.get("Diffusion", diffusion);
        config.get("Convection", convection);
        config.get("Reaction", reaction);
    }

    double GetConvection() const {
        return convection;
    }

    double GetDiffusion() const {
        return diffusion;
    }

    double GetReaction() const {
        return reaction;
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
59
    virtual ~IStochasticReactionProblem() = default;
60
61
62
63
64

    virtual Scalar Concentration(double t, const Point &x) const {
        return 0.0;
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
65
    virtual Scalar Load(const Point &x) const {
66
67
68
69
70
71
72
        return 0.0;
    }

    virtual VectorField Flux(int, const Point &) const {
        return zero;
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
73
    virtual VectorField Convection(const Point &x) const {
74
75
76
        return convection * zero;
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
77
78
    virtual VectorField CellConvection(const cell &c, const Point &x) const {
        return Convection(x);
79
80
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
81
82
83
84
85
    virtual double FaceConvection(const cell &c,
                                  int f,
                                  const VectorField &N,
                                  const Point &x) const {
        return Convection(x) * N;
86
87
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
88
    virtual Tensor Diffusion(const Point &x) const {
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
        return diffusion * Zero;
    }

    virtual Scalar Reaction(const Point &, double c) const {
        return reaction * 0.0;
    }

    virtual Scalar DerivativeReaction(const Point &, double c) const {
        return 0.0;
    }
};

inline Scalar u0(const Point &x) {
    double thickness = 1.0 / 8.0;
    double centerStart = 1 - 1.5 * thickness;
    if ((abs(centerStart - x[1]) < thickness / 2) &&
        (abs(0.5 - x[0]) < 3 * thickness))
        return 1;
    else return 0.0;
}

niklas.baumgarten's avatar
niklas.baumgarten committed
110
class ExponentialReaction2D : public IStochasticReactionProblem {
111
112
public:

niklas.baumgarten's avatar
niklas.baumgarten committed
113
    ExponentialReaction2D(Meshes &meshes) :
114
        IStochasticReactionProblem(meshes) {};
115
116
117
118
119

    Scalar Concentration(double t, const Point &x) const override {
        return u0(x);
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
120
    VectorField Convection(const Point &x) const override {
121
122
123
        return VectorField(0.0, -1.0);
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
124
    Tensor Diffusion(const Point &x) const override {
125
        return this->diffusion * One;
126
127
128
    }

    Scalar Reaction(const Point &, double c) const override {
129
        return this->reaction * c;
130
131
132
    }

    Scalar DerivativeReaction(const Point &, double c) const override {
133
        return this->reaction;
134
135
136
137
138
    }

    string Name() const override { return "ExponentialReaction2D"; }
};

niklas.baumgarten's avatar
niklas.baumgarten committed
139
class LogisticReaction2D : public IStochasticReactionProblem {
140
141
142
143
private:
    double r0 = 1.0;
    double r1 = 0.0;
public:
niklas.baumgarten's avatar
niklas.baumgarten committed
144
    LogisticReaction2D(Meshes &meshes) :
145
        IStochasticReactionProblem(meshes) {
146
147
148
149
150
151
152
153
        config.get("Reaction0", r0);
        config.get("Reaction1", r1);
    }

    Scalar Concentration(double t, const Point &x) const override {
        return u0(x);
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
154
    VectorField Convection(const Point &x) const override {
155
156
157
        return VectorField(0.0, -1.0);
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
158
    Tensor Diffusion(const Point &x) const override {
159
        return this->diffusion * One;
160
161
162
163
164
165
166
167
168
169
170
171
172
    }

    Scalar Reaction(const Point &x, double c) const override {
        return r0 * c - r1 * c * c;
    }

    Scalar DerivativeReaction(const Point &x, double c) const override {
        return r0 - 2 * r1 * c;
    }

    string Name() const override { return "LogisticReaction2D"; }
};

niklas.baumgarten's avatar
niklas.baumgarten committed
173
class HybridExponentialReaction2D : public ExponentialReaction2D,
niklas.baumgarten's avatar
niklas.baumgarten committed
174
                                    public StochasticHybridFlux {
175

niklas.baumgarten's avatar
niklas.baumgarten committed
176
    VectorField CellConvection(const cell &c, const Point &x) const override {
177
        return this->convection * hybrid->EvaluateCellFlux(*flux, c);
178
179
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
180
181
182
183
    double FaceConvection(const cell &c,
                          int f,
                          const VectorField &N,
                          const Point &x) const override {
184
        return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
185
186
187
188
189
    }

    string Name() const override { return "HybridExponentialReaction2D"; }
};

niklas.baumgarten's avatar
niklas.baumgarten committed
190
class HybridLogisticReaction2D : public LogisticReaction2D,
niklas.baumgarten's avatar
niklas.baumgarten committed
191
                                 public StochasticHybridFlux {
192

niklas.baumgarten's avatar
niklas.baumgarten committed
193
    VectorField CellConvection(const cell &c, const Point &x) const override {
194
        return this->convection * hybrid->EvaluateCellFlux(*flux, c);
195
196
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
197
198
199
200
    double FaceConvection(const cell &c,
                          int f,
                          const VectorField &N,
                          const Point &x) const override {
201
        return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
202
203
204
205
206
    }

    string Name() const override { return "HybridLogisticReaction2D"; }
};

niklas.baumgarten's avatar
niklas.baumgarten committed
207
#endif //STOCHASTICREACTIONPROBLEM_HPP