Carmen Code
 All Classes Files Functions Variables Macros Pages
Functions
IntermediaryStates.cpp File Reference
#include "Carmen.h"
Include dependency graph for IntermediaryStates.cpp:

Functions

Matrix stateUstar (const Vector &AvgL, const Vector &AvgR, const real prel, const real prer, real &slopeLeft, real &slopeRight, real &slopeM, real &slopeLeftStar, real &slopeRightStar, int AxisNo)
 Returns the intermediary states of HLLD numerical flux for MHD equations. More...
 

Function Documentation

Matrix stateUstar ( const Vector &  AvgL,
const Vector &  AvgR,
const real  prel,
const real  prer,
real &  slopeLeft,
real &  slopeRight,
real &  slopeM,
real &  slopeLeftStar,
real &  slopeRightStar,
int  AxisNo 
)

Returns the intermediary states of HLLD numerical flux for MHD equations.

Parameters
AvgLLeft average vector
AvgRRight average vector
prelLeft pressure
prerRight pressure
slopeLeftLeft slope
slopeRightRight slope
slopeMSlope value computed for HLLD
slopeLeftStarLeft star slope (intermediary states)
slopeRightStarRight star slope (intermediary states)
AxisNoAxis of interest
Returns
Matrix

variables U-star

variables U-star-star

13 {
14  real rhol=0.,rhor=0.;
15  real psil=0.,psir=0.;
16  real vxl=0.,vxr=0.,vyl=0.,vyr=0.,vzl=0.,vzr=0.;
17  real Bxr=0.,Bxl=0.,Byl=0.,Byr=0.,Bzl=0.,Bzr=0.;
18  real vl=0.,vr=0.,mf=0.;
19  real Bl=0.,Br=0.;
20  real pTl=0.,pTr=0.;
21  real vBl=0.,vBr=0.;
22  real el=0., er=0.;
23  real rhols=0.,rhors=0.;
24  real vxls=0.,vxrs=0.,vyls=0.,vyrs=0.,vzls=0.,vzrs=0.;
25  real Bxls=0.,Bxrs=0.,Byls=0.,Byrs=0.,Bzls=0.,Bzrs=0.;
26  real pTs=0.;
27  real vBls=0.,vBrs=0.;
28  real Els=0.,Ers=0.;
29  real rholss=0.,rhorss=0.;
30  real vxlss=0.,vxrss=0.,vylss=0.,vyrss=0.,vzlss=0.,vzrss=0.;
31  real Bxlss=0.,Bxrss=0.,Bylss=0.,Byrss=0.,Bzlss=0.,Bzrss=0.;
32  real vBlss=0.,vBrss=0.;
33  real Elss=0.,Erss=0.;
34  int Bsign=0;
35  real half=0.5;
36  real epsilon2=1e-16;
37  real auxl=0.,auxr=0., Sqrtrhols=0.;
38 
39  Matrix U(QuantityNb,4);
40 
41  //Variables
42  rhol = AvgL.value(1);
43  vxl = AvgL.value(2)/rhol;
44  vyl = AvgL.value(3)/rhol;
45  vzl = AvgL.value(4)/rhol;
46  el = AvgL.value(5);
47  psil = AvgL.value(6);
48  Bxl = AvgL.value(7);
49  Byl = AvgL.value(8);
50  Bzl = AvgL.value(9);
51 
52  rhor = AvgR.value(1);
53  vxr = AvgR.value(2)/rhor;
54  vyr = AvgR.value(3)/rhor;
55  vzr = AvgR.value(4)/rhor;
56  er = AvgR.value(5);
57  psir = AvgR.value(6);
58  Bxr = AvgR.value(7);
59  Byr = AvgR.value(8);
60  Bzr = AvgR.value(9);
61 
62  //v_x and v_y
63  vl = AvgL.value(AxisNo+1)/rhol;
64  vr = AvgR.value(AxisNo+1)/rhor;
65 
66  // Average B value
67  mf = (AvgL.value(AxisNo+6)+AvgR.value(AxisNo+6))/(double (2.0));
68 
69  if(AxisNo==1)
70  {
71  //|B|
72  Bl = sqrt( mf*mf + Byl*Byl + Bzl*Bzl );
73  Br = sqrt( mf*mf + Byr*Byr + Bzr*Bzr );
74  //Inner product v.B
75  vBl = vxl*mf + vyl*Byl + vzl*Bzl;
76  vBr = vxr*mf + vyr*Byr + vzr*Bzr;
77  }else if(AxisNo==2){
78  //|B|
79  Bl = sqrt( Bxl*Bxl + mf*mf + Bzl*Bzl );
80  Br = sqrt( Bxr*Bxr + mf*mf + Bzr*Bzr );
81  //Inner product v.B
82  vBl = vxl*Bxl + vyl*mf + vzl*Bzl;
83  vBr = vxr*Bxr + vyr*mf + vzr*Bzr;
84  }else{
85  //|B|
86  Bl = sqrt( Bxl*Bxl + Byl*Byl + mf*mf );
87  Br = sqrt( Bxr*Bxr + Byr*Byr + mf*mf );
88  //Inner product v.B
89  vBl = vxl*Bxl + vyl*Byl + vzl*mf;
90  vBr = vxr*Bxr + vyr*Byr + vzr*mf;
91  }
92 
93  //Total Pressure
94  pTl = prel + half*Bl*Bl;
95  pTr = prer + half*Br*Br;
96 
97  // S_M - Equation 38
98  slopeM = ((slopeRight - vr)*rhor*vr - (slopeLeft - vl)*rhol*vl - pTr +pTl)/
99  ((slopeRight - vr)*rhor - (slopeLeft - vl)*rhol);
100 
101  //Sign function
102  if(mf > 0) Bsign = 1;
103  else Bsign = -1;
104 
106 
107  //density - Equation 43
108  rhols = rhol*(slopeLeft-vl)/(slopeLeft-slopeM);
109  rhors = rhor*(slopeRight-vr)/(slopeRight-slopeM);
110 
111  if(AxisNo==1){
112  //velocities
113  //Equation 39
114  vxls = slopeM;
115  vxrs = vxls;
116  //B_x
117  Bxls = mf;
118  Bxrs = Bxls;
119 
120  auxl= (rhol*(slopeLeft - vl)*(slopeLeft - slopeM)-mf*mf);
121  auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);
122 
123  if( fabs(auxl)<epsilon2 ||
124  ( ( ( fabs(slopeM -vl) ) < epsilon2 ) &&
125  ( ( fabs(Byl) + fabs(Bzl) ) < epsilon2 ) &&
126  ( ( mf*mf ) > (Gamma*prel) ) ) ) {
127  vyls = vyl;
128  vzls = vzl;
129  Byls = Byl;
130  Bzls = Bzl;
131  }else{
132  //Equation 44
133  vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
134  //Equation 46
135  vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
136  //Equation 45
137  Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
138  //Equation 47
139  Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
140  }
141  if( fabs(auxr)<epsilon2 ||
142  ( ( ( fabs(slopeM -vr) ) < epsilon2 ) &&
143  ( ( fabs(Byr) + fabs(Bzr) ) < epsilon2 ) &&
144  ( ( mf*mf ) > (Gamma*prer) ) ) ) {
145  vyrs = vyr;
146  vzrs = vzr;
147  Byrs = Byr;
148  Bzrs = Bzr;
149  }else{
150  //Equation 44
151  vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
152  //Equation 46
153  vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
154  //Equation 45;
155  Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
156  //Equation 47
157  Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
158  }
159 
160  }else if(AxisNo==2){
161  //velocities
162  //Equation 39
163  vyls = slopeM;
164  vyrs = vyls;
165  //B_y
166  Byls = mf;
167  Byrs = Byls;
168 
169  auxl= (rhol*(slopeLeft - vl)*(slopeLeft - slopeM)-mf*mf);
170  auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);
171 
172  if( fabs(auxl)<epsilon2 ||
173  ( ( ( fabs(slopeM -vl) ) < epsilon2 ) &&
174  ( ( fabs(Bxl) + fabs(Bzl) ) < epsilon2 ) &&
175  ( ( mf*mf ) > (Gamma*prel) ) ) ) {
176  vxls = vxl;
177  vzls = vzl;
178  Bxls = Bxl;
179  Bzls = Bzl;
180  }else{
181  //Equation 44
182  vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
183  //Equation 46
184  vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
185  //Equation 45
186  Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
187  //Equation 47
188  Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
189  }
190 
191  if( fabs(auxr)<epsilon2 ||
192  ( ( ( fabs(slopeM -vr) ) < epsilon2 ) &&
193  ( ( fabs(Bxr) + fabs(Bzr) ) < epsilon2 ) &&
194  ( ( mf*mf ) > (Gamma*prer) ) ) ) {
195  vxrs = vxr;
196  vzrs = vzr;
197  Bxrs = Bxr;
198  Bzrs = Bzr;
199  }else{
200  //Equation 44
201  vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
202  //Equation 46
203  vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
204  //Equation 45
205  Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
206  //Equation 47
207  Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
208  }
209  }else{
210  //velocities
211  //Equation 39
212  vzls = slopeM;
213  vzrs = vzls;
214  //B_z
215  Bzls = mf;
216  Bzrs = Bzls;
217 
218  auxl= (rhol*(slopeLeft - vl)*(slopeLeft - slopeM)-mf*mf);
219  auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);
220 
221  if( fabs(auxl)<epsilon2 ||
222  ( ( ( fabs(slopeM -vl) ) < epsilon2 ) &&
223  ( ( fabs(Bxl) + fabs(Byl) ) < epsilon2 ) &&
224  ( ( mf*mf ) > (Gamma*prel) ) ) ) {
225  vxls = vxl;
226  vyls = vyl;
227  Bxls = Bxl;
228  Byls = Byl;
229  }else{
230  //Equation 44
231  vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
232  //Equation 46
233  vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
234  //Equation 45
235  Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
236  //Equation 47
237  Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
238  }
239 
240  if( fabs(auxr)<epsilon2 ||
241  ( ( ( fabs(slopeM -vr) ) < epsilon2 ) &&
242  ( ( fabs(Bxr) + fabs(Byr) ) < epsilon2 ) &&
243  ( ( mf*mf ) > (Gamma*prer) ) ) ) {
244  vxrs = vxr;
245  vyrs = vyr;
246  Bxrs = Bxr;
247  Byrs = Byr;
248  }else{
249  //Equation 44
250  vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
251  //Equation 46
252  vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
253  //Equation 45
254  Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
255  //Equation 47
256  Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
257  }
258  }
259 
260  //total pressure - Equation 41
261  pTs = pTl + rhol*(slopeLeft - vl)*(slopeM - vl);
262 
263  //inner product vs*Bs
264  vBls = vxls*Bxls + vyls*Byls + vzls*Bzls;
265  vBrs = vxrs*Bxrs + vyrs*Byrs + vzrs*Bzrs;
266 
267  //energy - Equation 48
268  Els = ((slopeLeft -vl)*el - pTl*vl+pTs*slopeM+mf*(vBl - vBls))/(slopeLeft - slopeM);
269  Ers = ((slopeRight-vr)*er - pTr*vr+pTs*slopeM+mf*(vBr - vBrs))/(slopeRight - slopeM);
270 
271  //U-star variables
272  //Left
273  U.setValue(1,1,rhols);
274  U.setValue(2,1,rhols*vxls);
275  U.setValue(3,1,rhols*vyls);
276  U.setValue(4,1,rhols*vzls);
277  U.setValue(5,1,Els);
278  U.setValue(6,1,psil);
279  U.setValue(7,1,Bxls);
280  U.setValue(8,1,Byls);
281  U.setValue(9,1,Bzls);
282  //Right
283  U.setValue(1,2,rhors);
284  U.setValue(2,2,rhors*vxrs);
285  U.setValue(3,2,rhors*vyrs);
286  U.setValue(4,2,rhors*vzrs);
287  U.setValue(5,2,Ers);
288  U.setValue(6,2,psir);
289  U.setValue(7,2,Bxrs);
290  U.setValue(8,2,Byrs);
291  U.setValue(9,2,Bzrs);
292 
293 
295 
296 
297  Sqrtrhols = (sqrt(rhols)+sqrt(rhors));
298 
299  if((mf*mf/min((Bl*Bl),(Br*Br))) < epsilon2){
300  for(int k=1;k<=QuantityNb;k++){
301  U.setValue(k,3, U.value(k,1));
302  U.setValue(k,4, U.value(k,2));
303  }
304  }else{
305  //density - Equation 49
306  rholss = rhols;
307  rhorss = rhors;
308 
309  if(AxisNo==1){
310  //velocities
311  //Equation 39
312  vxlss = slopeM;
313  vxrss = slopeM;
314  //Equation 59
315  vylss = (sqrt(rhols)*vyls + sqrt(rhors)*vyrs + (Byrs - Byls)*Bsign)/Sqrtrhols;
316  vyrss = vylss;
317  //Equation 60
318  vzlss = (sqrt(rhols)*vzls + sqrt(rhors)*vzrs + (Bzrs - Bzls)*Bsign)/Sqrtrhols;
319  vzrss = vzlss;
320 
321  //Magnetic Field components
322  //B_x
323  Bxlss = mf;
324  Bxrss = mf;
325  //Equation 61
326  Bylss = (sqrt(rhols)*Byrs + sqrt(rhors)*Byls + sqrt(rhols*rhors)*(vyrs - vyls)*Bsign)/Sqrtrhols;
327  Byrss = Bylss;
328  //Equation 62
329  Bzlss = (sqrt(rhols)*Bzrs + sqrt(rhors)*Bzls + sqrt(rhols*rhors)*(vzrs - vzls)*Bsign)/Sqrtrhols;
330  Bzrss = Bzlss;
331  }
332  else if(AxisNo==2){
333  //velocities
334  //Equation 59
335  vxlss = (sqrt(rhols)*vxls + sqrt(rhors)*vxrs + (Bxrs - Bxls)*Bsign)/Sqrtrhols;
336  vxrss = vxlss;
337  //Equation 39
338  vylss = slopeM;
339  vyrss = slopeM;
340  //Equation 60
341  vzlss = (sqrt(rhols)*vzls + sqrt(rhors)*vzrs + (Bzrs - Bzls)*Bsign)/Sqrtrhols;
342  vzrss = vzlss;
343 
344  //Magnetic Field components
345  //Equation 61
346  Bxlss = (sqrt(rhols)*Bxrs + sqrt(rhors)*Bxls + sqrt(rhols*rhors)*(vxrs - vxls)*Bsign)/Sqrtrhols;
347  Bxrss = Bxlss;
348  //B_y
349  Bylss = mf;
350  Byrss = mf;
351  //Equation 62
352  Bzlss = (sqrt(rhols)*Bzrs + sqrt(rhors)*Bzls + sqrt(rhols*rhors)*(vzrs - vzls)*Bsign)/Sqrtrhols;
353  Bzrss = Bzlss;
354 
355  }else{
356  //velocities
357  //Equation 59
358  vxlss = (sqrt(rhols)*vxls + sqrt(rhors)*vxrs + (Bxrs - Bxls)*Bsign)/Sqrtrhols;
359  vxrss = vxlss;
360  //Equation 60
361  vylss = (sqrt(rhols)*vyls + sqrt(rhors)*vyrs + (Byrs - Byls)*Bsign)/Sqrtrhols;
362  vyrss = vylss;
363  //Equation 39
364  vzlss = slopeM;
365  vzrss = slopeM;
366 
367  //Magnetic Field components
368  //Equation 61
369  Bxlss = (sqrt(rhols)*Bxrs + sqrt(rhors)*Bxls + sqrt(rhols*rhors)*(vxrs - vxls)*Bsign)/Sqrtrhols;
370  Bxrss = Bxlss;
371  //Equation 62
372  Bylss = (sqrt(rhols)*Byrs + sqrt(rhors)*Byls + sqrt(rhols*rhors)*(vyrs - vyls)*Bsign)/Sqrtrhols;
373  Byrss = Bylss;
374  //B_y
375  Bzlss = mf;
376  Bzrss = mf;
377 
378  }
379 
380 
381  //inner product vss*Bss
382  vBlss = vxlss*Bxlss + vylss*Bylss + vzlss*Bzlss;
383  vBrss = vxrss*Bxrss + vyrss*Byrss + vzrss*Bzrss;
384 
385  //Energy - Equation 63
386  Elss = Els - sqrt(rhols)*(vBls - vBlss)*Bsign;
387  Erss = Ers + sqrt(rhors)*(vBrs - vBrss)*Bsign;
388 
389 
390 
391 
392  //U-star-star variables
393  //Left
394  U.setValue(1,3,rholss);
395  U.setValue(2,3,rholss*vxlss);
396  U.setValue(3,3,rholss*vylss);
397  U.setValue(4,3,rholss*vzlss);
398  U.setValue(5,3,Elss);
399  U.setValue(6,3,psil);
400  U.setValue(7,3,Bxlss);
401  U.setValue(8,3,Bylss);
402  U.setValue(9,3,Bzlss);
403  //Right
404  U.setValue(1,4,rhorss);
405  U.setValue(2,4,rhorss*vxrss);
406  U.setValue(3,4,rhorss*vyrss);
407  U.setValue(4,4,rhorss*vzrss);
408  U.setValue(5,4,Erss);
409  U.setValue(6,4,psir);
410  U.setValue(7,4,Bxrss);
411  U.setValue(8,4,Byrss);
412  U.setValue(9,4,Bzrss);
413  }
414  //Equation 61 - S-star left and S-star right
415  slopeLeftStar = slopeM - Abs(mf)/sqrt(rhols);
416  slopeRightStar = slopeM + Abs(mf)/sqrt(rhors);
417 
418  return U;
419 }
QuantityNb
int QuantityNb
Definition: Parameters.cpp:171
Gamma
real Gamma
Definition: Parameters.cpp:109
Matrix
Standard class for every matrix in Carmen.
Definition: Matrix.h:28
Abs
#define Abs(x)
Definition: Carmen.h:78
Vector::value
real value(const int n) const
Returns the value of the component n.
Definition: Vector.h:565
real
#define real
Definition: PreProcessor.h:31

Here is the caller graph for this function:

Generated on Wed Sep 6 2017 14:34:41 for Carmen Code by   doxygen 1.8.6