HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimFifthOrderCompactUpwind.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <mathfunctions.h>
10 #include <interpolation.h>
11 #include <tridiagLU.h>
12 #include <mpivars.h>
13 #include <hypar.h>
14 
15 #ifdef with_omp
16 #include <omp.h>
17 #endif
18 
19 #undef _MINIMUM_GHOSTS_
20 
24 #define _MINIMUM_GHOSTS_ 3
25 
73  double *fI,
74  double *fC,
75  double *u,
76  double *x,
77  int upw,
78  int dir,
79  void *s,
80  void *m,
81  int uflag
82  )
83 {
84  HyPar *solver = (HyPar*) s;
85  MPIVariables *mpi = (MPIVariables*) m;
86  CompactScheme *compact= (CompactScheme*) solver->compact;
87  TridiagLU *lu = (TridiagLU*) solver->lusolver;
88  int sys,Nsys,d,v;
90 
91  int ghosts = solver->ghosts;
92  int ndims = solver->ndims;
93  int nvars = solver->nvars;
94  int *dim = solver->dim_local;
95  int *stride= solver->stride_with_ghosts;
96 
97  /* define some constants */
98  static const double one_third = 1.0/3.0,
99  thirteen_by_sixty = 13.0/60.0,
100  fortyseven_by_sixty = 47.0/60.0,
101  twentyseven_by_sixty = 27.0/60.0,
102  one_by_twenty = 1.0/20.0,
103  one_by_thirty = 1.0/30.0,
104  nineteen_by_thirty = 19.0/30.0,
105  three_by_ten = 3.0/10.0,
106  six_by_ten = 6.0/10.0,
107  one_by_ten = 1.0/10.0;
108 
109  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
110  dimension "dir" */
111  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
112  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
113  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
114  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
115 
116  /* calculate total number of tridiagonal systems to solve */
117  _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
118 
119  /* Allocate arrays for tridiagonal system */
120  double *A = compact->A;
121  double *B = compact->B;
122  double *C = compact->C;
123  double *R = compact->R;
124 
125 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
126  for (sys=0; sys < N_outer; sys++) {
127  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
128  _ArrayCopy1D_(index_outer,indexC,ndims);
129  _ArrayCopy1D_(index_outer,indexI,ndims);
130  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
131  int qm1,qm2,qm3,qp1,qp2,p;
132  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
133  if (upw > 0) {
134  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
135  qm3 = qm1 - 2*stride[dir];
136  qm2 = qm1 - stride[dir];
137  qp1 = qm1 + stride[dir];
138  qp2 = qm1 + 2*stride[dir];
139  } else {
140  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
141  qm3 = qm1 + 2*stride[dir];
142  qm2 = qm1 + stride[dir];
143  qp1 = qm1 - stride[dir];
144  qp2 = qm1 - 2*stride[dir];
145  }
146 
147  /* Defining stencil points */
148  double *fm3, *fm2, *fm1, *fp1, *fp2;
149  fm3 = fC+qm3*nvars;
150  fm2 = fC+qm2*nvars;
151  fm1 = fC+qm1*nvars;
152  fp1 = fC+qp1*nvars;
153  fp2 = fC+qp2*nvars;
154 
155  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
156  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
157 
158  /* Use 5th order upwind at the physical boundaries */
159  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
160  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,1.0)
161  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
162  for (v=0; v<nvars; v++) {
163  (R+Nsys*indexI[dir]+sys*nvars)[v] = one_by_thirty * fm3[v]
164  - thirteen_by_sixty * fm2[v]
165  + fortyseven_by_sixty * fm1[v]
166  + twentyseven_by_sixty * fp1[v]
167  - one_by_twenty * fp2[v];
168  }
169 
170  } else {
171 
172  /* Use 5th order upwind at the physical boundaries */
173  if (upw > 0) {
174  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
175  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten );
176  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten );
177  } else {
178  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
179  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten );
180  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten );
181  }
182  for (v=0; v<nvars; v++) {
183  (R+Nsys*indexI[dir]+sys*nvars)[v] = one_by_thirty * fm2[v]
184  + nineteen_by_thirty * fm1[v]
185  + one_third * fp1[v];
186  }
187  }
188  }
189  }
190 
191 #ifdef serial
192 
193  /* Solve the tridiagonal system */
194  IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
195 
196 #else
197 
198  /* Solve the tridiagonal system */
199  /* all processes except the last will solve without the last interface to avoid overlap */
200  if (mpi->ip[dir] != mpi->iproc[dir]-1) { IERR tridiagLU(A,B,C,R,dim[dir] ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
201  else { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
202 
203  /* Now get the solution to the last interface from the next proc */
204  double *sendbuf = compact->sendbuf;
205  double *recvbuf = compact->recvbuf;
206  MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
207  if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
208  if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
209  if (mpi->ip[dir]) MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
210  MPI_Status status_arr[2];
211  MPI_Waitall(2,&req[0],status_arr);
212  if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
213 
214 #endif
215 
216  /* save the solution to fI */
217 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
218  for (sys=0; sys < N_outer; sys++) {
219  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
220  _ArrayCopy1D_(index_outer,indexI,ndims);
221  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
222  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
223  _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
224  }
225  }
226 
227  return(0);
228 }
Header file for TridiagLU.
MPI_Comm * comm
#define _ArraySetValue_(x, size, value)
Definitions for the functions computing the interpolated value of the primitive at the cell interface...
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int tridiagLU(double *, double *, double *, double *, int, int, void *, void *)
Definition: tridiagLU.c:83
int * dim_local
Definition: hypar.h:37
MPI related function definitions.
double * recvbuf
Structure of variables/parameters needed by the compact schemes.
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int Interp1PrimFifthOrderCompactUpwind(double *, double *, double *, double *, int, int, void *, void *, int)
5th order compact upwind reconstruction (component-wise) on a uniform grid
void * compact
Definition: hypar.h:364
Contains function definitions for common mathematical functions.
int * stride_with_ghosts
Definition: hypar.h:414
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
#define CHECKERR(ierr)
Definition: basic.h:18
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 IERR
Definition: basic.h:16
#define _ArrayProduct1D_(x, size, p)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * sendbuf
void * lusolver
Definition: hypar.h:368
#define _DECLARE_IERR_
Definition: basic.h:17