HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TimeInitialize.c File Reference

Initialize time integration. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions_gpu.h>
#include <simulation_object.h>
#include <timeintegration.h>

Go to the source code of this file.

Functions

int TimeRHSFunctionExplicit (double *, double *, void *, void *, double)
 
int TimeInitialize (void *s, int nsims, int rank, int nproc, void *ts)
 

Detailed Description

Initialize time integration.

Author
Debojyoti Ghosh, Youngdae Kim

Definition in file TimeInitialize.c.

Function Documentation

int TimeRHSFunctionExplicit ( double *  rhs,
double *  u,
void *  s,
void *  m,
double  t 
)

This function computes the right-hand-side of the ODE given by

\begin{equation} \frac {{\bf u}}{dt} = {\bf F}\left({\bf u}\right) \end{equation}

for explicit time integration methods, i.e., where

\begin{equation} {\bf F}\left({\bf u}\right) = - {\bf F}_{\rm hyperbolic}\left({\bf u}\right) + {\bf F}_{\rm parabolic} \left({\bf u}\right) + {\bf F}_{\rm source} \left({\bf u}\right), \end{equation}

given the solution \({\bf u}\) and the current simulation time.

Parameters
rhsArray to hold the computed right-hand-side
uArray holding the solution
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time

Definition at line 30 of file TimeRHSFunctionExplicit.c.

37 {
38  HyPar *solver = (HyPar*) s;
39  MPIVariables *mpi = (MPIVariables*) m;
40  int d;
41 
42  int size = 1;
43  for (d=0; d<solver->ndims; d++) size *= (solver->dim_local[d]+2*solver->ghosts);
44 
45  /* apply boundary conditions and exchange data over MPI interfaces */
46  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
47  solver->ApplyIBConditions(solver,mpi,u,t);
48 
49  /* Evaluate hyperbolic, parabolic and source terms and the RHS */
50 #if defined(HAVE_CUDA)
51  if (solver->use_gpu) {
53  solver->nvars,
54  solver->gpu_dim_local,
55  solver->ghosts,
56  mpi,
57  u );
58  } else {
59 #endif
61  solver->nvars,
62  solver->dim_local,
63  solver->ghosts,
64  mpi,
65  u);
66 #if defined(HAVE_CUDA)
67  }
68 #endif
69 
70  solver->HyperbolicFunction( solver->hyp,
71  u,
72  solver,
73  mpi,
74  t,
75  1,
76  solver->FFunction,
77  solver->Upwind );
78  solver->ParabolicFunction(solver->par,u,solver,mpi,t);
79  solver->SourceFunction(solver->source,u,solver,mpi,t);
80 
81 #if defined(HAVE_CUDA)
82  if (solver->use_gpu) {
83  gpuArraySetValue(rhs, size*solver->nvars, 0.0);
84  gpuArrayAXPY(solver->hyp, -1.0, rhs, size*solver->nvars);
85  gpuArrayAXPY(solver->par, 1.0, rhs, size*solver->nvars);
86  gpuArrayAXPY(solver->source, 1.0, rhs, size*solver->nvars);
87  } else {
88 #endif
89  _ArraySetValue_(rhs,size*solver->nvars,0.0);
90  _ArrayAXPY_(solver->hyp ,-1.0,rhs,size*solver->nvars);
91  _ArrayAXPY_(solver->par , 1.0,rhs,size*solver->nvars);
92  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
93 #if defined(HAVE_CUDA)
94  }
95 #endif
96 
97  return(0);
98 }
double * hyp
Definition: hypar.h:119
void gpuArrayAXPY(const double *, double, double *, int)
#define _ArraySetValue_(x, size, value)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int * dim_local
Definition: hypar.h:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
void gpuArraySetValue(double *, int, double)
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
double * par
Definition: hypar.h:122
int * gpu_dim_local
Definition: hypar.h:455
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int nvars
Definition: hypar.h:29
int gpuMPIExchangeBoundariesnD(int, int, const int *, int, void *, double *)
int use_gpu
Definition: hypar.h:449
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int TimeInitialize ( void *  s,
int  nsims,
int  rank,
int  nproc,
void *  ts 
)

Initialize time integration: This function initializes all that is required for time integration.

  • It sets the number of iterations, time step size, simulation time, etc.
  • It allocates solution, right-hand-side, and stage solution arrays needed by specific time integration methods.
  • It calls the method-specific initialization functions.
Parameters
sArray of simulation objects of type SimulationObject
nsimsnumber of simulation objects
rankMPI rank of this process
nprocnumber of MPI processes
tsTime integration object of type TimeIntegration

Definition at line 28 of file TimeInitialize.c.

34 {
35  TimeIntegration* TS = (TimeIntegration*) ts;
37  int ns, d, i;
38 
39  if (!sim) return(1);
40 
41 
42  TS->simulation = sim;
43  TS->nsims = nsims;
44 
45  TS->rank = rank;
46  TS->nproc = nproc;
47 
48  TS->n_iter = sim[0].solver.n_iter;
49  TS->restart_iter = sim[0].solver.restart_iter;
50  TS->dt = sim[0].solver.dt;
51 
52  TS->waqt = (double) TS->restart_iter * TS->dt;
53  TS->max_cfl = 0.0;
54  TS->norm = 0.0;
55  TS->TimeIntegrate = sim[0].solver.TimeIntegrate;
56  TS->iter_wctime_total = 0.0;
57 
58  TS->u_offsets = (long*) calloc (nsims, sizeof(long));
59  TS->u_sizes = (long*) calloc (nsims, sizeof(long));
60 
61  TS->u_size_total = 0;
62  for (ns = 0; ns < nsims; ns++) {
63  TS->u_offsets[ns] = TS->u_size_total;
64  TS->u_sizes[ns] = sim[ns].solver.npoints_local_wghosts * sim[ns].solver.nvars ;
65  TS->u_size_total += TS->u_sizes[ns];
66  }
67 
68  TS->u = (double*) calloc (TS->u_size_total,sizeof(double));
69  TS->rhs = (double*) calloc (TS->u_size_total,sizeof(double));
70  _ArraySetValue_(TS->u ,TS->u_size_total,0.0);
71  _ArraySetValue_(TS->rhs,TS->u_size_total,0.0);
72 
73  /* initialize arrays to NULL, then allocate as necessary */
74  TS->U = NULL;
75  TS->Udot = NULL;
76  TS->BoundaryFlux = NULL;
77 
78  TS->bf_offsets = (int*) calloc (nsims, sizeof(int));
79  TS->bf_sizes = (int*) calloc (nsims, sizeof(int));
80  TS->bf_size_total = 0;
81  for (ns = 0; ns < nsims; ns++) {
82  TS->bf_offsets[ns] = TS->bf_size_total;
83  TS->bf_sizes[ns] = 2 * sim[ns].solver.ndims * sim[ns].solver.nvars;
84  TS->bf_size_total += TS->bf_sizes[ns];
85  }
86 
87 #if defined(HAVE_CUDA)
88  if (sim[0].solver.use_gpu) {
89 
90  if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
91 
92  /* explicit Runge-Kutta methods */
93  ExplicitRKParameters *params = (ExplicitRKParameters*) sim[0].solver.msti;
94  int nstages = params->nstages;
95 
96  gpuMalloc((void**)&TS->gpu_U, TS->u_size_total*nstages*sizeof(double));
97  gpuMalloc((void**)&TS->gpu_Udot, TS->u_size_total*nstages*sizeof(double));
98  gpuMalloc((void**)&TS->gpu_BoundaryFlux, TS->bf_size_total*nstages*sizeof(double));
99  gpuMemset(TS->gpu_U, 0, TS->u_size_total*nstages*sizeof(double));
100  gpuMemset(TS->gpu_Udot, 0, TS->u_size_total*nstages*sizeof(double));
101  gpuMemset(TS->gpu_BoundaryFlux, 0, TS->bf_size_total*nstages*sizeof(double));
102 
103  } else {
104 
105  fprintf(stderr,"ERROR in TimeInitialize(): %s is not yet implemented on GPUs.\n",
106  sim[0].solver.time_scheme );
107  return 1;
108 
109  }
110 
111  } else {
112 #endif
113 
114  if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
115 
116  /* explicit Runge-Kutta methods */
117  ExplicitRKParameters *params = (ExplicitRKParameters*) sim[0].solver.msti;
118  int nstages = params->nstages;
119  TS->U = (double**) calloc (nstages,sizeof(double*));
120  TS->Udot = (double**) calloc (nstages,sizeof(double*));
121  for (i = 0; i < nstages; i++) {
122  TS->U[i] = (double*) calloc (TS->u_size_total,sizeof(double));
123  TS->Udot[i] = (double*) calloc (TS->u_size_total,sizeof(double));
124  }
125 
126  TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
127  for (i=0; i<nstages; i++) {
128  TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
129  }
130 
131  } else if (!strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
132 
133  int nstages = 1;
134  TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
135  for (i=0; i<nstages; i++) {
136  TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
137  }
138 
139  } else if (!strcmp(sim[0].solver.time_scheme,_GLM_GEE_)) {
140 
141  /* General Linear Methods with Global Error Estimate */
142  GLMGEEParameters *params = (GLMGEEParameters*) sim[0].solver.msti;
143  int nstages = params->nstages;
144  int r = params->r;
145  TS->U = (double**) calloc (2*r-1 ,sizeof(double*));
146  TS->Udot = (double**) calloc (nstages,sizeof(double*));
147  for (i=0; i<2*r-1; i++) TS->U[i] = (double*) calloc (TS->u_size_total,sizeof(double));
148  for (i=0; i<nstages; i++) TS->Udot[i] = (double*) calloc (TS->u_size_total,sizeof(double));
149 
150  TS->BoundaryFlux = (double**) calloc (nstages,sizeof(double*));
151  for (i=0; i<nstages; i++) {
152  TS->BoundaryFlux[i] = (double*) calloc (TS->bf_size_total,sizeof(double));
153  }
154 
155 
156  if (!strcmp(params->ee_mode,_GLM_GEE_YYT_)) {
157  for (ns = 0; ns < nsims; ns++) {
158  for (i=0; i<r-1; i++) {
159  _ArrayCopy1D_( (sim[ns].solver.u),
160  (TS->U[r+i] + TS->u_offsets[ns]),
161  (TS->u_sizes[ns]) );
162  }
163  }
164  } else {
165  for (i=0; i<r-1; i++) {
166  _ArraySetValue_(TS->U[r+i],TS->u_size_total,0.0);
167  }
168  }
169  }
170 #if defined(HAVE_CUDA)
171  }
172 #endif
173 
174  /* set right-hand side function pointer */
176 
177  /* open files for writing */
178  if (!rank) {
179  if (sim[0].solver.write_residual) TS->ResidualFile = (void*) fopen("residual.out","w");
180  else TS->ResidualFile = NULL;
181  } else TS->ResidualFile = NULL;
182 
183  for (ns = 0; ns < nsims; ns++) {
184  sim[ns].solver.time_integrator = TS;
185  }
186 
187  return 0;
188 }
#define _FORWARD_EULER_
int(* RHSFunction)(double *, double *, void *, void *, double)
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
Structure of variables/parameters and function pointers for time integration.
double dt
Definition: hypar.h:67
#define _GLM_GEE_
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Structure containing the parameters for an explicit Runge-Kutta method.
#define _GLM_GEE_YYT_
#define _RK_
char ee_mode[_MAX_STRING_SIZE_]
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
void gpuMemset(void *, int, size_t)
void gpuMalloc(void **, size_t)
int n_iter
Definition: hypar.h:55
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* TimeIntegrate)(void *)
int(* TimeIntegrate)(void *)
Definition: hypar.h:220
int restart_iter
Definition: hypar.h:58
int TimeRHSFunctionExplicit(double *, double *, void *, void *, double)
void * time_integrator
Definition: hypar.h:165