1 /* -----------------------------------------------------------------------------
2 * Programmer(s): David J. Gardner @ LLNL
3 * -----------------------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2020, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 * -----------------------------------------------------------------------------
14 * This the implementation file for the IDA nonlinear solver interface.
15 * ---------------------------------------------------------------------------*/
16
17 #include "idas_impl.h"
18 #include "sundials/sundials_math.h"
19
20 /* constant macros */
21 #define PT0001 RCONST(0.0001) /* real 0.0001 */
22 #define ONE RCONST(1.0) /* real 1.0 */
23 #define TWENTY RCONST(20.0) /* real 20.0 */
24
25 /* nonlinear solver parameters */
26 #define MAXIT 4 /* default max number of nonlinear iterations */
27 #define RATEMAX RCONST(0.9) /* max convergence rate used in divergence check */
28
29 /* private functions passed to nonlinear solver */
30 static int idaNlsResidual(N_Vector ycor, N_Vector res, void* ida_mem);
31 static int idaNlsLSetup(booleantype jbad, booleantype* jcur, void* ida_mem);
32 static int idaNlsLSolve(N_Vector delta, void* ida_mem);
33 static int idaNlsConvTest(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector del,
34 realtype tol, N_Vector ewt, void* ida_mem);
35
36 /* -----------------------------------------------------------------------------
37 * Exported functions
38 * ---------------------------------------------------------------------------*/
39
IDASetNonlinearSolver(void * ida_mem,SUNNonlinearSolver NLS)40 int IDASetNonlinearSolver(void *ida_mem, SUNNonlinearSolver NLS)
41 {
42 IDAMem IDA_mem;
43 int retval;
44
45 /* return immediately if IDA memory is NULL */
46 if (ida_mem == NULL) {
47 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS",
48 "IDASetNonlinearSolver", MSG_NO_MEM);
49 return(IDA_MEM_NULL);
50 }
51 IDA_mem = (IDAMem) ida_mem;
52
53 /* return immediately if NLS memory is NULL */
54 if (NLS == NULL) {
55 IDAProcessError(NULL, IDA_ILL_INPUT, "IDAS",
56 "IDASetNonlinearSolver",
57 "NLS must be non-NULL");
58 return(IDA_ILL_INPUT);
59 }
60
61 /* check for required nonlinear solver functions */
62 if ( NLS->ops->gettype == NULL ||
63 NLS->ops->solve == NULL ||
64 NLS->ops->setsysfn == NULL ) {
65 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
66 "IDASetNonlinearSolver",
67 "NLS does not support required operations");
68 return(IDA_ILL_INPUT);
69 }
70
71 /* check for allowed nonlinear solver types */
72 if (SUNNonlinSolGetType(NLS) != SUNNONLINEARSOLVER_ROOTFIND) {
73 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
74 "IDASetNonlinearSolver",
75 "NLS type must be SUNNONLINEARSOLVER_ROOTFIND");
76 return(IDA_ILL_INPUT);
77 }
78
79 /* free any existing nonlinear solver */
80 if ((IDA_mem->NLS != NULL) && (IDA_mem->ownNLS))
81 retval = SUNNonlinSolFree(IDA_mem->NLS);
82
83 /* set SUNNonlinearSolver pointer */
84 IDA_mem->NLS = NLS;
85
86 /* Set NLS ownership flag. If this function was called to attach the default
87 NLS, IDA will set the flag to SUNTRUE after this function returns. */
88 IDA_mem->ownNLS = SUNFALSE;
89
90 /* set the nonlinear residual function */
91 retval = SUNNonlinSolSetSysFn(IDA_mem->NLS, idaNlsResidual);
92 if (retval != IDA_SUCCESS) {
93 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
94 "IDASetNonlinearSolver",
95 "Setting nonlinear system function failed");
96 return(IDA_ILL_INPUT);
97 }
98
99 /* set convergence test function */
100 retval = SUNNonlinSolSetConvTestFn(IDA_mem->NLS, idaNlsConvTest, ida_mem);
101 if (retval != IDA_SUCCESS) {
102 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
103 "IDASetNonlinearSolver",
104 "Setting convergence test function failed");
105 return(IDA_ILL_INPUT);
106 }
107
108 /* set max allowed nonlinear iterations */
109 retval = SUNNonlinSolSetMaxIters(IDA_mem->NLS, MAXIT);
110 if (retval != IDA_SUCCESS) {
111 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
112 "IDASetNonlinearSolver",
113 "Setting maximum number of nonlinear iterations failed");
114 return(IDA_ILL_INPUT);
115 }
116
117 return(IDA_SUCCESS);
118 }
119
120
121 /* -----------------------------------------------------------------------------
122 * Private functions
123 * ---------------------------------------------------------------------------*/
124
idaNlsInit(IDAMem IDA_mem)125 int idaNlsInit(IDAMem IDA_mem)
126 {
127 int retval;
128
129 /* set the linear solver setup wrapper function */
130 if (IDA_mem->ida_lsetup)
131 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLS, idaNlsLSetup);
132 else
133 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLS, NULL);
134
135 if (retval != IDA_SUCCESS) {
136 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
137 "Setting the linear solver setup function failed");
138 return(IDA_NLS_INIT_FAIL);
139 }
140
141 /* set the linear solver solve wrapper function */
142 if (IDA_mem->ida_lsolve)
143 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLS, idaNlsLSolve);
144 else
145 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLS, NULL);
146
147 if (retval != IDA_SUCCESS) {
148 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
149 "Setting linear solver solve function failed");
150 return(IDA_NLS_INIT_FAIL);
151 }
152
153 /* initialize nonlinear solver */
154 retval = SUNNonlinSolInitialize(IDA_mem->NLS);
155
156 if (retval != IDA_SUCCESS) {
157 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
158 MSG_NLS_INIT_FAIL);
159 return(IDA_NLS_INIT_FAIL);
160 }
161
162 return(IDA_SUCCESS);
163 }
164
165
idaNlsLSetup(booleantype jbad,booleantype * jcur,void * ida_mem)166 static int idaNlsLSetup(booleantype jbad, booleantype* jcur, void* ida_mem)
167 {
168 IDAMem IDA_mem;
169 int retval;
170
171 if (ida_mem == NULL) {
172 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsLSetup", MSG_NO_MEM);
173 return(IDA_MEM_NULL);
174 }
175 IDA_mem = (IDAMem) ida_mem;
176
177 IDA_mem->ida_nsetups++;
178 IDA_mem->ida_forceSetup = SUNFALSE;
179
180 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy, IDA_mem->ida_yp,
181 IDA_mem->ida_savres, IDA_mem->ida_tempv1,
182 IDA_mem->ida_tempv2, IDA_mem->ida_tempv3);
183
184 /* update Jacobian status */
185 *jcur = SUNTRUE;
186
187 /* update convergence test constants */
188 IDA_mem->ida_cjold = IDA_mem->ida_cj;
189 IDA_mem->ida_cjratio = ONE;
190 IDA_mem->ida_ss = TWENTY;
191 IDA_mem->ida_ssS = TWENTY;
192
193 if (retval < 0) return(IDA_LSETUP_FAIL);
194 if (retval > 0) return(IDA_LSETUP_RECVR);
195
196 return(IDA_SUCCESS);
197 }
198
199
idaNlsLSolve(N_Vector delta,void * ida_mem)200 static int idaNlsLSolve(N_Vector delta, void* ida_mem)
201 {
202 IDAMem IDA_mem;
203 int retval;
204
205 if (ida_mem == NULL) {
206 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsLSolve", MSG_NO_MEM);
207 return(IDA_MEM_NULL);
208 }
209 IDA_mem = (IDAMem) ida_mem;
210
211 retval = IDA_mem->ida_lsolve(IDA_mem, delta, IDA_mem->ida_ewt,
212 IDA_mem->ida_yy, IDA_mem->ida_yp,
213 IDA_mem->ida_savres);
214
215 if (retval < 0) return(IDA_LSOLVE_FAIL);
216 if (retval > 0) return(IDA_LSOLVE_RECVR);
217
218 return(IDA_SUCCESS);
219 }
220
221
idaNlsResidual(N_Vector ycor,N_Vector res,void * ida_mem)222 static int idaNlsResidual(N_Vector ycor, N_Vector res, void* ida_mem)
223 {
224 IDAMem IDA_mem;
225 int retval;
226
227 if (ida_mem == NULL) {
228 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsResidual", MSG_NO_MEM);
229 return(IDA_MEM_NULL);
230 }
231 IDA_mem = (IDAMem) ida_mem;
232
233 /* update yy and yp based on the current correction */
234 N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, ycor, IDA_mem->ida_yy);
235 N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, ycor, IDA_mem->ida_yp);
236
237 /* evaluate residual */
238 retval = IDA_mem->ida_res(IDA_mem->ida_tn, IDA_mem->ida_yy, IDA_mem->ida_yp,
239 res, IDA_mem->ida_user_data);
240
241 /* increment the number of residual evaluations */
242 IDA_mem->ida_nre++;
243
244 /* save a copy of the residual vector in savres */
245 N_VScale(ONE, res, IDA_mem->ida_savres);
246
247 if (retval < 0) return(IDA_RES_FAIL);
248 if (retval > 0) return(IDA_RES_RECVR);
249
250 return(IDA_SUCCESS);
251 }
252
253
idaNlsConvTest(SUNNonlinearSolver NLS,N_Vector ycor,N_Vector del,realtype tol,N_Vector ewt,void * ida_mem)254 static int idaNlsConvTest(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector del,
255 realtype tol, N_Vector ewt, void* ida_mem)
256 {
257 IDAMem IDA_mem;
258 int m, retval;
259 realtype delnrm;
260 realtype rate;
261
262 if (ida_mem == NULL) {
263 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsConvTest", MSG_NO_MEM);
264 return(IDA_MEM_NULL);
265 }
266 IDA_mem = (IDAMem) ida_mem;
267
268 /* compute the norm of the correction */
269 delnrm = N_VWrmsNorm(del, ewt);
270
271 /* get the current nonlinear solver iteration count */
272 retval = SUNNonlinSolGetCurIter(NLS, &m);
273 if (retval != IDA_SUCCESS) return(IDA_MEM_NULL);
274
275 /* test for convergence, first directly, then with rate estimate. */
276 if (m == 0){
277 IDA_mem->ida_oldnrm = delnrm;
278 if (delnrm <= PT0001 * IDA_mem->ida_toldel) return(SUN_NLS_SUCCESS);
279 } else {
280 rate = SUNRpowerR( delnrm/IDA_mem->ida_oldnrm, ONE/m );
281 if (rate > RATEMAX) return(SUN_NLS_CONV_RECVR);
282 IDA_mem->ida_ss = rate/(ONE - rate);
283 }
284
285 if (IDA_mem->ida_ss*delnrm <= tol) return(SUN_NLS_SUCCESS);
286
287 /* not yet converged */
288 return(SUN_NLS_CONTINUE);
289 }
290