Carmen Code
 All Classes Files Functions Variables Macros Pages
Functions
SchemeHLLD.cpp File Reference

Computes the HLLD Riemann solver. More...

#include "Carmen.h"
Include dependency graph for SchemeHLLD.cpp:

Functions

Vector SchemeHLLD (const Cell &Cell1, const Cell &Cell2, const Cell &Cell3, const Cell &Cell4, int AxisNo)
 Returns the HLLD numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides. More...
 

Detailed Description

Computes the HLLD Riemann solver.

Author
Anna Karina Fontes Gomes
Version
4.0
Date
July-2016

Function Documentation

Vector SchemeHLLD ( const Cell &  Cell1,
const Cell &  Cell2,
const Cell &  Cell3,
const Cell &  Cell4,
const int  AxisNo 
)

Returns the HLLD numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides.

Parameters
Cell1second neighbour on the left side
Cell2first neighbour on the left side
Cell3first neighbour on the right side
Cell4second neighbour on the right side
AxisNoAxis of interest.
Returns
Vector
13 {
14 
15  // General variables
16 
17  Vector LeftAverage(QuantityNb); //
18  Vector RightAverage(QuantityNb); // Conservative quantities
19  Vector Result(QuantityNb); // MHD numerical flux
20  int aux=0;
21 
22  // Variables for the HLL scheme
23  Vector FL(QuantityNb), FR(QuantityNb); //Left and right physical fluxes
24  Vector VL(3), VR(3); // Left and right velocities
25  Vector BL(3), BR(3); // Left and right velocities
26  real rhoL=0., rhoR=0.; // Left and right densities
27  real eL=0., eR=0.; // Left and right energies
28  real preL=0., preR=0.;
29  real bkL=0., bkR=0.;
30  real aL=0., aR=0.;
31  real bL=0., bR=0.;
32  real cfL=0., cfR=0.;
33  real SL=0., SR=0.;
34  real SLS=0., SRS=0.;
35  real SM=0.;
36  Matrix U(QuantityNb,4);
37  Matrix F(QuantityNb,2);
38  real dx=0.;
39  dx = Cell2.size(AxisNo);
40  real r, Limit, LeftSlope = 0., RightSlope = 0.; // Left and right slopes
41  int i;
42 
43 // --- Limiter function ---------------------------------------------------------
44 
45  for (i=1; i<=QuantityNb; i++)
46  {
47  // --- Compute left cell-average value ---
48 
49  if (Cell2.average(i) != Cell1.average(i))
50  {
51  RightSlope = Cell3.average(i)-Cell2.average(i);
52  LeftSlope = Cell2.average(i)-Cell1.average(i);
53  r = RightSlope/LeftSlope;
54  Limit = Limiter(r);
55  LeftAverage.setValue(i, Cell2.average(i) + 0.5*Limit*LeftSlope);
56  aux = 1;
57  }
58  else
59  LeftAverage.setValue(i, Cell2.average(i));
60 
61  // --- Compute right cell-average value ---
62 
63  if (Cell3.average(i) != Cell2.average(i))
64  {
65  RightSlope = Cell4.average(i)-Cell3.average(i);
66  LeftSlope = Cell3.average(i)-Cell2.average(i);
67  r = RightSlope/LeftSlope;
68  Limit = Limiter(r);
69  RightAverage.setValue(i, Cell3.average(i) - 0.5*Limit*LeftSlope);
70  aux = 1;
71  }
72  else
73  RightAverage.setValue(i, Cell3.average(i));
74  }
75 
76  // --- HLLD scheme -------------------------------------------------------------
77 
78  // --- Conservative variables ---
79 
80  // Left and right densities
81  rhoL = LeftAverage.value(1);
82  rhoR = RightAverage.value(1);
83 
84  // Left and right momentum and magnetic field
85  for (int i=1;i<=3;i++)
86  {
87  VL.setValue( i, LeftAverage.value (i+1));
88  VR.setValue( i, RightAverage.value(i+1));
89  BL.setValue( i, LeftAverage.value (i+6));
90  BR.setValue( i, RightAverage.value(i+6));
91  }
92 
93  // Left and right energies
94  eL = LeftAverage.value(5);
95  eR = RightAverage.value(5);
96 
97  // Left and right pressures
98  preL = (Gamma -1.)*(eL - 0.5*(VL*VL)/rhoL - 0.5*(BL*BL));
99  preR = (Gamma -1.)*(eR - 0.5*(VR*VR)/rhoR - 0.5*(BR*BR));
100 
101  // --- Magnetoacoustic waves computation --
102  bkL = power2(BL.value(AxisNo))/rhoL;
103  bkR = power2(BR.value(AxisNo))/rhoR;
104 
105  aL = Gamma*preL/rhoL;
106  aR = Gamma*preR/rhoR;
107 
108  bL = (BL*BL)/rhoL;
109  bR = (BR*BR)/rhoR;
110 
111  // Left and Right fast speeds
112  cfL = sqrt(0.5*(aL + bL + sqrt(power2(aL + bL) - 4.0*aL*bkL)));
113  cfR = sqrt(0.5*(aR + bR + sqrt(power2(aR + bR) - 4.0*aR*bkR)));
114 
115  // Left and Right slopes
116  SL = Min((VL.value(AxisNo))/rhoL, (VR.value(AxisNo))/rhoR) - Max(cfL,cfR);
117  SR = Max((VL.value(AxisNo))/rhoL, (VR.value(AxisNo))/rhoR) + Max(cfL,cfR);
118 
119  // --- Physical flux ---
120  if(AxisNo ==1){
121  EigenvalueX = Max(Max(Abs(SL),Abs(SR)),EigenvalueX);
122  FL = FluxX(LeftAverage);
123  FR = FluxX(RightAverage);
124  }else if(AxisNo ==2){
125  EigenvalueY = Max(Max(Abs(SL),Abs(SR)),EigenvalueY);
126  FL = FluxY(LeftAverage);
127  FR = FluxY(RightAverage);
128  }else{
129  EigenvalueZ = Max(Max(Abs(SL),Abs(SR)),EigenvalueZ);
130  FL = FluxZ(LeftAverage);
131  FR = FluxZ(RightAverage);
132  }
133 
134  // Intermediary states U* and U**
135  U = stateUstar(LeftAverage, RightAverage, preL, preR, SL, SR, SM, SLS, SRS,AxisNo);
136 
137  // --- HLLD Riemann Solver ---
138 
139  for(int i=1;i<=QuantityNb;i++)
140  {
141  //Flux Function - Equation 66
142  //F_L
143  if(SL>=0.)
144  Result.setValue(i, FL.value(i));
145  //F-star left // FL=FLstar
146  else if(SLS>=0. && SL<0.)
147  Result.setValue(i, FL.value(i) + SL*(U.value(i,1) - LeftAverage.value(i)));
148  //F-star-star left
149  else if(SM>=0. && SLS<0.)
150  Result.setValue(i, FL.value(i) + SLS*U.value(i,3) - (SLS - SL)*U.value(i,1) - SL*LeftAverage.value(i));
151  //F-star-star right
152  else if(SRS>=0. && SM<0.)
153  Result.setValue(i, FR.value(i) + SRS*U.value(i,4) - (SRS - SR)*U.value(i,2) - SR*RightAverage.value(i));
154  //F-star right
155  else if(SR>=0. && SRS<0.)
156  Result.setValue(i, FR.value(i) + SR*(U.value(i,2) - RightAverage.value(i)));
157  //F_R
158  else
159  Result.setValue(i, FR.value(i));
160  }
161  // Parabolic-Hyperbolic divergence Cleaning (Dedner, 2002)
162  //fluxCorrection(Result, Cell2.average(), Cell3.average(), AxisNo);
163  fluxCorrection(Result, LeftAverage, RightAverage, AxisNo);
164 
165  // Artificial diffusion terms
166  if(Diffusivity && aux==1) Result = Result - ArtificialViscosity(LeftAverage,RightAverage,dx,AxisNo);
167 
168  return Result;
169 
170 }
power2
#define power2(x)
Definition: Carmen.h:70
EigenvalueZ
real EigenvalueZ
Definition: Parameters.cpp:160
FluxY
Vector FluxY(const Vector &Avg)
Returns the physical flux of MHD equations in Y direction.
Definition: PhysicalFluxMHD.cpp:58
QuantityNb
int QuantityNb
Definition: Parameters.cpp:171
Vector
Standard class for every vector in Carmen.
Definition: Vector.h:29
ArtificialViscosity
Vector ArtificialViscosity(const Vector &Cell1, const Vector &Cell2, real dx, int AxisNo)
Returns the artificial diffusion source terms in the cell UserCell.
Definition: ArtificialViscosity.cpp:11
Min
#define Min(x, y)
Definition: Carmen.h:62
Gamma
real Gamma
Definition: Parameters.cpp:109
Limiter
Vector Limiter(const Vector u, const Vector v)
Returns the value of the slope limiter between the slopes u and v.
Definition: Limiter.cpp:56
Diffusivity
bool Diffusivity
Definition: Parameters.cpp:121
FluxX
Vector FluxX(const Vector &Avg)
Returns the physical flux of MHD equations in X direction.
Definition: PhysicalFluxMHD.cpp:8
Cell::size
real size(const int AxisNo) const
Returns the cell size in the direction AxisNo.
Definition: Cell.h:1095
fluxCorrection
void fluxCorrection(Vector &Flux, const Vector &AvgL, const Vector &AvgR, int AxisNo)
This function apply the divergence-free correction to the numerical flux.
Definition: FluxCorrection.cpp:9
Cell::average
real average(const int QuantityNo) const
Returns the component no. QuantityNo of the cell-average value.
Definition: Cell.h:1128
Matrix
Standard class for every matrix in Carmen.
Definition: Matrix.h:28
FluxZ
Vector FluxZ(const Vector &Avg)
Returns the physical flux of MHD equations in Z direction.
Definition: PhysicalFluxMHD.cpp:106
Abs
#define Abs(x)
Definition: Carmen.h:78
EigenvalueY
real EigenvalueY
Definition: Parameters.cpp:159
EigenvalueX
real EigenvalueX
Definition: Parameters.cpp:158
Max
#define Max(x, y)
Definition: Carmen.h:54
stateUstar
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.
Definition: IntermediaryStates.cpp:12
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