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  /* space-dependent refinement */
237  double x = x1*cos(x2), xc = 1.0/sqrt(2.0);
238  double y = x1*sin(x2), yc = 1.0/sqrt(2.0);
239  double r0 = 0.2;
240  double r = sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
241  if (r < r0) grad[k][j][i] = 1.0;
242  else grad[k][j][i] = 0.0;
243 
244  }
245 
246 /* --------------------------------------------------------------
247  5. Multiply U times volume again.
248  -------------------------------------------------------------- */
249 
250  #if GEOMETRY != CARTESIAN
251 
252  #if (CHOMBO_CONS_AM == YES) && ((CHOMBO_REF_VAR == iMPHI) || (CHOMBO_REF_VAR == -1))
253  #if ROTATING_FRAME == YES
254  for(BoxIterator bit(curBox); bit.ok(); ++bit) {
255  const IntVect& iv = bit();
256  UFab(iv,iMPHI) += UFab(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
257  UFab(iv,iMPHI) *= a_dV(iv,1);
258  }
259  #else
260  UFab.mult(a_dV,1,iMPHI);
261  #endif
262  #endif
263 
264  #if CHOMBO_REF_VAR >= 0
265  UFab.mult(a_dV,0,CHOMBO_REF_VAR);
266  #else
267  for (nv = 0; nv < NVAR; nv++) UFab.mult(a_dV,0,nv);
268  #endif
269 
270  #else
271 
272  #if CHOMBO_REF_VAR >= 0
273  if (g_stretch_fact != 1.) UFab.mult(g_stretch_fact,CHOMBO_REF_VAR);
274  #else
275  if (g_stretch_fact != 1.) UFab.mult(g_stretch_fact);
276  #endif
277 
278  #endif
279 
280 /* --------------------------------------------------------------
281  6. Free array
282  -------------------------------------------------------------- */
283 
284  FreeArrayBoxMap(grad, Gbox.kb, Gbox.ke, Gbox.jb, Gbox.je, Gbox.ib, Gbox.ie);
285 
286  #if CHOMBO_REF_VAR >= 0
287  FreeArrayBoxMap(q, Ubox.kb, Ubox.ke, Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
288  #else
289  for (nv = 0; nv < NVAR; nv++){
290  FreeArrayBoxMap(UU[nv], Ubox.kb, Ubox.ke,
291  Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
292  }
293  FreeArrayBox(q, Ubox.kb, Ubox.jb, Ubox.ib);
294  #endif
295 
296 }
297 
298 /* ********************************************************************* */
299 void computeRefVar(double ***UU[], double ***q, double dx, RBox *Ubox)
300 /*!
301  * Compute a user-defined array q(U) function of the conserved
302  * variables.
303  *
304  *
305  *********************************************************************** */
306 {
307  int nv, i, j, k;
308  double us[NVAR], vs[NVAR];
309 
310  BOX_LOOP(Ubox, k, j, i) {
311  VAR_LOOP(nv) us[nv] = UU[nv][k][j][i];
312  q[k][j][i] = us[RHO];
313  }
314 
315 }
316 #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
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
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
REF_CRIT
#define REF_CRIT
Definition: TagCells.cpp:22
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
GRID::x
double * x
Definition: structs.h:80
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
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
computeRefVar
static void computeRefVar(double ***UU[], double ***q, double, RBox *Ubox)
Definition: TagCells.cpp:299
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
CHOMBO_REF_VAR
#define CHOMBO_REF_VAR
Definition: TagCells.cpp:16