vec_pot_update.c

Go to the documentation of this file.

/* ///////////////////////////////////////////////////////////////////// */
/*!
  \file 
  \brief Function to handle vector potential.
 
  Update the cell-centered vector potential by integrating
  \f[
     \pd{\vec{A}}{t} + \vec{E} = 0
  \f] 
  Depending on the time-stepping algorithm, we advance the solution
  in time using the following increments
  \f[
   \begin{array}{ll}
    \rm{RK2:}\quad &\DS A^{n+1} = A^n + \frac{\Delta t^n}{2} (E^n + E^*) 
      \\ \noalign{\medskip}
    \rm{RK3:}\quad &\DS A^{n+1} = A^n + \frac{\Delta t^n}{6} (E^n + E^* + 4E^{**}) 
      \\ \noalign{\medskip}
    \rm{CTU:}\quad &\DS A^{n+1} = A^n + \Delta t^nE^{n+\HALF}
   \end{array}        
  \f]
  The cell-centered electric field \b E is computed by averaging values 
  available at faces or edges, depending on the MHD formulation:
  - staggered_mhd: we average the electric field values 
                   available at cell edges and previously computed
                   . In this
                   case the function does the job only 
                   if called from CT;
  - Cell-cent. mhd: we average the face values using upwind fluxes. 
                    In this case the function does the job only 
                    if called from time stepping algorithms.
                    The flux components of induction equation
                    contain the emf: 
                    \f$ F_{i+\HALF} = (   0,-E_z, E_y)_{i+\HALF}\f$,
                    \f$ F_{j+\HALF} = ( E_z,   0,-E_x)_{j+\HALF}\f$,
                    \f$ F_{k+\HALF} = (-E_y, E_x,   0)_{k+\HALF}\f$,
                    
 
  Note that the vector potential is NOT used in the rest
  of the code and serves solely as a diagnostic tool.
 
  \param [in,out] d   pointer to PLUTO Data structure
  \param [in] vp      pointer to an EMF structure, used only with CT
  \param [in] state   pointer to State_1D structure containing the 
                      electric field components (cell-centered MHD)
  \param [in] grid    pointer to an array of Grid structures
 
  \attention In cylindrical coordinates:
  PLUTO treats cylindrical (but not polar) coordinates as 
  Cartesian so that \f$ (x_1,x_2) = (R,z) \f$ .
  Nevertheless, the correct right-handed cylindrical system is 
  defined by \f$ (x_1,x_2,x_3) = (R,\phi,z)\f$. 
  The fact that the \f$ \phi \f$ coordinate is missing causes an 
  ambiguity when computing the third component of the EMF in 2.5 D:
  \f[
     E_3 = v_2B_1 - v_1B_2
  \f]
  which is (wrongly) called \f$ E_z \f$  but is, indeed, 
  \f$ -E_\phi \f$ .
  This is not a problem in the update of the magnetic field, since
  the contribution coming from \f$ E_3 \f$ is also reversed:
  \f[
  \Delta B_1 = -\Delta t dE_3/dx_2 \quad\Longrightarrow\quad
  \Delta B_r =  \Delta t dE_\phi/dz
  \f]
  as it should be.
  However, it is a problem in the definition of the vector potential
  which is consistently computed using its physical definition.
  Therefore, to correctly update the phi component of the vector
  potential we change sign to \f$ E_z \f$ (\f$ = -E_\phi \f$).
 
  \authors A. Mignone (mignone@ph.unito.it)\n
           P. Tzeferacos (petros.tzeferacos@ph.unito.it)
  \date   Sep 24, 2012
*/
/* ///////////////////////////////////////////////////////////////////// */
#include "pluto.h"

#if (PHYSICS == MHD || PHYSICS == RMHD) && (UPDATE_VECTOR_POTENTIAL == YES)
/* ********************************************************************* */
void VectorPotentialUpdate (const Data *d, const void *vp, 
                            const State_1D *state, const Grid *grid)
/*!
 *
 *********************************************************************** */
{
  int    i,j,k;
  double scrh, dt, sEz = 0.25;
  double Ex, Ey, Ez;

/* ---------------------------------------------------------
    When this function is called from Time stepping 
    routines and CT is used, do nothing.
   --------------------------------------------------------- */

  #ifdef STAGGERED_MHD
   if (vp == NULL) return;
  #endif

  #if GEOMETRY == CYLINDRICAL 
   sEz = -0.25;
  #endif

/* ------------------------------------------------------
     Determine dt according to the integration step
   ------------------------------------------------------ */

  dt = g_dt;
  #if TIME_STEPPING == RK2
   dt *= 0.5;
  #elif TIME_STEPPING == RK3
   dt /= 6.0;
   if (g_intStage == 3) dt *= 4.0;
  #endif

  #ifdef CTU  
   if (g_intStage == 1) return;
  #endif

  #ifdef STAGGERED_MHD
/* -----------------------------------------------------
     If staggered fields are used, we average the 
     EMF available at cell edges:

     A-----------B
     |           |
     |           |
     |     o     |    o = (A+B+C+D)/4
     |           |
     |           |
     D-----------C
   ----------------------------------------------------- */
  {
    EMF *emf;
    emf = (EMF *) vp;

    DOM_LOOP(k,j,i){
      Ez = sEz*(emf->ez[k][j][i]     + emf->ez[k][j-1][i] + 
                emf->ez[k][j-1][i-1] + emf->ez[k][j][i-1]);

      #if DIMENSIONS == 3
       Ex = 0.25*(emf->ex[k][j][i]     + emf->ex[k][j-1][i] + 
                  emf->ex[k-1][j-1][i] + emf->ex[k-1][j][i]);

       Ey = 0.25*(emf->ey[k][j][i]     + emf->ey[k][j][i-1] + 
                  emf->ey[k-1][j][i-1] + emf->ey[k-1][j][i]);
      #endif

      d->Ax3[k][j][i] -= dt*Ez; 
      #if DIMENSIONS == 3
       d->Ax1[k][j][i] -= dt*Ex; 
       d->Ax2[k][j][i] -= dt*Ey; 
      #endif

    }   
  }
  #else 
/* ------------------------------------------------------------
     If cell-center is used, we average the 
     EMF available at cell faces:

     +-----B-----+
     |           |
     |           |
     A     o     C    o = (A+B+C+D)/4
     |           |
     |           |
     +-----D-----+
     
     The components of the flux in the induction equation give 
     the emf:
     
     when g_dir == IDIR we have  F(i+1/2) = (0,-Ez, Ey)
     when g_dir == JDIR we have  F(j+1/2) = (Ez,0,-Ex)
     when g_dir == KDIR we have  F(k+1/2) = (-Ey,Ex,0)
     
   ------------------------------------------------------------ */
  {
    double **f;

    f = state->flux;
    if (g_dir == IDIR){

      IDOM_LOOP(i){
        d->Ax3[g_k][g_j][i]  -= -sEz*dt*(f[i][BX2] + f[i-1][BX2]); /* -Ez */
        #if DIMENSIONS == 3
         d->Ax2[g_k][g_j][i] -= 0.25*dt*(f[i][BX3] + f[i-1][BX3]); /* Ey */
        #endif
      }
    }else if (g_dir == JDIR){

      JDOM_LOOP(j) {
        #if DIMENSIONS == 3
         d->Ax1[g_k][j][g_i] -= -0.25*dt*(f[j][BX3] + f[j-1][BX3]); /* -Ex */
        #endif
        d->Ax3[g_k][j][g_i]  -= sEz*dt*(f[j][BX1] + f[j-1][BX1]);  /* Ez */
      }

    }else if (g_dir == KDIR){
      KDOM_LOOP(k){
        d->Ax1[k][g_j][g_i] -=  0.25*dt*(f[k][BX2] + f[k-1][BX2]);
        d->Ax2[k][g_j][g_i] -= -0.25*dt*(f[k][BX1] + f[k-1][BX1]);
      }
    }
  }
  #endif /* STAGGERED_MHD */
}
#endif  /* UPDATE_VECTOR_POTENTIAL == YES */