TagCells.cpp

Go to the documentation of this file.

#include <cstdio>
#include <string>
using std::string;

#include "PatchPluto.H"
#include "LoHiSide.H"

static void computeRefVar(double ***UU[], double ***q, double, RBox *Ubox);

#if (EOS != ISOTHERMAL) && (CHOMBO_EN_SWITCH == NO)
 #ifndef CHOMBO_REF_VAR  
  #define CHOMBO_REF_VAR ENG 
 #endif
#else
 #ifndef CHOMBO_REF_VAR  
  #define CHOMBO_REF_VAR RHO
 #endif
#endif

// #define CHOMBO_REF_VAR -1  /* means user-defined */

#define REF_CRIT 2   /* 1 == first derivative, 2 == second derivative */

/* ************************************************************************* */
void PatchPluto::computeRefGradient(FArrayBox& gFab, FArrayBox& UFab, 
                                    const FArrayBox& a_dV, const Box& b)
/*!
 * Tag zones for refinement using gradient of the conservative 
 * variables.
 * The gradient is computed by standard finite differences using
 *
 * - REF_CRIT equal to 1 --> compute (normalized) gradient using 1st 
 *                           derivative of the solution;
 * - REF_CRIT equal to 2 --> compute (normalized) gradient using 2nd 
 *                           derivative of the solution (default);
 *                           This approach is based on Lohner (1987).
 *
 * Zones will be flagged for refinement whenever grad[k][j][i] exceeds 
 * the threshold value specified by the 'Refine_thresh' parameter read in
 * pluto.ini.
 *
 * Derivatives are computed using the conserved variable
 * U[CHOMBO_REF_VAR] 
 * where CHOMBO_REF_VAR is taken to be energy density (default).
 * However, by setting CHOMBO_REF_VAR = -1, you can provide your own 
 * physical variable through the function computeRefVar().
 * 
 * \authors C. Zanni   (zanni@oato.inaf.it)\n
 *          A. Mignone (mignone@ph.unito.it)
 * \date    Oct 11, 2012
 *************************************************************************** */
{
  CH_assert(m_isDefined);

  int nv, i, j, k;
  double x1, dqx_p, dqx_m, dqx, d2qx, den_x;
  double x2, dqy_p, dqy_m, dqy, d2qy, den_y;
  double x3, dqz_p, dqz_m, dqz, d2qz, den_z;
  double gr1, gr2, eps = 0.01;
 #if CHOMBO_REF_VAR == -1
  double ***UU[NVAR];
 #endif
  double ***q, ***grad;
  RBox  Ubox, Gbox;

/* -- check ref criterion -- */

  #if REF_CRIT != 1 && REF_CRIT != 2
   print ("! TagCells.cpp: Refinement criterion not valid\n");
   QUIT_PLUTO(1);
  #endif

/* -----------------------------------------------------
   1. The solution array U is defined on the box 
      [Uib, Uie] x [Ujb, Uje] x [Ukb, Uke], which 
      differs from that of gFab ([Gib,...Gke]), 
      typically one point larger in each direction. 
   ----------------------------------------------------- */
    
  Ubox.jb = Ubox.je = Ubox.kb = Ubox.ke = 0;
  Gbox.jb = Gbox.je = Gbox.kb = Gbox.ke = 0;

  D_EXPAND(Ubox.ib = UFab.loVect()[IDIR]; Ubox.ie = UFab.hiVect()[IDIR]; ,
           Ubox.jb = UFab.loVect()[JDIR]; Ubox.je = UFab.hiVect()[JDIR]; ,
           Ubox.kb = UFab.loVect()[KDIR]; Ubox.ke = UFab.hiVect()[KDIR]; );

  D_EXPAND(Gbox.ib = gFab.loVect()[IDIR]; Gbox.ie = gFab.hiVect()[IDIR]; ,
           Gbox.jb = gFab.loVect()[JDIR]; Gbox.je = gFab.hiVect()[JDIR]; ,
           Gbox.kb = gFab.loVect()[KDIR]; Gbox.ke = gFab.hiVect()[KDIR]; );

/* --------------------------------------------------------
   2. Input solution array (UFab.dataPtr(nv)) is defined 
      as dV*U/dx, where U is an array of conservative 
      variables and dV is the zone volume. 
      To obtain U we must divide by volume.
   -------------------------------------------------------- */
 
  #if GEOMETRY != CARTESIAN

   #if CHOMBO_REF_VAR  >= 0
    UFab.divide(a_dV,0,CHOMBO_REF_VAR);
   #else
    for (nv = 0; nv < NVAR; nv++) UFab.divide(a_dV,0,nv); 
   #endif

   #if (CHOMBO_CONS_AM == YES) && \
       ((CHOMBO_REF_VAR == iMPHI) ||  (CHOMBO_REF_VAR == -1))
    #if ROTATING_FRAME == YES
     Box curBox = UFab.box();
     for(BoxIterator bit(curBox); bit.ok(); ++bit) {
       const IntVect& iv = bit();
       UFab(iv,iMPHI) /= a_dV(iv,1);
       UFab(iv,iMPHI) -= UFab(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
     }
    #else
     UFab.divide(a_dV,1,iMPHI);
    #endif
   #endif

  #else

   #if CHOMBO_REF_VAR  >= 0
    if (ch_vol != 1.) UFab.divide(ch_vol,CHOMBO_REF_VAR);
   #else
    if (ch_vol != 1.) UFab.divide(ch_vol);
   #endif

  #endif

/* ---------------------------------------------
   3. Set refinement variable
   --------------------------------------------- */

  #if CHOMBO_REF_VAR >= 0
   q = ArrayBoxMap(Ubox.kb, Ubox.ke,
                   Ubox.jb, Ubox.je,
                   Ubox.ib, Ubox.ie, UFab.dataPtr(CHOMBO_REF_VAR));
  #else
   for (nv = 0; nv < NVAR; nv++)
     UU[nv] = ArrayBoxMap(Ubox.kb, Ubox.ke,
                          Ubox.jb, Ubox.je,
                          Ubox.ib, Ubox.ie, UFab.dataPtr(nv));
   q = ArrayBox(Ubox.kb, Ubox.ke, Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
   computeRefVar(UU, q, m_dx, &Ubox);
  #endif
  
  grad = ArrayBoxMap(Gbox.kb, Gbox.ke, 
                     Gbox.jb, Gbox.je, 
                     Gbox.ib, Gbox.ie, gFab.dataPtr(0));

/* ----------------------------------------------------------------
   4. Main spatial loop for zone tagging based on 1st 
     (REF_CRIT = 1) or 2nd (REF_CRIT = 2) derivative error norm. 
   ---------------------------------------------------------------- */

  BOX_LOOP(&Gbox, k, j, i){
    x3 = (k + 0.5)*m_dx*zstretch + g_domBeg[KDIR];
    x2 = (j + 0.5)*m_dx*ystretch + g_domBeg[JDIR];
   #if CHOMBO_LOGR == NO
    x1 = (i + 0.5)*m_dx          + g_domBeg[IDIR];
   #else
    double xl = g_domBeg[IDIR] + i*m_dx;
    double xr = xl + m_dx;
    x1 = g_domBeg[IDIR]*0.5*(exp(xr)+exp(xl));
   #endif 

    D_EXPAND(dqx_p =    q[k][j][i+1] - q[k][j][i];
             dqx_m = - (q[k][j][i-1] - q[k][j][i]);  ,
             dqy_p =    q[k][j+1][i] - q[k][j][i];
             dqy_m = - (q[k][j-1][i] - q[k][j][i]);  ,
             dqz_p =    q[k+1][j][i] - q[k][j][i];
             dqz_m = - (q[k-1][j][i] - q[k][j][i]);)

  /* --------------------------------------------------------------
      Physical boundary values are not up to date and should be 
      excluded from gradient computation. 
      In this case, left and right derivatives are set equal to 
      each other. This will not trigger refinement in the leftmost 
      and rightmost internal zones (using 2nd derivative) but we 
      really don't care since buffer size will do the job.
     -------------------------------------------------------------- */
      
    D_EXPAND(if (i == 0) dqx_m = dqx_p;  ,
             if (j == 0) dqy_m = dqy_p;  ,
             if (k == 0) dqz_m = dqz_p;)

    D_EXPAND(if (i == m_domain.size(IDIR)-1) dqx_p = dqx_m;  ,
             if (j == m_domain.size(JDIR)-1) dqy_p = dqy_m;  ,
             if (k == m_domain.size(KDIR)-1) dqz_p = dqz_m;)

  /* -----------------------------------------------
         Compute gradient using 1st derivative 
      ---------------------------------------------- */

    #if REF_CRIT == 1
     D_EXPAND(dqx = dqx_p + dqx_m;  ,
              dqy = dqy_p + dqy_m;  ,
              dqz = dqz_p + dqz_m;)

     D_EXPAND(den_x = fabs(q[k][j][i+1]) + fabs(q[k][j][i-1]);  ,
              den_y = fabs(q[k][j+1][i]) + fabs(q[k][j-1][i]);  ,
              den_z = fabs(q[k+1][j][i]) + fabs(q[k-1][j][i]);)

     gr1  = D_EXPAND(dqx*dqx, + dqy*dqy, + dqz*dqz);
     gr1 /= D_EXPAND(den_x*den_x, + den_y*den_y, + den_z*den_z);

     grad[k][j][i] = sqrt(gr1);
    #endif

  /* -----------------------------------------------
         Compute gradient using 2nd derivative 
      ---------------------------------------------- */

    #if REF_CRIT == 2
     D_EXPAND(d2qx = dqx_p - dqx_m;  ,
              d2qy = dqy_p - dqy_m;  ,
              d2qz = dqz_p - dqz_m;)

     D_EXPAND(
       den_x = 2.0*fabs(q[k][j][i]) + fabs(q[k][j][i+1]) + fabs(q[k][j][i-1]);
       den_x = fabs(dqx_p) + fabs(dqx_m) + eps*den_x;    ,

       den_y = 2.0*fabs(q[k][j][i]) + fabs(q[k][j+1][i]) + fabs(q[k][j-1][i]);
       den_y = fabs(dqy_p) + fabs(dqy_m) + eps*den_y;    ,

       den_z = 2.0*fabs(q[k][j][i]) + fabs(q[k+1][j][i]) + fabs(q[k-1][j][i]);
       den_z = fabs(dqz_p) + fabs(dqz_m) + eps*den_z;
     )

     gr2  = D_EXPAND(d2qx*d2qx,   + d2qy*d2qy,   + d2qz*d2qz);
     gr2 /= D_EXPAND(den_x*den_x, + den_y*den_y, + den_z*den_z);

     grad[k][j][i] = sqrt(gr2);
    #endif
if (x1 < 0.5 && m_level <= 15) grad[k][j][i] = 1.0;
if (x1 > 1.8 && m_level <= 15) grad[k][j][i] = 1.0;

  }

/* --------------------------------------------------------------
   5. Multiply U times volume again.
   -------------------------------------------------------------- */

  #if GEOMETRY != CARTESIAN

   #if (CHOMBO_CONS_AM == YES) && ((CHOMBO_REF_VAR == iMPHI) || (CHOMBO_REF_VAR == -1))
    #if ROTATING_FRAME == YES
     for(BoxIterator bit(curBox); bit.ok(); ++bit) {
       const IntVect& iv = bit();
       UFab(iv,iMPHI) += UFab(iv,RHO)*a_dV(iv,1)*g_OmegaZ;
       UFab(iv,iMPHI) *= a_dV(iv,1);
     }
    #else
     UFab.mult(a_dV,1,iMPHI);
    #endif
   #endif

   #if CHOMBO_REF_VAR >= 0
    UFab.mult(a_dV,0,CHOMBO_REF_VAR);
   #else
    for (nv = 0; nv < NVAR; nv++) UFab.mult(a_dV,0,nv);
   #endif

  #else

   #if CHOMBO_REF_VAR  >= 0
    if (ch_vol != 1.) UFab.mult(ch_vol,CHOMBO_REF_VAR);
   #else
    if (ch_vol != 1.) UFab.mult(ch_vol);
   #endif

  #endif

/* --------------------------------------------------------------
   6. Free array
   -------------------------------------------------------------- */
   
  FreeArrayBoxMap(grad, Gbox.kb, Gbox.ke, Gbox.jb, Gbox.je, Gbox.ib, Gbox.ie);

  #if CHOMBO_REF_VAR >= 0
   FreeArrayBoxMap(q, Ubox.kb, Ubox.ke, Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
  #else
   for (nv = 0; nv < NVAR; nv++){
     FreeArrayBoxMap(UU[nv], Ubox.kb, Ubox.ke, 
                             Ubox.jb, Ubox.je, Ubox.ib, Ubox.ie);
   }
   FreeArrayBox(q, Ubox.kb, Ubox.jb, Ubox.ib);
  #endif

}

/* ********************************************************************* */
void computeRefVar(double ***UU[], double ***q, double dx, RBox *Ubox)
/*!
 * Compute a user-defined array q(U) function of the conserved
 * variables.
 *
 *
 *********************************************************************** */
{
  int nv, i, j, k;
  double us[NVAR], vs[NVAR];

  BOX_LOOP(Ubox, k, j, i) {
    VAR_LOOP(nv) us[nv] = UU[nv][k][j][i]; 
    q[k][j][i] = us[RHO];
  }

}
#undef REF_CRIT