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
TimeGLMGEEInitialize.c File Reference

Initialize the _GLM_GEE_ time integrator. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <matops.h>
#include <arrayfunctions.h>
#include <mpivars.h>
#include <timeintegration.h>

Go to the source code of this file.

Functions

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

Detailed Description

Initialize the _GLM_GEE_ time integrator.

Author
Debojyoti Ghosh

Definition in file TimeGLMGEEInitialize.c.

Function Documentation

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

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