-
Notifications
You must be signed in to change notification settings - Fork 0
/
RadMom1D_Recon.h
252 lines (215 loc) · 8.06 KB
/
RadMom1D_Recon.h
1
2
3
4
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
/*! \file RadMom1D_Recon.h
\brief Implementation of subroutines prototyped in RadMom1D_Recon.h file. */
#ifndef _RADMOM_1D_RECON_INCLUDED
#define _RADMOM_1D_RECON_INCLUDED
/* Include 1D RadMom solution header file. */
#ifndef _RADMOM1D_IO_INCLUDED
#include "RadMom1D_IO.h"
#endif // _RADMOM1D_IO_INCLUDED
/******************************************************//**
* Routine: Linear_Reconstruction_LeastSquares
*
* Peforms the reconstruction of a limited piecewise
* linear solution state within each cell of the
* computational mesh. A least squares approach is
* used in the evaluation of the unlimited solution
* gradients. Several slope limiters may be used.
*
********************************************************/
template<class cState, class pState>
void Linear_Reconstruction_LeastSquares(RadMom1D_UniformMesh<cState, pState> *Soln,
const int Limiter) {
int i, n, n2, n_pts, index[2];
double u0Min, u0Max, uQuad[2], phi;
double Dx, DxDx_ave;
pState DU, DUDx_ave;
int ICl, ICu;
ICl = Soln[0].ICl;
ICu = Soln[0].ICu;
/* Carry out the limited solution reconstruction in
*each cell. */
for ( i = ICl ; i <= ICu ; ++i ) {
n_pts = 2;
index[0] = i-1;
index[1] = i+1;
DUDx_ave.Vacuum();
DxDx_ave = ZERO;
for ( n2 = 0 ; n2 <= n_pts-1 ; ++n2 ) {
Dx = Soln[ index[n2] ].X.x - Soln[i].X.x;
DU = Soln[ index[n2] ].W - Soln[i].W;
DUDx_ave += DU*Dx;
DxDx_ave += Dx*Dx;
} /* endfor */
DUDx_ave = DUDx_ave/double(n_pts);
DxDx_ave = DxDx_ave/double(n_pts);
Soln[i].dWdx = DUDx_ave/DxDx_ave;
for ( n = 1 ; n <= DU.NumVar() ; ++n ) {
u0Min = Soln[i].W[n];
u0Max = u0Min;
for ( n2 = 0 ; n2 <= n_pts-1 ; ++n2 ) {
u0Min = min(u0Min, Soln[ index[n2] ].W[n]);
u0Max = max(u0Max, Soln[ index[n2] ].W[n]);
} /* endfor */
// Compute unlimited solution
uQuad[0] = Soln[i].W[n] - HALF*Soln[i].dWdx[n]*Soln[i].X.dx;
uQuad[1] = Soln[i].W[n] + HALF*Soln[i].dWdx[n]*Soln[i].X.dx;
switch(Limiter) {
case LIMITER_ZERO :
phi = ZERO;
break;
case LIMITER_ONE :
phi = ONE;
break;
case LIMITER_BARTH_JESPERSEN :
phi = Limiter_BarthJespersen(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VENKATAKRISHNAN :
phi = Limiter_Venkatakrishnan(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VENKATAKRISHNAN_CORRECTED :
phi = Limiter_Venkatakrishnan_Modified(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VANLEER :
phi = Limiter_VanLeer(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VANALBADA :
phi = Limiter_VanAlbada(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
default:
phi = Limiter_BarthJespersen(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
} /* endswitch */
Soln[i].phi[n] = phi;
} /* endfor */
} /* endfor */
// Check to make sure reconstructed solution at each of the cell's interface is realizable
// If not, we set the slope limiters back to zero
// pState W_recon;
// double dX;
// // East boundary
// dX = Soln[i].X.dx/2.0; // dX < 0 in this case;
// W_recon.Reconstruct(Soln[i].W, Soln[i].phi, Soln[i].dWdx, dX);
// if (W_recon.Unphysical_Properties()) {
// Soln[i].phi.Vacuum();
// }
//
// // West boundary
// dX = -Soln[i].X.dx/2.0; // dX < 0 in this case;
// W_recon.Reconstruct(Soln[i].W, Soln[i].phi, Soln[i].dWdx, dX);
// if (W_recon.Unphysical_Properties()) {
// Soln[i].phi.Vacuum();
// }
// Gradients and limiters at ghost cells
// Soln[0].dWdx = Soln[1].phi^Soln[1].dWdx;
// Soln[0].phi.Ones();
// Soln[Number_of_Cells+1].dWdx = Soln[Number_of_Cells].phi^Soln[Number_of_Cells].dWdx;
// Soln[Number_of_Cells+1].phi.Ones();
}
/******************************************************//**
* Routine: Linear_Reconstruction_GreenGauss
*
* Peforms the reconstruction of a limited piecewise
* linear solution state within each cell of the
* computational mesh. A Green-Gauss approach is used
* in the evaluation of the unlimited solution
* gradients. Several slope limiters may be used.
*
********************************************************/
template<class cState, class pState>
void Linear_Reconstruction_GreenGauss(RadMom1D_UniformMesh<cState, pState> *Soln,
const int Limiter) {
int i, n;
double u0Min, u0Max, uQuad[2], phi;
int ICl, ICu;
ICl = Soln[0].ICl;
ICu = Soln[0].ICu;
/* Carry out the limited solution reconstruction in
*each cell. */
for ( i = ICl ; i <= ICu ; ++i ) {
Soln[i].dWdx = HALF*(Soln[i+1].W - Soln[i-1].W)/Soln[i].X.dx;
for ( n = 1 ; n <= Soln[i].W.NumVar() ; ++n ) {
u0Min = min(Soln[i-1].W[n], Soln[i].W[n]);
u0Min = min(u0Min, Soln[i+1].W[n]);
u0Max = max(Soln[i-1].W[n], Soln[i].W[n]);
u0Max = max(u0Max, Soln[i+1].W[n]);
// Compute unlimited solution
uQuad[0] = Soln[i].W[n] - HALF*Soln[i].dWdx[n]*Soln[i].X.dx;
uQuad[1] = Soln[i].W[n] + HALF*Soln[i].dWdx[n]*Soln[i].X.dx;
switch(Limiter) {
case LIMITER_ZERO :
phi = ZERO;
break;
case LIMITER_ONE :
phi = ONE;
break;
case LIMITER_BARTH_JESPERSEN :
phi = Limiter_BarthJespersen(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VENKATAKRISHNAN :
phi = Limiter_Venkatakrishnan(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VENKATAKRISHNAN_CORRECTED :
phi = Limiter_Venkatakrishnan_Modified(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VANLEER :
phi = Limiter_VanLeer(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
case LIMITER_VANALBADA :
phi = Limiter_VanAlbada(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
default:
phi = Limiter_BarthJespersen(uQuad, Soln[i].W[n], u0Min, u0Max, 2);
break;
} /* endswitch */
Soln[i].phi[n] = phi;
} /* endfor */
} /* endfor */
// Check to make sure reconstructed solution at each of the cell's interface is realizable
// If not, we set the slope limiters back to zero
pState W_recon;
double dX;
// East boundary
dX = Soln[i].X.dx/2.0; // dX < 0 in this case;
W_recon.Reconstruct(Soln[i].W, Soln[i].phi, Soln[i].dWdx, dX);
if (W_recon.Unphysical_Properties()) {
Soln[i].phi.Vacuum();
}
// West boundary
dX = -Soln[i].X.dx/2.0; // dX < 0 in this case;
W_recon.Reconstruct(Soln[i].W, Soln[i].phi, Soln[i].dWdx, dX);
if (W_recon.Unphysical_Properties()) {
Soln[i].phi.Vacuum();
}
// Gradients and limiters at ghost cells
// Soln[0].dWdx = Soln[1].phi^Soln[1].dWdx;
// Soln[0].phi.Ones();
// Soln[Number_of_Cells+1].dWdx = Soln[Number_of_Cells].phi^Soln[Number_of_Cells].dWdx;
// Soln[Number_of_Cells+1].phi.Ones();
}
/******************************************************//**
* Routine: LimitedLinearReconstructionOverDomain
*
* Peforms the reconstruction of a limited piecewise
* linear solution state within each cell of the
* computational mesh. The input parameters object specifies
* all the parameters necessary to perform the reconstruction
* (e.g. method, limiter etc.).
*
********************************************************/
template<class cState, class pState>
void LimitedLinearReconstructionOverDomain(RadMom1D_UniformMesh<cState, pState> *Soln, const RadMom1D_Input_Parameters<cState, pState> &IP){
switch(IP.i_Reconstruction) {
case RECONSTRUCTION_GREEN_GAUSS :
Linear_Reconstruction_GreenGauss(Soln,
IP.i_Limiter);
break;
case RECONSTRUCTION_LEAST_SQUARES :
Linear_Reconstruction_LeastSquares(Soln,
IP.i_Limiter);
break;
default:
throw runtime_error("LimitedLinearReconstructionOverDomain() ERROR: Unknown reconstruction type");
break;
} /* endswitch */
}
#endif // _RADMOM_1D_RECON_INCLUDED