Line data Source code
1 : #ifndef SOLVER_LOCAL_RIEMANN_H
2 : #define SOLVER_LOCAL_RIEMANN_H
3 :
4 : #include "Set/Base.H"
5 : #include "Set/Set.H"
6 :
7 : namespace Solver
8 : {
9 : namespace Local
10 : {
11 : namespace Riemann
12 : {
13 : struct State {
14 : Set::Scalar rho = NAN;
15 : Set::Scalar M_normal = NAN;
16 : Set::Scalar M_tangent = NAN;
17 : Set::Scalar E = NAN;
18 : // Construtor for convenience
19 : State() { rho = 0.0; M_normal = 0.0, M_tangent = 0.0; E = 0.0;}
20 164505624 : State(Set::Scalar a_rho, Set::Scalar a_M_normal, Set::Scalar a_M_tangent, Set::Scalar a_E)
21 164505624 : : rho(a_rho), M_normal(a_M_normal), M_tangent(a_M_tangent), E(a_E) {}
22 109670400 : State(Set::Patch<const Set::Scalar> density_mf, Set::Patch<const Set::Scalar> momentum_mf, Set::Patch<const Set::Scalar> energy_mf, int i, int j, int k, int direction)
23 109670400 : {
24 109670400 : rho = density_mf(i,j,k);
25 109670400 : if (direction == 0)
26 : {
27 54835200 : M_normal = momentum_mf(i,j,k,0);
28 54835200 : M_tangent = momentum_mf(i,j,k,1);
29 : }
30 54835200 : else if (direction == 1)
31 : {
32 54835200 : M_normal = momentum_mf(i,j,k,1);
33 54835200 : M_tangent = momentum_mf(i,j,k,0);
34 : }
35 : else
36 : {
37 0 : Util::Abort(INFO, "Not supported yet");
38 : }
39 109670400 : E = energy_mf(i,j,k);
40 109670400 : }
41 :
42 :
43 0 : friend std::ostream &operator<<(std::ostream &os, const State &state)
44 : {
45 0 : os << "rho=" << state.rho << ", ";
46 0 : os << "Mn=" << state.M_normal << ", ";
47 0 : os << "Mt=" << state.M_tangent << ", ";
48 0 : os << "E=" << state.E << ", ";
49 : // do stuf
50 0 : return os;
51 : }
52 : void operator += (const State &a)
53 : {
54 : rho += a.rho; M_normal += a.M_normal; M_tangent += a.M_tangent; E += a.E;
55 : };
56 : void operator -= (const State &a)
57 : {
58 : rho -= a.rho; M_normal -= a.M_normal; M_tangent -= a.M_tangent; E -= a.E;
59 : };
60 : void operator *= (const Set::Scalar alpha)
61 : {
62 : rho *= alpha; M_normal *= alpha; M_tangent *= alpha; E *= alpha;
63 : };
64 : void operator /= (const Set::Scalar alpha)
65 : {
66 : rho /= alpha; M_normal /= alpha; M_tangent /= alpha; E /= alpha;
67 : };
68 : friend State operator + (const State &a, const State &b)
69 : {
70 : return State(a.rho+b.rho, a.M_normal+b.M_normal, a.M_tangent+b.M_tangent, a.E+b.E);
71 : }
72 54835200 : friend State operator - (const State &a, const State &b)
73 : {
74 54835200 : return State(a.rho-b.rho, a.M_normal-b.M_normal, a.M_tangent-b.M_tangent, a.E-b.E);
75 : }
76 54835200 : friend State operator * (const Set::Scalar alpha, const State &b)
77 : {
78 54835200 : return State(b.rho*alpha, b.M_normal*alpha, b.M_tangent*alpha, b.E*alpha);
79 : }
80 : friend State operator * (const State &b, const Set::Scalar alpha)
81 : {
82 : return State(b.rho*alpha, b.M_normal*alpha, b.M_tangent*alpha, b.E*alpha);
83 : }
84 54835200 : friend State operator / (const State &b, const Set::Scalar alpha)
85 : {
86 54835200 : return State(b.rho/alpha, b.M_normal/alpha, b.M_tangent/alpha, b.E/alpha);
87 : }
88 : };
89 :
90 : struct Flux {
91 : Set::Scalar mass;
92 : Set::Scalar momentum_normal;
93 : Set::Scalar momentum_tangent;
94 : Set::Scalar energy;
95 73113616 : Flux() : mass(0.0), momentum_normal(0.0), momentum_tangent(0.0), energy(0.0) {}
96 36556800 : Flux(Set::Scalar a_mass, Set::Scalar a_momentum_normal, Set::Scalar a_momentum_tangent, Set::Scalar a_energy) :
97 36556800 : mass(a_mass), momentum_normal(a_momentum_normal),
98 36556800 : momentum_tangent(a_momentum_tangent), energy(a_energy) {}
99 :
100 0 : Flux(State &in) :
101 0 : mass(in.rho),
102 0 : momentum_normal(in.M_normal),
103 0 : momentum_tangent(in.M_tangent),
104 0 : energy(in.E) {}
105 :
106 0 : friend std::ostream &operator<<(std::ostream &os, const Flux &flux)
107 : {
108 0 : os << "mass=" << flux.mass << ", ";
109 0 : os << "Mn=" << flux.momentum_normal << ", ";
110 0 : os << "Mt=" << flux.momentum_tangent << ", ";
111 0 : os << "E=" << flux.energy << ", ";
112 : // do stuff
113 0 : return os;
114 : }
115 : void operator += (const Flux &a)
116 : {
117 : mass += a.mass;
118 : momentum_normal += a.momentum_normal; momentum_tangent += a.momentum_tangent;
119 : energy += a.energy;
120 : }
121 : void operator -= (const Flux &a)
122 : {
123 : mass -= a.mass;
124 : momentum_normal -= a.momentum_normal; momentum_tangent -= a.momentum_tangent;
125 : energy -= a.energy;
126 : }
127 : void operator *= (const Set::Scalar alpha)
128 : {
129 : mass *= alpha;
130 : momentum_normal *= alpha;
131 : momentum_tangent *= alpha;
132 : energy *= alpha;
133 : }
134 : void operator /= (const Set::Scalar alpha)
135 : {
136 : mass /= alpha;
137 : momentum_normal /= alpha;
138 : momentum_tangent /= alpha;
139 : energy /= alpha;
140 : }
141 0 : friend Flux operator + (const Flux &a, const Flux &b)
142 : {
143 0 : return Flux(a.mass+b.mass, a.momentum_normal+b.momentum_normal,
144 0 : a.momentum_tangent+b.momentum_tangent, a.energy+b.energy);
145 : }
146 0 : friend Flux operator - (const Flux &a, const Flux &b)
147 : {
148 0 : return Flux(a.mass-b.mass, a.momentum_normal-b.momentum_normal,
149 0 : a.momentum_tangent-b.momentum_tangent, a.energy-b.energy);
150 : }
151 36556800 : friend Flux operator * (const Flux &a, const Set::Scalar beta)
152 : {
153 36556800 : return Flux(a.mass*beta, a.momentum_normal*beta,
154 36556800 : a.momentum_tangent*beta, a.energy*beta);
155 : }
156 0 : friend Flux operator * (const Set::Scalar beta,const Flux &a)
157 : {
158 0 : return Flux(a.mass*beta, a.momentum_normal*beta,
159 0 : a.momentum_tangent*beta, a.energy*beta);
160 : }
161 0 : friend Flux operator / (const Flux &a, const Set::Scalar beta)
162 : {
163 0 : return Flux(a.mass/beta, a.momentum_normal/beta,
164 0 : a.momentum_tangent/beta, a.energy/beta);
165 : }
166 :
167 : };
168 :
169 :
170 : class Riemann
171 : {
172 : public:
173 : virtual Flux Solve(State lo, State hi, Set::Scalar gamma, Set::Scalar p_ref, Set::Scalar small) = 0;
174 8 : virtual ~Riemann() = default;
175 : };
176 :
177 :
178 : }
179 : }
180 : }
181 :
182 :
183 :
184 :
185 : #endif
|
