Interp1PrimFifthOrderUpwind.c

Go to the documentation of this file.

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <interpolation.h>
#include <mpivars.h>
#include <hypar.h>

#ifdef with_omp
#include <omp.h>
#endif

#undef  _MINIMUM_GHOSTS_

#define _MINIMUM_GHOSTS_ 3

int Interp1PrimFifthOrderUpwind(
                                double *fI,  
                                double *fC,  
                                double *u,   
                                double *x,   
                                int    upw,  
                                int    dir,  
                                void   *s,   
                                void   *m,   
                                int    uflag 
                               )
{
  HyPar           *solver = (HyPar*)          s;

  int ghosts = solver->ghosts;
  int ndims  = solver->ndims;
  int nvars  = solver->nvars;
  int *dim   = solver->dim_local;
  int *stride= solver->stride_with_ghosts;

  /* define some constants */
  static const double one_by_thirty        = 1.0/30.0,
                      thirteen_by_sixty    = 13.0/60.0,
                      fortyseven_by_sixty  = 47.0/60.0,
                      twentyseven_by_sixty = 27.0/60.0,
                      one_by_twenty        = 1.0/20.0;

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);

  int i;
#pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
  for (i=0; i<N_outer; i++) {
    _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
    _ArrayCopy1D_(index_outer,indexC,ndims);
    _ArrayCopy1D_(index_outer,indexI,ndims);
    for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
      int qm1,qm2,qm3,qp1,qp2,p;
      _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
      if (upw > 0) {
        indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
        qm3 = qm1 - 2*stride[dir];
        qm2 = qm1 -   stride[dir];
        qp1 = qm1 +   stride[dir];
        qp2 = qm1 + 2*stride[dir];
      } else {
        indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
        qm3 = qm1 + 2*stride[dir];
        qm2 = qm1 +   stride[dir];
        qp1 = qm1 -   stride[dir];
        qp2 = qm1 - 2*stride[dir];
      }

      /* Defining stencil points */
      double *fm3, *fm2, *fm1, *fp1, *fp2;
      fm3 = (fC+qm3*nvars);
      fm2 = (fC+qm2*nvars);
      fm1 = (fC+qm1*nvars);
      fp1 = (fC+qp1*nvars);
      fp2 = (fC+qp2*nvars);

      int v;
      for (v=0; v<nvars; v++) {
        (fI+p*nvars)[v] =   one_by_thirty         * fm3[v]
                          - thirteen_by_sixty     * fm2[v]
                          + fortyseven_by_sixty   * fm1[v]
                          + twentyseven_by_sixty  * fp1[v]
                          - one_by_twenty         * fp2[v];
      }
    }
  }

  return(0);
}