1C ---------------------------------------------------------------- 2C Programmer(s): Daniel R. Reynolds @ SMU 3C ---------------------------------------------------------------- 4C SUNDIALS Copyright Start 5C Copyright (c) 2002-2020, Lawrence Livermore National Security 6C and Southern Methodist University. 7C All rights reserved. 8C 9C See the top-level LICENSE and NOTICE files for details. 10C 11C SPDX-License-Identifier: BSD-3-Clause 12C SUNDIALS Copyright End 13C ---------------------------------------------------------------- 14C Diagonal ODE example. Nonstiff case: alpha = 10/NEQ. 15C ---------------------------------------------------------------- 16C 17C Include MPI-Fortran header file for MPI_COMM_WORLD, MPI types. 18C 19 IMPLICIT NONE 20C 21 INCLUDE "mpif.h" 22 INTEGER*8 NLOCAL 23 PARAMETER (NLOCAL=2) 24C 25 INTEGER*4 IER, MYPE, NPES, NOUT, LNST, LNST_ATT, LNFE, LNFI, LNNI 26 INTEGER*4 LNCF, LNETF, METH, IATOL, ITASK, JOUT 27 INTEGER*8 NEQ, I, NST, NST_ATT, NFE, NFI, NNI, NCFN, NETF 28 INTEGER*8 IOUT(35), IPAR(2) 29 DOUBLE PRECISION Y(128), ROUT(6), RPAR(1) 30 DOUBLE PRECISION ATOL, RTOL, DTOUT, T, ALPHA, TOUT 31 DOUBLE PRECISION ERMAX, ERRI, GERMAX 32C 33 DATA ATOL/1.0D-10/, RTOL/1.0D-5/, DTOUT/0.1D0/, NOUT/10/ 34 DATA LNST/3/, LNST_ATT/6/, LNFE/7/, LNFI/8/, LNNI/11/, LNCF/12/, 35 1 LNETF/10/ 36C 37C Get NPES and MYPE. Requires initialization of MPI. 38 CALL MPI_INIT(IER) 39 IF (IER .NE. 0) THEN 40 WRITE(6,5) IER 41 5 FORMAT(///' MPI_ERROR: MPI_INIT returned IER = ', I5) 42 STOP 43 ENDIF 44 CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPES, IER) 45 IF (IER .NE. 0) THEN 46 WRITE(6,6) IER 47 6 FORMAT(///' MPI_ERROR: MPI_COMM_SIZE returned IER = ', I5) 48 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 49 STOP 50 ENDIF 51 CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYPE, IER) 52 IF (IER .NE. 0) THEN 53 WRITE(6,7) IER 54 7 FORMAT(///' MPI_ERROR: MPI_COMM_RANK returned IER = ', I5) 55 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 56 STOP 57 ENDIF 58C 59C Set input arguments. 60 NEQ = NPES * NLOCAL 61 T = 0.0D0 62 METH = 1 63 IATOL = 1 64 ITASK = 1 65c Set parameter ALPHA 66 ALPHA = 10.0D0 / NEQ 67C 68C Load IPAR and RPAR 69 IPAR(1) = NLOCAL 70 IPAR(2) = MYPE 71 RPAR(1) = ALPHA 72C 73 DO 10 I = 1, NLOCAL 74 10 Y(I) = 1.0D0 75C 76 IF (MYPE .EQ. 0) THEN 77 WRITE(6,11) NEQ, NLOCAL, ALPHA 78 11 FORMAT('Diagonal test problem:'//' NEQ = ', I3, / 79 1 ' NLOCAL = ', I3, / 80 2 ' parameter alpha = ', F8.3) 81 WRITE(6,12) 82 12 FORMAT(' ydot_i = -alpha*i * y_i (i = 1,...,NEQ)') 83 WRITE(6,13) RTOL, ATOL 84 13 FORMAT(' RTOL, ATOL = ', 2E10.1) 85 WRITE(6,14) 86 14 FORMAT(' Method is ERK') 87 WRITE(6,15) NPES 88 15 FORMAT(' Number of processors = ', I3//) 89 ENDIF 90C 91 CALL FNVINITP(MPI_COMM_WORLD, 4, NLOCAL, NEQ, IER) 92C 93 IF (IER .NE. 0) THEN 94 WRITE(6,20) IER 95 20 FORMAT(///' SUNDIALS_ERROR: FNVINITP returned IER = ', I5) 96 CALL MPI_FINALIZE(IER) 97 STOP 98 ENDIF 99C 100 CALL FARKMALLOC(T, Y, METH, IATOL, RTOL, ATOL, 101 1 IOUT, ROUT, IPAR, RPAR, IER) 102C 103 IF (IER .NE. 0) THEN 104 WRITE(6,30) IER 105 30 FORMAT(///' SUNDIALS_ERROR: FARKMALLOC returned IER = ', I5) 106 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 107 STOP 108 ENDIF 109C 110C Loop through tout values, call solver, print output, test for failure. 111 TOUT = DTOUT 112 DO 70 JOUT = 1, NOUT 113C 114 CALL FARKODE(TOUT, T, Y, ITASK, IER) 115C 116 IF (MYPE .EQ. 0) WRITE(6,40) T, IOUT(LNST), IOUT(LNST_ATT), 117 & IOUT(LNFE), IOUT(LNFI) 118 40 FORMAT(' t = ', D10.2, 5X, 'steps = ', I5, 119 & ' (attempted = ', I5, '), fe = ', I5, 120 & ' fi = ', I5) 121C 122 IF (IER .NE. 0) THEN 123 WRITE(6,60) IER, IOUT(16) 124 60 FORMAT(///' SUNDIALS_ERROR: FARKODE returned IER = ', I5, /, 125 & ' Linear Solver returned IER = ', I5) 126 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 127 STOP 128 ENDIF 129C 130 TOUT = TOUT + DTOUT 131 70 CONTINUE 132C 133C Get max. absolute error in the local vector. 134 ERMAX = 0.0D0 135 DO 75 I = 1, NLOCAL 136 ERRI = Y(I) - EXP(-ALPHA * (MYPE * NLOCAL + I) * T) 137 ERMAX = MAX(ERMAX, ABS(ERRI)) 138 75 CONTINUE 139C Get global max. error from MPI_REDUCE call. 140 CALL MPI_REDUCE(ERMAX, GERMAX, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 141 1 0, MPI_COMM_WORLD, IER) 142 IF (IER .NE. 0) THEN 143 WRITE(6,80) IER 144 80 FORMAT(///' MPI_ERROR: MPI_REDUCE returned IER = ', I5) 145 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 146 STOP 147 ENDIF 148 IF (MYPE .EQ. 0) WRITE(6,85) GERMAX 149 85 FORMAT(/'Max. absolute error is ', E10.2/) 150C 151C Print final statistics. 152 NST = IOUT(LNST) 153 NST_ATT = IOUT(LNST_ATT) 154 NFE = IOUT(LNFE) 155 NFI = IOUT(LNFI) 156 NNI = IOUT(LNNI) 157 NCFN = IOUT(LNCF) 158 NETF = IOUT(LNETF) 159 IF (MYPE .EQ. 0) WRITE (6,90) NST, NST_ATT, NFE, NFI, NNI, NCFN, 160 & NETF 161 90 FORMAT(/'Final statistics:'// 162 & ' number of steps = ', I5/ 163 & ' number of step attempts = ', I5/ 164 & ' number of fe evals. = ', I5/ 165 & ' number of fi evals. = ', I5/ 166 & ' number of nonlinear iters. = ', I5/ 167 & ' number of nonlinear conv. failures = ', I3/ 168 & ' number of error test failures = ', I3) 169C 170C Free the memory and finalize MPI. 171 CALL FARKFREE 172 CALL MPI_FINALIZE(IER) 173 IF (IER .NE. 0) THEN 174 WRITE(6,95) IER 175 95 FORMAT(///' MPI_ERROR: MPI_FINALIZE returned IER = ', I5) 176 STOP 177 ENDIF 178C 179 STOP 180 END 181C 182C ------------------------------------------------------------------------ 183C 184 SUBROUTINE FARKIFUN(T, Y, YDOT, IPAR, RPAR, IER) 185C Routine for right-hand side function fi 186C 187 IMPLICIT NONE 188C 189 INTEGER*4 IER 190 INTEGER*8 IPAR(*) 191 DOUBLE PRECISION T, Y(*), YDOT(*), RPAR(*) 192C 193 INTEGER*8 NLOCAL, I 194C 195 NLOCAL = IPAR(1) 196C 197 DO I = 1, NLOCAL 198 YDOT(I) = 0.D0 199 ENDDO 200C 201 IER = 0 202C 203 RETURN 204 END 205C ------------------------------------------------------------------------ 206C 207 SUBROUTINE FARKEFUN(T, Y, YDOT, IPAR, RPAR, IER) 208C Routine for right-hand side function fe 209C 210 IMPLICIT NONE 211C 212 INTEGER*4 IER 213 INTEGER*8 IPAR(*) 214 DOUBLE PRECISION T, Y(*), YDOT(*), RPAR(*) 215C 216 INTEGER*8 MYPE 217 INTEGER*8 NLOCAL, I 218 DOUBLE PRECISION ALPHA 219C 220 NLOCAL = IPAR(1) 221 MYPE = IPAR(2) 222 ALPHA = RPAR(1) 223C 224 DO I = 1, NLOCAL 225 YDOT(I) = -ALPHA * (MYPE * NLOCAL + I) * Y(I) 226 ENDDO 227C 228 IER = 0 229C 230 RETURN 231 END 232