The FORTRAN example problem programs supplied with the IDA package are all written in standard FORTRAN77 and use double precision arithmetic. However, when the FORTRAN examples are built, the source code is automatically modified according to the configure options supplied by the user and the system type. Integer variables are declared as INTEGER*n, where n denotes the number of bytes in the corresponding C type (long int or int). Floating-point variable declarations remain unchanged if double precision is used, but are changed to REAL*n, where n denotes the number of bytes in the SUNDIALS type realtype, if using single-precision. Also, if using single-precision, then declarations of floating-point constants are appropriately modified; e.g. 0.5D-4 is changed to 0.5E-4.
The fidadenx example is a FORTRAN equivalent of the idadenx problem. The source program fidadenx.f is listed in Appendix F.
The main program begins with declarations and initializations. It calls the routines FNVINITS, FIDAMALLOC, FIDAROOTINIT, FIDADENSE, and FIDADENSESETJAC, to initialize the NVECTOR-SERIAL module, the main solver memory, the rootfinding module, and the IDADENSE module, and to specify user-supplied Jacobian routine, respectively. It calls FIDASOLVE in a loop over TOUT values, with printing of the solution values and performance data (current order and step count from the IOUT array, and current step size from the ROUT array). In the case of a root return, an extra line is printed with the root information from FIDAROOTINFO. At the end, it prints a number of performance counters, and frees memory with calls to FIDAROOTFREE and FIDAFREE.
In fidadenx.f, the FIDARESFUN routine is a straghtforward implementation of Eqns. (1). In FIDADJAC, the 3×3 system Jacobian is supplied. The FIDAROOTFN routine defines the two root functions, which are set to determine the points at which y1 = 10-4 or y3 = .01. The final two routines are for the printing of a header and of the final run statistics.
The following is sample output from fidadenx. The performance of FIDA here is similar to that of IDA on the idadenx problem, with somewhat lower cost counters owing to the larger absolute error tolerances.
fidadenx: Robertson kinetics DAE serial exampleproblem for IDA
Three equation chemicalkinetics problem.
Tolerance parameters: rtol = 0.10E-03 atol = 0.10E-05 0.10E-09 0.10E-05
Initial conditions y0 = ( 0.10E+01 0.00E+00 0.00E+00)
t y1 y2 y3 nst k h
0.2640E+00 0.9900E+00 0.3471E-04 0.1000E-01 75 2 0.5716E-01
Above is a root, INFO() = 0 1
0.4000E+00 0.9852E+00 0.3386E-04 0.1480E-01 77 3 0.1143E+00
0.4000E+01 0.9055E+00 0.2240E-04 0.9447E-01 91 4 0.3704E+00
0.4000E+02 0.7158E+00 0.9185E-05 0.2842E+00 127 4 0.2963E+01
0.4000E+03 0.4505E+00 0.3223E-05 0.5495E+00 177 3 0.1241E+02
0.4000E+04 0.1832E+00 0.8940E-06 0.8168E+00 228 3 0.2765E+03
0.4000E+05 0.3899E-01 0.1622E-06 0.9610E+00 278 5 0.2614E+04
0.4000E+06 0.4939E-02 0.1985E-07 0.9951E+00 324 5 0.2770E+05
0.4000E+07 0.5176E-03 0.2072E-08 0.9995E+00 355 4 0.3979E+06
0.2075E+08 0.1000E-03 0.4000E-09 0.9999E+00 374 4 0.1592E+07
Above is a root, INFO() = 1 0
0.4000E+08 0.5191E-04 0.2076E-09 0.9999E+00 380 3 0.6366E+07
0.4000E+09 0.5882E-05 0.2353E-10 0.1000E+01 394 1 0.9167E+08
0.4000E+10 0.7054E-06 0.2822E-11 0.1000E+01 402 1 0.1467E+10
0.4000E+11 -0.7300E-06 -0.2920E-11 0.1000E+01 407 1 0.2347E+11
Final Run Statistics:
Number of steps = 407
Number of residual evaluations = 557
Number of Jacobian evaluations = 65
Number of nonlinear iterations = 557
Number of error test failures = 6
Number of nonlinear conv. failures = 0
Number of root function evals. = 437
This example, fidakryx_bbd_p, is the FORTRAN equivalent of the idakryx1_bbd_p example. The heat equation problem is described under the idakryx example above, but here it is solved in parallel, using the IDABBDPRE (band-block-diagonal) preconditioner module. The decomposition of the problem onto a processor array is identical to that in the idakryx1_p example above. The source file, fidakryx_bbd_p.f, is listed in Appendix G.
The problem is solved twice -- once with half-bandwidths of 5 in the difference-quotient banded preconditioner blocks, and once with half-bandwidths of 1 (which results in lumping of Jacobian values). In both cases, the retained banded blocks are tridiagonal, even though the true Jacobian is not.
The main program begins with initializations, including MPI calls, a call to FNVINITP to initialize NVECTOR-PARALLEL, and a call to SETINITPROFILE to initialize the UU, UP, ID, and CONSTR arrays (containing the solution vector, solution derivative vector, the differential/algebraic bit vector, and the contraint specification vector, respectively). A call to FIDASETIIN and two calls to FIDASETVIN are made to suppress error control on the algebraic variables, and to supply the ID array and constraints array (making the computed solution non-negative). The call to FIDAMALLOC initializes the FIDA main memory, and the calls to FIDABBDINIT and FIDABBDSPGMR initialize the FIDABBD module.
In the first loop over TOUT values, the main program calls FIDASOLVE and prints the max-norm of the solution and selected counters. When finished, it calls PRNTFINALSTATS to print a few more counters.
The second solution is initialized by a second call to SETINITPROFILE, and calls to FIDAREINIT and FIDABBDREINIT. After completing the second solution, the program frees memory and terminates MPI.
The FIDARESFUN routine simply calls two other routines: FIDACOMMFN, to communicate boundary needed data from U to an extension of it called UEXT; and FIDAGLOCFN, to compute the residuals in terms of UEXT and UP.
The following is a sample output from fidakryx_bbd_p, with a 10×10 mesh and NPES = 4 processors. The performance is similar for the two solutions. The second case requires more linear iterations, as expected, but their cost is offset by the much cheaper preconditioner evaluations.
fidakryx_bbd_p: Heat equation, parallel example problem for FIDA
Discretized heat equation on 2D unit square.
Zero boundary conditions, polynomial conditions.
Mesh dimensions: 10 x 10 Total system size: 100
Subgrid dimensions: 5 x 5 Processor array: 2 x 2
Tolerance parameters: rtol = 0.00E+00 atol = 0.10E-02
Constraints set to force all solution components >= 0.
SUPPRESSALG = TRUE to remove boundary components from the error test.
Linear solver: SPGMR. Preconditioner: BBDPRE - Banded-block-diagonal.
Case 1
Difference quotient half-bandwidths = 5
Retained matrix half-bandwidths = 1
Output Summary
umax = max-norm of solution
nre = nre + nreLS (total number of RES evals.)
time umax k nst nni nli nre nge h npe nps
--------------------------------------------------------------------------
0.1000E-01 0.82411E+00 2 12 14 7 14+ 7 96 0.26E-02 8 21
0.2000E-01 0.68812E+00 3 15 18 12 18+12 96 0.51E-02 8 30
0.4000E-01 0.47075E+00 3 18 24 22 24+22 108 0.66E-02 9 46
0.8000E-01 0.21660E+00 3 22 29 30 29+30 108 0.13E-01 9 59
0.1600E+00 0.45659E-01 4 28 37 43 37+43 120 0.26E-01 10 80
0.3200E+00 0.21096E-02 4 35 45 59 45+59 120 0.24E-01 10 104
0.6400E+00 0.55368E-04 1 40 54 71 54+71 156 0.19E+00 13 125
0.1280E+01 0.15597E-18 1 42 56 71 56+71 180 0.76E+00 15 127
0.2560E+01 0.33865E-20 1 43 57 71 57+71 192 0.15E+01 16 128
0.5120E+01 0.86074E-20 1 44 58 71 58+71 204 0.30E+01 17 129
0.1024E+02 0.16630E-19 1 45 59 71 59+71 216 0.61E+01 18 130
Error test failures = 1
Nonlinear convergence failures = 0
Linear convergence failures = 0
Case 2
Difference quotient half-bandwidths = 1
Retained matrix half-bandwidths = 1
Output Summary
umax = max-norm of solution
nre = nre + nreLS (total number of RES evals.)
time umax k nst nni nli nre nge h npe nps
--------------------------------------------------------------------------
0.1000E-01 0.82411E+00 2 12 14 7 14+ 7 32 0.26E-02 8 21
0.2000E-01 0.68812E+00 3 15 18 12 18+12 32 0.51E-02 8 30
0.4000E-01 0.47093E+00 3 19 23 20 23+20 36 0.10E-01 9 43
0.8000E-01 0.21655E+00 3 23 27 32 27+32 36 0.10E-01 9 59
0.1600E+00 0.45225E-01 4 27 33 44 33+44 40 0.20E-01 10 77
0.3200E+00 0.21868E-02 3 34 41 67 41+67 44 0.41E-01 11 108
0.6400E+00 0.48847E-18 1 39 49 86 49+86 52 0.16E+00 13 135
0.1280E+01 0.53982E-18 1 41 51 86 51+86 60 0.66E+00 15 137
0.2560E+01 0.74194E-17 1 42 52 86 52+86 64 0.13E+01 16 138
0.5120E+01 0.61081E-16 1 43 53 86 53+86 68 0.26E+01 17 139
0.1024E+02 0.40536E-15 1 44 54 86 54+86 72 0.52E+01 18 140
Error test failures = 0
Nonlinear convergence failures = 0
Linear convergence failures = 0