1 /* -----------------------------------------------------------------------------
2 * Programmer(s): David J. Gardner @ LLNL
3 * -----------------------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2021, 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 IDAGetNonlinearSystemData:
123
124 This routine provides access to the relevant data needed to
125 compute the nonlinear system function.
126 ---------------------------------------------------------------*/
IDAGetNonlinearSystemData(void * ida_mem,realtype * tcur,N_Vector * yypred,N_Vector * yppred,N_Vector * yyn,N_Vector * ypn,N_Vector * res,realtype * cj,void ** user_data)127 int IDAGetNonlinearSystemData(void *ida_mem, realtype *tcur, N_Vector *yypred,
128 N_Vector *yppred, N_Vector *yyn, N_Vector *ypn,
129 N_Vector *res, realtype *cj, void **user_data)
130 {
131 IDAMem IDA_mem;
132
133 if (ida_mem == NULL) {
134 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "IDAGetNonlinearSystemData",
135 MSG_NO_MEM);
136 return(IDA_MEM_NULL);
137 }
138
139 IDA_mem = (IDAMem) ida_mem;
140
141 *tcur = IDA_mem->ida_tn;
142 *yypred = IDA_mem->ida_yypredict;
143 *yppred = IDA_mem->ida_yppredict;
144 *yyn = IDA_mem->ida_yy;
145 *ypn = IDA_mem->ida_yp;
146 *res = IDA_mem->ida_savres;
147 *cj = IDA_mem->ida_cj;
148 *user_data = IDA_mem->ida_user_data;
149
150 return(IDA_SUCCESS);
151 }
152
153
154 /* -----------------------------------------------------------------------------
155 * Private functions
156 * ---------------------------------------------------------------------------*/
157
idaNlsInit(IDAMem IDA_mem)158 int idaNlsInit(IDAMem IDA_mem)
159 {
160 int retval;
161
162 /* set the linear solver setup wrapper function */
163 if (IDA_mem->ida_lsetup)
164 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLS, idaNlsLSetup);
165 else
166 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLS, NULL);
167
168 if (retval != IDA_SUCCESS) {
169 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
170 "Setting the linear solver setup function failed");
171 return(IDA_NLS_INIT_FAIL);
172 }
173
174 /* set the linear solver solve wrapper function */
175 if (IDA_mem->ida_lsolve)
176 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLS, idaNlsLSolve);
177 else
178 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLS, NULL);
179
180 if (retval != IDA_SUCCESS) {
181 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
182 "Setting linear solver solve function failed");
183 return(IDA_NLS_INIT_FAIL);
184 }
185
186 /* initialize nonlinear solver */
187 retval = SUNNonlinSolInitialize(IDA_mem->NLS);
188
189 if (retval != IDA_SUCCESS) {
190 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInit",
191 MSG_NLS_INIT_FAIL);
192 return(IDA_NLS_INIT_FAIL);
193 }
194
195 return(IDA_SUCCESS);
196 }
197
198
idaNlsLSetup(booleantype jbad,booleantype * jcur,void * ida_mem)199 static int idaNlsLSetup(booleantype jbad, booleantype* jcur, void* ida_mem)
200 {
201 IDAMem IDA_mem;
202 int retval;
203
204 if (ida_mem == NULL) {
205 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsLSetup", MSG_NO_MEM);
206 return(IDA_MEM_NULL);
207 }
208 IDA_mem = (IDAMem) ida_mem;
209
210 IDA_mem->ida_nsetups++;
211 IDA_mem->ida_forceSetup = SUNFALSE;
212
213 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy, IDA_mem->ida_yp,
214 IDA_mem->ida_savres, IDA_mem->ida_tempv1,
215 IDA_mem->ida_tempv2, IDA_mem->ida_tempv3);
216
217 /* update Jacobian status */
218 *jcur = SUNTRUE;
219
220 /* update convergence test constants */
221 IDA_mem->ida_cjold = IDA_mem->ida_cj;
222 IDA_mem->ida_cjratio = ONE;
223 IDA_mem->ida_ss = TWENTY;
224 IDA_mem->ida_ssS = TWENTY;
225
226 if (retval < 0) return(IDA_LSETUP_FAIL);
227 if (retval > 0) return(IDA_LSETUP_RECVR);
228
229 return(IDA_SUCCESS);
230 }
231
232
idaNlsLSolve(N_Vector delta,void * ida_mem)233 static int idaNlsLSolve(N_Vector delta, void* ida_mem)
234 {
235 IDAMem IDA_mem;
236 int retval;
237
238 if (ida_mem == NULL) {
239 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsLSolve", MSG_NO_MEM);
240 return(IDA_MEM_NULL);
241 }
242 IDA_mem = (IDAMem) ida_mem;
243
244 retval = IDA_mem->ida_lsolve(IDA_mem, delta, IDA_mem->ida_ewt,
245 IDA_mem->ida_yy, IDA_mem->ida_yp,
246 IDA_mem->ida_savres);
247
248 if (retval < 0) return(IDA_LSOLVE_FAIL);
249 if (retval > 0) return(IDA_LSOLVE_RECVR);
250
251 return(IDA_SUCCESS);
252 }
253
254
idaNlsResidual(N_Vector ycor,N_Vector res,void * ida_mem)255 static int idaNlsResidual(N_Vector ycor, N_Vector res, void* ida_mem)
256 {
257 IDAMem IDA_mem;
258 int retval;
259
260 if (ida_mem == NULL) {
261 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsResidual", MSG_NO_MEM);
262 return(IDA_MEM_NULL);
263 }
264 IDA_mem = (IDAMem) ida_mem;
265
266 /* update yy and yp based on the current correction */
267 N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, ycor, IDA_mem->ida_yy);
268 N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, ycor, IDA_mem->ida_yp);
269
270 /* evaluate residual */
271 retval = IDA_mem->ida_res(IDA_mem->ida_tn, IDA_mem->ida_yy, IDA_mem->ida_yp,
272 res, IDA_mem->ida_user_data);
273
274 /* increment the number of residual evaluations */
275 IDA_mem->ida_nre++;
276
277 /* save a copy of the residual vector in savres */
278 N_VScale(ONE, res, IDA_mem->ida_savres);
279
280 if (retval < 0) return(IDA_RES_FAIL);
281 if (retval > 0) return(IDA_RES_RECVR);
282
283 return(IDA_SUCCESS);
284 }
285
286
idaNlsConvTest(SUNNonlinearSolver NLS,N_Vector ycor,N_Vector del,realtype tol,N_Vector ewt,void * ida_mem)287 static int idaNlsConvTest(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector del,
288 realtype tol, N_Vector ewt, void* ida_mem)
289 {
290 IDAMem IDA_mem;
291 int m, retval;
292 realtype delnrm;
293 realtype rate;
294
295 if (ida_mem == NULL) {
296 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsConvTest", MSG_NO_MEM);
297 return(IDA_MEM_NULL);
298 }
299 IDA_mem = (IDAMem) ida_mem;
300
301 /* compute the norm of the correction */
302 delnrm = N_VWrmsNorm(del, ewt);
303
304 /* get the current nonlinear solver iteration count */
305 retval = SUNNonlinSolGetCurIter(NLS, &m);
306 if (retval != IDA_SUCCESS) return(IDA_MEM_NULL);
307
308 /* test for convergence, first directly, then with rate estimate. */
309 if (m == 0){
310 IDA_mem->ida_oldnrm = delnrm;
311 if (delnrm <= PT0001 * IDA_mem->ida_toldel) return(SUN_NLS_SUCCESS);
312 } else {
313 rate = SUNRpowerR( delnrm/IDA_mem->ida_oldnrm, ONE/m );
314 if (rate > RATEMAX) return(SUN_NLS_CONV_RECVR);
315 IDA_mem->ida_ss = rate/(ONE - rate);
316 }
317
318 if (IDA_mem->ida_ss*delnrm <= tol) return(SUN_NLS_SUCCESS);
319
320 /* not yet converged */
321 return(SUN_NLS_CONTINUE);
322 }
323