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

Contains function declarations for time integration. More...

#include <timeintegration_struct.h>

Go to the source code of this file.

Functions

int TimeExplicitRKInitialize (char *, char *, void *, void *)
 
int TimeExplicitRKCleanup (void *)
 
int TimeGLMGEEInitialize (char *, char *, void *, void *)
 
int TimeGLMGEECleanup (void *)
 
int TimeInitialize (void *, int, int, int, void *)
 
int TimeCleanup (void *)
 
int TimePreStep (void *)
 
int TimeStep (void *)
 
int TimePostStep (void *)
 
int TimePrintStep (void *)
 
int TimeError (void *, void *, double *)
 
int TimeGetAuxSolutions (int *, double **, void *, int, int)
 
int TimeForwardEuler (void *)
 
int TimeRK (void *)
 
int TimeGLMGEE (void *)
 

Detailed Description

Contains function declarations for time integration.

Author
Debojyoti Ghosh

Definition in file timeintegration.h.

Function Documentation

int TimeExplicitRKInitialize ( char *  class,
char *  type,
void *  s,
void *  m 
)

Initialize the explicit Runge-Kutta time-integration method

Initialize the explicit Runge-Kutta time integrator: Depending on the specific explicit Runge-Kutta (ERK) method chosen, this function allocates memory for the Butcher tableaux and sets their coefficients.

Parameters
className of time integrator class; must match _RK_
typeType of explicit Runge-Kutta method
sObject of type ExplicitRKParameters
mUnused argument

Definition at line 17 of file TimeExplicitRKInitialize.c.

23 {
24  ExplicitRKParameters *params = (ExplicitRKParameters*) s;
25 
26  if (!strcmp(class,_RK_)) {
27  if (!strcmp(type,_RK_1FE_)) {
28  params->nstages = 1;
29  params->A = (double*) calloc (params->nstages*params->nstages,sizeof(double));
30  params->b = (double*) calloc (params->nstages ,sizeof(double));
31  params->c = (double*) calloc (params->nstages ,sizeof(double));
32  _ArraySetValue_(params->A,params->nstages*params->nstages,0.0);
33  _ArraySetValue_(params->b,params->nstages ,0.0);
34  _ArraySetValue_(params->c,params->nstages ,0.0);
35  params->b[0] = 1.0;
36  } else if (!strcmp(type,_RK_22_)) {
37  params->nstages = 2;
38  params->A = (double*) calloc (params->nstages*params->nstages,sizeof(double));
39  params->b = (double*) calloc (params->nstages ,sizeof(double));
40  params->c = (double*) calloc (params->nstages ,sizeof(double));
41  _ArraySetValue_(params->A,params->nstages*params->nstages,0.0);
42  _ArraySetValue_(params->b,params->nstages ,0.0);
43  _ArraySetValue_(params->c,params->nstages ,0.0);
44  params->A[2] = 1.0;;
45  params->c[1] = 1.0;;
46  params->b[0] = params->b[1] = 0.5;
47  } else if (!strcmp(type,_RK_33_)) {
48  params->nstages = 3;
49  params->A = (double*) calloc (params->nstages*params->nstages,sizeof(double));
50  params->b = (double*) calloc (params->nstages ,sizeof(double));
51  params->c = (double*) calloc (params->nstages ,sizeof(double));
52  _ArraySetValue_(params->A,params->nstages*params->nstages,0.0);
53  _ArraySetValue_(params->b,params->nstages ,0.0);
54  _ArraySetValue_(params->c,params->nstages ,0.0);
55  params->A[3] = 2.0/3.0; params->A[6] = 2.0/3.0-1.0/4.0; params->A[7] = 1.0/4.0;
56  params->c[1] = 2.0/3.0; params->c[2] = 2.0/3.0;
57  params->b[0] = 1.0/4.0; params->b[1] = -1.0/4.0; params->b[2] = 1.0;
58  } else if (!strcmp(type,_RK_44_)) {
59  params->nstages = 4;
60  params->A = (double*) calloc (params->nstages*params->nstages,sizeof(double));
61  params->b = (double*) calloc (params->nstages ,sizeof(double));
62  params->c = (double*) calloc (params->nstages ,sizeof(double));
63  _ArraySetValue_(params->A,params->nstages*params->nstages,0.0);
64  _ArraySetValue_(params->b,params->nstages ,0.0);
65  _ArraySetValue_(params->c,params->nstages ,0.0);
66  params->A[4] = 0.5; params->A[9] = 0.5; params->A[14] = 1.0;
67  params->c[1] = params->c[2] = 0.5; params->c[3] = 1.0;
68  params->b[0] = 1.0/6.0; params->b[1] = 1.0/3.0; params->b[2] = 1.0/3.0; params->b[3] = 1.0/6.0;
69  } else if ((!strcmp(type,_RK_SSP3_)) || (!strcmp(type,_RK_TVD3_))) {
70  params->nstages = 3;
71  params->A = (double*) calloc (params->nstages*params->nstages,sizeof(double));
72  params->b = (double*) calloc (params->nstages ,sizeof(double));
73  params->c = (double*) calloc (params->nstages ,sizeof(double));
74  _ArraySetValue_(params->A,params->nstages*params->nstages,0.0);
75  _ArraySetValue_(params->b,params->nstages ,0.0);
76  _ArraySetValue_(params->c,params->nstages ,0.0);
77  params->A[3] = 1.0; params->A[6] = 0.25; params->A[7] = 0.25;
78  params->c[1] = 1.0; params->c[2] = 0.5;
79  params->b[0] = params->b[1] = 1.0/6.0; params->b[2] = 2.0/3.0;
80  } else {
81  fprintf(stderr,"Error in TimeExplicitRKInitialize(): %s is not a supported ",type);
82  fprintf(stderr,"multi-stage time integration scheme of class %s.\n",class);
83  return(1);
84  }
85  } else {
86  fprintf(stderr,"Error in TimeExplicitRKInitialize(): Code should not have ");
87  fprintf(stderr,"reached here for %s class of time-integrators. This is a ",class);
88  fprintf(stderr,"coding mistake.\n");
89  return(1);
90  }
91  return(0);
92 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_RK_1FE_
#define _RK_1FE_
Definition: timeintegration_struct.h:128
_RK_TVD3_
#define _RK_TVD3_
Definition: timeintegration_struct.h:166
_RK_SSP3_
#define _RK_SSP3_
Definition: timeintegration_struct.h:165
ExplicitRKParameters
Structure containing the parameters for an explicit Runge-Kutta method.
Definition: timeintegration_struct.h:177
_RK_44_
#define _RK_44_
Definition: timeintegration_struct.h:152
_RK_
#define _RK_
Definition: timeintegration_struct.h:16
_RK_22_
#define _RK_22_
Definition: timeintegration_struct.h:136
_RK_33_
#define _RK_33_
Definition: timeintegration_struct.h:144
int TimeExplicitRKCleanup ( void *  s)

Clean up variables related to the explicit Runge-Kutta time-integration method

Clean up allocations related to explicit Runge-Kutta time integration: This function frees up the arrays for the Butcher tableaux.

Parameters
sObject of type ExplicitRKParameters

Definition at line 12 of file TimeExplicitRKCleanup.c.

13 {
14  ExplicitRKParameters *params = (ExplicitRKParameters*) s;
15  if (params->A) free(params->A);
16  if (params->b) free(params->b);
17  if (params->c) free(params->c);
18  return(0);
19 }
ExplicitRKParameters
Structure containing the parameters for an explicit Runge-Kutta method.
Definition: timeintegration_struct.h:177
int TimeGLMGEEInitialize ( char *  class,
char *  type,
void *  s,
void *  m 
)

Initialize the General Linear Methods with Global Error Estimators (GLM-GEE)

Initialize the GLM-GEE (_GLM_GEE_) time integation method: This function allocates the arrays to store the Butcher tableaux, and sets their coefficients, as well as other parameters for the GLM-GEE methods.

Reference:

Parameters
classClass of time integration method; must match _GLM_GEE_
typeName of the _GLM_GEE_ method to use
sObject of type GLMGEEParameters
mMPI object of type MPIVariables

Definition at line 24 of file TimeGLMGEEInitialize.c.

30 {
31  GLMGEEParameters *params = (GLMGEEParameters*) s;
32  MPIVariables *mpi = (MPIVariables*) m;
33  int i,j;
34 
35  if (!strcmp(class,_GLM_GEE_)) {
36 
37  if (!strcmp(type,_GLM_GEE_23_)) {
38  params->nstages = 3;
39  params->r = 2;
40  params->gamma = 0.0;
41  } else if (!strcmp(type,_GLM_GEE_24_)) {
42  params->nstages = 4;
43  params->r = 2;
44  params->gamma = 0.0;
45  } else if (!strcmp(type,_GLM_GEE_25I_)) {
46  params->nstages = 5;
47  params->r = 2;
48  params->gamma = 0.0;
49  } else if (!strcmp(type,_GLM_GEE_35_)) {
50  params->nstages = 5;
51  params->r = 2;
52  params->gamma = 0.0;
53  } else if (!strcmp(type,_GLM_GEE_EXRK2A_)) {
54  params->nstages = 6;
55  params->r = 2;
56  params->gamma = 0.25;
57  } else if (!strcmp(type,_GLM_GEE_RK32G1_)) {
58  params->nstages = 8;
59  params->r = 2;
60  params->gamma = 0.0;
61  } else if (!strcmp(type,_GLM_GEE_RK285EX_)) {
62  params->nstages = 9;
63  params->r = 2;
64  params->gamma = 0.25;
65  } else {
66  fprintf(stderr,"Error in TimeGLMGEEInitialize(): %s is not a supported ",type);
67  fprintf(stderr,"multi-stage time integration scheme of class %s.\n",class);
68  return(1);
69  }
70 
71  int s = params->nstages;
72  int r = params->r;
73 
74  params->A_yyt = (double*) calloc (s*s,sizeof(double));
75  params->B_yyt = (double*) calloc (s*r,sizeof(double));
76  params->C_yyt = (double*) calloc (s*r,sizeof(double));
77  params->D_yyt = (double*) calloc (r*r,sizeof(double));
78  params->c_yyt = (double*) calloc (s ,sizeof(double));
79  _ArraySetValue_(params->A_yyt,s*s,0.0);
80  _ArraySetValue_(params->B_yyt,s*r,0.0);
81  _ArraySetValue_(params->C_yyt,s*r,0.0);
82  _ArraySetValue_(params->D_yyt,r*r,0.0);
83  _ArraySetValue_(params->c_yyt,s ,0.0);
84 
85  params->A_yeps = (double*) calloc (s*s,sizeof(double));
86  params->B_yeps = (double*) calloc (s*r,sizeof(double));
87  params->C_yeps = (double*) calloc (s*r,sizeof(double));
88  params->D_yeps = (double*) calloc (r*r,sizeof(double));
89  params->c_yeps = (double*) calloc (s ,sizeof(double));
90  _ArraySetValue_(params->A_yeps,s*s,0.0);
91  _ArraySetValue_(params->B_yeps,s*r,0.0);
92  _ArraySetValue_(params->C_yeps,s*r,0.0);
93  _ArraySetValue_(params->D_yeps,r*r,0.0);
94  _ArraySetValue_(params->c_yeps,s ,0.0);
95 
96  if (!strcmp(type,_GLM_GEE_23_)) {
97 
98  params->A_yeps[1*s+0] = 1.0;
99  params->A_yeps[2*s+0] = params->A_yeps[2*s+1] = 0.25;
100 
101  params->B_yeps[0*s+0] = params->B_yeps[0*s+1] = 1.0/12.0; params->B_yeps[0*s+2] = 5.0/6.0;
102  params->B_yeps[1*s+0] = params->B_yeps[1*s+1] = 1.0/12.0; params->B_yeps[1*s+2] = -1.0/6.0;
103 
104  params->C_yeps[0*r+0] = 1.0;
105  params->C_yeps[1*r+0] = 1.0; params->C_yeps[1*r+1] = 10.0;
106  params->C_yeps[2*r+0] = 1.0; params->C_yeps[2*r+1] = -1.0;
107 
108  params->D_yeps[0*r+0] = 1.0;
109  params->D_yeps[1*r+1] = 1.0;
110 
111  params->A_yyt[1*s+0] = 1.0;
112  params->A_yyt[2*s+0] = params->A_yyt[2*s+1] = 0.25;
113 
114  params->B_yyt[0*s+0] = params->B_yyt[0*s+1] = 1.0/12.0; params->B_yyt[0*s+2] = 5.0/6.0;
115  params->B_yyt[1*s+0] = params->B_yyt[1*s+1] = 1.0/6.0; params->B_yyt[1*s+2] = 2.0/3.0;
116 
117  params->C_yyt[0*r+0] = 1.0;
118  params->C_yyt[1*r+0] = -9.0; params->C_yyt[1*r+1] = 10.0;
119  params->C_yyt[2*r+0] = 2.0; params->C_yyt[2*r+1] = -1.0;
120 
121  params->D_yyt[0*r+0] = 1.0;
122  params->D_yyt[1*r+1] = 1.0;
123 
124  } else if (!strcmp(type,_GLM_GEE_24_)) {
125 
126  params->A_yyt[1*s+0] = 0.75;
127  params->A_yyt[2*s+0] = 0.25;
128  params->A_yyt[2*s+1] = 29.0/60.0;
129  params->A_yyt[3*s+0] = -21.0/44.0;
130  params->A_yyt[3*s+1] = 145.0/44.0;
131  params->A_yyt[3*s+2] = -20.0/11.0;
132 
133  params->B_yyt[0*s+0] = 109.0/275.0;
134  params->B_yyt[0*s+1] = 58.0/75.0;
135  params->B_yyt[0*s+2] = -37.0/110.0;
136  params->B_yyt[0*s+3] = 1.0/6.0;
137  params->B_yyt[1*s+0] = 3.0/11.0;
138  params->B_yyt[1*s+2] = 75.0/88.0;
139  params->B_yyt[1*s+3] = -1.0/8.0;
140 
141  params->C_yyt[0*r+1] = 1.0;
142  params->C_yyt[1*r+0] = 75.0/58.0;
143  params->C_yyt[1*r+1] = -17.0/58.0;
144  params->C_yyt[2*r+1] = params->C_yyt[3*r+1] = 1.0;
145 
146  params->D_yyt[0*r+0] = 1.0;
147  params->D_yyt[1*r+1] = 1.0;
148 
149  double T[r*r],Tinv[r*r];
150  T[0*r+0] = 1.0;
151  T[0*r+1] = 0.0;
152  T[1*r+0] = 1.0;
153  T[1*r+1] = 1.0-params->gamma;
154  _MatrixInvert_(T,Tinv,r);
155 
156  _ArrayCopy1D_(params->A_yyt,params->A_yeps,(s*s));
157  _MatrixMultiplyNonSquare_(Tinv,params->B_yyt,params->B_yeps,r,r,s);
158  _MatrixMultiplyNonSquare_(params->C_yyt,T,params->C_yeps,s,r,r);
159  _ArrayCopy1D_(params->D_yyt,params->D_yeps,(r*r));
160 
161  } else if (!strcmp(type,_GLM_GEE_25I_)) {
162 
163  params->A_yyt[1*s+0]=-0.94079244066783383269;
164  params->A_yyt[2*s+0]= 0.64228187778301907108;
165  params->A_yyt[2*s+1]= 0.10915356933958500042;
166  params->A_yyt[3*s+0]=-0.51764297742287450812;
167  params->A_yyt[3*s+1]= 0.74414270351096040738;
168  params->A_yyt[3*s+2]=-0.71404164927824538121;
169  params->A_yyt[4*s+0]=-0.44696561556825969206;
170  params->A_yyt[4*s+1]=-0.76768425657590196518;
171  params->A_yyt[4*s+2]= 0.20111608138142987881;
172  params->A_yyt[4*s+3]= 0.93828186737840469796;
173 
174  params->B_yyt[0*s+0]=-0.029309178948150356153;
175  params->B_yyt[0*s+1]=-0.49671981884013874923;
176  params->B_yyt[0*s+2]= 0.34275801517650053274;
177  params->B_yyt[0*s+3]= 0.32941112623949194988;
178  params->B_yyt[0*s+4]= 0.85385985637229662276;
179  params->B_yyt[1*s+0]= 0.78133219686062535272;
180  params->B_yyt[1*s+1]= 0.074238691892675897635;
181  params->B_yyt[1*s+2]= 0.57957363498384957966;
182  params->B_yyt[1*s+3]=-0.24638502829674959968;
183  params->B_yyt[1*s+4]=-0.18875949544040123033;
184 
185  params->C_yyt[0*r+0]= 0.16911424754448327735;
186  params->C_yyt[0*r+1]= 0.83088575245551672265;
187  params->C_yyt[1*r+0]= 0.53638465733199574340;
188  params->C_yyt[1*r+1]= 0.46361534266800425660;
189  params->C_yyt[2*r+0]= 0.39901579167169582526;
190  params->C_yyt[2*r+1]= 0.60098420832830417474;
191  params->C_yyt[3*r+0]= 0.87689005530618575480;
192  params->C_yyt[3*r+1]= 0.12310994469381424520;
193  params->C_yyt[4*r+0]= 0.99056100455550913009;
194  params->C_yyt[4*r+1]= 0.0094389954444908699092;
195 
196  params->D_yyt[0*r+0] = 1.0;
197  params->D_yyt[1*r+1] = 1.0;
198 
199  double T[r*r],Tinv[r*r];
200  T[0*r+0] = 1.0;
201  T[0*r+1] = 0.0;
202  T[1*r+0] = 1.0;
203  T[1*r+1] = 1.0-params->gamma;
204  _MatrixInvert_(T,Tinv,r);
205 
206  _ArrayCopy1D_(params->A_yyt,params->A_yeps,(s*s));
207  _MatrixMultiplyNonSquare_(Tinv,params->B_yyt,params->B_yeps,r,r,s);
208  _MatrixMultiplyNonSquare_(params->C_yyt,T,params->C_yeps,s,r,r);
209  _ArrayCopy1D_(params->D_yyt,params->D_yeps,(r*r));
210 
211  } else if (!strcmp(type,_GLM_GEE_35_)) {
212 
213  params->A_yyt[1*s+0] = - 2169604947363702313.0 / 24313474998937147335.0;
214  params->A_yyt[2*s+0] = 46526746497697123895.0 / 94116917485856474137.0;
215  params->A_yyt[2*s+1] = -10297879244026594958.0 / 49199457603717988219.0;
216  params->A_yyt[3*s+0] = 23364788935845982499.0 / 87425311444725389446.0;
217  params->A_yyt[3*s+1] = -79205144337496116638.0 / 148994349441340815519.0;
218  params->A_yyt[3*s+2] = 40051189859317443782.0 / 36487615018004984309.0;
219  params->A_yyt[4*s+0] = 42089522664062539205.0 / 124911313006412840286.0;
220  params->A_yyt[4*s+1] = -15074384760342762939.0 / 137927286865289746282.0;
221  params->A_yyt[4*s+2] = -62274678522253371016.0 / 125918573676298591413.0;
222  params->A_yyt[4*s+3] = 13755475729852471739.0 / 79257927066651693390.0;
223 
224  params->B_yyt[0*s+0] = 61546696837458703723.0 / 56982519523786160813.0;
225  params->B_yyt[0*s+1] = -55810892792806293355.0 / 206957624151308356511.0;
226  params->B_yyt[0*s+2] = 24061048952676379087.0 / 158739347956038723465.0;
227  params->B_yyt[0*s+3] = 3577972206874351339.0 / 7599733370677197135.0;
228  params->B_yyt[0*s+4] = -59449832954780563947.0 / 137360038685338563670.0;
229  params->B_yyt[1*s+0] = - 9738262186984159168.0 / 99299082461487742983.0;
230  params->B_yyt[1*s+1] = -32797097931948613195.0 / 61521565616362163366.0;
231  params->B_yyt[1*s+2] = 42895514606418420631.0 / 71714201188501437336.0;
232  params->B_yyt[1*s+3] = 22608567633166065068.0 / 55371917805607957003.0;
233  params->B_yyt[1*s+4] = 94655809487476459565.0 / 151517167160302729021.0;
234 
235  params->C_yyt[0*r+0] = 70820309139834661559.0 / 80863923579509469826.0;
236  params->C_yyt[0*r+1] = 10043614439674808267.0 / 80863923579509469826.0;
237  params->C_yyt[1*r+0] = 161694774978034105510.0 / 106187653640211060371.0;
238  params->C_yyt[1*r+1] = -55507121337823045139.0 / 106187653640211060371.0;
239  params->C_yyt[2*r+0] = 78486094644566264568.0 / 88171030896733822981.0;
240  params->C_yyt[2*r+1] = 9684936252167558413.0 / 88171030896733822981.0;
241  params->C_yyt[3*r+0] = 65394922146334854435.0 / 84570853840405479554.0;
242  params->C_yyt[3*r+1] = 19175931694070625119.0 / 84570853840405479554.0;
243  params->C_yyt[4*r+0] = 8607282770183754108.0 / 108658046436496925911.0;
244  params->C_yyt[4*r+1] = 100050763666313171803.0 / 108658046436496925911.0;
245 
246  params->D_yyt[0*r+0] = 1.0;
247  params->D_yyt[1*r+1] = 1.0;
248 
249  double T[r*r],Tinv[r*r];
250  T[0*r+0] = 1.0;
251  T[0*r+1] = 0.0;
252  T[1*r+0] = 1.0;
253  T[1*r+1] = 1.0-params->gamma;
254  _MatrixInvert_(T,Tinv,r);
255 
256  _ArrayCopy1D_(params->A_yyt,params->A_yeps,(s*s));
257  _MatrixMultiplyNonSquare_(Tinv,params->B_yyt,params->B_yeps,r,r,s);
258  _MatrixMultiplyNonSquare_(params->C_yyt,T,params->C_yeps,s,r,r);
259  _ArrayCopy1D_(params->D_yyt,params->D_yeps,(r*r));
260 
261  } else if (!strcmp(type,_GLM_GEE_EXRK2A_)) {
262 
263  params->A_yeps[1*s+0] = 1.0;
264  params->A_yeps[3*s+2] = 0.5;
265  params->A_yeps[4*s+2] = 0.25;
266  params->A_yeps[4*s+3] = 0.25;
267  params->A_yeps[5*s+2] = 0.25;
268  params->A_yeps[5*s+3] = 0.25;
269  params->A_yeps[5*s+4] = 0.5;
270 
271  params->B_yeps[0*s+0] = params->B_yeps[0*s+1] = 0.5;
272  params->B_yeps[1*s+0] = params->B_yeps[1*s+1] = -2.0/3.0;
273  params->B_yeps[1*s+2] = params->B_yeps[1*s+3] = params->B_yeps[1*s+4]
274  = params->B_yeps[1*s+5] = 1.0/3.0;
275 
276  params->C_yeps[0*r+0] = 1.0;
277  params->C_yeps[1*r+0] = 1.0;
278  params->C_yeps[2*r+0] = 1.0; params->C_yeps[2*r+1] = 0.75;
279  params->C_yeps[3*r+0] = 1.0; params->C_yeps[3*r+1] = 0.75;
280  params->C_yeps[4*r+0] = 1.0; params->C_yeps[4*r+1] = 0.75;
281  params->C_yeps[5*r+0] = 1.0; params->C_yeps[5*r+1] = 0.75;
282 
283  params->D_yeps[0*r+0] = 1.0;
284  params->D_yeps[1*r+1] = 1.0;
285 
286  params->A_yyt[1*s+0] = 1.0;
287  params->A_yyt[3*s+2] = 0.5;
288  params->A_yyt[4*s+2] = 0.25;
289  params->A_yyt[4*s+3] = 0.25;
290  params->A_yyt[5*s+2] = 0.25;
291  params->A_yyt[5*s+3] = 0.25;
292  params->A_yyt[5*s+4] = 0.5;
293 
294  params->B_yyt[0*s+0] = params->B_yyt[0*s+1] = 0.5;
295  params->B_yyt[1*s+2] = params->B_yyt[1*s+3] = params->B_yyt[1*s+4]
296  = params->B_yyt[1*s+5] = 0.25;
297 
298  params->C_yyt[0*r+0] = 1.0;
299  params->C_yyt[1*r+0] = 1.0;
300  params->C_yyt[2*r+1] = 1.0;
301  params->C_yyt[3*r+1] = 1.0;
302  params->C_yyt[4*r+1] = 1.0;
303  params->C_yyt[5*r+1] = 1.0;
304 
305  params->D_yyt[0*r+0] = 1.0;
306  params->D_yyt[1*r+1] = 1.0;
307 
308  } else if (!strcmp(type,_GLM_GEE_RK32G1_)) {
309 
310  params->A_yeps[1*s+0] = 0.5;
311 
312  params->A_yeps[2*s+0] = -1.0;
313  params->A_yeps[2*s+1] = 2.0;
314 
315  params->A_yeps[3*s+0] = 1.0/6.0;
316  params->A_yeps[3*s+1] = 2.0/3.0;
317  params->A_yeps[3*s+2] = 1.0/6.0;
318 
319  params->A_yeps[5*s+0] = -7.0/24.0;
320  params->A_yeps[5*s+1] = 1.0/3.0;
321  params->A_yeps[5*s+2] = 1.0/12.0;
322  params->A_yeps[5*s+3] = -1.0/8.0;
323  params->A_yeps[5*s+4] = 0.5;
324 
325  params->A_yeps[6*s+0] = 7.0/6.0;
326  params->A_yeps[6*s+1] = -4.0/3.0;
327  params->A_yeps[6*s+2] = -1.0/3.0;
328  params->A_yeps[6*s+3] = 0.5;
329  params->A_yeps[6*s+4] = -1.0;
330  params->A_yeps[6*s+5] = 2.0;
331 
332  params->A_yeps[7*s+4] = 1.0/6.0;
333  params->A_yeps[7*s+5] = 2.0/3.0;
334  params->A_yeps[7*s+6] = 1.0/6.0;
335 
336  params->B_yeps[0*s+0] = 1.0/6.0;
337  params->B_yeps[0*s+1] = 2.0/3.0;
338  params->B_yeps[0*s+2] = 1.0/6.0;
339 
340  params->B_yeps[1*s+0] = -1.0/6.0;
341  params->B_yeps[1*s+1] = -2.0/3.0;
342  params->B_yeps[1*s+2] = -1.0/6.0;
343  params->B_yeps[1*s+4] = 1.0/6.0;
344  params->B_yeps[1*s+5] = 2.0/3.0;
345  params->B_yeps[1*s+6] = 1.0/6.0;
346 
347  params->C_yeps[0*r+0] = 1.0;
348  params->C_yeps[1*r+0] = 1.0;
349  params->C_yeps[2*r+0] = 1.0;
350  params->C_yeps[3*r+0] = 1.0;
351  params->C_yeps[4*r+0] = 1.0;
352  params->C_yeps[5*r+0] = 1.0;
353  params->C_yeps[6*r+0] = 1.0;
354  params->C_yeps[7*r+0] = 1.0;
355 
356  params->C_yeps[4*r+1] = 1.0;
357  params->C_yeps[5*r+1] = 1.0;
358  params->C_yeps[6*r+1] = 1.0;
359  params->C_yeps[7*r+1] = 1.0;
360 
361  params->D_yeps[0*r+0] = 1.0;
362  params->D_yeps[1*r+1] = 1.0;
363 
364  double T[r*r],Tinv[r*r];
365  T[0*r+0] = 1.0;
366  T[0*r+1] = 0.0;
367  T[1*r+0] = 1.0;
368  T[1*r+1] = 1.0-params->gamma;
369  _MatrixInvert_(T,Tinv,r);
370 
371  _ArrayCopy1D_(params->A_yeps,params->A_yyt,(s*s));
372  _MatrixMultiplyNonSquare_(T,params->B_yeps,params->B_yyt,r,r,s);
373  _MatrixMultiplyNonSquare_(params->C_yeps,Tinv,params->C_yyt,s,r,r);
374  _ArrayCopy1D_(params->D_yeps,params->D_yyt,(r*r));
375 
376  } else if (!strcmp(type,_GLM_GEE_RK285EX_)) {
377 
378  params->A_yeps[1*s+0] = 0.585786437626904966;
379  params->A_yeps[2*s+0] = 0.149999999999999994;
380  params->A_yeps[2*s+1] = 0.849999999999999978;
381  params->A_yeps[4*s+3] = 0.292893218813452483;
382  params->A_yeps[5*s+3] = 0.0749999999999999972;
383  params->A_yeps[5*s+4] = 0.424999999999999989;
384  params->A_yeps[6*s+3] = 0.176776695296636893;
385  params->A_yeps[6*s+4] = 0.176776695296636893;
386  params->A_yeps[6*s+5] = 0.146446609406726241;
387  params->A_yeps[7*s+3] = 0.176776695296636893;
388  params->A_yeps[7*s+4] = 0.176776695296636893;
389  params->A_yeps[7*s+5] = 0.146446609406726241;
390  params->A_yeps[7*s+6] = 0.292893218813452483;
391  params->A_yeps[8*s+3] = 0.176776695296636893;
392  params->A_yeps[8*s+4] = 0.176776695296636893;
393  params->A_yeps[8*s+5] = 0.146446609406726241;
394  params->A_yeps[8*s+6] = 0.0749999999999999972;
395  params->A_yeps[8*s+7] = 0.424999999999999989;
396 
397  params->B_yeps[0*s+0] = 0.353553390593273786;
398  params->B_yeps[0*s+1] = 0.353553390593273786;
399  params->B_yeps[0*s+2] = 0.292893218813452483;
400 
401  params->B_yeps[1*s+0] =-0.471404520791031678;
402  params->B_yeps[1*s+1] =-0.471404520791031678;
403  params->B_yeps[1*s+2] =-0.390524291751269959;
404  params->B_yeps[1*s+3] = 0.235702260395515839;
405  params->B_yeps[1*s+4] = 0.235702260395515839;
406  params->B_yeps[1*s+5] = 0.195262145875634979;
407  params->B_yeps[1*s+6] = 0.235702260395515839;
408  params->B_yeps[1*s+7] = 0.235702260395515839;
409  params->B_yeps[1*s+8] = 0.195262145875634979;
410 
411  params->C_yeps[0*r+0] = 1.0;
412  params->C_yeps[1*r+0] = 1.0;
413  params->C_yeps[2*r+0] = 1.0;
414  params->C_yeps[3*r+0] = 1.0; params->C_yeps[3*r+1] = 0.75;
415  params->C_yeps[4*r+0] = 1.0; params->C_yeps[4*r+1] = 0.75;
416  params->C_yeps[5*r+0] = 1.0; params->C_yeps[5*r+1] = 0.75;
417  params->C_yeps[6*r+0] = 1.0; params->C_yeps[6*r+1] = 0.75;
418  params->C_yeps[7*r+0] = 1.0; params->C_yeps[7*r+1] = 0.75;
419  params->C_yeps[8*r+0] = 1.0; params->C_yeps[8*r+1] = 0.75;
420 
421  params->D_yeps[0*r+0] = 1.0;
422  params->D_yeps[1*r+1] = 1.0;
423 
424  double T[r*r],Tinv[r*r];
425  T[0*r+0] = 1.0;
426  T[0*r+1] = 0.0;
427  T[1*r+0] = 1.0;
428  T[1*r+1] = 1.0-params->gamma;
429  _MatrixInvert_(T,Tinv,r);
430 
431  _ArrayCopy1D_(params->A_yeps,params->A_yyt,(s*s));
432  _MatrixMultiplyNonSquare_(T,params->B_yeps,params->B_yyt,r,r,s);
433  _MatrixMultiplyNonSquare_(params->C_yeps,Tinv,params->C_yyt,s,r,r);
434  _ArrayCopy1D_(params->D_yeps,params->D_yyt,(r*r));
435 
436  }
437 
438  for (i=0; i<s; i++) {
439  for (j=0; j<s; j++) {
440  params->c_yyt[i] += params->A_yyt [i*s+j];
441  params->c_yeps[i] += params->A_yeps[i*s+j];
442  }
443  }
444 
445  if (!mpi->rank) {
446  FILE *in;
447  int ferr;
448  in = fopen("glm_gee.inp","r");
449  strcat(params->ee_mode,_GLM_GEE_YEPS_);
450  if (in) {
451  printf("Reading GLM-GEE method parameters from glm_gee.inp.\n");
452  char word[_MAX_STRING_SIZE_];
453  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
454  if (!strcmp(word,"begin")) {
455  while (strcmp(word,"end")) {
456  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
457  if (!strcmp(word,"ee_mode")) {
458  ferr = fscanf(in,"%s",params->ee_mode);
459  if (ferr != 1) return(1);
460  } else if (strcmp(word,"end")) {
461  char useless[_MAX_STRING_SIZE_];
462  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
463  printf("Warning: keyword %s in file \"glm_gee.inp\" with value %s not ",word,useless);
464  printf("recognized or extraneous. Ignoring.\n");
465  }
466  }
467  } else {
468  fprintf(stderr,"Error: Illegal format in file \"glm_gee.inp\".\n");
469  return(1);
470  }
471  fclose(in);
472  if (strcmp(params->ee_mode,_GLM_GEE_YEPS_) && strcmp(params->ee_mode,_GLM_GEE_YYT_)) {
473  fprintf(stderr,"Error in TimeGLMGEEInitialize(): %s is not a valid value for ",params->ee_mode);
474  fprintf(stderr,"ee_mode. Acceptable inputs are %s or %s.\n",_GLM_GEE_YEPS_,_GLM_GEE_YYT_);
475  strcat(params->ee_mode,_GLM_GEE_YEPS_);
476  }
477  }
478  printf("GLM-GEE time integration error estimation mode: %s\n",params->ee_mode);
479  }
480  IERR MPIBroadcast_character(params->ee_mode,_MAX_STRING_SIZE_,0,&mpi->world);
481 
482  if (!strcmp(params->ee_mode,_GLM_GEE_YYT_)) {
483  params->A = params->A_yyt;
484  params->B = params->B_yyt;
485  params->C = params->C_yyt;
486  params->D = params->D_yyt;
487  params->c = params->c_yyt;
488  } else {
489  params->A = params->A_yeps;
490  params->B = params->B_yeps;
491  params->C = params->C_yeps;
492  params->D = params->D_yeps;
493  params->c = params->c_yeps;
494  }
495 
496  } else {
497  fprintf(stderr,"Error in TimeGLMGEEInitialize(): Code should not have ");
498  fprintf(stderr,"reached here for %s class of time-integrators. This is a ",class);
499  fprintf(stderr,"coding mistake.\n");
500  return(1);
501  }
502  return(0);
503 }
_GLM_GEE_RK285EX_
#define _GLM_GEE_RK285EX_
Definition: timeintegration_struct.h:257
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_GLM_GEE_25I_
#define _GLM_GEE_25I_
Definition: timeintegration_struct.h:225
_MatrixMultiplyNonSquare_
#define _MatrixMultiplyNonSquare_(A, B, C, NRowA, NColA, NColB)
Definition: matops.h:40
_GLM_GEE_35_
#define _GLM_GEE_35_
Definition: timeintegration_struct.h:233
_GLM_GEE_23_
#define _GLM_GEE_23_
Definition: timeintegration_struct.h:209
_GLM_GEE_24_
#define _GLM_GEE_24_
Definition: timeintegration_struct.h:217
_GLM_GEE_RK32G1_
#define _GLM_GEE_RK32G1_
Definition: timeintegration_struct.h:249
_GLM_GEE_
#define _GLM_GEE_
Definition: timeintegration_struct.h:18
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
_GLM_GEE_YYT_
#define _GLM_GEE_YYT_
Definition: timeintegration_struct.h:193
GLMGEEParameters::ee_mode
char ee_mode[_MAX_STRING_SIZE_]
Definition: timeintegration_struct.h:270
_GLM_GEE_YEPS_
#define _GLM_GEE_YEPS_
Definition: timeintegration_struct.h:201
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
MPIBroadcast_character
int MPIBroadcast_character(char *, int, int, void *)
Definition: MPIBroadcast.c:37
_MatrixInvert_
#define _MatrixInvert_(A, B, N)
Definition: matops.h:93
IERR
#define IERR
Definition: basic.h:16
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_GLM_GEE_EXRK2A_
#define _GLM_GEE_EXRK2A_
Definition: timeintegration_struct.h:241
int TimeGLMGEECleanup ( void *  s)

Clean up variables relted to General Linear Methods with Global Error Estimators (GLM-GEE)

Clean up allocations for the GLM-GEE (_GLM_GEE_) time integration method: This function frees the arrays used to store the Butcher tableaux

Parameters
sObject of type GLMGEEParameters

Definition at line 14 of file TimeGLMGEECleanup.c.

15 {
16  GLMGEEParameters *params = (GLMGEEParameters*) s;
17  if (params->A_yyt ) free(params->A_yyt);
18  if (params->B_yyt ) free(params->B_yyt);
19  if (params->C_yyt ) free(params->C_yyt);
20  if (params->D_yyt ) free(params->D_yyt);
21  if (params->c_yyt ) free(params->c_yyt);
22  if (params->A_yeps) free(params->A_yeps);
23  if (params->B_yeps) free(params->B_yeps);
24  if (params->C_yeps) free(params->C_yeps);
25  if (params->D_yeps) free(params->D_yeps);
26  if (params->c_yeps) free(params->c_yeps);
27  return(0);
28 }
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
int TimeInitialize ( void *  s,
int  nsims,
int  rank,
int  nproc,
void *  ts 
)

Initialize the time integration

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;
36  SimulationObject* sim = (SimulationObject*) s;
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 */
175  TS->RHSFunction = TimeRHSFunctionExplicit;
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 }
_FORWARD_EULER_
#define _FORWARD_EULER_
Definition: timeintegration_struct.h:14
TimeIntegration::RHSFunction
int(* RHSFunction)(double *, double *, void *, void *, double)
Definition: timeintegration_struct.h:98
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::dt
double dt
Definition: hypar.h:67
_GLM_GEE_
#define _GLM_GEE_
Definition: timeintegration_struct.h:18
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
ExplicitRKParameters
Structure containing the parameters for an explicit Runge-Kutta method.
Definition: timeintegration_struct.h:177
_GLM_GEE_YYT_
#define _GLM_GEE_YYT_
Definition: timeintegration_struct.h:193
_RK_
#define _RK_
Definition: timeintegration_struct.h:16
GLMGEEParameters::ee_mode
char ee_mode[_MAX_STRING_SIZE_]
Definition: timeintegration_struct.h:270
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
gpuMemset
void gpuMemset(void *, int, size_t)
Definition: ArrayImplementations_GPU.cu:131
gpuMalloc
void gpuMalloc(void **, size_t)
Definition: ArrayImplementations_GPU.cu:118
HyPar::n_iter
int n_iter
Definition: hypar.h:55
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::ndims
int ndims
Definition: hypar.h:26
TimeIntegration::TimeIntegrate
int(* TimeIntegrate)(void *)
Definition: timeintegration_struct.h:96
HyPar::TimeIntegrate
int(* TimeIntegrate)(void *)
Definition: hypar.h:220
HyPar::restart_iter
int restart_iter
Definition: hypar.h:58
TimeRHSFunctionExplicit
int TimeRHSFunctionExplicit(double *, double *, void *, void *, double)
Definition: TimeRHSFunctionExplicit.c:30
HyPar::time_integrator
void * time_integrator
Definition: hypar.h:165
int TimeCleanup ( void *  ts)

Clean up variables related to time integration

Clean up all allocations related to time integration

Parameters
tsObject of type TimeIntegration

Definition at line 18 of file TimeCleanup.c.

19 {
20  TimeIntegration* TS = (TimeIntegration*) ts;
21  SimulationObject* sim = (SimulationObject*) TS->simulation;
22  int ns, nsims = TS->nsims;
23 
24  /* close files opened for writing */
25  if (!TS->rank) if (sim[0].solver.write_residual) fclose((FILE*)TS->ResidualFile);
26 
27 #if defined(HAVE_CUDA)
28  if (sim[0].solver.use_gpu) {
29  if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
30  gpuFree(TS->gpu_U);
31  gpuFree(TS->gpu_Udot);
32  gpuFree(TS->gpu_BoundaryFlux);
33  }
34  } else {
35 #endif
36  if (!strcmp(sim[0].solver.time_scheme,_RK_)) {
37  int i;
38  ExplicitRKParameters *params = (ExplicitRKParameters*) sim[0].solver.msti;
39  for (i=0; i<params->nstages; i++) free(TS->U[i]); free(TS->U);
40  for (i=0; i<params->nstages; i++) free(TS->Udot[i]); free(TS->Udot);
41  for (i=0; i<params->nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
42  } else if (!strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
43  int nstages = 1, i;
44  for (i=0; i<nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
45  } else if (!strcmp(sim[0].solver.time_scheme,_GLM_GEE_)) {
46  int i;
47  GLMGEEParameters *params = (GLMGEEParameters*) sim[0].solver.msti;
48  for (i=0; i<2*params->r-1 ; i++) free(TS->U[i]); free(TS->U);
49  for (i=0; i<params->nstages; i++) free(TS->Udot[i]); free(TS->Udot);
50  for (i=0; i<params->nstages; i++) free(TS->BoundaryFlux[i]); free(TS->BoundaryFlux);
51  }
52 #if defined(HAVE_CUDA)
53  }
54 #endif
55 
56  /* deallocate arrays */
57  free(TS->u_offsets);
58  free(TS->u_sizes);
59  free(TS->bf_offsets);
60  free(TS->bf_sizes);
61  free(TS->u );
62  free(TS->rhs);
63  for (ns = 0; ns < nsims; ns++) {
64  sim[ns].solver.time_integrator = NULL;
65  }
66  return(0);
67 }
_FORWARD_EULER_
#define _FORWARD_EULER_
Definition: timeintegration_struct.h:14
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
_GLM_GEE_
#define _GLM_GEE_
Definition: timeintegration_struct.h:18
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
ExplicitRKParameters
Structure containing the parameters for an explicit Runge-Kutta method.
Definition: timeintegration_struct.h:177
_RK_
#define _RK_
Definition: timeintegration_struct.h:16
gpuFree
void gpuFree(void *)
Definition: ArrayImplementations_GPU.cu:140
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::time_integrator
void * time_integrator
Definition: hypar.h:165
int TimePreStep ( void *  ts)

Function called at the beginning of a time step

Pre-time-step function: This function is called before each time step. Some notable things this does are:

  • Computes CFL and diffusion numbers.
  • Call the physics-specific pre-time-step function, if defined.
Parameters
tsObject of type TimeIntegration

Definition at line 22 of file TimePreStep.c.

23 {
24  TimeIntegration* TS = (TimeIntegration*) ts;
25  _DECLARE_IERR_;
26 
27  SimulationObject* sim = (SimulationObject*) TS->simulation;
28  int ns, nsims = TS->nsims;
29 
30  gettimeofday(&TS->iter_start_time,NULL);
31 
32  for (ns = 0; ns < nsims; ns++) {
33 
34  HyPar* solver = &(sim[ns].solver);
35  MPIVariables* mpi = &(sim[ns].mpi);
36 
37  double *u = NULL;
38 #if defined(HAVE_CUDA)
39  if (solver->use_gpu) {
40  u = solver->gpu_u;
41  } else{
42 #endif
43  u = solver->u;
44 #if defined(HAVE_CUDA)
45  }
46 #endif
47 
48  /* apply boundary conditions and exchange data over MPI interfaces */
49 
50  solver->ApplyBoundaryConditions( solver,
51  mpi,
52  u,
53  NULL,
54  TS->waqt );
55 
56  solver->ApplyIBConditions( solver,
57  mpi,
58  u,
59  TS->waqt );
60 
61 #if defined(HAVE_CUDA)
62  if (solver->use_gpu) {
63  gpuMPIExchangeBoundariesnD( solver->ndims,
64  solver->nvars,
65  solver->gpu_dim_local,
66  solver->ghosts,
67  mpi,
68  u );
69  } else {
70 #endif
71  MPIExchangeBoundariesnD( solver->ndims,
72  solver->nvars,
73  solver->dim_local,
74  solver->ghosts,
75  mpi,
76  u );
77 #if defined(HAVE_CUDA)
78  }
79 #endif
80 
81  if ((TS->iter+1)%solver->screen_op_iter == 0) {
82 
83 #if defined(HAVE_CUDA)
84  if (!solver->use_gpu) {
85 #endif
86  _ArrayCopy1D_( solver->u,
87  (TS->u + TS->u_offsets[ns]),
88  (solver->npoints_local_wghosts*solver->nvars) );
89 #if defined(HAVE_CUDA)
90  }
91 #endif
92 
93  /* compute max CFL and diffusion number over the domain */
94  if (solver->ComputeCFL) {
95  double local_max_cfl = -1.0;
96  local_max_cfl = solver->ComputeCFL (solver,mpi,TS->dt,TS->waqt);
97  MPIMax_double(&TS->max_cfl ,&local_max_cfl ,1,&mpi->world);
98  } else {
99  TS->max_cfl = -1;
100  }
101  if (solver->ComputeDiffNumber) {
102  double local_max_diff = -1.0;
103  local_max_diff = solver->ComputeDiffNumber (solver,mpi,TS->dt,TS->waqt);
104  MPIMax_double(&TS->max_diff,&local_max_diff,1,&mpi->world);
105  } else {
106  TS->max_diff = -1;
107  }
108 
109  }
110 
111  /* set the step boundary flux integral value to zero */
112 #if defined(HAVE_CUDA)
113  if (solver->use_gpu) {
114  gpuArraySetValue(solver->StepBoundaryIntegral,2*solver->ndims*solver->nvars,0.0);
115  } else {
116 #endif
117  _ArraySetValue_(solver->StepBoundaryIntegral,2*solver->ndims*solver->nvars,0.0);
118 #if defined(HAVE_CUDA)
119  }
120 #endif
121 
122  if (solver->PreStep) {
123  solver->PreStep(u,solver,mpi,TS->waqt);
124  }
125 
126  }
127 
128  return 0;
129 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::ComputeCFL
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
HyPar::u
double * u
Definition: hypar.h:116
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
TimeIntegration::iter_start_time
struct timeval iter_start_time
Definition: timeintegration_struct.h:101
gpuArraySetValue
void gpuArraySetValue(double *, int, double)
Definition: ArrayImplementations_GPU.cu:161
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::screen_op_iter
int screen_op_iter
Definition: hypar.h:168
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
HyPar::StepBoundaryIntegral
double * StepBoundaryIntegral
Definition: hypar.h:384
HyPar::gpu_u
double * gpu_u
Definition: hypar.h:459
HyPar::ComputeDiffNumber
double(* ComputeDiffNumber)(void *, void *, double, double)
Definition: hypar.h:272
HyPar::gpu_dim_local
int * gpu_dim_local
Definition: hypar.h:455
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
gpuMPIExchangeBoundariesnD
int gpuMPIExchangeBoundariesnD(int, int, const int *, int, void *, double *)
HyPar::use_gpu
int use_gpu
Definition: hypar.h:449
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::ApplyBoundaryConditions
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
MPIExchangeBoundariesnD
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Definition: MPIExchangeBoundariesnD.c:42
HyPar::PreStep
int(* PreStep)(double *, void *, void *, double)
Definition: hypar.h:339
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::ApplyIBConditions
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
int TimeStep ( void *  ts)

Take one step in time

Advance one time step.

Parameters
tsObject of type TimeIntegration

Definition at line 13 of file TimeStep.c.

14 {
15  TimeIntegration *TS = (TimeIntegration*) ts;
16  SimulationObject* sim = (SimulationObject*) TS->simulation;
17 
18  if (TS->TimeIntegrate) { TS->TimeIntegrate(TS); }
19 
20  return(0);
21 }
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
TimeIntegration::TimeIntegrate
int(* TimeIntegrate)(void *)
Definition: timeintegration_struct.h:96
int TimePostStep ( void *  ts)

Function called at the end of a time step

Post-time-step function: this function is called at the end of each time step.

  • It updates the current simulation time.
  • It calls functions to print information and to write transient solution to file.
  • It will also call any physics-specific post-time-step function, if defined.
Parameters
tsObject of type TimeIntegration

Definition at line 28 of file TimePostStep.c.

29 {
30  TimeIntegration* TS = (TimeIntegration*) ts;
31  SimulationObject* sim = (SimulationObject*) TS->simulation;
32  int ns, nsims = TS->nsims;
33 
34  /* update current time */
35  TS->waqt += TS->dt;
36 
37  if ((TS->iter+1)%sim[0].solver.screen_op_iter == 0) {
38 
39 #if defined(HAVE_CUDA)
40  if (sim[0].solver.use_gpu) {
41  TS->norm = -1;
42  } else {
43 #endif
44  /* Calculate norm for this time step */
45  double sum = 0.0;
46  double npts = 0;
47  for (ns = 0; ns < nsims; ns++) {
48  _ArrayAXPY_(sim[ns].solver.u,-1.0,(TS->u+TS->u_offsets[ns]),TS->u_sizes[ns]);
49  sum += ArraySumSquarenD( sim[ns].solver.nvars,
50  sim[ns].solver.ndims,
51  sim[ns].solver.dim_local,
52  sim[ns].solver.ghosts,
53  sim[ns].solver.index,
54  (TS->u+TS->u_offsets[ns]) );
55  npts += (double)sim[ns].solver.npoints_global;
56  }
57 
58  double global_sum = 0;
59  MPISum_double( &global_sum,
60  &sum,1,
61  &(sim[0].mpi.world) );
62 
63  TS->norm = sqrt(global_sum/npts);
64 #if defined(HAVE_CUDA)
65  }
66 #endif
67 
68  /* write to file */
69  if (TS->ResidualFile) {
70  fprintf((FILE*)TS->ResidualFile,"%10d\t%E\t%E\n",TS->iter+1,TS->waqt,TS->norm);
71  }
72 
73  }
74 
75 
76 #if defined(HAVE_CUDA)
77  if (!sim[0].solver.use_gpu) {
78 #endif
79  for (ns = 0; ns < nsims; ns++) {
80 
81  if (!strcmp(sim[ns].solver.ConservationCheck,"yes")) {
82  /* calculate volume integral of the solution at this time step */
83  IERR sim[ns].solver.VolumeIntegralFunction( sim[ns].solver.VolumeIntegral,
84  sim[ns].solver.u,
85  &(sim[ns].solver),
86  &(sim[ns].mpi) ); CHECKERR(ierr);
87  /* calculate surface integral of the flux at this time step */
88  IERR sim[ns].solver.BoundaryIntegralFunction( &(sim[ns].solver),
89  &(sim[ns].mpi)); CHECKERR(ierr);
90  /* calculate the conservation error at this time step */
91  IERR sim[ns].solver.CalculateConservationError( &(sim[ns].solver),
92  &(sim[ns].mpi)); CHECKERR(ierr);
93  }
94 
95  if (sim[ns].solver.PostStep) {
96  sim[ns].solver.PostStep(sim[ns].solver.u,&(sim[ns].solver),&(sim[ns].mpi),TS->waqt,TS->iter);
97  }
98 
99  }
100 #if defined(HAVE_CUDA)
101  }
102 #endif
103 
104  gettimeofday(&TS->iter_end_time,NULL);
105  long long walltime;
106  walltime = ( (TS->iter_end_time.tv_sec * 1000000 + TS->iter_end_time.tv_usec)
107  - (TS->iter_start_time.tv_sec * 1000000 + TS->iter_start_time.tv_usec));
108  TS->iter_wctime = (double) walltime / 1000000.0;
109  TS->iter_wctime_total += TS->iter_wctime;
110 
111  double global_total = 0, global_wctime = 0, global_mpi_total = 0, global_mpi_wctime = 0;
112 
113  MPIMax_double(&global_wctime, &TS->iter_wctime, 1, &(sim[0].mpi.world));
114  MPIMax_double(&global_total, &TS->iter_wctime_total, 1, &(sim[0].mpi.world));
115 
116 #if defined(HAVE_CUDA)
117  MPIMax_double(&global_mpi_wctime, &sim[0].mpi.wctime, 1, &(sim[0].mpi.world));
118  MPIMax_double(&global_mpi_total, &sim[0].mpi.wctime_total, 1, &(sim[0].mpi.world));
119 #endif
120 
121  return(0);
122 }
HyPar::u
double * u
Definition: hypar.h:116
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
TimeIntegration::iter_start_time
struct timeval iter_start_time
Definition: timeintegration_struct.h:101
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::screen_op_iter
int screen_op_iter
Definition: hypar.h:168
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
HyPar::BoundaryIntegralFunction
int(* BoundaryIntegralFunction)(void *, void *)
Definition: hypar.h:390
TimeIntegration::iter_end_time
struct timeval iter_end_time
Definition: timeintegration_struct.h:103
HyPar::VolumeIntegralFunction
int(* VolumeIntegralFunction)(double *, double *, void *, void *)
Definition: hypar.h:388
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::index
int * index
Definition: hypar.h:102
IERR
#define IERR
Definition: basic.h:16
HyPar::CalculateConservationError
int(* CalculateConservationError)(void *, void *)
Definition: hypar.h:392
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar::PostStep
int(* PostStep)(double *, void *, void *, double, int)
Definition: hypar.h:341
int TimePrintStep ( void *  ts)

Print time integration related information

Print information to screen (also calls any physics-specific printing function, if defined).

Parameters
tsObject of type TimeIntegration

Definition at line 16 of file TimePrintStep.c.

17 {
18  TimeIntegration* TS = (TimeIntegration*) ts;
19  SimulationObject* sim = (SimulationObject*) TS->simulation;
20  int ns, nsims = TS->nsims;
21 
22  if ((!TS->rank) && ((TS->iter+1)%sim[0].solver.screen_op_iter == 0)) {
23  if (nsims > 1) {
24  printf("--\n");
25  printf("iter=%7d, t=%1.3e\n", TS->iter+1, TS->waqt);
26  if (TS->max_cfl >= 0) printf(" CFL=%1.3E\n", TS->max_cfl);
27  if (TS->norm >= 0) printf(" norm=%1.4E\n", TS->norm);
28  printf(" wctime=%1.1E (s)\n",TS->iter_wctime);
29  } else {
30  printf("iter=%7d ",TS->iter+1 );
31  printf("t=%1.3E ",TS->waqt );
32  if (TS->max_cfl >= 0) printf("CFL=%1.3E ",TS->max_cfl );
33  if (TS->norm >= 0) printf("norm=%1.4E ",TS->norm );
34  printf("wctime: %1.1E (s) ",TS->iter_wctime);
35  }
36 
37  /* calculate and print conservation error */
38  if (!strcmp(sim[0].solver.ConservationCheck,"yes")) {
39 
40  double error = 0;
41  for (ns = 0; ns < nsims; ns++) {
42  int v;
43  for (v=0; v<sim[ns].solver.nvars; v++) {
44  error += (sim[ns].solver.ConservationError[v] * sim[ns].solver.ConservationError[v]);
45  }
46  }
47  error = sqrt(error);
48  printf(" cons_err=%1.4E\n", error);
49 
50  } else {
51 
52  if (nsims == 1) printf("\n");
53 
54  }
55 
56  /* print physics-specific info, if available */
57  for (ns = 0; ns < nsims; ns++) {
58  if (sim[ns].solver.PrintStep) {
59  if (nsims > 1) printf("Physics-specific output for domain %d:\n", ns);
60  sim[ns].solver.PrintStep( &(sim[ns].solver),
61  &(sim[ns].mpi),
62  TS->waqt );
63  }
64  }
65  if (nsims > 1) {
66  printf("--\n");
67  printf("\n");
68  }
69  }
70 
71  return(0);
72 }
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::nvars
int nvars
Definition: hypar.h:29
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::PrintStep
int(* PrintStep)(void *, void *, double)
Definition: hypar.h:344
HyPar::ConservationError
double * ConservationError
Definition: hypar.h:374
int TimeError ( void *  s,
void *  m,
double *  uex 
)

Compute/estimate error in solution

Computes the time integration error in the solution: If the time integration method chosen has a mechanism to compute the error in the numerical integration, it is computed in this function. In addition, if an exact solution is available (see function argument uex, the error of the computed solution (stored in HyPar::u) with respect to this exact solution is also computed. These errors are written to a text file whose name depends on the time integration method being used.

Time integration method with error estimation currently implemented:

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
uexExact solution (stored with the same array layout as HyPar::u (may be NULL)

Definition at line 27 of file TimeError.c.

34 {
35  HyPar *solver = (HyPar*) s;
36  MPIVariables *mpi = (MPIVariables*) m;
37  TimeIntegration *TS = (TimeIntegration*) solver->time_integrator;
38  int size = solver->npoints_local_wghosts * solver->nvars;
39  double sum = 0.0, global_sum = 0.0;
40  static const double tolerance = 1e-15;
41  if (!TS) return(0);
42 
43  if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
44  /* For GLM-GEE methods, calculate the norm of the estimated global error */
45  GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
46  double error[6] = {0,0,0,0,0,0}, *Uerr;
47  if (!strcmp(params->ee_mode,_GLM_GEE_YEPS_)) Uerr = TS->U[params->r];
48  else {
49  Uerr = TS->U[0];
50  _ArraySubtract1D_(Uerr,solver->u,TS->U[params->r],size);
51  _ArrayScale1D_(Uerr,(1.0/(1.0-params->gamma)),size);
52  }
53 
54  /* calculate solution norm for relative errors */
55  double sol_norm[3] = {0.0,0.0,0.0};
56  /* L1 */
57  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
58  solver->ghosts,solver->index,solver->u);
59  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
60  sol_norm[0] = global_sum/((double)solver->npoints_global);
61  /* L2 */
62  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
63  solver->ghosts,solver->index,solver->u);
64  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
65  sol_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
66  /* Linf */
67  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
68  solver->ghosts,solver->index,solver->u);
69  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
70  sol_norm[2] = global_sum;
71 
72  /* calculate L1 norm of error */
73  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
74  solver->ghosts,solver->index,Uerr);
75  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
76  error[0] = global_sum/((double)solver->npoints_global);
77  /* calculate L2 norm of error */
78  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
79  solver->ghosts,solver->index,Uerr);
80  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
81  error[1] = sqrt(global_sum/((double)solver->npoints_global));
82  /* calculate Linf norm of error */
83  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
84  solver->ghosts,solver->index,Uerr);
85  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
86  error[2] = global_sum;
87 
88  if (uex) {
89  _ArrayAXBY_(TS->Udot[0],1.0,solver->u,-1.0,uex,size);
90  _ArrayAXPY_(Uerr,-1.0,TS->Udot[0],size);
91  /* calculate L1 norm of error */
92  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
93  solver->ghosts,solver->index,TS->Udot[0]);
94  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
95  error[3] = global_sum/((double)solver->npoints_global);
96  /* calculate L2 norm of error */
97  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
98  solver->ghosts,solver->index,TS->Udot[0]);
99  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
100  error[4] = sqrt(global_sum/((double)solver->npoints_global));
101  /* calculate Linf norm of error */
102  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
103  solver->ghosts,solver->index,TS->Udot[0]);
104  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
105  error[5] = global_sum;
106  } else error[3] = error[4] = error[5] = -1;
107 
108  if ( (sol_norm[0] > tolerance)
109  && (sol_norm[1] > tolerance)
110  && (sol_norm[2] > tolerance) ) {
111  error[0] /= sol_norm[0];
112  error[1] /= sol_norm[1];
113  error[2] /= sol_norm[2];
114  if (uex) {
115  error[3] /= sol_norm[0];
116  error[4] /= sol_norm[1];
117  error[5] /= sol_norm[2];
118  }
119  }
120 
121  /* write to file */
122  if (!mpi->rank) {
123  FILE *out;
124  out = fopen("glm_err.dat","w");
125  fprintf(out,"%1.16E %1.16E %1.16E %1.16E ",TS->dt,error[0],error[1],error[2]);
126  fprintf(out,"%1.16E %1.16E %1.16E\n",error[3],error[4],error[5]);
127  fclose(out);
128  printf("Estimated time integration errors (GLM-GEE time-integration):-\n");
129  printf(" L1 Error : %1.16E\n",error[0]);
130  printf(" L2 Error : %1.16E\n",error[1]);
131  printf(" Linfinity Error : %1.16E\n",error[2]);
132  }
133  }
134 
135  return(0);
136 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::u
double * u
Definition: hypar.h:116
_ArrayAXBY_
#define _ArrayAXBY_(z, a, x, b, y, size)
Definition: arrayfunctions.h:283
HyPar::npoints_global
int npoints_global
Definition: hypar.h:40
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
_ArrayScale1D_
#define _ArrayScale1D_(x, a, size)
Definition: arrayfunctions.h:190
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_GLM_GEE_
#define _GLM_GEE_
Definition: timeintegration_struct.h:18
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
HyPar::time_scheme
char time_scheme[_MAX_STRING_SIZE_]
Definition: hypar.h:78
_GLM_GEE_YEPS_
#define _GLM_GEE_YEPS_
Definition: timeintegration_struct.h:201
HyPar::nvars
int nvars
Definition: hypar.h:29
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
HyPar::msti
void * msti
Definition: hypar.h:366
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::index
int * index
Definition: hypar.h:102
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
HyPar::time_integrator
void * time_integrator
Definition: hypar.h:165
int TimeGetAuxSolutions ( int *  N,
double **  uaux,
void *  s,
int  nu,
int  ns 
)

Function to get auxiliary solutions if available (for example, in GLM-GEE methods)

Return auxiliary solution: Some time integrators may have the concept of auxiliary solutions that they evolve along with the main solution HyPar::u (these may be used for error estimation, for example). This function returns a pointer to such an auxiliary solution. Note that the auxiliary solution has the same dimensions and array layout as the main solution.

  • Call with the final argument n less than 0 to get the total number of auxiliary solutions (this value will be stored in the argument N at exit).
  • Call with the final argument n less than or equal to zero to get the n-th auxiliary solution (the argument uaux will point to this at exit).

Note that auxiliary solutions are numbered in the C convention: 0,1,...,N-1.

Time integration methods which use auxiliary solutions currently implemented:

Parameters
NNumber of auxiliary solutions
uauxPointer to the array holding the auxiliary solution
sSolver object of type HyPar
nuIndex of the auxiliary solution to return
nsIndex of the simulation domain of which the auxiliary solution to return

Definition at line 29 of file TimeGetAuxSolutions.c.

36 {
37  HyPar *solver = (HyPar*) s;
38  TimeIntegration *TS = (TimeIntegration*) solver->time_integrator;
39 
40  if (nu >= 0) {
41  if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
42  GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
43  *uaux = (TS->U[params->r+nu] + TS->u_offsets[ns]);
44  }
45  } else {
46  if (!TS) *N = 0;
47  else {
48  if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
49  GLMGEEParameters *params = (GLMGEEParameters*) solver->msti;
50  *N = params->r - 1;
51  } else *N = 0;
52  }
53  }
54 
55  return(0);
56 }
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
_GLM_GEE_
#define _GLM_GEE_
Definition: timeintegration_struct.h:18
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
HyPar::time_scheme
char time_scheme[_MAX_STRING_SIZE_]
Definition: hypar.h:78
HyPar::msti
void * msti
Definition: hypar.h:366
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
HyPar::time_integrator
void * time_integrator
Definition: hypar.h:165
int TimeForwardEuler ( void *  ts)

Take a step in time using the Forward Euler method

Advance the ODE given by

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

by one time step of size HyPar::dt using the forward Euler method given by

\begin{equation} {\bf u}^{n+1} = {\bf u}^n + \Delta t {\bf F}\left( {\bf u}^n \right) \end{equation}

where the superscript represents the time level, \(\Delta t\) is the time step size HyPar::dt, and \({\bf F}\left({\bf u}\right)\) is computed by TimeIntegration::RHSFunction.

Parameters
tsTime integrator object of type TimeIntegration

Definition at line 25 of file TimeForwardEuler.c.

28 {
29  TimeIntegration* TS = (TimeIntegration*) ts;
30  SimulationObject* sim = (SimulationObject*) TS->simulation;
31  int ns, nsims = TS->nsims;
32  _DECLARE_IERR_;
33 
34  for (ns = 0; ns < nsims; ns++) {
35 
36  HyPar* solver = &(sim[ns].solver);
37  MPIVariables* mpi = &(sim[ns].mpi);
38 
39  double* rhs = TS->rhs + TS->u_offsets[ns];
40 
41  /* Evaluate right-hand side */
42  IERR TS->RHSFunction( rhs,
43  solver->u,
44  solver,
45  mpi,
46  TS->waqt ); CHECKERR(ierr);
47 
48  /* update solution */
49  _ArrayAXPY_( rhs,
50  TS->dt,
51  solver->u,
52  TS->u_sizes[ns] );
53 
54  _ArrayScaleCopy1D_( solver->StageBoundaryIntegral,
55  TS->dt,
56  solver->StepBoundaryIntegral,
57  TS->bf_sizes[ns] );
58  }
59 
60  return(0);
61 }
TimeIntegration::RHSFunction
int(* RHSFunction)(double *, double *, void *, void *, double)
Definition: timeintegration_struct.h:98
HyPar::u
double * u
Definition: hypar.h:116
_ArrayScaleCopy1D_
#define _ArrayScaleCopy1D_(x, a, y, size)
Definition: arrayfunctions.h:356
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::StepBoundaryIntegral
double * StepBoundaryIntegral
Definition: hypar.h:384
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
HyPar::StageBoundaryIntegral
double * StageBoundaryIntegral
Definition: hypar.h:382
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
IERR
#define IERR
Definition: basic.h:16
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
int TimeRK ( void *  ts)

Take a step in time using the explicit Runge-Kutta method

Advance the ODE given by

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

by one time step of size HyPar::dt using the forward Euler method given by

\begin{align} {\bf U}^{\left(i\right)} &= {\bf u}_n + \Delta t \sum_{j=1}^{i-1} a_{ij} {\bf F}\left({\bf U}^{\left(j\right)}\right), \\ {\bf u}_{n+1} &= {\bf u}_n + \Delta t \sum_{i=1}^s b_{i} {\bf F}\left({\bf U}^{\left(i\right)}\right), \end{align}

where the subscript represents the time level, the superscripts represent the stages, \(\Delta t\) is the time step size HyPar::dt, and \({\bf F}\left({\bf u}\right)\) is computed by TimeIntegration::RHSFunction. The Butcher tableaux coefficients are \(a_{ij}\) (ExplicitRKParameters::A) and \(b_i\) (ExplicitRKParameters::b).

Note: In the code TimeIntegration::Udot is equivalent to \({\bf F}\left({\bf u}\right)\).

Parameters
tsObject of type TimeIntegration

Definition at line 35 of file TimeRK.c.

36 {
37  TimeIntegration* TS = (TimeIntegration*) ts;
38  SimulationObject* sim = (SimulationObject*) TS->simulation;
39  ExplicitRKParameters *params = (ExplicitRKParameters*) sim[0].solver.msti;
40  int ns, stage, i, nsims = TS->nsims;
41 
42 #if defined(HAVE_CUDA)
43  if (sim[0].solver.use_gpu) {
44 
45  /* Calculate stage values */
46  for (stage = 0; stage < params->nstages; stage++) {
47 
48  double stagetime = TS->waqt + params->c[stage]*TS->dt;
49 
50  for (ns = 0; ns < nsims; ns++) {
51  gpuArrayCopy1D( sim[ns].solver.gpu_u,
52  (TS->gpu_U + stage*TS->u_size_total + TS->u_offsets[ns]),
53  (TS->u_sizes[ns]) );
54 
55  }
56 
57  for (i = 0; i < stage; i++) {
58  gpuArrayAXPY( TS->gpu_Udot + i*TS->u_size_total,
59  (TS->dt * params->A[stage*params->nstages+i]),
60  TS->gpu_U + stage*TS->u_size_total,
61  TS->u_size_total );
62 
63  }
64 
65  for (ns = 0; ns < nsims; ns++) {
66  if (sim[ns].solver.PreStage) {
67  fprintf(stderr,"ERROR in TimeRK(): Call to solver->PreStage() commented out!\n");
68  return 1;
69 // sim[ns].solver.PreStage( stage,
70 // (TS->gpu_U),
71 // &(sim[ns].solver),
72 // &(sim[ns].mpi),
73 // stagetime ); CHECKERR(ierr);
74  }
75  }
76 
77  for (ns = 0; ns < nsims; ns++) {
78  if (sim[ns].solver.PostStage) {
79  sim[ns].solver.PostStage( (TS->gpu_U + stage*TS->u_size_total + TS->u_offsets[ns]),
80  &(sim[ns].solver),
81  &(sim[ns].mpi),
82  stagetime);
83  }
84  }
85 
86  for (ns = 0; ns < nsims; ns++) {
87  TS->RHSFunction( (TS->gpu_Udot + stage*TS->u_size_total + TS->u_offsets[ns]),
88  (TS->gpu_U + stage*TS->u_size_total + TS->u_offsets[ns]),
89  &(sim[ns].solver),
90  &(sim[ns].mpi),
91  stagetime );
92  }
93 
94  gpuArraySetValue(TS->gpu_BoundaryFlux + stage*TS->bf_size_total, TS->bf_size_total, 0.0);
95 
96  for (ns = 0; ns < nsims; ns++) {
97  gpuArrayCopy1D( sim[ns].solver.StageBoundaryIntegral,
98  (TS->gpu_BoundaryFlux + stage*TS->bf_size_total + TS->bf_offsets[ns]),
99  TS->bf_sizes[ns] );
100 
101  }
102 
103  }
104 
105  /* Step completion */
106  for (stage = 0; stage < params->nstages; stage++) {
107 
108  for (ns = 0; ns < nsims; ns++) {
109  gpuArrayAXPY( (TS->gpu_Udot + stage*TS->u_size_total + TS->u_offsets[ns]),
110  (TS->dt * params->b[stage]),
111  (sim[ns].solver.gpu_u),
112  (TS->u_sizes[ns]) );
113  gpuArrayAXPY( (TS->gpu_BoundaryFlux + stage*TS->bf_size_total + TS->bf_offsets[ns]),
114  (TS->dt * params->b[stage]),
115  (sim[ns].solver.StepBoundaryIntegral),
116  (TS->bf_sizes[ns]) );
117 
118  }
119 
120  }
121 
122  } else {
123 #endif
124 
125  /* Calculate stage values */
126  for (stage = 0; stage < params->nstages; stage++) {
127 
128  double stagetime = TS->waqt + params->c[stage]*TS->dt;
129 
130  for (ns = 0; ns < nsims; ns++) {
131  _ArrayCopy1D_( sim[ns].solver.u,
132  (TS->U[stage] + TS->u_offsets[ns]),
133  (TS->u_sizes[ns]) );
134  }
135 
136  for (i = 0; i < stage; i++) {
137  _ArrayAXPY_( TS->Udot[i],
138  (TS->dt * params->A[stage*params->nstages+i]),
139  TS->U[stage],
140  TS->u_size_total );
141  }
142 
143  for (ns = 0; ns < nsims; ns++) {
144  if (sim[ns].solver.PreStage) {
145  fprintf(stderr,"ERROR in TimeRK(): Call to solver->PreStage() commented out!\n");
146  return 1;
147  // sim[ns].solver.PreStage( stage,
148  // (TS->U),
149  // &(sim[ns].solver),
150  // &(sim[ns].mpi),
151  // stagetime ); CHECKERR(ierr);
152  }
153  }
154 
155  for (ns = 0; ns < nsims; ns++) {
156  if (sim[ns].solver.PostStage) {
157  sim[ns].solver.PostStage( (TS->U[stage] + TS->u_offsets[ns]),
158  &(sim[ns].solver),
159  &(sim[ns].mpi),
160  stagetime); CHECKERR(ierr);
161  }
162  }
163 
164  for (ns = 0; ns < nsims; ns++) {
165  TS->RHSFunction( (TS->Udot[stage] + TS->u_offsets[ns]),
166  (TS->U[stage] + TS->u_offsets[ns]),
167  &(sim[ns].solver),
168  &(sim[ns].mpi),
169  stagetime);
170  }
171 
172  _ArraySetValue_(TS->BoundaryFlux[stage], TS->bf_size_total, 0.0);
173  for (ns = 0; ns < nsims; ns++) {
174  _ArrayCopy1D_( sim[ns].solver.StageBoundaryIntegral,
175  (TS->BoundaryFlux[stage] + TS->bf_offsets[ns]),
176  TS->bf_sizes[ns] );
177  }
178 
179  }
180 
181  /* Step completion */
182  for (stage = 0; stage < params->nstages; stage++) {
183 
184  for (ns = 0; ns < nsims; ns++) {
185  _ArrayAXPY_( (TS->Udot[stage] + TS->u_offsets[ns]),
186  (TS->dt * params->b[stage]),
187  (sim[ns].solver.u),
188  (TS->u_sizes[ns]) );
189  _ArrayAXPY_( (TS->BoundaryFlux[stage] + TS->bf_offsets[ns]),
190  (TS->dt * params->b[stage]),
191  (sim[ns].solver.StepBoundaryIntegral),
192  (TS->bf_sizes[ns]) );
193  }
194 
195  }
196 
197 #if defined(HAVE_CUDA)
198  }
199 #endif
200 
201  return 0;
202 }
TimeIntegration::RHSFunction
int(* RHSFunction)(double *, double *, void *, void *, double)
Definition: timeintegration_struct.h:98
gpuArrayAXPY
void gpuArrayAXPY(const double *, double, double *, int)
Definition: ArrayImplementations_GPU.cu:177
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::u
double * u
Definition: hypar.h:116
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
gpuArraySetValue
void gpuArraySetValue(double *, int, double)
Definition: ArrayImplementations_GPU.cu:161
HyPar::PostStage
int(* PostStage)(double *, void *, void *, double)
Definition: hypar.h:336
HyPar::StepBoundaryIntegral
double * StepBoundaryIntegral
Definition: hypar.h:384
ExplicitRKParameters
Structure containing the parameters for an explicit Runge-Kutta method.
Definition: timeintegration_struct.h:177
HyPar::gpu_u
double * gpu_u
Definition: hypar.h:459
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
HyPar::msti
void * msti
Definition: hypar.h:366
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
gpuArrayCopy1D
void gpuArrayCopy1D(const double *, double *, int)
Definition: ArrayImplementations_GPU.cu:148
int TimeGLMGEE ( void *  ts)

Take a step in time using the General Linear Methods with Global Error Estimators

Advance the ODE given by

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

by one time step of size HyPar::dt using the \(s\)-stage General Linear Method with Global Error Estimation (GLM-GEE), given by

\begin{align} {\bf U}^{\left(i\right)} &= c_{i0}{\bf u}_n + \sum_{j=1}^{r-1} c_{ij} \tilde{\bf u}_n^j + \Delta t \sum_{j=1}^{i-1} a_{ij} {\bf F}\left({\bf U}^{\left(j\right)}\right), i=1,\cdots,s\\ {\bf u}_{n+1} &= d_{00} {\bf u}_n + \sum_{j=1}^{r-1} d_{0j} \tilde{\bf u}_n^j + \Delta t \sum_{j=1}^s b_{0j} {\bf F}\left({\bf U}^{\left(j\right)}\right), \\ \tilde{\bf u}_{n+1}^i &= d_{i0} {\bf u}_n + \sum_{j=1}^{r-1} d_{ij} \tilde{\bf u}_n^j + \Delta t \sum_{j=1}^s b_{ij} {\bf F}\left({\bf U}^{\left(j\right)}\right), i=1,\cdots,r-1 \end{align}

where the superscripts in parentheses represent stages, the subscripts represent the time level, the superscripts without parentheses represent auxiliary solutions, \(\Delta t\) is the time step size HyPar::dt, \(\tilde{\bf u}^i, i=1,\cdots,r-1\) are the auxiliary solutions, and \({\bf F}\left({\bf u}\right)\) is computed by TimeIntegration::RHSFunction. The coefficients defining this methods are:

where \(s\) is the number of stages (GLMGEEParameters::nstages) and \(r\) is the number of auxiliary solutions propagated with the solution (GLMGEEParameters::r).

Note: In the code TimeIntegration::Udot is equivalent to \({\bf F}\left({\bf u}\right)\).

References:

Parameters
tsObject of type TimeIntegration

Definition at line 45 of file TimeGLMGEE.c.

46 {
47  TimeIntegration* TS = (TimeIntegration*) ts;
48  SimulationObject* sim = (SimulationObject*) TS->simulation;
49  GLMGEEParameters* params = (GLMGEEParameters*) sim[0].solver.msti;
50  int ns, stage, i, j, nsims = TS->nsims;
51  _DECLARE_IERR_;
52 
53  int s = params->nstages;
54  int r = params->r;
55  double dt = TS->dt;
56  double *A =params->A,
57  *B =params->B,
58  *C=params->C,
59  *D=params->D,
60  *c=params->c,
61  **U = TS->U,
62  **Udot = TS->Udot,
63  **Uaux = &TS->U[r];
64 
65  /* Calculate stage values */
66  for (j=0; j<s; j++) {
67 
68  double stagetime = TS->waqt + c[j]*dt;
69 
70  for (ns = 0; ns < nsims; ns++) {
71  _ArrayScaleCopy1D_( sim[ns].solver.u,
72  C[j*r+0],
73  (U[0] + TS->u_offsets[ns]),
74  (TS->u_sizes[ns]) );
75  }
76 
77  for (i=1;i<r;i++) _ArrayAXPY_(Uaux[i-1], C[j*r+i] , U[0], TS->u_size_total);
78  for (i=0;i<j;i++) _ArrayAXPY_(Udot[i] , dt*A[j*s+i], U[0], TS->u_size_total);
79 
80  for (ns = 0; ns < nsims; ns++) {
81  if (sim[ns].solver.PreStage) {
82  fprintf(stderr,"Call to solver->PreStage() commented out in TimeGLMGEE()!\n");
83  return 1;
84 // IERR sim[ns].solver.PreStage( stage,
85 // (U),
86 // &(sim[ns].solver),
87 // &(sim[ns].mpi),
88 // stagetime ); CHECKERR(ierr);
89  }
90  }
91 
92  for (ns = 0; ns < nsims; ns++) {
93  IERR TS->RHSFunction( (Udot[j] + TS->u_offsets[ns]),
94  (U[0] + TS->u_offsets[ns]),
95  &(sim[ns].solver),
96  &(sim[ns].mpi),
97  stagetime);
98  }
99 
100  for (ns = 0; ns < nsims; ns++) {
101  if (sim[ns].solver.PostStage) {
102  IERR sim[ns].solver.PostStage( (U[j] + TS->u_offsets[ns]),
103  &(sim[ns].solver),
104  &(sim[ns].mpi),
105  stagetime); CHECKERR(ierr);
106  }
107  }
108 
109  _ArraySetValue_(TS->BoundaryFlux[j], TS->bf_size_total, 0.0);
110  for (ns = 0; ns < nsims; ns++) {
111  _ArrayCopy1D_( sim[ns].solver.StageBoundaryIntegral,
112  (TS->BoundaryFlux[j] + TS->bf_offsets[ns]),
113  TS->bf_sizes[ns] );
114  }
115  }
116 
117  /* Step completion */
118  for (j=0; j<r; j++) {
119 
120  for (ns = 0; ns < nsims; ns++) {
121 
122  _ArrayScaleCopy1D_( sim[ns].solver.u,
123  D[j*r+0],
124  (U[j]+TS->u_offsets[ns]),
125  (TS->u_sizes[ns]) );
126 
127  }
128 
129  for (i=1; i<r; i++) _ArrayAXPY_(Uaux[i-1], D[j*r+i] , U[j], TS->u_size_total);
130  for (i=0; i<s; i++) _ArrayAXPY_(Udot[i] , dt*B[j*s+i], U[j], TS->u_size_total);
131 
132  }
133 
134  for (ns = 0; ns < nsims; ns++) {
135  for (i=0; i<s; i++) {
136  _ArrayAXPY_( (TS->BoundaryFlux[i] + TS->bf_offsets[ns]),
137  dt*B[0*s+i],
138  sim[ns].solver.StepBoundaryIntegral,
139  TS->bf_sizes[ns] );
140  }
141 
142  }
143 
144  for (ns = 0; ns < nsims; ns++) {
145  _ArrayCopy1D_( (U[0] + TS->u_offsets[ns]),
146  sim[ns].solver.u,
147  TS->u_sizes[ns] );
148  }
149  for (i=1; i<r; i++) {
150  _ArrayCopy1D_(U[i],Uaux[i-1],TS->u_size_total);
151  }
152 
153  return(0);
154 }
TimeIntegration::RHSFunction
int(* RHSFunction)(double *, double *, void *, void *, double)
Definition: timeintegration_struct.h:98
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::u
double * u
Definition: hypar.h:116
_ArrayScaleCopy1D_
#define _ArrayScaleCopy1D_(x, a, y, size)
Definition: arrayfunctions.h:356
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::PostStage
int(* PostStage)(double *, void *, void *, double)
Definition: hypar.h:336
HyPar::StepBoundaryIntegral
double * StepBoundaryIntegral
Definition: hypar.h:384
GLMGEEParameters
Structure containing the parameters for the General Linear Methods with Global Error Estimators (GLM-...
Definition: timeintegration_struct.h:267
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
HyPar::msti
void * msti
Definition: hypar.h:366
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
IERR
#define IERR
Definition: basic.h:16
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17