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 #include "sundials/sundials_nvector_senswrapper.h"
20
21 /* constant macros */
22 #define PT0001 RCONST(0.0001) /* real 0.0001 */
23 #define ONE RCONST(1.0) /* real 1.0 */
24 #define TWENTY RCONST(20.0) /* real 20.0 */
25
26 /* nonlinear solver parameters */
27 #define MAXIT 4 /* default max number of nonlinear iterations */
28 #define RATEMAX RCONST(0.9) /* max convergence rate used in divergence check */
29
30 /* private functions passed to nonlinear solver */
31 static int idaNlsResidualSensSim(N_Vector ycor, N_Vector res, void* ida_mem);
32 static int idaNlsLSetupSensSim(booleantype jbad, booleantype* jcur,
33 void* ida_mem);
34 static int idaNlsLSolveSensSim(N_Vector delta, void* ida_mem);
35 static int idaNlsConvTestSensSim(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector del,
36 realtype tol, N_Vector ewt, void* ida_mem);
37
38 /* -----------------------------------------------------------------------------
39 * Exported functions
40 * ---------------------------------------------------------------------------*/
41
IDASetNonlinearSolverSensSim(void * ida_mem,SUNNonlinearSolver NLS)42 int IDASetNonlinearSolverSensSim(void *ida_mem, SUNNonlinearSolver NLS)
43 {
44 IDAMem IDA_mem;
45 int retval, is;
46
47 /* return immediately if IDA memory is NULL */
48 if (ida_mem == NULL) {
49 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS",
50 "IDASetNonlinearSolverSensSim", MSG_NO_MEM);
51 return(IDA_MEM_NULL);
52 }
53 IDA_mem = (IDAMem) ida_mem;
54
55 /* return immediately if NLS memory is NULL */
56 if (NLS == NULL) {
57 IDAProcessError(NULL, IDA_ILL_INPUT, "IDAS",
58 "IDASetNonlinearSolverSensSim",
59 "NLS must be non-NULL");
60 return(IDA_ILL_INPUT);
61 }
62
63 /* check for required nonlinear solver functions */
64 if ( NLS->ops->gettype == NULL ||
65 NLS->ops->solve == NULL ||
66 NLS->ops->setsysfn == NULL ) {
67 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
68 "IDASetNonlinearSolverSensSim",
69 "NLS does not support required operations");
70 return(IDA_ILL_INPUT);
71 }
72
73 /* check for allowed nonlinear solver types */
74 if (SUNNonlinSolGetType(NLS) != SUNNONLINEARSOLVER_ROOTFIND) {
75 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
76 "IDASetNonlinearSolverSensSim",
77 "NLS type must be SUNNONLINEARSOLVER_ROOTFIND");
78 return(IDA_ILL_INPUT);
79 }
80
81 /* check that sensitivities were initialized */
82 if (!(IDA_mem->ida_sensi)) {
83 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
84 "IDASetNonlinearSolverSensSim",
85 MSG_NO_SENSI);
86 return(IDA_ILL_INPUT);
87 }
88
89 /* check that the simultaneous corrector was selected */
90 if (IDA_mem->ida_ism != IDA_SIMULTANEOUS) {
91 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
92 "IDASetNonlinearSolverSensSim",
93 "Sensitivity solution method is not IDA_SIMULTANEOUS");
94 return(IDA_ILL_INPUT);
95 }
96
97 /* free any existing nonlinear solver */
98 if ((IDA_mem->NLSsim != NULL) && (IDA_mem->ownNLSsim))
99 retval = SUNNonlinSolFree(IDA_mem->NLSsim);
100
101 /* set SUNNonlinearSolver pointer */
102 IDA_mem->NLSsim = NLS;
103
104 /* Set NLS ownership flag. If this function was called to attach the default
105 NLS, IDA will set the flag to SUNTRUE after this function returns. */
106 IDA_mem->ownNLSsim = SUNFALSE;
107
108 /* set the nonlinear residual function */
109 retval = SUNNonlinSolSetSysFn(IDA_mem->NLSsim, idaNlsResidualSensSim);
110 if (retval != IDA_SUCCESS) {
111 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
112 "IDASetNonlinearSolverSensSim",
113 "Setting nonlinear system function failed");
114 return(IDA_ILL_INPUT);
115 }
116
117 /* set convergence test function */
118 retval = SUNNonlinSolSetConvTestFn(IDA_mem->NLSsim, idaNlsConvTestSensSim,
119 ida_mem);
120 if (retval != IDA_SUCCESS) {
121 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
122 "IDASetNonlinearSolverSensSim",
123 "Setting convergence test function failed");
124 return(IDA_ILL_INPUT);
125 }
126
127 /* set max allowed nonlinear iterations */
128 retval = SUNNonlinSolSetMaxIters(IDA_mem->NLSsim, MAXIT);
129 if (retval != IDA_SUCCESS) {
130 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS",
131 "IDASetNonlinearSolverSensSim",
132 "Setting maximum number of nonlinear iterations failed");
133 return(IDA_ILL_INPUT);
134 }
135
136 /* create vector wrappers if necessary */
137 if (IDA_mem->simMallocDone == SUNFALSE) {
138
139 IDA_mem->ypredictSim = N_VNewEmpty_SensWrapper(IDA_mem->ida_Ns+1);
140 if (IDA_mem->ypredictSim == NULL) {
141 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAS",
142 "IDASetNonlinearSolverSensSim", MSG_MEM_FAIL);
143 return(IDA_MEM_FAIL);
144 }
145
146 IDA_mem->ycorSim = N_VNewEmpty_SensWrapper(IDA_mem->ida_Ns+1);
147 if (IDA_mem->ycorSim == NULL) {
148 N_VDestroy(IDA_mem->ypredictSim);
149 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAS",
150 "IDASetNonlinearSolverSensSim", MSG_MEM_FAIL);
151 return(IDA_MEM_FAIL);
152 }
153
154 IDA_mem->ewtSim = N_VNewEmpty_SensWrapper(IDA_mem->ida_Ns+1);
155 if (IDA_mem->ewtSim == NULL) {
156 N_VDestroy(IDA_mem->ypredictSim);
157 N_VDestroy(IDA_mem->ycorSim);
158 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAS",
159 "IDASetNonlinearSolverSensSim", MSG_MEM_FAIL);
160 return(IDA_MEM_FAIL);
161 }
162
163 IDA_mem->simMallocDone = SUNTRUE;
164 }
165
166 /* attach vectors to vector wrappers */
167 NV_VEC_SW(IDA_mem->ypredictSim, 0) = IDA_mem->ida_yypredict;
168 NV_VEC_SW(IDA_mem->ycorSim, 0) = IDA_mem->ida_ee;
169 NV_VEC_SW(IDA_mem->ewtSim, 0) = IDA_mem->ida_ewt;
170
171 for (is=0; is < IDA_mem->ida_Ns; is++) {
172 NV_VEC_SW(IDA_mem->ypredictSim, is+1) = IDA_mem->ida_yySpredict[is];
173 NV_VEC_SW(IDA_mem->ycorSim, is+1) = IDA_mem->ida_eeS[is];
174 NV_VEC_SW(IDA_mem->ewtSim, is+1) = IDA_mem->ida_ewtS[is];
175 }
176
177 return(IDA_SUCCESS);
178 }
179
180
181 /* -----------------------------------------------------------------------------
182 * Private functions
183 * ---------------------------------------------------------------------------*/
184
idaNlsInitSensSim(IDAMem IDA_mem)185 int idaNlsInitSensSim(IDAMem IDA_mem)
186 {
187 int retval;
188
189 /* set the linear solver setup wrapper function */
190 if (IDA_mem->ida_lsetup)
191 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLSsim, idaNlsLSetupSensSim);
192 else
193 retval = SUNNonlinSolSetLSetupFn(IDA_mem->NLSsim, NULL);
194
195 if (retval != IDA_SUCCESS) {
196 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInitSnesSim",
197 "Setting the linear solver setup function failed");
198 return(IDA_NLS_INIT_FAIL);
199 }
200
201 /* set the linear solver solve wrapper function */
202 if (IDA_mem->ida_lsolve)
203 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLSsim, idaNlsLSolveSensSim);
204 else
205 retval = SUNNonlinSolSetLSolveFn(IDA_mem->NLSsim, NULL);
206
207 if (retval != IDA_SUCCESS) {
208 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInitSnesSim",
209 "Setting linear solver solve function failed");
210 return(IDA_NLS_INIT_FAIL);
211 }
212
213 /* initialize nonlinear solver */
214 retval = SUNNonlinSolInitialize(IDA_mem->NLSsim);
215
216 if (retval != IDA_SUCCESS) {
217 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "idaNlsInitSnesSim",
218 MSG_NLS_INIT_FAIL);
219 return(IDA_NLS_INIT_FAIL);
220 }
221
222 return(IDA_SUCCESS);
223 }
224
225
idaNlsLSetupSensSim(booleantype jbad,booleantype * jcur,void * ida_mem)226 static int idaNlsLSetupSensSim(booleantype jbad, booleantype* jcur,
227 void* ida_mem)
228 {
229 IDAMem IDA_mem;
230 int retval;
231
232 if (ida_mem == NULL) {
233 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS",
234 "idaNlsLSetupSensSim", MSG_NO_MEM);
235 return(IDA_MEM_NULL);
236 }
237 IDA_mem = (IDAMem) ida_mem;
238
239 IDA_mem->ida_nsetups++;
240 IDA_mem->ida_forceSetup = SUNFALSE;
241
242 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy, IDA_mem->ida_yp,
243 IDA_mem->ida_savres, IDA_mem->ida_tempv1,
244 IDA_mem->ida_tempv2, IDA_mem->ida_tempv3);
245
246 /* update Jacobian status */
247 *jcur = SUNTRUE;
248
249 /* update convergence test constants */
250 IDA_mem->ida_cjold = IDA_mem->ida_cj;
251 IDA_mem->ida_cjratio = ONE;
252 IDA_mem->ida_ss = TWENTY;
253 IDA_mem->ida_ssS = TWENTY;
254
255 if (retval < 0) return(IDA_LSETUP_FAIL);
256 if (retval > 0) return(IDA_LSETUP_RECVR);
257
258 return(IDA_SUCCESS);
259 }
260
261
idaNlsLSolveSensSim(N_Vector deltaSim,void * ida_mem)262 static int idaNlsLSolveSensSim(N_Vector deltaSim, void* ida_mem)
263 {
264 IDAMem IDA_mem;
265 int retval, is;
266 N_Vector delta;
267 N_Vector *deltaS;
268
269 if (ida_mem == NULL) {
270 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS",
271 "idaNlsLSolveSensSim", MSG_NO_MEM);
272 return(IDA_MEM_NULL);
273 }
274 IDA_mem = (IDAMem) ida_mem;
275
276 /* extract state update vector from the vector wrapper */
277 delta = NV_VEC_SW(deltaSim,0);
278
279 /* solve the state linear system */
280 retval = IDA_mem->ida_lsolve(IDA_mem, delta, IDA_mem->ida_ewt,
281 IDA_mem->ida_yy, IDA_mem->ida_yp,
282 IDA_mem->ida_savres);
283
284 if (retval < 0) return(IDA_LSOLVE_FAIL);
285 if (retval > 0) return(IDA_LSOLVE_RECVR);
286
287 /* extract sensitivity deltas from the vector wrapper */
288 deltaS = NV_VECS_SW(deltaSim)+1;
289
290 /* solve the sensitivity linear systems */
291 for(is=0; is<IDA_mem->ida_Ns; is++) {
292 retval = IDA_mem->ida_lsolve(IDA_mem, deltaS[is], IDA_mem->ida_ewtS[is],
293 IDA_mem->ida_yy, IDA_mem->ida_yp,
294 IDA_mem->ida_savres);
295
296 if (retval < 0) return(IDA_LSOLVE_FAIL);
297 if (retval > 0) return(IDA_LSOLVE_RECVR);
298 }
299
300 return(IDA_SUCCESS);
301 }
302
303
idaNlsResidualSensSim(N_Vector ycorSim,N_Vector resSim,void * ida_mem)304 static int idaNlsResidualSensSim(N_Vector ycorSim, N_Vector resSim, void* ida_mem)
305 {
306 IDAMem IDA_mem;
307 int retval;
308 N_Vector ycor, res;
309 N_Vector *ycorS, *resS;
310
311 if (ida_mem == NULL) {
312 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS",
313 "idaNlsResidualSensSim", MSG_NO_MEM);
314 return(IDA_MEM_NULL);
315 }
316 IDA_mem = (IDAMem) ida_mem;
317
318 /* extract state and residual vectors from the vector wrapper */
319 ycor = NV_VEC_SW(ycorSim,0);
320 res = NV_VEC_SW(resSim,0);
321
322 /* update yy and yp based on the current correction */
323 N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, ycor, IDA_mem->ida_yy);
324 N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, ycor, IDA_mem->ida_yp);
325
326 /* evaluate residual */
327 retval = IDA_mem->ida_res(IDA_mem->ida_tn, IDA_mem->ida_yy, IDA_mem->ida_yp,
328 res, IDA_mem->ida_user_data);
329
330 /* increment the number of residual evaluations */
331 IDA_mem->ida_nre++;
332
333 /* save a copy of the residual vector in savres */
334 N_VScale(ONE, res, IDA_mem->ida_savres);
335
336 if (retval < 0) return(IDA_RES_FAIL);
337 if (retval > 0) return(IDA_RES_RECVR);
338
339 /* extract sensitivity and residual vectors from the vector wrapper */
340 ycorS = NV_VECS_SW(ycorSim)+1;
341 resS = NV_VECS_SW(resSim)+1;
342
343 /* update yS and ypS based on the current correction */
344 N_VLinearSumVectorArray(IDA_mem->ida_Ns,
345 ONE, IDA_mem->ida_yySpredict,
346 ONE, ycorS, IDA_mem->ida_yyS);
347 N_VLinearSumVectorArray(IDA_mem->ida_Ns,
348 ONE, IDA_mem->ida_ypSpredict,
349 IDA_mem->ida_cj, ycorS, IDA_mem->ida_ypS);
350
351 /* evaluate sens residual */
352 retval = IDA_mem->ida_resS(IDA_mem->ida_Ns, IDA_mem->ida_tn,
353 IDA_mem->ida_yy, IDA_mem->ida_yp, res,
354 IDA_mem->ida_yyS, IDA_mem->ida_ypS, resS,
355 IDA_mem->ida_user_dataS, IDA_mem->ida_tmpS1,
356 IDA_mem->ida_tmpS2, IDA_mem->ida_tmpS3);
357
358 /* increment the number of sens residual evaluations */
359 IDA_mem->ida_nrSe++;
360
361 if (retval < 0) return(IDA_SRES_FAIL);
362 if (retval > 0) return(IDA_SRES_RECVR);
363
364 return(IDA_SUCCESS);
365 }
366
367
idaNlsConvTestSensSim(SUNNonlinearSolver NLS,N_Vector ycor,N_Vector del,realtype tol,N_Vector ewt,void * ida_mem)368 static int idaNlsConvTestSensSim(SUNNonlinearSolver NLS, N_Vector ycor, N_Vector del,
369 realtype tol, N_Vector ewt, void* ida_mem)
370 {
371 IDAMem IDA_mem;
372 int m, retval;
373 realtype delnrm;
374 realtype rate;
375
376 if (ida_mem == NULL) {
377 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "idaNlsConvTestSensSim", MSG_NO_MEM);
378 return(IDA_MEM_NULL);
379 }
380 IDA_mem = (IDAMem) ida_mem;
381
382 /* compute the norm of the correction */
383 delnrm = N_VWrmsNorm(del, ewt);
384
385 /* get the current nonlinear solver iteration count */
386 retval = SUNNonlinSolGetCurIter(NLS, &m);
387 if (retval != IDA_SUCCESS) return(IDA_MEM_NULL);
388
389 /* test for convergence, first directly, then with rate estimate. */
390 if (m == 0){
391 IDA_mem->ida_oldnrm = delnrm;
392 if (delnrm <= PT0001 * IDA_mem->ida_toldel) return(SUN_NLS_SUCCESS);
393 } else {
394 rate = SUNRpowerR( delnrm/IDA_mem->ida_oldnrm, ONE/m );
395 if (rate > RATEMAX) return(SUN_NLS_CONV_RECVR);
396 IDA_mem->ida_ss = rate/(ONE - rate);
397 }
398
399 if (IDA_mem->ida_ss*delnrm <= tol) return(SUN_NLS_SUCCESS);
400
401 /* not yet converged */
402 return(SUN_NLS_CONTINUE);
403 }
404