@@ -306,13 +306,13 @@ void lsrkStep_Free(ARKodeMem ark_mem)
306
306
{
307
307
free (step_mem -> cvals );
308
308
step_mem -> cvals = NULL ;
309
- ark_mem -> lrw -= 5 ;
309
+ ark_mem -> lrw -= step_mem -> nfusedopvecs ;
310
310
}
311
311
if (step_mem -> Xvecs != NULL )
312
312
{
313
313
free (step_mem -> Xvecs );
314
314
step_mem -> Xvecs = NULL ;
315
- ark_mem -> liw -= 5 ;
315
+ ark_mem -> liw -= step_mem -> nfusedopvecs ;
316
316
}
317
317
318
318
/* free the time stepper module itself */
@@ -429,15 +429,15 @@ int lsrkStep_Init(ARKodeMem ark_mem, int init_type)
429
429
/* Allocate reusable arrays for fused vector interface */
430
430
if (step_mem -> cvals == NULL )
431
431
{
432
- step_mem -> cvals = (sunrealtype * )calloc (5 , sizeof (sunrealtype ));
432
+ step_mem -> cvals = (sunrealtype * )calloc (step_mem -> nfusedopvecs , sizeof (sunrealtype ));
433
433
if (step_mem -> cvals == NULL ) { return (ARK_MEM_FAIL ); }
434
- ark_mem -> lrw += 5 ;
434
+ ark_mem -> lrw += step_mem -> nfusedopvecs ;
435
435
}
436
436
if (step_mem -> Xvecs == NULL )
437
437
{
438
- step_mem -> Xvecs = (N_Vector * )calloc (5 , sizeof (N_Vector ));
438
+ step_mem -> Xvecs = (N_Vector * )calloc (step_mem -> nfusedopvecs , sizeof (N_Vector ));
439
439
if (step_mem -> Xvecs == NULL ) { return (ARK_MEM_FAIL ); }
440
- ark_mem -> liw += 5 ; /* pointers */
440
+ ark_mem -> liw += step_mem -> nfusedopvecs ; /* pointers */
441
441
}
442
442
443
443
/* Signal to shared arkode module that full RHS evaluations are required */
@@ -564,7 +564,7 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
564
564
sunrealtype * cvals ;
565
565
sunrealtype w0 , w1 , temp1 , temp2 , arg , bjm1 , bjm2 , mus , thjm1 , thjm2 , zjm1 ,
566
566
zjm2 , dzjm1 , dzjm2 , d2zjm1 , d2zjm2 , zj , dzj , d2zj , bj , ajm1 , mu , nu , thj ;
567
- sunrealtype onep54 = 1.54 , c13 = 13.0 , p8 = 0.8 , p4 = 0.4 ;
567
+ static sunrealtype onep54 = SUN_RCONST ( 1.54 ) , c13 = SUN_RCONST ( 13.0 ) , p8 = SUN_RCONST ( 0.8 ) , p4 = SUN_RCONST ( 0.4 ) ;
568
568
N_Vector * Xvecs ;
569
569
ARKodeLSRKStepMem step_mem ;
570
570
@@ -673,7 +673,7 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
673
673
cvals [4 ] = - mus * ajm1 * ark_mem -> h ;
674
674
Xvecs [4 ] = ark_mem -> fn ;
675
675
676
- retval = N_VLinearCombination (5 , cvals , Xvecs , ark_mem -> ycur );
676
+ retval = N_VLinearCombination (step_mem -> nfusedopvecs , cvals , Xvecs , ark_mem -> ycur );
677
677
if (retval != 0 ) { return (ARK_VECTOROP_ERR ); }
678
678
679
679
thj = mu * thjm1 + nu * thjm2 + mus * (ONE - ajm1 );
@@ -757,7 +757,7 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
757
757
int retval ;
758
758
sunrealtype * cvals ;
759
759
sunrealtype w1 , bjm1 , bjm2 , mus , bj , ajm1 , cjm1 , temj , cj , mu , nu ;
760
- sunrealtype p8 = 0.8 , p4 = 0.4 ;
760
+ static sunrealtype p8 = SUN_RCONST ( 0.8 ) , p4 = SUN_RCONST ( 0.4 ) ;
761
761
N_Vector * Xvecs ;
762
762
ARKodeLSRKStepMem step_mem ;
763
763
@@ -850,7 +850,7 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
850
850
cvals [4 ] = - mus * ajm1 * ark_mem -> h ;
851
851
Xvecs [4 ] = ark_mem -> fn ;
852
852
853
- retval = N_VLinearCombination (5 , cvals , Xvecs , ark_mem -> ycur );
853
+ retval = N_VLinearCombination (step_mem -> nfusedopvecs , cvals , Xvecs , ark_mem -> ycur );
854
854
if (retval != 0 ) { return (ARK_VECTOROP_ERR ); }
855
855
856
856
/* Shift the data for the next stage */
@@ -938,10 +938,11 @@ int lsrkStep_TakeStepSSPs2(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr
938
938
sunrealtype sm1inv = ONE / (rs - ONE );
939
939
sunrealtype bt1 , bt2 , bt3 ;
940
940
941
+ /* Embedding coefficients differ when reqstages == 2 */
941
942
if (step_mem -> reqstages == 2 )
942
943
{
943
- bt1 = 0.694021459207626 ;
944
- bt3 = 1 - 0.694021459207626 ;
944
+ bt1 = SUN_RCONST ( 0.694021459207626 ) ;
945
+ bt3 = ONE - bt1 ;
945
946
}
946
947
else
947
948
{
@@ -1183,7 +1184,9 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
1183
1184
* nflagPtr = ARK_SUCCESS ;
1184
1185
* dsmPtr = ZERO ;
1185
1186
1186
- sunrealtype onesixth = 1.0 / 6.0 ;
1187
+ static sunrealtype onesixth = ONE / SIX , onefifth = ONE / FIVE , p1 = SUN_RCONST (0.1 ), p3 = SUN_RCONST (0.3 ),
1188
+ p6 = SUN_RCONST (0.6 ), onetwntyfth = SUN_RCONST (1.0 /25.0 ), onetwntynth = SUN_RCONST (9.0 /25.0 ),
1189
+ fifteen = SUN_RCONST (15.0 ), negfive = SUN_RCONST (-5.0 );
1187
1190
1188
1191
/* access ARKodeLSRKStepMem structure */
1189
1192
retval = lsrkStep_AccessStepMem (ark_mem , __func__ , & step_mem );
@@ -1211,7 +1214,7 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
1211
1214
1212
1215
N_VLinearSum (ONE , ark_mem -> yn , onesixth * ark_mem -> h , ark_mem -> fn ,
1213
1216
ark_mem -> ycur );
1214
- N_VLinearSum (ONE , ark_mem -> yn , 1.0 / 5.0 * ark_mem -> h , ark_mem -> fn ,
1217
+ N_VLinearSum (ONE , ark_mem -> yn , onefifth * ark_mem -> h , ark_mem -> fn ,
1215
1218
ark_mem -> tempv1 );
1216
1219
1217
1220
/* Evaluate stages j = 2,...,step_mem->reqstages */
@@ -1226,17 +1229,17 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
1226
1229
ark_mem -> ycur );
1227
1230
if (j == 4 )
1228
1231
{
1229
- N_VLinearSum (ONE , ark_mem -> tempv1 , 3.0 / 10.0 * ark_mem -> h , ark_mem -> fn ,
1232
+ N_VLinearSum (ONE , ark_mem -> tempv1 , p3 * ark_mem -> h , ark_mem -> fn ,
1230
1233
ark_mem -> tempv1 );
1231
1234
}
1232
1235
}
1233
- N_VLinearSum (1.0 / 25.0 , ark_mem -> tempv2 , 9.0 / 25.0 , ark_mem -> ycur ,
1236
+ N_VLinearSum (onetwntyfth , ark_mem -> tempv2 , onetwntynth , ark_mem -> ycur ,
1234
1237
ark_mem -> tempv2 );
1235
- N_VLinearSum (15 , ark_mem -> tempv2 , -5 , ark_mem -> ycur , ark_mem -> ycur );
1238
+ N_VLinearSum (fifteen , ark_mem -> tempv2 , negfive , ark_mem -> ycur , ark_mem -> ycur );
1236
1239
for (int j = 6 ; j <= 9 ; j ++ )
1237
1240
{
1238
1241
retval = step_mem -> fe (ark_mem -> tcur +
1239
- ((sunrealtype )j - 4.0 ) * onesixth * ark_mem -> h ,
1242
+ ((sunrealtype )j - FOUR ) * onesixth * ark_mem -> h ,
1240
1243
ark_mem -> ycur , ark_mem -> fn , ark_mem -> user_data );
1241
1244
if (retval != ARK_SUCCESS ) { return (ARK_RHSFUNC_FAIL ); }
1242
1245
@@ -1245,12 +1248,12 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
1245
1248
1246
1249
if (j == 7 )
1247
1250
{
1248
- N_VLinearSum (ONE , ark_mem -> tempv1 , 1.0 / 5.0 * ark_mem -> h , ark_mem -> fn ,
1251
+ N_VLinearSum (ONE , ark_mem -> tempv1 , onefifth * ark_mem -> h , ark_mem -> fn ,
1249
1252
ark_mem -> tempv1 );
1250
1253
}
1251
1254
if (j == 9 )
1252
1255
{
1253
- N_VLinearSum (ONE , ark_mem -> tempv1 , 3.0 / 10.0 * ark_mem -> h , ark_mem -> fn ,
1256
+ N_VLinearSum (ONE , ark_mem -> tempv1 , p3 * ark_mem -> h , ark_mem -> fn ,
1254
1257
ark_mem -> tempv1 );
1255
1258
}
1256
1259
}
@@ -1259,8 +1262,8 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
1259
1262
ark_mem -> user_data );
1260
1263
if (retval != ARK_SUCCESS ) { return (ARK_RHSFUNC_FAIL ); }
1261
1264
1262
- N_VLinearSum (ONE , ark_mem -> tempv2 , 3.0 / 5.0 , ark_mem -> ycur , ark_mem -> ycur );
1263
- N_VLinearSum (ONE , ark_mem -> ycur , 1.0 / 10.0 * ark_mem -> h , ark_mem -> fn ,
1265
+ N_VLinearSum (ONE , ark_mem -> tempv2 , p6 , ark_mem -> ycur , ark_mem -> ycur );
1266
+ N_VLinearSum (ONE , ark_mem -> ycur , p1 * ark_mem -> h , ark_mem -> fn ,
1264
1267
ark_mem -> ycur );
1265
1268
1266
1269
step_mem -> nfe += 9 ;
0 commit comments