HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
FirstDerivativeFourthOrder.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <basic.h>
9 #include <arrayfunctions.h>
10 #include <firstderivative.h>
11 
12 #include <mpivars.h>
13 #include <hypar.h>
16 
17 #ifdef with_omp
18 #include <omp.h>
19 #endif
20 
37  double *Df,
38  double *f,
40  int dir,
41  int bias,
43  void *s,
44  void *m
45  )
46 {
47  SolverContext *solver = (SolverContext*) s;
48  int i, j, v;
49 
50  int ghosts = solver->ghosts;
51  int ndims = solver->ndims;
52  int nvars = solver->nvars;
53  int *dim = solver->dim_local;
54 
55 
56  if ((!Df) || (!f)) {
57  fprintf(stderr, "Error in FirstDerivativeFourthOrder(): input arrays not allocated.\n");
58  return(1);
59  }
60 
61  static double one_twelve = 1.0/12.0;
62 
63  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
64  dimension "dir" */
65  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
66  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
67  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
68 
69 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
70  for (j=0; j<N_outer; j++) {
71  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
72  _ArrayCopy1D_(index_outer,indexC,ndims);
73  /* left boundary */
74  for (i = -ghosts; i < -ghosts+1; i++) {
75  int qC, qp1, qp2, qp3, qp4;
76  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
77  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
78  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
79  indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
80  indexC[dir] = i+4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp4);
81  for (v=0; v<nvars; v++)
82  Df[qC*nvars+v] = (-25*f[qC*nvars+v]+48*f[qp1*nvars+v]-36*f[qp2*nvars+v]+16*f[qp3*nvars+v]-3*f[qp4*nvars+v])*one_twelve;
83  }
84  for (i = -ghosts+1; i < -ghosts+2; i++) {
85  int qC, qm1, qp1, qp2, qp3;
86  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
87  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
88  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
89  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
90  indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
91  for (v=0; v<nvars; v++)
92  Df[qC*nvars+v] = (-3*f[qm1*nvars+v]-10*f[qC*nvars+v]+18*f[qp1*nvars+v]-6*f[qp2*nvars+v]+f[qp3*nvars+v])*one_twelve;
93  }
94  /* interior */
95  for (i = -ghosts+2; i < dim[dir]+ghosts-2; i++) {
96  int qC, qm1, qm2, qp1, qp2;
97  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
98  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
99  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
100  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
101  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
102  for (v=0; v<nvars; v++)
103  Df[qC*nvars+v] = (f[qm2*nvars+v]-8*f[qm1*nvars+v]+8*f[qp1*nvars+v]-f[qp2*nvars+v])*one_twelve;
104  }
105  /* right boundary */
106  for (i = dim[dir]+ghosts-2; i < dim[dir]+ghosts-1; i++) {
107  int qC, qm3, qm2, qm1, qp1;
108  indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
109  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
110  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
111  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
112  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
113  for (v=0; v<nvars; v++)
114  Df[qC*nvars+v] = (-f[qm3*nvars+v]+6*f[qm2*nvars+v]-18*f[qm1*nvars+v]+10*f[qC*nvars+v]+3*f[qp1*nvars+v])*one_twelve;
115  }
116  for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
117  int qC, qm4, qm3, qm2, qm1;
118  indexC[dir] = i-4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm4);
119  indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
120  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
121  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
122  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
123  for (v=0; v<nvars; v++)
124  Df[qC*nvars+v] = (3*f[qm4*nvars+v]-16*f[qm3*nvars+v]+36*f[qm2*nvars+v]-48*f[qm1*nvars+v]+25*f[qC*nvars+v])*one_twelve;
125  }
126  }
127 
128  return(0);
129 }
HyPar SolverContext
Definitions for the functions computing the first derivative.
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int FirstDerivativeFourthOrderCentral(double *, double *, int, int, void *, void *)
int * dim_local
Definition: hypar.h:37
MPI related function definitions.
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
Contains structure definition for hypar.
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
#define _ArrayProduct1D_(x, size, p)
Structure of MPI-related variables.
MPIVariables MPIContext
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23