Interp1PrimSecondOrderMUSCL.c

Go to the documentation of this file.

#include <stdio.h>
#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_ 2

int Interp1PrimSecondOrderMUSCL (
                                  double *fI,  
                                  double *fC,  
                                  double *u,   
                                  double *x,   
                                  int    upw,  
                                  int    dir,  
                                  void   *s,   
                                  void   *m,   
                                  int    uflag 
                                )
{
  HyPar           *solver = (HyPar*)          s;
  MUSCLParameters *muscl  = (MUSCLParameters*) solver->interp;

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

  /* define some constants */
  double half = 1.0/2.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;
  if (upw > 0) {
#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 p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
        int qm2,qm1,qp1,v;
        indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
        indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
        indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
        for (v=0; v<nvars; v++)  {
          /* Defining stencil points */
          double m2, m1, p1;
          m2 = fC[qm2*nvars+v];
          m1 = fC[qm1*nvars+v];
          p1 = fC[qp1*nvars+v];

          double slope_ratio = (m1 - m2) / ((p1 - m1) + 1e-40);
          double phi = muscl->LimiterFunction(slope_ratio);
          fI[p*nvars+v] =  m1 + 0.5 * phi * (p1-m1);
        }
      }
    }
  } else {
#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 p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
        int qm1,qp1,qp2,v;
        indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
        indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
        indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
        for (v=0; v<nvars; v++)  {
          /* Defining stencil points */
          double m1, p1, p2;
          m1 = fC[qm1*nvars+v];
          p1 = fC[qp1*nvars+v];
          p2 = fC[qp2*nvars+v];

          double slope_ratio = (p1 - m1) / ((p2 - p1) + 1e-40);
          double phi = muscl->LimiterFunction(slope_ratio);
          fI[p*nvars+v] = p1 + 0.5 * phi * (p1-p2);
        }
      }
    }
  }
  
  return(0);
}