PLUTO
Functions
mappers.c File Reference

Convert between primitive and conservative variables. More...

#include "pluto.h"
Include dependency graph for mappers.c:

Go to the source code of this file.

Functions

void PrimToCons (double *uprim[], double *ucons[], int ibeg, int iend)
 
int ConsToPrim (double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
 

Detailed Description

Convert between primitive and conservative variables.

The PrimToCons() converts an array of primitive quantities to an array of conservative variables for the RHD equations.

The ConsToPrim() converts an array of conservative quantities to an array of primitive quantities. During the conversion, pressure is normally recovered from total energy unless zone has been tagged with FLAG_ENTROPY. In this case we recover pressure from conserved entropy: 

if (FLAG_ENTROPY is TRUE)  --> p = p(S)
else                       --> p = p(E)
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
June 25, 2015

Definition in file mappers.c.

Function Documentation

int ConsToPrim ( double **  ucons,
double **  uprim,
int  ibeg,
int  iend,
unsigned char *  flag 
)

Convert from conservative to primitive variables.

Parameters
[in]uconsarray of conservative variables
[out]uprimarray of primitive variables
[in]begstarting index of computation
[in]endfinal index of computation
[in,out]flagarray of flags tagging, in input, zones where entropy must be used to recover pressure and, on output, zones where conversion was not successful.
Returns
Return 0 if conversion was successful in all zones in the range [ibeg,iend]. Return 1 if one or more zones could not be converted correctly and either pressure, density or energy took non-physical values.

Definition at line 102 of file mappers.c.

122 {
123  int nv, i, err, ifail;
124  int use_entropy, use_energy=1;
125  double scrh, m, g;
126  double *u, *v;
127  Map_param par;
128 
129  ifail = 0;
130  for (i = ibeg; i <= iend; i++) {
131 
132  u = ucons[i];
133  v = uprim[i];
134 
135  /* -----------------------------------------------------------
136  Define the input parameters of the parameter structure
137  ----------------------------------------------------------- */
138 
139  par.D = u[RHO];
140  par.E = u[ENG];
141  par.m2 = EXPAND(u[MX1]*u[MX1], + u[MX2]*u[MX2], + u[MX3]*u[MX3]);
142 
143  /* -------------------------------------------
144  Check density and energy positivity
145  ------------------------------------------- */
146 
147  if (u[RHO] < 0.0) {
148  print("! ConsToPrim: negative density (%8.2e), ", u[RHO]);
149  Where (i, NULL);
150  u[RHO] = g_smallDensity;
151  flag[i] |= FLAG_CONS2PRIM_FAIL;
152  ifail = 1;
153  }
154 
155  if (u[ENG] < 0.0) {
156  WARNING(
157  print("! ConsToPrim: negative energy (%8.2e), ", u[ENG]);
158  Where (i, NULL);
159  )
160  u[ENG] = sqrt(1.e-8 + par.m2 + u[RHO]*u[RHO]);
161  flag[i] |= FLAG_CONS2PRIM_FAIL;
162  ifail = 1;
163  }
164 
165  /* --------------------------------
166  recover pressure by inverting
167  the energy or entropy equation.
168  -------------------------------- */
169 
170 #if ENTROPY_SWITCH
171  use_entropy = (flag[i] & FLAG_ENTROPY);
172  use_energy = !use_entropy;
173  par.sigma_c = u[ENTR];
174  if (use_entropy) {
175  err = EntropySolve(&par);
176  if (err) {
177  WARNING(Where (i, NULL);)
178  err = PressureFix(&par);
179  if (err){
180  Where(i,NULL);
181  QUIT_PLUTO(1);
182  }
183  flag[i] |= FLAG_CONS2PRIM_FAIL;
184  ifail = 1;
185  }
186  u[ENG] = par.E; /* Redefine energy */
187  }
188 #endif
189 
190  if (use_energy){
191  err = EnergySolve(&par);
192  if (err){
193  WARNING(Where(i,NULL);)
194  err = PressureFix(&par);
195  if (err){
196  Where(i,NULL);
197  QUIT_PLUTO(1);
198  }
199  u[ENG] = par.E;
200  flag[i] |= FLAG_CONS2PRIM_FAIL;
201  ifail = 1;
202  }
203 #if ENTROPY_SWITCH
204  u[ENTR] = par.sigma_c; /* Redefine entropy */
205 #endif
206  }
207 
208  /* ------------------------------------------
209  complete conversion
210  ------------------------------------------ */
211 
212  v[PRS] = par.prs;
213  scrh = 1.0/(u[ENG] + v[PRS]); /* = 1 / W */
214 
215  EXPAND(v[VX1] = u[MX1]*scrh; ,
216  v[VX2] = u[MX2]*scrh; ,
217  v[VX3] = u[MX3]*scrh;)
218 
219  g = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
220  g = 1.0/sqrt(1.0 - g);
221  v[RHO] = u[RHO]/g;
222 
223 #if NSCL > 0
224  NSCL_LOOP(nv) v[nv] = u[nv]/u[RHO];
225 #endif
226 
227  }
228 
229  return ifail;
230 }
MX3
#define MX3
Definition: mod_defs.h:22
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
MX1
#define MX1
Definition: mod_defs.h:20
MAP_PARAM::m2
double m2
Square of total momentum (input).
Definition: mod_defs.h:89
VX2
#define VX2
Definition: mod_defs.h:29
MAP_PARAM::D
double D
Lab density (input).
Definition: mod_defs.h:86
g_smallDensity
double g_smallDensity
Small value for density fix.
Definition: globals.h:109
PressureFix
int PressureFix(Map_param *)
Definition: pressure_fix.c:22
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
WARNING
#define WARNING(a)
Definition: macros.h:217
func_param::v
double v[NVAR]
Definition: eos.h:106
EnergySolve
int EnergySolve(Map_param *par)
Definition: energy_solve.c:23
VX1
#define VX1
Definition: mod_defs.h:28
FLAG_CONS2PRIM_FAIL
#define FLAG_CONS2PRIM_FAIL
Definition: pluto.h:189
FLAG_ENTROPY
#define FLAG_ENTROPY
Update pressure using entropy equation.
Definition: pluto.h:182
EntropySolve
int EntropySolve(Map_param *par)
Definition: entropy_solve.c:21
MX2
#define MX2
Definition: mod_defs.h:21
Where
void Where(int, Grid *)
Definition: tools.c:235
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MAP_PARAM
Definition: mod_defs.h:85
i
int i
Definition: analysis.c:2
MAP_PARAM::prs
double prs
Thermal pressure (output).
Definition: mod_defs.h:94
VX3
#define VX3
Definition: mod_defs.h:30
MAP_PARAM::sigma_c
double sigma_c
Conserved entropy (input).
Definition: mod_defs.h:87
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
MAP_PARAM::E
double E
Total energy (input).
Definition: mod_defs.h:88

Here is the call graph for this function:

void PrimToCons ( double *  uprim[],
double *  ucons[],
int  ibeg,
int  iend 
)

Convert primitive variables to conservative variables.

Parameters
[in]uprimarray of primitive variables
[out]uconsarray of conservative variables
[in]begstarting index of computation
[in]endfinal index of computation

Definition at line 25 of file mappers.c.

36 {
37  int nv, ii;
38  double rhoh_g2, scrh, g;
39  double beta_fix = 0.9999;
40  double *u, *v;
41  static double *h;
42 
43  if (h == NULL){
44  h = ARRAY_1D(NMAX_POINT, double);
45  }
46 
47  Enthalpy (uprim, h, ibeg, iend);
48 
49  for (ii = ibeg; ii <= iend; ii++) {
50 
51  u = ucons[ii];
52  v = uprim[ii];
53 
54  g = EXPAND(v[VX1]*v[VX1], +v[VX2]*v[VX2], +v[VX3]*v[VX3]);
55 
56  if (g >= 1.0){
57  WARNING(
58  print ("! u^2 > 1 (%f) in PrimToCons\n", scrh);
59  Where (ii, NULL);
60  )
61 
62  g = beta_fix/sqrt(g);
63  EXPAND(v[VX1] *= g; ,
64  v[VX2] *= g; ,
65  v[VX3] *= g;)
66 
67  g = beta_fix*beta_fix;
68  }
69  g = 1.0/(1.0 - g);
70  scrh = rhoh_g2 = v[RHO]*h[ii]*g;
71  g = sqrt(g);
72 
73  u[RHO] = v[RHO]*g;
74  EXPAND(u[MX1] = scrh*v[VX1]; ,
75  u[MX2] = scrh*v[VX2]; ,
76  u[MX3] = scrh*v[VX3];)
77 
78  ucons[ii][ENG] = rhoh_g2 - v[PRS];
79 
80 #if NSCL > 0
81  NSCL_LOOP(nv) u[nv] = v[nv]*u[RHO];
82 #endif
83 
84  #if CHECK_CONSERVATIVE_VAR == YES
85  m2 = EXPAND(u[MX1]*u[MX1], + u[MX2]*u[MX2], + u[MX3]*u[MX3]);
86  g = u[ENG] - sqrt(m2);
87  if (g <= 0.0){
88  printf ("! E - m < 0 in PrimToCons (%12.6e)\n",g);
89  QUIT_PLUTO(1);
90  }
91  g = u[RHO]/u[ENG] - sqrt(1.0 - m2/u[ENG]/u[ENG]);
92  if (g >=0){
93  print ("! g > 0.0 in PrimToCons(2) (%12.6e)\n",g);
94  Show(ucons,ii);
95  QUIT_PLUTO(1);
96  }
97  #endif
98  }
99 }
MX3
#define MX3
Definition: mod_defs.h:22
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
MX1
#define MX1
Definition: mod_defs.h:20
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
WARNING
#define WARNING(a)
Definition: macros.h:217
Enthalpy
void Enthalpy(double **v, real *h, int beg, int end)
Definition: eos.c:48
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
MX2
#define MX2
Definition: mod_defs.h:21
Where
void Where(int, Grid *)
Definition: tools.c:235
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
Show
void Show(double **, int)
Definition: tools.c:132
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
VX3
#define VX3
Definition: mod_defs.h:30
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function: