HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
WENOInitialize.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9 #if defined(HAVE_CUDA)
10 #include <arrayfunctions_gpu.h>
11 #else
12 #include <arrayfunctions.h>
13 #endif
14 #include <interpolation.h>
15 #include <mpivars.h>
16 #include <hypar.h>
17 
18 int WENOFifthOrderInitializeWeights (double* const,double* const,double* const,
19  const int* const, int,void*,void*);
20 int WENOFifthOrderCalculateWeights (double*,double*,double*,int,void*,void*);
21 int WENOFifthOrderCalculateWeightsChar(double*,double*,double*,int,void*,void*);
22 
23 #if defined(HAVE_CUDA)
24 int gpuWENOFifthOrderCalculateWeights (double*,double*,double*,int,void*,void*);
25 #endif
26 
35  void *s,
36  void *m,
37  char *scheme,
38  char *type
39  )
40 {
41  HyPar *solver = (HyPar*) s;
42  MPIVariables *mpi = (MPIVariables*) m;
43  WENOParameters *weno = (WENOParameters*) solver->interp;
44 
45  static int count = 0;
46 
47  int nvars = solver->nvars;
48  int ndims = solver->ndims;
49 
50  /* default parameters */
51  weno->mapped = 0;
52  weno->borges = 0;
53  weno->yc = 0;
54  weno->no_limiting = 0;
55  weno->eps = 1e-6;
56  weno->p = 2.0;
57 
58  weno->rc = 0.3;
59  weno->xi = 0.001;
60  weno->tol = 1e-16;
61 
62  if (!mpi->rank) {
63  FILE *in;
64  int ferr;
65  in = fopen("weno.inp","r");
66  if (!in) printf("Warning: File weno.inp not found. Using default parameters for WENO5/CRWENO5/HCWENO5 scheme.\n");
67  else {
68  if (!count) printf("Reading WENO parameters from weno.inp.\n");
69  char word[_MAX_STRING_SIZE_];
70  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
71  if (!strcmp(word, "begin")){
72  while (strcmp(word, "end")){
73  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
74  if (!strcmp(word,"mapped" )) { ferr = fscanf(in,"%d" ,&weno->mapped ); if (ferr != 1) return(1); }
75  else if (!strcmp(word,"borges" )) { ferr = fscanf(in,"%d" ,&weno->borges ); if (ferr != 1) return(1); }
76  else if (!strcmp(word,"yc" )) { ferr = fscanf(in,"%d" ,&weno->yc ); if (ferr != 1) return(1); }
77  else if (!strcmp(word,"no_limiting")) { ferr = fscanf(in,"%d" ,&weno->no_limiting); if (ferr != 1) return(1); }
78  else if (!strcmp(word,"epsilon" )) { ferr = fscanf(in,"%lf",&weno->eps ); if (ferr != 1) return(1); }
79  else if (!strcmp(word,"p" )) { ferr = fscanf(in,"%lf",&weno->p ); if (ferr != 1) return(1); }
80  else if (!strcmp(word,"rc" )) { ferr = fscanf(in,"%lf",&weno->rc ); if (ferr != 1) return(1); }
81  else if (!strcmp(word,"xi" )) { ferr = fscanf(in,"%lf",&weno->xi ); if (ferr != 1) return(1); }
82  else if (!strcmp(word,"tol" )) { ferr = fscanf(in,"%lf",&weno->tol ); if (ferr != 1) return(1); }
83  else if (strcmp(word,"end")) {
84  char useless[_MAX_STRING_SIZE_];
85  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
86  printf("Warning: keyword %s in file \"weno.inp\" with value %s not ",word,useless);
87  printf("recognized or extraneous. Ignoring.\n");
88  }
89  }
90  } else {
91  fprintf(stderr,"Error: Illegal format in file \"weno.inp\".\n");
92  return(1);
93  }
94  fclose(in);
95  }
96  }
97 
98  int integer_data[4];
99  double real_data[5];
100  if (!mpi->rank) {
101  integer_data[0] = weno->mapped;
102  integer_data[1] = weno->borges;
103  integer_data[2] = weno->yc;
104  integer_data[3] = weno->no_limiting;
105  real_data[0] = weno->eps;
106  real_data[1] = weno->p;
107  real_data[2] = weno->rc;
108  real_data[3] = weno->xi;
109  real_data[4] = weno->tol;
110  }
111  MPIBroadcast_integer(integer_data,4,0,&mpi->world);
112  MPIBroadcast_double (real_data ,5,0,&mpi->world);
113 
114  weno->mapped = integer_data[0];
115  weno->borges = integer_data[1];
116  weno->yc = integer_data[2];
117  weno->no_limiting = integer_data[3];
118  weno->eps = real_data [0];
119  weno->p = real_data [1];
120  weno->rc = real_data [2];
121  weno->xi = real_data [3];
122  weno->tol = real_data [4];
123 
124  /* WENO weight calculation is hard-coded for p=2, so return error if p != 2 in
125  * user input file, so that there's no confusion */
126  if (weno->p != 2.0) {
127  if (!mpi->rank && !count) printf("Warning from WENOInitialize(): \"p\" parameter is 2.0. Any other value will be ignored!\n");
128  }
129 
130  weno->offset = NULL;
131  weno->w1 = NULL;
132  weno->w2 = NULL;
133  weno->w3 = NULL;
134 
135  weno->offset = (int*) calloc (ndims,sizeof(int));
136  int dir,d;
137  for (dir=0; dir<ndims; dir++) {
138  weno->offset[dir] = 0;
139  for (d=0; d<dir; d++) {
140  int size = nvars, i;
141  for (i=0; i<ndims; i++)
142  size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
143  weno->offset[dir] += size;
144  }
145  }
146 
147  int total_size = 0;
148  for (d=0; d<ndims; d++) {
149  int size = nvars, i;
150  for (i=0; i<ndims; i++)
151  size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
152  total_size += size;
153  }
154  weno->size = total_size;
155 
156  if ((!strcmp(type,_CHARACTERISTIC_)) && (nvars > 1))
158  else {
159 #if defined(HAVE_CUDA)
160  if (solver->use_gpu) {
162  } else {
163 #endif
165 #if defined(HAVE_CUDA)
166  }
167 #endif
168  }
169 
170  /* initialize WENO weights to their optimal values */
171  double* tmp_w1 = (double*) calloc (4*total_size,sizeof(double));
172  double* tmp_w2 = (double*) calloc (4*total_size,sizeof(double));
173  double* tmp_w3 = (double*) calloc (4*total_size,sizeof(double));
174  for (d=0; d<ndims; d++) WENOFifthOrderInitializeWeights( tmp_w1,
175  tmp_w2,
176  tmp_w3,
177  weno->offset,
178  d,
179  solver,
180  mpi);
181  count++;
182 
183 #if defined(HAVE_CUDA)
184  if (solver->use_gpu) {
185  //gpuMalloc((void**)&weno->gpu_offset, solver->ndims*sizeof(int));
186  gpuMalloc((void**)&weno->w1, 4*total_size*sizeof(double));
187  gpuMalloc((void**)&weno->w2, 4*total_size*sizeof(double));
188  gpuMalloc((void**)&weno->w3, 4*total_size*sizeof(double));
189 
190  //gpuMemcpy(weno->gpu_offset, weno->offset, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
191  gpuMemcpy(weno->w1, tmp_w1, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
192  gpuMemcpy(weno->w2, tmp_w2, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
193  gpuMemcpy(weno->w3, tmp_w3, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
194  } else {
195 #endif
196  weno->w1 = (double*) calloc (4*total_size,sizeof(double));
197  weno->w2 = (double*) calloc (4*total_size,sizeof(double));
198  weno->w3 = (double*) calloc (4*total_size,sizeof(double));
199 
200  _ArrayCopy1D_(tmp_w1, weno->w1, 4*total_size);
201  _ArrayCopy1D_(tmp_w2, weno->w2, 4*total_size);
202  _ArrayCopy1D_(tmp_w3, weno->w3, 4*total_size);
203 #if defined(HAVE_CUDA)
204  }
205 #endif
206 
207  free(tmp_w1);
208  free(tmp_w2);
209  free(tmp_w3);
210 
211  return 0;
212 }
void * interp
Definition: hypar.h:362
int gpuWENOFifthOrderCalculateWeights(double *fC, double *uC, double *x, int dir, void *s, void *m)
#define _CHARACTERISTIC_
Definition: interpolation.h:33
int WENOFifthOrderCalculateWeightsChar(double *fC, double *uC, double *x, int dir, void *s, void *m)
Definitions for the functions computing the interpolated value of the primitive at the cell interface...
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
int WENOFifthOrderCalculateWeights(double *fC, double *uC, double *x, int dir, void *s, void *m)
int * dim_local
Definition: hypar.h:37
void gpuMemcpy(void *, const void *, size_t, enum gpuMemcpyKind)
MPI related function definitions.
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)
Structure of variables/parameters needed by the WENO-type scheme.
int(* SetInterpLimiterVar)(double *, double *, double *, int, void *, void *)
Definition: hypar.h:234
MPI_Comm world
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int WENOInitialize(void *, void *, char *, char *)
Contains function definitions for common array operations on GPU.
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
void gpuMalloc(void **, size_t)
Contains structure definition for hypar.
int use_gpu
Definition: hypar.h:449
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23