TYCHO  1.3.0
 All Data Structures Files Functions Variables Enumerations Enumerator
/home/kapf/tycho_docu/remap.c
Go to the documentation of this file.
1 /*
2  * remap.c
3  *
4  * Author: Wolfgang Kapferer
5  */
6 
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <math.h>
10 
11 #include "variables_global.h"
12 #include "prototypes.h"
13 
17 inline double max_remap(double a, double b) {
18  double tmp;
19 
20  if (a < b) tmp = b;
21  if (a == b)tmp = b;
22  if (a > b) tmp = a;
23 
24  return tmp;
25 }
26 
30 int remap(int nmin, int nmax, int flag, double *a_coef, double *dx,
31  double *ai_coef, double *b_coef, double *bi_coef, double *c_coef,
32  double *d_x, double **para, double *ci_coef, double *dr, double *r6,
33  double *rl, double *diffa, double *da, double *ar, double *flat,
34  double *scratch1, double *scratch2, double *scratch3, double *du,
35  double *u6, double *ul, double *dv, double *v6, double *vl, double *w6,
36  double *wl, double *dq, double *q6, double *ql, double *de, double *e6,
37  double *el, double *xa, double *xa0, double *delta, double *fluxr,
38  double *fluxu, double *fluxv, double *fluxw, double *dw, double *fluxe,
39  double *fluxq, double *dm, double *rho_1D, double *dvol, double *dvol0,
40  double *dm0, double * vx_1D, double *vy_1D, double *vz_1D, double *eng_1D,
41  double *e_int_1D, double *pre_1D) {
42 
43  int n, nn;
44 
45  double fractn, fractn2, ekin;
46  double deltx;
47  double fourthd = 4.0 / 3.0;
48 
49  //Generate interpolation functions, saving da, al for
50  //constructing left and right total energy states.
51  para_coef(nmin - 1, nmax + 1, 1, a_coef, dx,
52  ai_coef, b_coef, bi_coef, c_coef,
53  d_x, para, ci_coef);
54  parabola(nmin - 1, nmax + 1, rho_1D, dr, r6, rl, 3, diffa,
55  da, para, ar, flat, scratch1, scratch2, scratch3);
56  parabola(nmin - 1, nmax + 1, vx_1D, du, u6, ul, 4, diffa,
57  da, para, ar, flat, scratch1, scratch2, scratch3);
58  parabola(nmin - 1, nmax + 1, vy_1D, dv, v6, vl, 5, diffa,
59  da, para, ar, flat, scratch1, scratch2, scratch3);
60  parabola(nmin - 1, nmax + 1, vz_1D, dw, w6, wl, 6, diffa,
61  da, para, ar, flat, scratch1, scratch2, scratch3);
62  parabola(nmin - 1, nmax + 1, e_int_1D, dq, q6, ql, 7, diffa,
63  da, para, ar, flat, scratch1, scratch2, scratch3);
64  parabola(nmin - 1, nmax + 1, eng_1D, de, e6, el, 8, diffa,
65  da, para, ar, flat, scratch1, scratch2, scratch3);
66 
67  for (n = nmin; n <= nmax + 1; n++) {
68  delta[n] = xa[n] - xa0[n];
69  }
70 
71 
72  //Calculate the total mass (fluxr), momentum (fluxu), and energy (fluxe)
73  //in the subshell created by the overlap of the Lagrangian and Eulerian grids.
74  //If the zone face has moved to the left (deltx > 0), use the integral from the
75  //left side of zone n (fluxrr). If the zone face has moved to the right
76  //(deltx < 0), use the integral from the right side of zone nn=n-1 (fluxrl).
77  for (n = nmin; n <= nmax + 1; n++) {
78  deltx = xa[n] - xa0[n];
79  if (deltx >= 0.0) {
80  nn = n - 1;
81  fractn = 0.5 * deltx / dx[nn];
82  fractn2 = 1.0 - fourthd*fractn;
83  fluxr[n] = (rl[nn] + dr[nn] - fractn * (dr[nn] - fractn2 * r6[nn])) * delta[n];
84  fluxu[n] = (ul[nn] + du[nn] - fractn * (du[nn] - fractn2 * u6[nn])) * fluxr[n];
85  fluxv[n] = (vl[nn] + dv[nn] - fractn * (dv[nn] - fractn2 * v6[nn])) * fluxr[n];
86  fluxw[n] = (wl[nn] + dw[nn] - fractn * (dw[nn] - fractn2 * w6[nn])) * fluxr[n];
87  fluxe[n] = (el[nn] + de[nn] - fractn * (de[nn] - fractn2 * e6[nn])) * fluxr[n];
88  fluxq[n] = (ql[nn] + dq[nn] - fractn * (dq[nn] - fractn2 * q6[nn])) * fluxr[n];
89  } else {
90  fractn = 0.5 * deltx / dx[n];
91  fractn2 = 1.0 + fourthd*fractn;
92  fluxr[n] = (rl[n] - fractn * (dr[n] + fractn2 * r6[n])) * delta[n];
93  fluxu[n] = (ul[n] - fractn * (du[n] + fractn2 * u6[n])) * fluxr[n];
94  fluxv[n] = (vl[n] - fractn * (dv[n] + fractn2 * v6[n])) * fluxr[n];
95  fluxw[n] = (wl[n] - fractn * (dw[n] + fractn2 * w6[n])) * fluxr[n];
96  fluxe[n] = (el[n] - fractn * (de[n] + fractn2 * e6[n])) * fluxr[n];
97  fluxq[n] = (ql[n] - fractn * (dq[n] + fractn2 * q6[n])) * fluxr[n];
98  }
99  }
100 
101  //Advect mass, momentum, and energy by moving the sub-shell quantities
102  //into the appropriate Eulerian zone.
103  for (n = nmin; n <= nmax; n++) {
104  dm[n] = rho_1D[n] * dvol[n];
105  dm0[n] = dm[n] + fluxr[n] - fluxr[n + 1];
106  rho_1D[n] = dm0[n] / dvol0[n];
107  rho_1D[n] = max_remap(smallr, rho_1D[n]);
108  dm0[n] = 1.0 / (rho_1D[n] * dvol0[n]);
109  vx_1D[n] = (vx_1D[n] * dm[n] + fluxu[n] - fluxu[n + 1]) * dm0[n];
110  vy_1D[n] = (vy_1D[n] * dm[n] + fluxv[n] - fluxv[n + 1]) * dm0[n];
111  vz_1D[n] = (vz_1D[n] * dm[n] + fluxw[n] - fluxw[n + 1]) * dm0[n];
112  eng_1D[n] = (eng_1D[n] * dm[n] + fluxe[n] - fluxe[n + 1]) * dm0[n];
113  e_int_1D[n] = (e_int_1D[n] * dm[n] + fluxq[n] - fluxq[n + 1]) * dm0[n];
114  //If flow is highly supersonic remap on internal energy, else on total energy
115  ekin = 0.5 * (pow(vx_1D[n], 2) + pow(vy_1D[n], 2) + pow(vz_1D[n], 2));
116  if (ekin / e_int_1D[n] < 100.0) e_int_1D[n] = eng_1D[n] - ekin;
117  pre_1D[n] = Gamma * rho_1D[n] * e_int_1D[n];
118  pre_1D[n] = max_remap(smallp, pre_1D[n]);
119  }
120 
121  return 0;
122 }