PLUTO
TagCells.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <string>
3 using std::string;
4 
5 #include "PatchPluto.H"
6 #include "LoHiSide.H"
7 
8 static void computeRefVar(double ***UU[], double ***q, double, RBox *Ubox);
9 
10 #if (EOS != ISOTHERMAL) && (CHOMBO_EN_SWITCH == NO)
11  #ifndef CHOMBO_REF_VAR
12  #define CHOMBO_REF_VAR ENG
13  #endif
14 #else
15  #ifndef CHOMBO_REF_VAR
16  #define CHOMBO_REF_VAR RHO
17  #endif
18 #endif
19 
20 // #define CHOMBO_REF_VAR -1 /* means user-defined */
21 
22 #define REF_CRIT 2 /* 1 == first derivative, 2 == second derivative */
23 
24 /* ************************************************************************* */
25 void PatchPluto::computeRefGradient(FArrayBox& gFab, FArrayBox& UFab,
26  const FArrayBox& a_dV, const Box& b)
27 /*!
28  * Tag zones for refinement using gradient of the conservative
29  * variables.
30  * The gradient is computed by standard finite differences using
31  *
32  * - REF_CRIT equal to 1 --> compute (normalized) gradient using 1st
33  * derivative of the solution;
34  * - REF_CRIT equal to 2 --> compute (normalized) gradient using 2nd
35  * derivative of the solution (default);
36  * This approach is based on Lohner (1987).
37  *
38  * Zones will be flagged for refinement whenever grad[k][j][i] exceeds
39  * the threshold value specified by the 'Refine_thresh' parameter read in
40  * pluto.ini.
41  *
42  * Derivatives are computed using the conserved variable
43  * U[CHOMBO_REF_VAR]
44  * where CHOMBO_REF_VAR is taken to be energy density (default).
45  * However, by setting CHOMBO_REF_VAR = -1, you can provide your own
46  * physical variable through the function computeRefVar().
47  *
48  * \authors C. Zanni (zanni@oato.inaf.it)\n
49  * A. Mignone (mignone@ph.unito.it)
50  * \date Oct 11, 2012
51  *************************************************************************** */
52 {
53  CH_assert(m_isDefined);
54 
55  int nv, i, j, k;
56  double x1, dqx_p, dqx_m, dqx, d2qx, den_x;
57  double x2, dqy_p, dqy_m, dqy, d2qy, den_y;
58  double x3, dqz_p, dqz_m, dqz, d2qz, den_z;
59  double gr1, gr2, eps = 0.01;
60  #if CHOMBO_REF_VAR == -1
61  double ***UU[NVAR];
62  #endif
63  double ***q, ***grad;
64  RBox Ubox, Gbox;
65 
66 /* -- check ref criterion -- */
67 
68  #if REF_CRIT != 1 && REF_CRIT != 2
69  print ("! TagCells.cpp: Refinement criterion not valid\n");
70  QUIT_PLUTO(1);
71  #endif
72 
73 /* -----------------------------------------------------
74  1. The solution array U is defined on the box
75  [Uib, Uie] x [Ujb, Uje] x [Ukb, Uke], which
76  differs from that of gFab ([Gib,...Gke]),
77  typically one point larger in each direction.
78  ----------------------------------------------------- */
79 
80  Ubox.jb = Ubox.je = Ubox.kb = Ubox.ke = 0;
81  Gbox.jb = Gbox.je = Gbox.kb = Gbox.ke = 0;
82 
83  D_EXPAND(Ubox.ib = UFab.loVect()[IDIR]; Ubox.ie = UFab.hiVect()[IDIR]; ,
84  Ubox.jb = UFab.loVect()[JDIR]; Ubox.je = UFab.hiVect()[JDIR]; ,
85  Ubox.kb = UFab.loVect()[KDIR]; Ubox.ke = UFab.hiVect()[KDIR]; );
86 
87  D_EXPAND(Gbox.ib = gFab.loVect()[IDIR]; Gbox.ie = gFab.hiVect()[IDIR]; ,
88  Gbox.jb = gFab.loVect()[JDIR]; Gbox.je = gFab.hiVect()[JDIR]; ,
89  Gbox.kb = gFab.loVect()[KDIR]; Gbox.ke = gFab.hiVect()[KDIR]; );
90 
91 /* --------------------------------------------------------
92  2. Input solution array (UFab.dataPtr(nv)) is defined
93  as dV*U/dx, where U is an array of conservative
94  variables and dV is the zone volume.
95  To obtain U we must divide by volume.
96  -------------------------------------------------------- */
97 
98  #if GEOMETRY != CARTESIAN
99 
100  #if CHOMBO_REF_VAR >= 0
101  UFab.divide(a_dV,0,CHOMBO_REF_VAR);
102  #else
103  for (nv = 0; nv < NVAR; nv++) UFab.divide(a_dV,0,nv);
104  #endif
105 
106  #if (CHOMBO_CONS_AM == YES) && \
107  ((CHOMBO_REF_VAR == iMPHI) || (CHOMBO_REF_VAR == -1))
108  #if ROTATING_FRAME == YES
109  Box curBox = UFab.box();
110  for(BoxIterator bit(curBox); bit.ok(); ++bit) {
111  const IntVect& iv = bit();
112  UFab(iv,iMPHI) /= a_dV(iv,1);
113  UFab(iv,iMPHI) -= UFab(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
114  }
115  #else
116  UFab.divide(a_dV,1,iMPHI);
117  #endif
118  #endif
119 
120  #else
121 
122  #if CHOMBO_REF_VAR >= 0
123  if (g_stretch_fact != 1.) UFab.divide(g_stretch_fact,CHOMBO_REF_VAR);
124  #else
125  if (g_stretch_fact != 1.) UFab.divide(g_stretch_fact);
126  #endif
127 
128  #endif
129 
130 /* ---------------------------------------------
131  3. Set refinement variable
132  --------------------------------------------- */
133 
134  #if CHOMBO_REF_VAR >= 0
135  q = ArrayBoxMap(Ubox.kb, Ubox.ke,
136  Ubox.jb, Ubox.je,
137  Ubox.ib, Ubox.ie, UFab.dataPtr(CHOMBO_REF_VAR));
138  #else
139  for (nv = 0; nv < NVAR; nv++)
140  UU[nv] = ArrayBoxMap(Ubox.kb, Ubox.ke,
141  Ubox.jb, Ubox.je,
142  Ubox.ib, Ubox.ie, UFab.dataPtr(nv));
143  q = ArrayBox(Ubox.kb, Ubox.ke, Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
144  computeRefVar(UU, q, m_dx, &Ubox);
145  #endif
146 
147  grad = ArrayBoxMap(Gbox.kb, Gbox.ke,
148  Gbox.jb, Gbox.je,
149  Gbox.ib, Gbox.ie, gFab.dataPtr(0));
150 
151 /* ----------------------------------------------------------------
152  4. Main spatial loop for zone tagging based on 1st
153  (REF_CRIT = 1) or 2nd (REF_CRIT = 2) derivative error norm.
154  ---------------------------------------------------------------- */
155 
156  BOX_LOOP(&Gbox, k, j, i){
157  x3 = (k + 0.5)*m_dx*g_x3stretch + g_domBeg[KDIR];
158  x2 = (j + 0.5)*m_dx*g_x2stretch + g_domBeg[JDIR];
159  #if CHOMBO_LOGR == NO
160  x1 = (i + 0.5)*m_dx + g_domBeg[IDIR];
161  #else
162  double xl = g_domBeg[IDIR] + i*m_dx;
163  double xr = xl + m_dx;
164  x1 = g_domBeg[IDIR]*0.5*(exp(xr)+exp(xl));
165  #endif
166 
167  D_EXPAND(dqx_p = q[k][j][i+1] - q[k][j][i];
168  dqx_m = - (q[k][j][i-1] - q[k][j][i]); ,
169  dqy_p = q[k][j+1][i] - q[k][j][i];
170  dqy_m = - (q[k][j-1][i] - q[k][j][i]); ,
171  dqz_p = q[k+1][j][i] - q[k][j][i];
172  dqz_m = - (q[k-1][j][i] - q[k][j][i]);)
173 
174  /* --------------------------------------------------------------
175  Physical boundary values are not up to date and should be
176  excluded from gradient computation.
177  In this case, left and right derivatives are set equal to
178  each other. This will not trigger refinement in the leftmost
179  and rightmost internal zones (using 2nd derivative) but we
180  really don't care since buffer size will do the job.
181  -------------------------------------------------------------- */
182 
183  D_EXPAND(if (i == 0) dqx_m = dqx_p; ,
184  if (j == 0) dqy_m = dqy_p; ,
185  if (k == 0) dqz_m = dqz_p;)
186 
187  D_EXPAND(if (i == m_domain.size(IDIR)-1) dqx_p = dqx_m; ,
188  if (j == m_domain.size(JDIR)-1) dqy_p = dqy_m; ,
189  if (k == m_domain.size(KDIR)-1) dqz_p = dqz_m;)
190 
191  /* -----------------------------------------------
192  Compute gradient using 1st derivative
193  ---------------------------------------------- */
194 
195  #if REF_CRIT == 1
196  D_EXPAND(dqx = dqx_p + dqx_m; ,
197  dqy = dqy_p + dqy_m; ,
198  dqz = dqz_p + dqz_m;)
199 
200  D_EXPAND(den_x = fabs(q[k][j][i+1]) + fabs(q[k][j][i-1]); ,
201  den_y = fabs(q[k][j+1][i]) + fabs(q[k][j-1][i]); ,
202  den_z = fabs(q[k+1][j][i]) + fabs(q[k-1][j][i]);)
203 
204  gr1 = D_EXPAND(dqx*dqx, + dqy*dqy, + dqz*dqz);
205  gr1 /= D_EXPAND(den_x*den_x, + den_y*den_y, + den_z*den_z);
206 
207  grad[k][j][i] = sqrt(gr1);
208  #endif
209 
210  /* -----------------------------------------------
211  Compute gradient using 2nd derivative
212  ---------------------------------------------- */
213 
214  #if REF_CRIT == 2
215  D_EXPAND(d2qx = dqx_p - dqx_m; ,
216  d2qy = dqy_p - dqy_m; ,
217  d2qz = dqz_p - dqz_m;)
218 
219  D_EXPAND(
220  den_x = 2.0*fabs(q[k][j][i]) + fabs(q[k][j][i+1]) + fabs(q[k][j][i-1]);
221  den_x = fabs(dqx_p) + fabs(dqx_m) + eps*den_x; ,
222 
223  den_y = 2.0*fabs(q[k][j][i]) + fabs(q[k][j+1][i]) + fabs(q[k][j-1][i]);
224  den_y = fabs(dqy_p) + fabs(dqy_m) + eps*den_y; ,
225 
226  den_z = 2.0*fabs(q[k][j][i]) + fabs(q[k+1][j][i]) + fabs(q[k-1][j][i]);
227  den_z = fabs(dqz_p) + fabs(dqz_m) + eps*den_z;
228  )
229 
230  gr2 = D_EXPAND(d2qx*d2qx, + d2qy*d2qy, + d2qz*d2qz);
231  gr2 /= D_EXPAND(den_x*den_x, + den_y*den_y, + den_z*den_z);
232 
233  grad[k][j][i] = sqrt(gr2);
234  #endif
235 
236  #if GEOMETRY == CYLINDRICAL
237  double R= sqrt(x1*x1 + x2*x2);
238  if (R <= 1.0) grad[k][j][i] = 1.e3;
239  if (i == 0) grad[k][j][i] = 1.e3;
240  #elif GEOMETRY == CARTESIAN
241  double R= sqrt(x1*x1 + x2*x2 + x3*x3);
242  if (R <= 1.0) grad[k][j][i] = 1.e3;
243  #endif
244 
245  }
246 
247 /* --------------------------------------------------------------
248  5. Multiply U times volume again.
249  -------------------------------------------------------------- */
250 
251  #if GEOMETRY != CARTESIAN
252 
253  #if (CHOMBO_CONS_AM == YES) && ((CHOMBO_REF_VAR == iMPHI) || (CHOMBO_REF_VAR == -1))
254  #if ROTATING_FRAME == YES
255  for(BoxIterator bit(curBox); bit.ok(); ++bit) {
256  const IntVect& iv = bit();
257  UFab(iv,iMPHI) += UFab(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
258  UFab(iv,iMPHI) *= a_dV(iv,1);
259  }
260  #else
261  UFab.mult(a_dV,1,iMPHI);
262  #endif
263  #endif
264 
265  #if CHOMBO_REF_VAR >= 0
266  UFab.mult(a_dV,0,CHOMBO_REF_VAR);
267  #else
268  for (nv = 0; nv < NVAR; nv++) UFab.mult(a_dV,0,nv);
269  #endif
270 
271  #else
272 
273  #if CHOMBO_REF_VAR >= 0
274  if (g_stretch_fact != 1.) UFab.mult(g_stretch_fact,CHOMBO_REF_VAR);
275  #else
276  if (g_stretch_fact != 1.) UFab.mult(g_stretch_fact);
277  #endif
278 
279  #endif
280 
281 /* --------------------------------------------------------------
282  6. Free array
283  -------------------------------------------------------------- */
284 
285  FreeArrayBoxMap(grad, Gbox.kb, Gbox.ke, Gbox.jb, Gbox.je, Gbox.ib, Gbox.ie);
286 
287  #if CHOMBO_REF_VAR >= 0
288  FreeArrayBoxMap(q, Ubox.kb, Ubox.ke, Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
289  #else
290  for (nv = 0; nv < NVAR; nv++){
291  FreeArrayBoxMap(UU[nv], Ubox.kb, Ubox.ke,
292  Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
293  }
294  FreeArrayBox(q, Ubox.kb, Ubox.jb, Ubox.ib);
295  #endif
296 
297 }
298 
299 /* ********************************************************************* */
300 void computeRefVar(double ***UU[], double ***q, double dx, RBox *Ubox)
301 /*!
302  * Compute a user-defined array q(U) function of the conserved
303  * variables.
304  *
305  *
306  *********************************************************************** */
307 {
308  int nv, i, j, k;
309  double us[NVAR], vs[NVAR];
310 
311  BOX_LOOP(Ubox, k, j, i) {
312  VAR_LOOP(nv) us[nv] = UU[nv][k][j][i];
313  q[k][j][i] = us[RHO];
314  }
315 
316 }
317 #undef REF_CRIT
iMPHI
#define iMPHI
Definition: mod_defs.h:71
RBOX::jb
int jb
Lower corner index in the x2 direction.
Definition: structs.h:349
CHOMBO_REF_VAR
#define CHOMBO_REF_VAR
Definition: TagCells.cpp:16
GRID::xr
double * xr
Definition: structs.h:81
RHO
#define RHO
Definition: mod_defs.h:19
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
RBOX::kb
int kb
Lower corner index in the x3 direction.
Definition: structs.h:351
FreeArrayBox
void FreeArrayBox(double ***t, long nrl, long ncl, long ndl)
Definition: arrays.c:403
g_OmegaZ
#define g_OmegaZ
Definition: init.c:64
ArrayBox
double *** ArrayBox(long int nrl, long int nrh, long int ncl, long int nch, long int ndl, long int ndh)
Definition: arrays.c:341
REF_CRIT
#define REF_CRIT
Definition: TagCells.cpp:22
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
RBOX::ib
int ib
Lower corner index in the x1 direction.
Definition: structs.h:347
GRID::xl
double * xl
Definition: structs.h:82
eps
#define eps
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
g_domBeg
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
computeRefVar
static void computeRefVar(double ***UU[], double ***q, double, RBox *Ubox)
Definition: TagCells.cpp:300
i
int i
Definition: analysis.c:2
RBOX::ie
int ie
Upper corner index in the x1 direction.
Definition: structs.h:348
RBOX::ke
int ke
Upper corner index in the x3 direction.
Definition: structs.h:352
RBOX::je
int je
Upper corner index in the x2 direction.
Definition: structs.h:350
ArrayBoxMap
double *** ArrayBoxMap(int nrl, int nrh, int ncl, int nch, int ndl, int ndh, double *uptr)
Definition: arrays.c:537
FreeArrayBoxMap
void FreeArrayBoxMap(double ***t, int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
Definition: arrays.c:586
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
RBOX
Definition: structs.h:346
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125