HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
InitialSolution.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9 #include <math.h>
10 #include <basic.h>
11 #include <common.h>
12 #if defined(HAVE_CUDA)
13 #include <arrayfunctions_gpu.h>
14 #else
15 #include <arrayfunctions.h>
16 #endif
17 #include <io.h>
18 #include <mpivars.h>
19 #include <simulation_object.h>
20 
21 int VolumeIntegral(double*,double*,void*,void*);
22 
25 int InitialSolution ( void *s,
26  int nsims
27  )
28 {
29  SimulationObject* simobj = (SimulationObject*) s;
30  int n, flag, d, i, offset, ierr;
31 
32  for (n = 0; n < nsims; n++) {
33 
34  int ghosts = simobj[n].solver.ghosts;
35 
36  char fname_root[_MAX_STRING_SIZE_] = "initial";
37  if (nsims > 1) {
38  char index[_MAX_STRING_SIZE_];
39  GetStringFromInteger(n, index, (int)log10(nsims)+1);
40  strcat(fname_root, "_");
41  strcat(fname_root, index);
42  }
43 
44  ierr = ReadArray( simobj[n].solver.ndims,
45  simobj[n].solver.nvars,
46  simobj[n].solver.dim_global,
47  simobj[n].solver.dim_local,
48  simobj[n].solver.ghosts,
49  &(simobj[n].solver),
50  &(simobj[n].mpi),
51  simobj[n].solver.x,
52  simobj[n].solver.u,
53  fname_root,
54  &flag );
55  if (ierr) {
56  fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
57  simobj[n].mpi.rank);
58  return ierr;
59  }
60  if (!flag) {
61  fprintf(stderr,"Error: initial solution file not found.\n");
62  return(1);
63  }
64  CHECKERR(ierr);
65 
66  /* exchange MPI-boundary values of u between processors */
67  MPIExchangeBoundariesnD( simobj[n].solver.ndims,
68  simobj[n].solver.nvars,
69  simobj[n].solver.dim_local,
70  simobj[n].solver.ghosts,
71  &(simobj[n].mpi),
72  simobj[n].solver.u );
73 
74  /* calculate dxinv */
75  offset = 0;
76  for (d = 0; d < simobj[n].solver.ndims; d++) {
77  for (i = 0; i < simobj[n].solver.dim_local[d]; i++) {
78  simobj[n].solver.dxinv[i+offset+ghosts]
79  = 2.0 / (simobj[n].solver.x[i+1+offset+ghosts]-simobj[n].solver.x[i-1+offset+ghosts]);
80  }
81  offset += (simobj[n].solver.dim_local[d] + 2*ghosts);
82  }
83 
84  /* exchange MPI-boundary values of dxinv between processors */
85  offset = 0;
86  for (d = 0; d < simobj[n].solver.ndims; d++) {
87  ierr = MPIExchangeBoundaries1D( &(simobj[n].mpi),
88  &(simobj[n].solver.dxinv[offset]),
89  simobj[n].solver.dim_local[d],
90  ghosts,
91  d,
92  simobj[n].solver.ndims ); CHECKERR(ierr);
93  if (ierr) {
94  fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
95  simobj[n].mpi.rank);
96  return ierr;
97  }
98  offset += (simobj[n].solver.dim_local[d] + 2*ghosts);
99  }
100 
101  /* fill in ghost values of dxinv at physical boundaries by extrapolation */
102  offset = 0;
103  for (d = 0; d < simobj[n].solver.ndims; d++) {
104  double *dxinv = &(simobj[n].solver.dxinv[offset]);
105  int *dim = simobj[n].solver.dim_local;
106  if (simobj[n].mpi.ip[d] == 0) {
107  /* fill left boundary along this dimension */
108  for (i = 0; i < ghosts; i++) dxinv[i] = dxinv[ghosts];
109  }
110  if (simobj[n].mpi.ip[d] == simobj[n].mpi.iproc[d]-1) {
111  /* fill right boundary along this dimension */
112  for (i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) dxinv[i] = dxinv[dim[d]+ghosts-1];
113  }
114  offset += (dim[d] + 2*ghosts);
115  }
116 
117  /* calculate volume integral of the initial solution */
118  ierr = VolumeIntegral( simobj[n].solver.VolumeIntegralInitial,
119  simobj[n].solver.u,
120  &(simobj[n].solver),
121  &(simobj[n].mpi) ); CHECKERR(ierr);
122  if (ierr) {
123  fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
124  simobj[n].mpi.rank);
125  return ierr;
126  }
127  if (!simobj[n].mpi.rank) {
128  if (nsims > 1) printf("Volume integral of the initial solution on domain %d:\n", n);
129  else printf("Volume integral of the initial solution:\n");
130  for (d=0; d<simobj[n].solver.nvars; d++) {
131  printf("%2d: %1.16E\n",d,simobj[n].solver.VolumeIntegralInitial[d]);
132  }
133  }
134  /* Set initial total boundary flux integral to zero */
135  _ArraySetValue_(simobj[n].solver.TotalBoundaryIntegral,simobj[n].solver.nvars,0);
136 
137  }
138 
139 #if defined(HAVE_CUDA)
140  if (simobj[0].solver.use_gpu) {
141  for (int n = 0; n < nsims; n++) {
142  int npoints_local_wghosts = simobj[n].solver.npoints_local_wghosts;
143  int nvars = simobj[n].solver.nvars;
144  int size_x = simobj[n].solver.size_x;
145 
146  gpuMemcpy(simobj[n].solver.gpu_x,
147  simobj[n].solver.x, simobj[n].solver.size_x*sizeof(double),
149  gpuMemcpy(simobj[n].solver.gpu_dxinv, simobj[n].solver.dxinv,
150  simobj[n].solver.size_x*sizeof(double),
152 
153 #ifdef CUDA_VAR_ORDERDING_AOS
154  gpuMemcpy(simobj[n].solver.gpu_u,
155  simobj[n].solver.u,
156  simobj[n].solver.ndof_cells_wghosts*sizeof(double),
158 #else
159  double *h_u = (double *) malloc(simobj[n].solver.ndof_cells_wghosts*sizeof(double));
160  for (int i=0; i<npoints_local_wghosts; i++) {
161  for (int v=0; v<nvars; v++) {
162  h_u[i+v*npoints_local_wghosts] = simobj[n].solver.u[i*nvars+v];
163  }
164  }
165  gpuMemcpy(simobj[n].solver.gpu_u,
166  h_u,
167  simobj[n].solver.ndof_cells_wghosts*sizeof(double),
169  free(h_u);
170 #endif
171  }
172  }
173 #endif
174 
175  return 0;
176 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int * dim_global
Definition: hypar.h:33
double * u
Definition: hypar.h:116
void GetStringFromInteger(int, char *, int)
int MPIExchangeBoundaries1D(void *, double *, int, int, int, int)
int * dim_local
Definition: hypar.h:37
void gpuMemcpy(void *, const void *, size_t, enum gpuMemcpyKind)
MPI related function definitions.
Function declarations for file I/O functions.
int ghosts
Definition: hypar.h:52
int InitialSolution(void *, int)
int ReadArray(int, int, int *, int *, int, void *, void *, double *, double *, char *, int *)
Definition: ReadArray.c:25
int size_x
Definition: hypar.h:48
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int ndof_cells_wghosts
Definition: hypar.h:47
Contains function definitions for common array operations on GPU.
Simulation object.
int nvars
Definition: hypar.h:29
#define CHECKERR(ierr)
Definition: basic.h:18
double * dxinv
Definition: hypar.h:110
Structure defining a simulation.
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
double * x
Definition: hypar.h:107
int VolumeIntegral(double *VolumeIntegral, double *u, void *s, void *m)
Some common functions used here and there.