PLUTO
Functions
source.c File Reference

Compute source terms for Powell formulation. More...

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

Go to the source code of this file.

Functions

void Roe_DivBSource (const State_1D *state, int is, int ie, Grid *grid)
 
void HLL_DivBSource (const State_1D *state, double **Uhll, int beg, int end, Grid *grid)
 

Detailed Description

Compute source terms for Powell formulation.

This file contains two implementations of Powell's source term in the 8-wave formulation for MHD:

\[ Q\nabla\cdot\vec{B} \approx Q_i\, \frac{A_{i\HALF}B_{i+\HALF} - A_{i-\HALF}B_{i-\HALF}}{\Delta{\cal V}_i} \]

where B is the magnetic field in the normal direction (g_dir), Q is a fluid quantity, A is the area and the denominator is the cell volume. The first function, Roe_DivBSource() is called by Roe_Solver() and TVDLF_Solver() and computes the normal component using arithmetic average of the left and right states: $ B_{i+\HALF} = (B^+_i + B^-_{i+1})/2 $. The second implementation contained in HLL_DivBSource() computes the term using upwinding:

\[ B_{i+\HALF} = \left\{\begin{array}{ll} B^+_i & \qquad{\rm if}\quad F_{\rho,i+\HALF} > 0\\ \noalign{\medskip} B^-_{i+1} & \qquad{\rm otherwise} \end{array}\right. \]

Reference:

  • "A positive conservative method for MHD based on HLL and Roe methods" P. Janhunen, JCP (2000), 160, 649.
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
June 8, 2007

Definition in file source.c.

Function Documentation

void HLL_DivBSource ( const State_1D *  state,
double **  Uhll,
int  beg,
int  end,
Grid *  grid 
)

Include div.B source term to momentum, induction and energy equation. Used in conjunction with an HLL-type Riemann solver.

Definition at line 179 of file source.c.

186 {
187  int i, nv;
188  double vc[NVAR], *A, *src, *vm;
189  double r, s, vB;
190  static double *divB, *vp;
191  Grid *GG;
192 
193  if (divB == NULL){
194  divB = ARRAY_1D(NMAX_POINT, double);
195  vp = ARRAY_1D(NMAX_POINT, double);
196  }
197 
198  GG = grid + g_dir;
199  vm = vp - 1;
200  A = grid[g_dir].A;
201 
202 /* -------------------------------------------------------
203  Compute normal component of the field using upwinding
204  ------------------------------------------------------- */
205 
206  for (i = beg - 1; i <= end; i++) {
207  vp[i] = state->flux[i][RHO] < 0.0 ? state->vR[i][BXn]: state->vL[i][BXn];
208  }
209 
210 /* --------------------------------------------
211  Compute div.B contribution from the normal
212  direction (1) in different geometries
213  -------------------------------------------- */
214 
215  #if GEOMETRY == CARTESIAN
216 
217  for (i = beg; i <= end; i++) {
218  divB[i] = (vp[i] - vm[i])/GG->dx[i];
219  }
220 
221  #elif GEOMETRY == CYLINDRICAL
222 
223  if (g_dir == IDIR){ /* -- r -- */
224  for (i = beg; i <= end; i++) {
225  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
226  }
227  }else if (g_dir == JDIR){ /* -- z -- */
228  for (i = beg; i <= end; i++) {
229  divB[i] = (vp[i] - vm[i])/GG->dx[i];
230  }
231  }
232 
233  #elif GEOMETRY == POLAR
234 
235  if (g_dir == IDIR){ /* -- r -- */
236  for (i = beg; i <= end; i++) {
237  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
238  }
239  }else if (g_dir == JDIR){ /* -- phi -- */
240  r = grid[IDIR].x[g_i];
241  for (i = beg; i <= end; i++) {
242  divB[i] = (vp[i] - vm[i])/(r*GG->dx[i]);
243  }
244  }else if (g_dir == KDIR){ /* -- z -- */
245  for (i = beg; i <= end; i++) {
246  divB[i] = (vp[i] - vm[i])/GG->dx[i];
247  }
248  }
249 
250  #elif GEOMETRY == SPHERICAL
251 
252  if (g_dir == IDIR){ /* -- r -- */
253  for (i = beg; i <= end; i++) {
254  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
255  }
256  }else if (g_dir == JDIR){ /* -- theta -- */
257  r = grid[IDIR].x[g_i];
258  for (i = beg; i <= end; i++) {
259  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/(r*GG->dV[i]);
260  }
261  }else if (g_dir == KDIR){ /* -- phi -- */
262  r = grid[IDIR].x[g_i];
263  s = sin(grid[JDIR].x[g_j]);
264  for (i = beg; i <= end; i++) {
265  divB[i] = (vp[i] - vm[i])/(r*s*GG->dx[i]);
266  }
267  }
268 
269  #endif
270 
271  /* -----------------------------------------
272  compute total source terms
273  ----------------------------------------- */
274 
275  for (i = beg; i <= end; i++) {
276 
277  src = state->src[i];
278  for (nv = NFLX; nv--; ) vc[nv] = state->vh[i][nv];
279 
280  vB = EXPAND(vc[VX1]*vc[BX1], + vc[VX2]*vc[BX2], + vc[VX3]*vc[BX3]);
281  src[RHO] = 0.0;
282  EXPAND(src[MX1] = -vc[BX1]*divB[i]; ,
283  src[MX2] = -vc[BX2]*divB[i]; ,
284  src[MX3] = -vc[BX3]*divB[i];)
285 
286  EXPAND(src[BX1] = -vc[VX1]*divB[i]; ,
287  src[BX2] = -vc[VX2]*divB[i]; ,
288  src[BX3] = -vc[VX3]*divB[i];)
289 
290  #if HAVE_ENERGY
291  src[ENG] = -vB*divB[i];
292  #endif
293 
294  }
295 }
MX3
#define MX3
Definition: mod_defs.h:22
MX1
#define MX1
Definition: mod_defs.h:20
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
VX2
#define VX2
Definition: mod_defs.h:29
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
RHO
#define RHO
Definition: mod_defs.h:19
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
GRID::dx
double * dx
Definition: structs.h:83
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
BXn
int BXn
Definition: globals.h:75
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
STATE_1D::src
double ** src
Definition: structs.h:160
NFLX
#define NFLX
Definition: mod_defs.h:32
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
GRID
Definition: structs.h:78
MX2
#define MX2
Definition: mod_defs.h:21
divB
double *** divB
Definition: analysis.c:4
GRID::x
double * x
Definition: structs.h:80
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395

Here is the caller graph for this function:

void Roe_DivBSource ( const State_1D *  state,
int  is,
int  ie,
Grid *  grid 
)

Include Powell div.B source term to momentum, induction and energy equation for Roe and TVDLF solvers.

Definition at line 39 of file source.c.

45 {
46  int i;
47  double btx, bty, btz, bx, by, bz, vx, vy, vz;
48  double r, s;
49  double *Ar, *Ath;
50  double *vm, **bgf;
51  double *src, *v;
52  static double *divB, *vp;
53  Grid *GG;
54 
55  if (divB == NULL){
56  divB = ARRAY_1D(NMAX_POINT, double);
57  vp = ARRAY_1D(NMAX_POINT, double);
58  }
59 
60 /* ----------------------- ---------------------
61  compute magnetic field normal component
62  interface value by arithmetic averaging
63  -------------------------------------------- */
64 
65  for (i = is - 1; i <= ie; i++) {
66  vp[i] = 0.5*(state->vL[i][BXn] + state->vR[i][BXn]);
67  }
68  vm = vp - 1;
69  GG = grid + g_dir;
70 
71 /* --------------------------------------------
72  Compute div.B contribution from the normal
73  direction (1) in different geometries
74  -------------------------------------------- */
75 
76  #if GEOMETRY == CARTESIAN
77 
78  for (i = is; i <= ie; i++) {
79  divB[i] = (vp[i] - vm[i])/GG->dx[i];
80  }
81 
82  #elif GEOMETRY == CYLINDRICAL
83 
84  if (g_dir == IDIR){ /* -- r -- */
85  Ar = grid[IDIR].A;
86  for (i = is; i <= ie; i++) {
87  divB[i] = (vp[i]*Ar[i] - vm[i]*Ar[i - 1])/GG->dV[i];
88  }
89  }else if (g_dir == JDIR){ /* -- z -- */
90  for (i = is; i <= ie; i++) {
91  divB[i] = (vp[i] - vm[i])/GG->dx[i];
92  }
93  }
94 
95  #elif GEOMETRY == POLAR
96 
97  if (g_dir == IDIR){ /* -- r -- */
98  Ar = grid[IDIR].A;
99  for (i = is; i <= ie; i++) {
100  divB[i] = (vp[i]*Ar[i] - vm[i]*Ar[i - 1])/GG->dV[i];
101  }
102  }else if (g_dir == JDIR){ /* -- phi -- */
103  r = grid[IDIR].x[g_i];
104  for (i = is; i <= ie; i++) {
105  divB[i] = (vp[i] - vm[i])/(r*GG->dx[i]);
106  }
107  }else if (g_dir == KDIR){ /* -- z -- */
108  for (i = is; i <= ie; i++) {
109  divB[i] = (vp[i] - vm[i])/GG->dx[i];
110  }
111  }
112 
113  #elif GEOMETRY == SPHERICAL
114 
115  if (g_dir == IDIR){ /* -- r -- */
116  Ar = grid[IDIR].A;
117  for (i = is; i <= ie; i++) {
118  divB[i] = (vp[i]*Ar[i] - vm[i]*Ar[i - 1])/GG->dV[i];
119  }
120  }else if (g_dir == JDIR){ /* -- theta -- */
121  Ath = grid[JDIR].A;
122  r = grid[IDIR].x[g_i];
123  for (i = is; i <= ie; i++) {
124  divB[i] = (vp[i]*Ath[i] - vm[i]*Ath[i - 1]) /
125  (r*GG->dV[i]);
126  }
127  }else if (g_dir == KDIR){ /* -- phi -- */
128  r = grid[IDIR].x[g_i];
129  s = sin(grid[JDIR].x[g_j]);
130  for (i = is; i <= ie; i++) {
131  divB[i] = (vp[i] - vm[i])/(r*s*GG->dx[i]);
132  }
133  }
134 
135  #endif
136 
137  #if BACKGROUND_FIELD == YES
138  bgf = GetBackgroundField (is - 1, ie, CELL_CENTER, grid);
139  #endif
140 
141 /* -------------------------------------------
142  Compute Powell's source term
143  ------------------------------------------- */
144 
145  for (i = is; i <= ie; i++) {
146 
147  v = state->vh[i];
148  src = state->src[i];
149 
150  EXPAND (vx = v[VX1]; ,
151  vy = v[VX2]; ,
152  vz = v[VX3];)
153 
154  EXPAND (bx = btx = v[BX1]; ,
155  by = bty = v[BX2]; ,
156  bz = btz = v[BX3];)
157 
158  #if BACKGROUND_FIELD == YES
159  btx += bgf[i][BX1];
160  bty += bgf[i][BX2];
161  btz += bgf[i][BX3];
162  #endif
163 
164  src[RHO] = 0.0;
165  EXPAND (src[MX1] = -divB[i]*btx; ,
166  src[MX2] = -divB[i]*bty; ,
167  src[MX3] = -divB[i]*btz;)
168 
169  #if HAVE_ENERGY
170  src[ENG] = -divB[i]*(EXPAND(vx*bx, +vy*by, +vz*bz));
171  #endif
172  EXPAND (src[BX1] = -divB[i]*vx; ,
173  src[BX2] = -divB[i]*vy; ,
174  src[BX3] = -divB[i]*vz;)
175  }
176 }
MX3
#define MX3
Definition: mod_defs.h:22
MX1
#define MX1
Definition: mod_defs.h:20
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
YES
#define YES
Definition: pluto.h:25
func_param::v
double v[NVAR]
Definition: eos.h:106
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
GRID::dx
double * dx
Definition: structs.h:83
GetBackgroundField
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
Definition: bckgrnd_field.c:5
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
BXn
int BXn
Definition: globals.h:75
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
STATE_1D::src
double ** src
Definition: structs.h:160
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
GRID
Definition: structs.h:78
MX2
#define MX2
Definition: mod_defs.h:21
divB
double *** divB
Definition: analysis.c:4
GRID::x
double * x
Definition: structs.h:80
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
BACKGROUND_FIELD
void BACKGROUND_FIELD(real x1, real x2, real x3, real *B0)
Definition: init.c:95
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395

Here is the call graph for this function:

Here is the caller graph for this function: