WENOInitialize.c File Reference

Initializes the WENO-type schemes. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <arrayfunctions_gpu.h>
#include <interpolation.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
int	WENOFifthOrderInitializeWeights (double *const, double *const, double *const, const int *const, int, void *, void *)
int	WENOFifthOrderCalculateWeights (double *, double *, double *, int, void *, void *)
int	WENOFifthOrderCalculateWeightsChar (double *, double *, double *, int, void *, void *)
int	gpuWENOFifthOrderCalculateWeights (double *, double *, double *, int, void *, void *)
int	WENOInitialize (void *s, void *m, char *scheme, char *type)

Detailed Description

Initializes the WENO-type schemes.

Author

Debojyoti Ghosh, Youngdae Kim

Definition in file WENOInitialize.c.

Function Documentation

int WENOFifthOrderInitializeWeights	(	double *const	a_w1,
		double *const	a_w2,
		double *const	a_w3,
		const int *const	a_offset,
		int	dir,
		void *	s,
		void *	m
	)

This function initializes the arrays to store the nonlinear weights for fifth order WENO-type schemes to their optimal values (i.e. values of the weights for a perfectly smooth solution). 5th order methods need three weights at each grid interface.

Notes:

• The three weights at all the grid interfaces along all the dimensions are stored in three big arrays (WENOParameters::w1, WENOParameters::w2, WENOParameters::w3), one for each weight.
• This function initializes the values for the spatial dimension given by the input dir.
• offset indicates the location in WENOParameters::w1,WENOParameters::w2, and WENOParameters::w3 from where the chunk of memory with the weights for the spatial dimension dir starts.
• This chunk of memory (1D array) represents a multidimensional array of the following size: 4 X (number of interfaces) X (number of solution components HyPar::nvars). The factor 4 comes from the fact that the array stores weights for the interpolation of left-biased flux function, left-biased solution, right-biased flux function, and right-biased solution, one after the other, in this order.
• The weights are initialized to their optimal values so that, if no limiting is specified (WENOParameters::no_limiting = 1), then no further computation of weights are required.

Parameters

a_w1	Weight array
a_w2	Weight array
a_w3	Weight array
a_offset	Offset array
dir	Spatial dimension
s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 33 of file WENOFifthOrderInitializeWeights.c.

{
  HyPar           *solver = (HyPar*)          s;
  WENOParameters  *weno   = (WENOParameters*) solver->interp;
  MPIVariables    *mpi    = (MPIVariables*)   m;
  int             done;
  double          *ww1, *ww2, *ww3;


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

  /* calculate dimension offset */
  int offset = a_offset[dir];

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int 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;

  /* calculate weights for a left-biased interpolation */
  ww1 = a_w1 + offset;
  ww2 = a_w2 + offset;
  ww3 = a_w3 + offset;
  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexI,ndims);
    for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
      int p, v;
      _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
      for (v=0; v<nvars; v++)  {
        /* optimal weights*/
        double c1, c2, c3;
        if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
          if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
              || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
            /* Use WENO5 at the physical boundaries */
            c1 = _WENO_OPTIMAL_WEIGHT_1_;
            c2 = _WENO_OPTIMAL_WEIGHT_2_;
            c3 = _WENO_OPTIMAL_WEIGHT_3_;
          } else {
            /* CRWENO5 at the interior points */
            c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
            c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
            c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
          }
        } else {
          /* WENO5 and HCWENO5 */
          c1 = _WENO_OPTIMAL_WEIGHT_1_;
          c2 = _WENO_OPTIMAL_WEIGHT_2_;
          c3 = _WENO_OPTIMAL_WEIGHT_3_;
        }

        /* save the weights */
        *(ww1+p*nvars+v) = c1;
        *(ww2+p*nvars+v) = c2;
        *(ww3+p*nvars+v) = c3;
      }
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  ww1 = a_w1 + weno->size + offset;
  ww2 = a_w2 + weno->size + offset;
  ww3 = a_w3 + weno->size + offset;
  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexI,ndims);
    for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
      int p, v;
      _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
      for (v=0; v<nvars; v++)  {

        /* optimal weights*/
        double c1, c2, c3;
        if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
          if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
              || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
            /* Use WENO5 at the physical boundaries */
            c1 = _WENO_OPTIMAL_WEIGHT_1_;
            c2 = _WENO_OPTIMAL_WEIGHT_2_;
            c3 = _WENO_OPTIMAL_WEIGHT_3_;
          } else {
            /* CRWENO5 at the interior points */
            c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
            c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
            c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
          }
        } else {
          /* WENO5 and HCWENO5 */
          c1 = _WENO_OPTIMAL_WEIGHT_1_;
          c2 = _WENO_OPTIMAL_WEIGHT_2_;
          c3 = _WENO_OPTIMAL_WEIGHT_3_;
        }

        /* save the weights */
        *(ww1+p*nvars+v) = c1;
        *(ww2+p*nvars+v) = c2;
        *(ww3+p*nvars+v) = c3;
      }
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  /* calculate weights for a right-biased interpolation */
  ww1 = a_w1 + 2*weno->size + offset;
  ww2 = a_w2 + 2*weno->size + offset;
  ww3 = a_w3 + 2*weno->size + offset;
  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexI,ndims);
    for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
      int p, v;
      _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
      for (v=0; v<nvars; v++)  {

        /* optimal weights*/
        double c1, c2, c3;
        if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
          if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
              || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
            /* Use WENO5 at the physical boundaries */
            c1 = _WENO_OPTIMAL_WEIGHT_1_;
            c2 = _WENO_OPTIMAL_WEIGHT_2_;
            c3 = _WENO_OPTIMAL_WEIGHT_3_;
          } else {
            /* CRWENO5 at the interior points */
            c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
            c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
            c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
          }
        } else {
          /* WENO5 and HCWENO5 */
          c1 = _WENO_OPTIMAL_WEIGHT_1_;
          c2 = _WENO_OPTIMAL_WEIGHT_2_;
          c3 = _WENO_OPTIMAL_WEIGHT_3_;
        }

        /* save the weights */
        *(ww1+p*nvars+v) = c1;
        *(ww2+p*nvars+v) = c2;
        *(ww3+p*nvars+v) = c3;
      }
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  ww1 = a_w1 + 2*weno->size + weno->size + offset;
  ww2 = a_w2 + 2*weno->size + weno->size + offset;
  ww3 = a_w3 + 2*weno->size + weno->size + offset;
  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexI,ndims);
    for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
      int p, v;
      _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
      for (v=0; v<nvars; v++)  {

        /* optimal weights*/
        double c1, c2, c3;
        if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
          if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
              || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
            /* Use WENO5 at the physical boundaries */
            c1 = _WENO_OPTIMAL_WEIGHT_1_;
            c2 = _WENO_OPTIMAL_WEIGHT_2_;
            c3 = _WENO_OPTIMAL_WEIGHT_3_;
          } else {
            /* CRWENO5 at the interior points */
            c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
            c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
            c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
          }
        } else {
          /* WENO5 and HCWENO5 */
          c1 = _WENO_OPTIMAL_WEIGHT_1_;
          c2 = _WENO_OPTIMAL_WEIGHT_2_;
          c3 = _WENO_OPTIMAL_WEIGHT_3_;
        }

        /* save the weights */
        *(ww1+p*nvars+v) = c1;
        *(ww2+p*nvars+v) = c2;
        *(ww3+p*nvars+v) = c3;
      }
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int WENOFifthOrderCalculateWeights	(	double *	fC,
		double *	uC,
		double *	x,
		int	dir,
		void *	s,
		void *	m
	)

Compute the nonlinear weights for 5th order WENO-type schemes. This function is a wrapper that calls the appropriate function, depending on the type of WENO weights.

Parameters

fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
uC	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables

Definition at line 40 of file WENOFifthOrderCalculateWeights.c.

{
  HyPar           *solver = (HyPar*)          s;
  WENOParameters  *weno   = (WENOParameters*) solver->interp;
  MPIVariables    *mpi    = (MPIVariables*)   m;

  int ret;

  if (weno->yc)           ret = WENOFifthOrderCalculateWeightsYC (fC,uC,x,dir,solver,mpi);
  else if (weno->borges)  ret = WENOFifthOrderCalculateWeightsZ  (fC,uC,x,dir,solver,mpi);
  else if (weno->mapped)  ret = WENOFifthOrderCalculateWeightsM  (fC,uC,x,dir,solver,mpi);
  else                    ret = WENOFifthOrderCalculateWeightsJS (fC,uC,x,dir,solver,mpi);

  return ret;
}

int WENOFifthOrderCalculateWeightsChar	(	double *	fC,
		double *	uC,
		double *	x,
		int	dir,
		void *	s,
		void *	m
	)

Compute the nonlinear weights for 5th order WENO-type schemes. This function is a wrapper that calls the appropriate function, depending on the type of WENO weights.

Parameters

fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
uC	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables

Definition at line 66 of file WENOFifthOrderCalculateWeights.c.

{
  HyPar           *solver = (HyPar*)          s;
  WENOParameters  *weno   = (WENOParameters*) solver->interp;
  MPIVariables    *mpi    = (MPIVariables*)   m;

  if (weno->yc)           return(WENOFifthOrderCalculateWeightsCharYC (fC,uC,x,dir,solver,mpi));
  else if (weno->borges)  return(WENOFifthOrderCalculateWeightsCharZ  (fC,uC,x,dir,solver,mpi));
  else if (weno->mapped)  return(WENOFifthOrderCalculateWeightsCharM  (fC,uC,x,dir,solver,mpi));
  else                    return(WENOFifthOrderCalculateWeightsCharJS (fC,uC,x,dir,solver,mpi));
}

int gpuWENOFifthOrderCalculateWeights	(	double *	fC,
		double *	uC,
		double *	x,
		int	dir,
		void *	s,
		void *	m
	)

Compute the nonlinear weights for 5th order WENO-type schemes. This function is a wrapper that calls the appropriate function, depending on the type of WENO weights.

Parameters

fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
uC	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables

Definition at line 30 of file WENOFifthOrderCalculateWeights_GPU.cu.

{
  HyPar           *solver = (HyPar*)          s;
  WENOParameters  *weno   = (WENOParameters*) solver->interp;
  MPIVariables    *mpi    = (MPIVariables*)   m;

  if (weno->yc)           return(gpuWENOFifthOrderCalculateWeightsYC (fC,uC,x,dir,solver,mpi));
  else if (weno->mapped)  return(gpuWENOFifthOrderCalculateWeightsM (fC,uC,x,dir,solver,mpi));
  else {
    printf("ERROR! WENO functions other than yc or mapped have not been implemented on GPUs.\n");
    return 1;
  }
}

int WENOInitialize	(	void *	s,
		void *	m,
		char *	scheme,
		char *	type
	)

This function initializes the WENO-type methods.

• Sets the parameters to default values.
• Reads in the parameters from optional input file "weno.inp", if available.
• Allocates memory for and initializes the nonlinear weights used by WENO-type schemes.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
scheme	Name of scheme
type	Type of interpolation

Definition at line 34 of file WENOInitialize.c.

{
  HyPar           *solver = (HyPar*) s;
  MPIVariables    *mpi    = (MPIVariables*) m;
  WENOParameters  *weno   = (WENOParameters*) solver->interp;

  static int count = 0;

  int nvars = solver->nvars;
  int ndims = solver->ndims;

  /* default parameters */
  weno->mapped      = 0;
  weno->borges      = 0;
  weno->yc          = 0;
  weno->no_limiting = 0;
  weno->eps         = 1e-6;
  weno->p           = 2.0;

  weno->rc          = 0.3;
  weno->xi          = 0.001;
  weno->tol         = 1e-16;

  if (!mpi->rank) {
    FILE *in;
    int ferr;
    in = fopen("weno.inp","r");
    if (!in) printf("Warning: File weno.inp not found. Using default parameters for WENO5/CRWENO5/HCWENO5 scheme.\n");
    else {
      if (!count) printf("Reading WENO parameters from weno.inp.\n");
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if      (!strcmp(word,"mapped"     )) { ferr = fscanf(in,"%d" ,&weno->mapped     ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"borges"     )) { ferr = fscanf(in,"%d" ,&weno->borges     ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"yc"         )) { ferr = fscanf(in,"%d" ,&weno->yc         ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"no_limiting")) { ferr = fscanf(in,"%d" ,&weno->no_limiting); if (ferr != 1) return(1); }
          else if (!strcmp(word,"epsilon"    )) { ferr = fscanf(in,"%lf",&weno->eps        ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"p"          )) { ferr = fscanf(in,"%lf",&weno->p          ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"rc"         )) { ferr = fscanf(in,"%lf",&weno->rc         ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"xi"         )) { ferr = fscanf(in,"%lf",&weno->xi         ); if (ferr != 1) return(1); }
          else if (!strcmp(word,"tol"        )) { ferr = fscanf(in,"%lf",&weno->tol        ); if (ferr != 1) return(1); }
          else if (strcmp(word,"end")) {
            char useless[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
            printf("Warning: keyword %s in file \"weno.inp\" with value %s not ",word,useless);
            printf("recognized or extraneous. Ignoring.\n");
          }
        }
        } else {
          fprintf(stderr,"Error: Illegal format in file \"weno.inp\".\n");
        return(1);
        }
      fclose(in);
    }
  }

  int     integer_data[4];
  double  real_data[5];
  if (!mpi->rank) {
    integer_data[0] = weno->mapped;
    integer_data[1] = weno->borges;
    integer_data[2] = weno->yc;
    integer_data[3] = weno->no_limiting;
    real_data[0]    = weno->eps;
    real_data[1]    = weno->p;
    real_data[2]    = weno->rc;
    real_data[3]    = weno->xi;
    real_data[4]    = weno->tol;
  }
  MPIBroadcast_integer(integer_data,4,0,&mpi->world);
  MPIBroadcast_double (real_data   ,5,0,&mpi->world);

  weno->mapped      = integer_data[0];
  weno->borges      = integer_data[1];
  weno->yc          = integer_data[2];
  weno->no_limiting = integer_data[3];
  weno->eps         = real_data   [0];
  weno->p           = real_data   [1];
  weno->rc          = real_data   [2];
  weno->xi          = real_data   [3];
  weno->tol         = real_data   [4];

  /* WENO weight calculation is hard-coded for p=2, so return error if p != 2 in
   * user input file, so that there's no confusion */
  if (weno->p != 2.0) {
    if (!mpi->rank && !count) printf("Warning from WENOInitialize(): \"p\" parameter is 2.0. Any other value will be ignored!\n");
  }

  weno->offset = NULL;
  weno->w1 = NULL;
  weno->w2 = NULL;
  weno->w3 = NULL;

  weno->offset = (int*) calloc (ndims,sizeof(int));
  int dir,d;
  for (dir=0; dir<ndims; dir++) {
    weno->offset[dir] = 0;
    for (d=0; d<dir; d++) {
      int size = nvars, i;
      for (i=0; i<ndims; i++)
        size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
      weno->offset[dir] += size;
    }
  }

  int total_size = 0;
  for (d=0; d<ndims; d++) {
    int size = nvars, i;
    for (i=0; i<ndims; i++)
      size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
    total_size += size;
  }
  weno->size = total_size;

  if ((!strcmp(type,_CHARACTERISTIC_)) && (nvars > 1))
    solver->SetInterpLimiterVar = WENOFifthOrderCalculateWeightsChar;
  else {
#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      solver->SetInterpLimiterVar = gpuWENOFifthOrderCalculateWeights;
    } else {
#endif
      solver->SetInterpLimiterVar = WENOFifthOrderCalculateWeights;
#if defined(HAVE_CUDA)
    }
#endif
  }

  /* initialize WENO weights to their optimal values */
  double* tmp_w1 = (double*) calloc (4*total_size,sizeof(double));
  double* tmp_w2 = (double*) calloc (4*total_size,sizeof(double));
  double* tmp_w3 = (double*) calloc (4*total_size,sizeof(double));
  for (d=0; d<ndims; d++) WENOFifthOrderInitializeWeights(  tmp_w1, 
                                                            tmp_w2, 
                                                            tmp_w3, 
                                                            weno->offset, 
                                                            d,
                                                            solver,
                                                            mpi);
  count++;

#if defined(HAVE_CUDA)
  if (solver->use_gpu) {
    //gpuMalloc((void**)&weno->gpu_offset, solver->ndims*sizeof(int));
    gpuMalloc((void**)&weno->w1, 4*total_size*sizeof(double));
    gpuMalloc((void**)&weno->w2, 4*total_size*sizeof(double));
    gpuMalloc((void**)&weno->w3, 4*total_size*sizeof(double));

    //gpuMemcpy(weno->gpu_offset, weno->offset, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
    gpuMemcpy(weno->w1, tmp_w1, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
    gpuMemcpy(weno->w2, tmp_w2, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
    gpuMemcpy(weno->w3, tmp_w3, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
  } else {
#endif
    weno->w1 = (double*) calloc (4*total_size,sizeof(double));
    weno->w2 = (double*) calloc (4*total_size,sizeof(double));
    weno->w3 = (double*) calloc (4*total_size,sizeof(double));

    _ArrayCopy1D_(tmp_w1, weno->w1, 4*total_size);
    _ArrayCopy1D_(tmp_w2, weno->w2, 4*total_size);
    _ArrayCopy1D_(tmp_w3, weno->w3, 4*total_size);
#if defined(HAVE_CUDA)
  }
#endif

  free(tmp_w1);
  free(tmp_w2);
  free(tmp_w3);

  return 0;
}