00001 C Work performed under the auspices of the U.S. Department of Energy 00002 C by Lawrence Livermore National Laboratory under contract number 00003 C W-7405-Eng-48. 00004 C 00005 SUBROUTINE DDASIC (X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL, 00006 * H, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E, YIC, YPIC, 00007 * PWK, WM, IWM, HMIN, UROUND, EPLI, SQRTN, RSQRTN, EPCONI, 00008 * STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC) 00009 C 00010 C***BEGIN PROLOGUE DDASIC 00011 C***REFER TO DDASPK 00012 C***DATE WRITTEN 940628 (YYMMDD) 00013 C***REVISION DATE 941206 (YYMMDD) 00014 C***REVISION DATE 950714 (YYMMDD) 00015 C 00016 C----------------------------------------------------------------------- 00017 C***DESCRIPTION 00018 C 00019 C DDASIC is a driver routine to compute consistent initial values 00020 C for Y and YPRIME. There are two different options: 00021 C Denoting the differential variables in Y by Y_d, and 00022 C the algebraic variables by Y_a, the problem solved is either: 00023 C 1. Given Y_d, calculate Y_a and Y_d', or 00024 C 2. Given Y', calculate Y. 00025 C In either case, initial values for the given components 00026 C are input, and initial guesses for the unknown components 00027 C must also be provided as input. 00028 C 00029 C The external routine NLSIC solves the resulting nonlinear system. 00030 C 00031 C The parameters represent 00032 C 00033 C X -- Independent variable. 00034 C Y -- Solution vector at X. 00035 C YPRIME -- Derivative of solution vector. 00036 C NEQ -- Number of equations to be integrated. 00037 C ICOPT -- Flag indicating initial condition option chosen. 00038 C ICOPT = 1 for option 1 above. 00039 C ICOPT = 2 for option 2. 00040 C ID -- Array of dimension NEQ, which must be initialized 00041 C if option 1 is chosen. 00042 C ID(i) = +1 if Y_i is a differential variable, 00043 C ID(i) = -1 if Y_i is an algebraic variable. 00044 C RES -- External user-supplied subroutine to evaluate the 00045 C residual. See RES description in DDASPK prologue. 00046 C JAC -- External user-supplied routine to update Jacobian 00047 C or preconditioner information in the nonlinear solver 00048 C (optional). See JAC description in DDASPK prologue. 00049 C PSOL -- External user-supplied routine to solve 00050 C a linear system using preconditioning. 00051 C See PSOL in DDASPK prologue. 00052 C H -- Scaling factor in iteration matrix. DDASIC may 00053 C reduce H to achieve convergence. 00054 C WT -- Vector of weights for error criterion. 00055 C NIC -- Input number of initial condition calculation call 00056 C (= 1 or 2). 00057 C IDID -- Completion code. See IDID in DDASPK prologue. 00058 C RPAR,IPAR -- Real and integer parameter arrays that 00059 C are used for communication between the 00060 C calling program and external user routines. 00061 C They are not altered by DNSK 00062 C PHI -- Work space for DDASIC of length at least 2*NEQ. 00063 C SAVR -- Work vector for DDASIC of length NEQ. 00064 C DELTA -- Work vector for DDASIC of length NEQ. 00065 C E -- Work vector for DDASIC of length NEQ. 00066 C YIC,YPIC -- Work vectors for DDASIC, each of length NEQ. 00067 C PWK -- Work vector for DDASIC of length NEQ. 00068 C WM,IWM -- Real and integer arrays storing 00069 C information required by the linear solver. 00070 C EPCONI -- Test constant for Newton iteration convergence. 00071 C ICNFLG -- Flag showing whether constraints on Y are to apply. 00072 C ICNSTR -- Integer array of length NEQ with constraint types. 00073 C 00074 C The other parameters are for use internally by DDASIC. 00075 C 00076 C----------------------------------------------------------------------- 00077 C***ROUTINES CALLED 00078 C DCOPY, NLSIC 00079 C 00080 C***END PROLOGUE DDASIC 00081 C 00082 C 00083 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00084 DIMENSION Y(*),YPRIME(*),ID(*),WT(*),PHI(NEQ,*) 00085 DIMENSION SAVR(*),DELTA(*),E(*),YIC(*),YPIC(*),PWK(*) 00086 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*), ICNSTR(*) 00087 EXTERNAL RES, JAC, PSOL, NLSIC 00088 C 00089 PARAMETER (LCFN=15) 00090 PARAMETER (LMXNH=34) 00091 C 00092 C The following parameters are data-loaded here: 00093 C RHCUT = factor by which H is reduced on retry of Newton solve. 00094 C RATEMX = maximum convergence rate for which Newton iteration 00095 C is considered converging. 00096 C 00097 SAVE RHCUT, RATEMX 00098 DATA RHCUT/0.1D0/, RATEMX/0.8D0/ 00099 C 00100 C 00101 C----------------------------------------------------------------------- 00102 C BLOCK 1. 00103 C Initializations. 00104 C JSKIP is a flag set to 1 when NIC = 2 and NH = 1, to signal that 00105 C the initial call to the JAC routine is to be skipped then. 00106 C Save Y and YPRIME in PHI. Initialize IDID, NH, and CJ. 00107 C----------------------------------------------------------------------- 00108 C 00109 MXNH = IWM(LMXNH) 00110 IDID = 1 00111 NH = 1 00112 JSKIP = 0 00113 IF (NIC .EQ. 2) JSKIP = 1 00114 CALL DCOPY (NEQ, Y, 1, PHI(1,1), 1) 00115 CALL DCOPY (NEQ, YPRIME, 1, PHI(1,2), 1) 00116 C 00117 IF (ICOPT .EQ. 2) THEN 00118 CJ = 0.0D0 00119 ELSE 00120 CJ = 1.0D0/H 00121 ENDIF 00122 C 00123 C----------------------------------------------------------------------- 00124 C BLOCK 2 00125 C Call the nonlinear system solver to obtain 00126 C consistent initial values for Y and YPRIME. 00127 C----------------------------------------------------------------------- 00128 C 00129 200 CONTINUE 00130 CALL NLSIC(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JAC,PSOL,H,WT,JSKIP, 00131 * RPAR,IPAR,SAVR,DELTA,E,YIC,YPIC,PWK,WM,IWM,CJ,UROUND, 00132 * EPLI,SQRTN,RSQRTN,EPCONI,RATEMX,STPTOL,JFLG,ICNFLG,ICNSTR, 00133 * IERNLS) 00134 C 00135 IF (IERNLS .EQ. 0) RETURN 00136 C 00137 C----------------------------------------------------------------------- 00138 C BLOCK 3 00139 C The nonlinear solver was unsuccessful. Increment NCFN. 00140 C Return with IDID = -12 if either 00141 C IERNLS = -1: error is considered unrecoverable, 00142 C ICOPT = 2: we are doing initialization problem type 2, or 00143 C NH = MXNH: the maximum number of H values has been tried. 00144 C Otherwise (problem 1 with IERNLS .GE. 1), reduce H and try again. 00145 C If IERNLS > 1, restore Y and YPRIME to their original values. 00146 C----------------------------------------------------------------------- 00147 C 00148 IWM(LCFN) = IWM(LCFN) + 1 00149 JSKIP = 0 00150 C 00151 IF (IERNLS .EQ. -1) GO TO 350 00152 IF (ICOPT .EQ. 2) GO TO 350 00153 IF (NH .EQ. MXNH) GO TO 350 00154 C 00155 NH = NH + 1 00156 H = H*RHCUT 00157 CJ = 1.0D0/H 00158 C 00159 IF (IERNLS .EQ. 1) GO TO 200 00160 C 00161 CALL DCOPY (NEQ, PHI(1,1), 1, Y, 1) 00162 CALL DCOPY (NEQ, PHI(1,2), 1, YPRIME, 1) 00163 GO TO 200 00164 C 00165 350 IDID = -12 00166 RETURN 00167 C 00168 C------END OF SUBROUTINE DDASIC----------------------------------------- 00169 END
1.7.1
