1 /*
2 * -----------------------------------------------------------------
3 * $Revision$
4 * $Date$
5 * -----------------------------------------------------------------
6 * Programmers: Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * SUNDIALS Copyright Start
9 * Copyright (c) 2002-2019, Lawrence Livermore National Security
10 * and Southern Methodist University.
11 * All rights reserved.
12 *
13 * See the top-level LICENSE and NOTICE files for details.
14 *
15 * SPDX-License-Identifier: BSD-3-Clause
16 * SUNDIALS Copyright End
17 * -----------------------------------------------------------------
18 * This is the implementation file for the IC calculation for IDAS.
19 * It is independent of the linear solver in use.
20 * -----------------------------------------------------------------
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "idas_impl.h"
27 #include <sundials/sundials_math.h>
28
29 /*
30 * =================================================================
31 * IDA Constants
32 * =================================================================
33 */
34
35 /* Private Constants */
36
37 #define ZERO RCONST(0.0) /* real 0.0 */
38 #define HALF RCONST(0.5) /* real 0.5 */
39 #define ONE RCONST(1.0) /* real 1.0 */
40 #define TWO RCONST(2.0) /* real 2.0 */
41 #define PT99 RCONST(0.99) /* real 0.99 */
42 #define PT1 RCONST(0.1) /* real 0.1 */
43 #define PT001 RCONST(0.001) /* real 0.001 */
44
45 /* IDACalcIC control constants */
46
47 #define ICRATEMAX RCONST(0.9) /* max. Newton conv. rate */
48 #define ALPHALS RCONST(0.0001) /* alpha in linesearch conv. test */
49
50 /* Return values for lower level routines used by IDACalcIC */
51
52 #define IC_FAIL_RECOV 1
53 #define IC_CONSTR_FAILED 2
54 #define IC_LINESRCH_FAILED 3
55 #define IC_CONV_FAIL 4
56 #define IC_SLOW_CONVRG 5
57
58 /*
59 * =================================================================
60 * Private Helper Functions Prototypes
61 * =================================================================
62 */
63
64 extern int IDAInitialSetup(IDAMem IDA_mem);
65 extern realtype IDAWrmsNorm(IDAMem IDA_mem, N_Vector x,
66 N_Vector w, booleantype mask);
67 extern realtype IDASensWrmsNorm(IDAMem IDA_mem, N_Vector *xS,
68 N_Vector *wS, booleantype mask);
69 extern realtype IDASensWrmsNormUpdate(IDAMem IDA_mem, realtype old_nrm,
70 N_Vector *xS, N_Vector *wS,
71 booleantype mask);
72
73 extern int IDASensEwtSet(IDAMem IDA_mem, N_Vector *yScur, N_Vector *weightS);
74
75 static int IDANlsIC(IDAMem IDA_mem);
76
77 static int IDANewtonIC(IDAMem IDA_mem);
78 static int IDALineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm);
79 static int IDAfnorm(IDAMem IDA_mem, realtype *fnorm);
80 static int IDANewyyp(IDAMem IDA_mem, realtype lambda);
81 static int IDANewy(IDAMem IDA_mem);
82
83 static int IDASensNewtonIC(IDAMem IDA_mem);
84 static int IDASensLineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm);
85 static int IDASensNewyyp(IDAMem IDA_mem, realtype lambda);
86 static int IDASensfnorm(IDAMem IDA_mem, realtype *fnorm);
87 static int IDASensNlsIC(IDAMem IDA_mem);
88
89 static int IDAICFailFlag(IDAMem IDA_mem, int retval);
90
91
92 /*
93 * =================================================================
94 * EXPORTED FUNCTIONS IMPLEMENTATION
95 * =================================================================
96 */
97
98 /*
99 * -----------------------------------------------------------------
100 * IDACalcIC
101 * -----------------------------------------------------------------
102 * IDACalcIC computes consistent initial conditions, given the
103 * user's initial guess for unknown components of yy0 and/or yp0.
104 *
105 * The return value is IDA_SUCCESS = 0 if no error occurred.
106 *
107 * The error return values (fully described in ida.h) are:
108 * IDA_MEM_NULL ida_mem is NULL
109 * IDA_NO_MALLOC ida_mem was not allocated
110 * IDA_ILL_INPUT bad value for icopt, tout1, or id
111 * IDA_LINIT_FAIL the linear solver linit routine failed
112 * IDA_BAD_EWT zero value of some component of ewt
113 * IDA_RES_FAIL res had a non-recoverable error
114 * IDA_FIRST_RES_FAIL res failed recoverably on the first call
115 * IDA_LSETUP_FAIL lsetup had a non-recoverable error
116 * IDA_LSOLVE_FAIL lsolve had a non-recoverable error
117 * IDA_NO_RECOVERY res, lsetup, or lsolve had a recoverable
118 * error, but IDACalcIC could not recover
119 * IDA_CONSTR_FAIL the inequality constraints could not be met
120 * IDA_LINESEARCH_FAIL if the linesearch failed (either on steptol test
121 * or on the maxbacks test)
122 * IDA_CONV_FAIL the Newton iterations failed to converge
123 * -----------------------------------------------------------------
124 */
125
IDACalcIC(void * ida_mem,int icopt,realtype tout1)126 int IDACalcIC(void *ida_mem, int icopt, realtype tout1)
127 {
128 int ewtsetOK;
129 int ier, nwt, nh, mxnh, icret, retval=0;
130 int is;
131 realtype tdist, troundoff, minid, hic, ypnorm;
132 IDAMem IDA_mem;
133 booleantype sensi_stg, sensi_sim;
134
135 /* Check if IDA memory exists */
136
137 if(ida_mem == NULL) {
138 IDAProcessError(NULL, IDA_MEM_NULL, "IDAS", "IDACalcIC", MSG_NO_MEM);
139 return(IDA_MEM_NULL);
140 }
141 IDA_mem = (IDAMem) ida_mem;
142
143 /* Check if problem was malloc'ed */
144
145 if(IDA_mem->ida_MallocDone == SUNFALSE) {
146 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDAS", "IDACalcIC", MSG_NO_MALLOC);
147 return(IDA_NO_MALLOC);
148 }
149
150 /* Check inputs to IDA for correctness and consistency */
151
152 ier = IDAInitialSetup(IDA_mem);
153 if(ier != IDA_SUCCESS) return(IDA_ILL_INPUT);
154 IDA_mem->ida_SetupDone = SUNTRUE;
155
156 /* Check legality of input arguments, and set IDA memory copies. */
157
158 if(icopt != IDA_YA_YDP_INIT && icopt != IDA_Y_INIT) {
159 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "IDACalcIC", MSG_IC_BAD_ICOPT);
160 return(IDA_ILL_INPUT);
161 }
162 IDA_mem->ida_icopt = icopt;
163
164 if(icopt == IDA_YA_YDP_INIT && (IDA_mem->ida_id == NULL)) {
165 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "IDACalcIC", MSG_IC_MISSING_ID);
166 return(IDA_ILL_INPUT);
167 }
168
169 tdist = SUNRabs(tout1 - IDA_mem->ida_tn);
170 troundoff = TWO * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(tout1));
171 if(tdist < troundoff) {
172 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "IDACalcIC", MSG_IC_TOO_CLOSE);
173 return(IDA_ILL_INPUT);
174 }
175
176 /* Are we computing sensitivities? */
177 sensi_stg = (IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_STAGGERED));
178 sensi_sim = (IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS));
179
180 /* Allocate space and initialize temporary vectors */
181
182 IDA_mem->ida_yy0 = N_VClone(IDA_mem->ida_ee);
183 IDA_mem->ida_yp0 = N_VClone(IDA_mem->ida_ee);
184 IDA_mem->ida_t0 = IDA_mem->ida_tn;
185 N_VScale(ONE, IDA_mem->ida_phi[0], IDA_mem->ida_yy0);
186 N_VScale(ONE, IDA_mem->ida_phi[1], IDA_mem->ida_yp0);
187
188 if (IDA_mem->ida_sensi) {
189
190 /* Allocate temporary space required for sensitivity IC: yyS0 and ypS0. */
191 IDA_mem->ida_yyS0 = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_ee);
192 IDA_mem->ida_ypS0 = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_ee);
193
194 /* Initialize sensitivity vector. */
195 for (is=0; is<IDA_mem->ida_Ns; is++) {
196 N_VScale(ONE, IDA_mem->ida_phiS[0][is], IDA_mem->ida_yyS0[is]);
197 N_VScale(ONE, IDA_mem->ida_phiS[1][is], IDA_mem->ida_ypS0[is]);
198 }
199
200 /* Initialize work space vectors needed for sensitivities. */
201 IDA_mem->ida_savresS = IDA_mem->ida_phiS[2];
202 IDA_mem->ida_delnewS = IDA_mem->ida_phiS[3];
203 IDA_mem->ida_yyS0new = IDA_mem->ida_phiS[4];
204 IDA_mem->ida_ypS0new = IDA_mem->ida_eeS;
205 }
206
207 /* For use in the IDA_YA_YP_INIT case, set sysindex and tscale. */
208
209 IDA_mem->ida_sysindex = 1;
210 IDA_mem->ida_tscale = tdist;
211 if(icopt == IDA_YA_YDP_INIT) {
212 minid = N_VMin(IDA_mem->ida_id);
213 if(minid < ZERO) {
214 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAS", "IDACalcIC", MSG_IC_BAD_ID);
215 return(IDA_ILL_INPUT);
216 }
217 if(minid > HALF) IDA_mem->ida_sysindex = 0;
218 }
219
220 /* Set the test constant in the Newton convergence test */
221
222 IDA_mem->ida_epsNewt = IDA_mem->ida_epiccon;
223
224 /* Initializations:
225 cjratio = 1 (for use in direct linear solvers);
226 set nbacktr = 0; */
227
228 IDA_mem->ida_cjratio = ONE;
229 IDA_mem->ida_nbacktr = 0;
230
231 /* Set hic, hh, cj, and mxnh. */
232
233 hic = PT001*tdist;
234 ypnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_yp0, IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
235
236 if (sensi_sim)
237 ypnorm = IDASensWrmsNormUpdate(IDA_mem, ypnorm, IDA_mem->ida_ypS0, IDA_mem->ida_ewtS, SUNFALSE);
238
239 if(ypnorm > HALF/hic) hic = HALF/ypnorm;
240 if(tout1 < IDA_mem->ida_tn) hic = -hic;
241 IDA_mem->ida_hh = hic;
242 if(icopt == IDA_YA_YDP_INIT) {
243 IDA_mem->ida_cj = ONE/hic;
244 mxnh = IDA_mem->ida_maxnh;
245 }
246 else {
247 IDA_mem->ida_cj = ZERO;
248 mxnh = 1;
249 }
250
251 /* Loop over nwt = number of evaluations of ewt vector. */
252
253 for(nwt = 1; nwt <= 2; nwt++) {
254
255 /* Loop over nh = number of h values. */
256 for(nh = 1; nh <= mxnh; nh++) {
257
258 /* Call the IC nonlinear solver function. */
259 retval = IDANlsIC(IDA_mem);
260
261 /* Cut h and loop on recoverable IDA_YA_YDP_INIT failure; else break. */
262 if(retval == IDA_SUCCESS) break;
263 IDA_mem->ida_ncfn++;
264 if(retval < 0) break;
265 if(nh == mxnh) break;
266
267 /* If looping to try again, reset yy0 and yp0 if not converging. */
268 if(retval != IC_SLOW_CONVRG) {
269 N_VScale(ONE, IDA_mem->ida_phi[0], IDA_mem->ida_yy0);
270 N_VScale(ONE, IDA_mem->ida_phi[1], IDA_mem->ida_yp0);
271 if (sensi_sim) {
272
273 /* Reset yyS0 and ypS0. */
274 /* Copy phiS[0] and phiS[1] into yyS0 and ypS0. */
275 for (is=0; is<IDA_mem->ida_Ns; is++) {
276 N_VScale(ONE, IDA_mem->ida_phiS[0][is], IDA_mem->ida_yyS0[is]);
277 N_VScale(ONE, IDA_mem->ida_phiS[1][is], IDA_mem->ida_ypS0[is]);
278 }
279 }
280 }
281 hic *= PT1;
282 IDA_mem->ida_cj = ONE/hic;
283 IDA_mem->ida_hh = hic;
284 } /* End of nh loop */
285
286 /* Break on failure */
287 if(retval != IDA_SUCCESS) break;
288
289 /* Reset ewt, save yy0, yp0 in phi, and loop. */
290 ewtsetOK = IDA_mem->ida_efun(IDA_mem->ida_yy0, IDA_mem->ida_ewt, IDA_mem->ida_edata);
291 if(ewtsetOK != 0) {
292 retval = IDA_BAD_EWT;
293 break;
294 }
295 N_VScale(ONE, IDA_mem->ida_yy0, IDA_mem->ida_phi[0]);
296 N_VScale(ONE, IDA_mem->ida_yp0, IDA_mem->ida_phi[1]);
297
298 if (sensi_sim) {
299
300 /* Reevaluate ewtS. */
301 ewtsetOK = IDASensEwtSet(IDA_mem, IDA_mem->ida_yyS0, IDA_mem->ida_ewtS);
302 if(ewtsetOK != 0) {
303 retval = IDA_BAD_EWT;
304 break;
305 }
306
307 /* Save yyS0 and ypS0. */
308 for (is=0; is<IDA_mem->ida_Ns; is++) {
309 N_VScale(ONE, IDA_mem->ida_yyS0[is], IDA_mem->ida_phiS[0][is]);
310 N_VScale(ONE, IDA_mem->ida_ypS0[is], IDA_mem->ida_phiS[1][is]);
311 }
312 }
313
314 } /* End of nwt loop */
315
316 /* Load the optional outputs. */
317
318 if(icopt == IDA_YA_YDP_INIT) IDA_mem->ida_hused = hic;
319
320 /* On any failure, free memory, print error message and return */
321
322 if(retval != IDA_SUCCESS) {
323 N_VDestroy(IDA_mem->ida_yy0);
324 N_VDestroy(IDA_mem->ida_yp0);
325
326 if(IDA_mem->ida_sensi) {
327 N_VDestroyVectorArray(IDA_mem->ida_yyS0, IDA_mem->ida_Ns);
328 N_VDestroyVectorArray(IDA_mem->ida_ypS0, IDA_mem->ida_Ns);
329 }
330
331 icret = IDAICFailFlag(IDA_mem, retval);
332 return(icret);
333 }
334
335 /* Unless using the STAGGERED approach for sensitivities, return now */
336
337 if (!sensi_stg) {
338
339 N_VDestroy(IDA_mem->ida_yy0);
340 N_VDestroy(IDA_mem->ida_yp0);
341
342 if(IDA_mem->ida_sensi) {
343 N_VDestroyVectorArray(IDA_mem->ida_yyS0, IDA_mem->ida_Ns);
344 N_VDestroyVectorArray(IDA_mem->ida_ypS0, IDA_mem->ida_Ns);
345 }
346
347 return(IDA_SUCCESS);
348 }
349
350 /* Find consistent I.C. for sensitivities using a staggered approach */
351
352
353 /* Evaluate res at converged y, needed for future evaluations of sens. RHS
354 If res() fails recoverably, treat it as a convergence failure and
355 attempt the step again */
356
357 retval = IDA_mem->ida_res(IDA_mem->ida_t0, IDA_mem->ida_yy0,
358 IDA_mem->ida_yp0, IDA_mem->ida_delta,
359 IDA_mem->ida_user_data);
360 IDA_mem->ida_nre++;
361 if(retval < 0)
362 /* res function failed unrecoverably. */
363 return(IDA_RES_FAIL);
364
365 if(retval > 0)
366 /* res function failed recoverably but no recovery possible. */
367 return(IDA_FIRST_RES_FAIL);
368
369 /* Loop over nwt = number of evaluations of ewt vector. */
370 for(nwt = 1; nwt <= 2; nwt++) {
371
372 /* Loop over nh = number of h values. */
373 for(nh = 1; nh <= mxnh; nh++) {
374
375 retval = IDASensNlsIC(IDA_mem);
376 if(retval == IDA_SUCCESS) break;
377
378 /* Increment the number of the sensitivity related corrector convergence failures. */
379 IDA_mem->ida_ncfnS++;
380
381 if(retval < 0) break;
382 if(nh == mxnh) break;
383
384 /* If looping to try again, reset yyS0 and ypS0 if not converging. */
385 if(retval != IC_SLOW_CONVRG) {
386 for (is=0; is<IDA_mem->ida_Ns; is++) {
387 N_VScale(ONE, IDA_mem->ida_phiS[0][is], IDA_mem->ida_yyS0[is]);
388 N_VScale(ONE, IDA_mem->ida_phiS[1][is], IDA_mem->ida_ypS0[is]);
389 }
390 }
391 hic *= PT1;
392 IDA_mem->ida_cj = ONE/hic;
393 IDA_mem->ida_hh = hic;
394
395 } /* End of nh loop */
396
397 /* Break on failure */
398 if(retval != IDA_SUCCESS) break;
399
400 /* Since it was successful, reevaluate ewtS with the new values of yyS0, save
401 yyS0 and ypS0 in phiS[0] and phiS[1] and loop one more time to check and
402 maybe correct the new sensitivities IC with respect to the new weights. */
403
404 /* Reevaluate ewtS. */
405 ewtsetOK = IDASensEwtSet(IDA_mem, IDA_mem->ida_yyS0, IDA_mem->ida_ewtS);
406 if(ewtsetOK != 0) {
407 retval = IDA_BAD_EWT;
408 break;
409 }
410
411 /* Save yyS0 and ypS0. */
412 for (is=0; is<IDA_mem->ida_Ns; is++) {
413 N_VScale(ONE, IDA_mem->ida_yyS0[is], IDA_mem->ida_phiS[0][is]);
414 N_VScale(ONE, IDA_mem->ida_ypS0[is], IDA_mem->ida_phiS[1][is]);
415 }
416
417 } /* End of nwt loop */
418
419
420 /* Load the optional outputs. */
421 if(icopt == IDA_YA_YDP_INIT) IDA_mem->ida_hused = hic;
422
423 /* Free temporary space */
424 N_VDestroy(IDA_mem->ida_yy0);
425 N_VDestroy(IDA_mem->ida_yp0);
426
427 /* Here sensi is SUNTRUE, so deallocate sensitivity temporary vectors. */
428 N_VDestroyVectorArray(IDA_mem->ida_yyS0, IDA_mem->ida_Ns);
429 N_VDestroyVectorArray(IDA_mem->ida_ypS0, IDA_mem->ida_Ns);
430
431
432 /* On any failure, print message and return proper flag. */
433 if(retval != IDA_SUCCESS) {
434 icret = IDAICFailFlag(IDA_mem, retval);
435 return(icret);
436 }
437
438 /* Otherwise return success flag. */
439
440 return(IDA_SUCCESS);
441
442 }
443
444 /*
445 * =================================================================
446 * PRIVATE FUNCTIONS IMPLEMENTATION
447 * =================================================================
448 */
449
450 /*
451 * -----------------------------------------------------------------
452 * IDANlsIC
453 * -----------------------------------------------------------------
454 * IDANlsIC solves a nonlinear system for consistent initial
455 * conditions. It calls IDANewtonIC to do most of the work.
456 *
457 * The return value is IDA_SUCCESS = 0 if no error occurred.
458 * The error return values (positive) considered recoverable are:
459 * IC_FAIL_RECOV if res, lsetup, or lsolve failed recoverably
460 * IC_CONSTR_FAILED if the constraints could not be met
461 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
462 * or on maxbacks test)
463 * IC_CONV_FAIL if the Newton iterations failed to converge
464 * IC_SLOW_CONVRG if the iterations are converging slowly
465 * (failed the convergence test, but showed
466 * norm reduction or convergence rate < 1)
467 * The error return values (negative) considered non-recoverable are:
468 * IDA_RES_FAIL if res had a non-recoverable error
469 * IDA_FIRST_RES_FAIL if res failed recoverably on the first call
470 * IDA_LSETUP_FAIL if lsetup had a non-recoverable error
471 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
472 * -----------------------------------------------------------------
473 */
474
IDANlsIC(IDAMem IDA_mem)475 static int IDANlsIC(IDAMem IDA_mem)
476 {
477 int retval, nj, is;
478 N_Vector tv1, tv2, tv3;
479 booleantype sensi_sim;
480
481 /* Are we computing sensitivities with the IDA_SIMULTANEOUS approach? */
482 sensi_sim = (IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS));
483
484 tv1 = IDA_mem->ida_ee;
485 tv2 = IDA_mem->ida_tempv2;
486 tv3 = IDA_mem->ida_phi[2];
487
488 /* Evaluate RHS. */
489 retval = IDA_mem->ida_res(IDA_mem->ida_t0, IDA_mem->ida_yy0, IDA_mem->ida_yp0,
490 IDA_mem->ida_delta, IDA_mem->ida_user_data);
491 IDA_mem->ida_nre++;
492 if(retval < 0) return(IDA_RES_FAIL);
493 if(retval > 0) return(IDA_FIRST_RES_FAIL);
494
495 /* Save the residual. */
496 N_VScale(ONE, IDA_mem->ida_delta, IDA_mem->ida_savres);
497
498 if(sensi_sim) {
499
500 /*Evaluate sensitivity RHS and save it in savresS. */
501 retval = IDA_mem->ida_resS(IDA_mem->ida_Ns, IDA_mem->ida_t0,
502 IDA_mem->ida_yy0, IDA_mem->ida_yp0,
503 IDA_mem->ida_delta,
504 IDA_mem->ida_yyS0, IDA_mem->ida_ypS0,
505 IDA_mem->ida_deltaS,
506 IDA_mem->ida_user_dataS,
507 IDA_mem->ida_tmpS1, IDA_mem->ida_tmpS2,
508 IDA_mem->ida_tmpS3);
509 IDA_mem->ida_nrSe++;
510 if(retval < 0) return(IDA_RES_FAIL);
511 if(retval > 0) return(IDA_FIRST_RES_FAIL);
512
513 for(is=0; is<IDA_mem->ida_Ns; is++)
514 N_VScale(ONE, IDA_mem->ida_deltaS[is], IDA_mem->ida_savresS[is]);
515 }
516
517 /* Loop over nj = number of linear solve Jacobian setups. */
518 for(nj = 1; nj <= IDA_mem->ida_maxnj; nj++) {
519
520 /* If there is a setup routine, call it. */
521 if(IDA_mem->ida_lsetup) {
522 IDA_mem->ida_nsetups++;
523 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy0,
524 IDA_mem->ida_yp0, IDA_mem->ida_delta,
525 tv1, tv2, tv3);
526 if(retval < 0) return(IDA_LSETUP_FAIL);
527 if(retval > 0) return(IC_FAIL_RECOV);
528 }
529
530 /* Call the Newton iteration routine, and return if successful. */
531 retval = IDANewtonIC(IDA_mem);
532 if(retval == IDA_SUCCESS) return(IDA_SUCCESS);
533
534 /* If converging slowly and lsetup is nontrivial, retry. */
535 if(retval == IC_SLOW_CONVRG && IDA_mem->ida_lsetup) {
536 N_VScale(ONE, IDA_mem->ida_savres, IDA_mem->ida_delta);
537
538 if(sensi_sim)
539 for(is=0; is<IDA_mem->ida_Ns; is++)
540 N_VScale(ONE, IDA_mem->ida_savresS[is], IDA_mem->ida_deltaS[is]);
541
542 continue;
543 } else {
544 return(retval);
545 }
546
547 } /* End of nj loop */
548
549 /* No convergence after maxnj tries; return with retval=IC_SLOW_CONVRG */
550 return(retval);
551
552 }
553
554 /*
555 * -----------------------------------------------------------------
556 * IDANewtonIC
557 * -----------------------------------------------------------------
558 * IDANewtonIC performs the Newton iteration to solve for consistent
559 * initial conditions. It calls IDALineSrch within each iteration.
560 * On return, savres contains the current residual vector.
561 *
562 * The return value is IDA_SUCCESS = 0 if no error occurred.
563 * The error return values (positive) considered recoverable are:
564 * IC_FAIL_RECOV if res or lsolve failed recoverably
565 * IC_CONSTR_FAILED if the constraints could not be met
566 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
567 * or on maxbacks test)
568 * IC_CONV_FAIL if the Newton iterations failed to converge
569 * IC_SLOW_CONVRG if the iterations appear to be converging slowly.
570 * They failed the convergence test, but showed
571 * an overall norm reduction (by a factor of < 0.1)
572 * or a convergence rate <= ICRATEMAX).
573 * The error return values (negative) considered non-recoverable are:
574 * IDA_RES_FAIL if res had a non-recoverable error
575 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
576 * -----------------------------------------------------------------
577 */
578
IDANewtonIC(IDAMem IDA_mem)579 static int IDANewtonIC(IDAMem IDA_mem)
580 {
581 int retval, mnewt, is;
582 realtype delnorm, fnorm, fnorm0, oldfnrm, rate;
583 booleantype sensi_sim;
584
585 /* Are we computing sensitivities with the IDA_SIMULTANEOUS approach? */
586 sensi_sim = (IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS));
587
588 /* Set pointer for vector delnew */
589 IDA_mem->ida_delnew = IDA_mem->ida_phi[2];
590
591 /* Call the linear solve function to get the Newton step, delta. */
592 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delta,
593 IDA_mem->ida_ewt, IDA_mem->ida_yy0,
594 IDA_mem->ida_yp0, IDA_mem->ida_savres);
595 if(retval < 0) return(IDA_LSOLVE_FAIL);
596 if(retval > 0) return(IC_FAIL_RECOV);
597
598 /* Compute the norm of the step. */
599 fnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delta, IDA_mem->ida_ewt, SUNFALSE);
600
601 /* Call the lsolve function to get correction vectors deltaS. */
602 if (sensi_sim) {
603 for(is=0;is<IDA_mem->ida_Ns;is++) {
604 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_deltaS[is],
605 IDA_mem->ida_ewtS[is], IDA_mem->ida_yy0,
606 IDA_mem->ida_yp0, IDA_mem->ida_savres);
607 if(retval < 0) return(IDA_LSOLVE_FAIL);
608 if(retval > 0) return(IC_FAIL_RECOV);
609 }
610 /* Update the norm of delta. */
611 fnorm = IDASensWrmsNormUpdate(IDA_mem, fnorm, IDA_mem->ida_deltaS,
612 IDA_mem->ida_ewtS, SUNFALSE);
613 }
614
615 /* Test for convergence. Return now if the norm is small. */
616 if(IDA_mem->ida_sysindex == 0)
617 fnorm *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
618 if(fnorm <= IDA_mem->ida_epsNewt) return(IDA_SUCCESS);
619 fnorm0 = fnorm;
620
621 /* Initialize rate to avoid compiler warning message */
622 rate = ZERO;
623
624 /* Newton iteration loop */
625
626 for(mnewt = 0; mnewt < IDA_mem->ida_maxnit; mnewt++) {
627
628 IDA_mem->ida_nni++;
629 delnorm = fnorm;
630 oldfnrm = fnorm;
631
632 /* Call the Linesearch function and return if it failed. */
633 retval = IDALineSrch(IDA_mem, &delnorm, &fnorm);
634 if(retval != IDA_SUCCESS) return(retval);
635
636 /* Set the observed convergence rate and test for convergence. */
637 rate = fnorm/oldfnrm;
638 if(fnorm <= IDA_mem->ida_epsNewt) return(IDA_SUCCESS);
639
640 /* If not converged, copy new step vector, and loop. */
641 N_VScale(ONE, IDA_mem->ida_delnew, IDA_mem->ida_delta);
642
643 if(sensi_sim) {
644 /* Update the iteration's step for sensitivities. */
645 for(is=0; is<IDA_mem->ida_Ns; is++)
646 N_VScale(ONE, IDA_mem->ida_delnewS[is], IDA_mem->ida_deltaS[is]);
647 }
648
649 } /* End of Newton iteration loop */
650
651 /* Return either IC_SLOW_CONVRG or recoverable fail flag. */
652 if(rate <= ICRATEMAX || fnorm < PT1*fnorm0) return(IC_SLOW_CONVRG);
653 return(IC_CONV_FAIL);
654 }
655
656 /*
657 * -----------------------------------------------------------------
658 * IDALineSrch
659 * -----------------------------------------------------------------
660 * IDALineSrch performs the Linesearch algorithm with the
661 * calculation of consistent initial conditions.
662 *
663 * On entry, yy0 and yp0 are the current values of y and y', the
664 * Newton step is delta, the current residual vector F is savres,
665 * delnorm is WRMS-norm(delta), and fnorm is the norm of the vector
666 * J-inverse F.
667 *
668 * On a successful return, yy0, yp0, and savres have been updated,
669 * delnew contains the current value of J-inverse F, and fnorm is
670 * WRMS-norm(delnew).
671 *
672 * The return value is IDA_SUCCESS = 0 if no error occurred.
673 * The error return values (positive) considered recoverable are:
674 * IC_FAIL_RECOV if res or lsolve failed recoverably
675 * IC_CONSTR_FAILED if the constraints could not be met
676 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
677 * or on maxbacks test)
678 * The error return values (negative) considered non-recoverable are:
679 * IDA_RES_FAIL if res had a non-recoverable error
680 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
681 * -----------------------------------------------------------------
682 */
683
IDALineSrch(IDAMem IDA_mem,realtype * delnorm,realtype * fnorm)684 static int IDALineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm)
685 {
686 booleantype conOK;
687 int retval, is, nbacks;
688 realtype f1norm, fnormp, f1normp, ratio, lambda, minlam, slpi;
689 N_Vector mc;
690 booleantype sensi_sim;
691
692 /* Initialize work space pointers, f1norm, ratio.
693 (Use of mc in constraint check does not conflict with ypnew.) */
694 mc = IDA_mem->ida_ee;
695 IDA_mem->ida_dtemp = IDA_mem->ida_phi[3];
696 IDA_mem->ida_ynew = IDA_mem->ida_tempv2;
697 IDA_mem->ida_ypnew = IDA_mem->ida_ee;
698 f1norm = (*fnorm)*(*fnorm)*HALF;
699 ratio = ONE;
700
701 /* If there are constraints, check and reduce step if necessary. */
702 if(IDA_mem->ida_constraintsSet) {
703
704 /* Update y and check constraints. */
705 IDANewy(IDA_mem);
706 conOK = N_VConstrMask(IDA_mem->ida_constraints,
707 IDA_mem->ida_ynew, mc);
708
709 if(!conOK) {
710 /* Not satisfied. Compute scaled step to satisfy constraints. */
711 N_VProd(mc, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
712 ratio = PT99*N_VMinQuotient(IDA_mem->ida_yy0, IDA_mem->ida_dtemp);
713 (*delnorm) *= ratio;
714 if((*delnorm) <= IDA_mem->ida_steptol)
715 return(IC_CONSTR_FAILED);
716 N_VScale(ratio, IDA_mem->ida_delta, IDA_mem->ida_delta);
717 }
718
719 } /* End of constraints check */
720
721 slpi = -TWO*f1norm*ratio;
722 minlam = IDA_mem->ida_steptol / (*delnorm);
723 lambda = ONE;
724 nbacks = 0;
725
726 /* Are we computing sensitivities with the IDA_SIMULTANEOUS approach? */
727 sensi_sim = (IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS));
728
729 /* In IDA_Y_INIT case, set ypnew = yp0 (fixed) for linesearch. */
730 if(IDA_mem->ida_icopt == IDA_Y_INIT) {
731 N_VScale(ONE, IDA_mem->ida_yp0, IDA_mem->ida_ypnew);
732
733 /* do the same for sensitivities. */
734 if(sensi_sim) {
735 for(is=0; is<IDA_mem->ida_Ns; is++)
736 N_VScale(ONE, IDA_mem->ida_ypS0[is], IDA_mem->ida_ypS0new[is]);
737 }
738 }
739
740 /* Loop on linesearch variable lambda. */
741
742 for(;;) {
743
744 if (nbacks == IDA_mem->ida_maxbacks)
745 return(IC_LINESRCH_FAILED);
746 /* Get new (y,y') = (ynew,ypnew) and norm of new function value. */
747 IDANewyyp(IDA_mem, lambda);
748 retval = IDAfnorm(IDA_mem, &fnormp);
749 if(retval != IDA_SUCCESS) return(retval);
750
751 /* If lsoff option is on, break out. */
752 if(IDA_mem->ida_lsoff) break;
753
754 /* Do alpha-condition test. */
755 f1normp = fnormp*fnormp*HALF;
756 if(f1normp <= f1norm + ALPHALS*slpi*lambda) break;
757 if(lambda < minlam) return(IC_LINESRCH_FAILED);
758 lambda /= TWO;
759 IDA_mem->ida_nbacktr++; nbacks++;
760
761 } /* End of breakout linesearch loop */
762
763 /* Update yy0, yp0. */
764 N_VScale(ONE, IDA_mem->ida_ynew, IDA_mem->ida_yy0);
765
766 if(sensi_sim) {
767 /* Update yyS0 and ypS0. */
768 for(is=0; is<IDA_mem->ida_Ns; is++)
769 N_VScale(ONE, IDA_mem->ida_yyS0new[is], IDA_mem->ida_yyS0[is]);
770 }
771
772 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
773 N_VScale(ONE, IDA_mem->ida_ypnew, IDA_mem->ida_yp0);
774
775 if(sensi_sim)
776 for(is=0; is<IDA_mem->ida_Ns; is++)
777 N_VScale(ONE, IDA_mem->ida_ypS0new[is], IDA_mem->ida_ypS0[is]);
778
779 }
780 /* Update fnorm, then return. */
781 *fnorm = fnormp;
782 return(IDA_SUCCESS);
783
784 }
785
786 /*
787 * -----------------------------------------------------------------
788 * IDAfnorm
789 * -----------------------------------------------------------------
790 * IDAfnorm computes the norm of the current function value, by
791 * evaluating the DAE residual function, calling the linear
792 * system solver, and computing a WRMS-norm.
793 *
794 * On return, savres contains the current residual vector F, and
795 * delnew contains J-inverse F.
796 *
797 * The return value is IDA_SUCCESS = 0 if no error occurred, or
798 * IC_FAIL_RECOV if res or lsolve failed recoverably, or
799 * IDA_RES_FAIL if res had a non-recoverable error, or
800 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error.
801 * -----------------------------------------------------------------
802 */
803
IDAfnorm(IDAMem IDA_mem,realtype * fnorm)804 static int IDAfnorm(IDAMem IDA_mem, realtype *fnorm)
805 {
806 int retval, is;
807
808 /* Get residual vector F, return if failed, and save F in savres. */
809 retval = IDA_mem->ida_res(IDA_mem->ida_t0, IDA_mem->ida_ynew,
810 IDA_mem->ida_ypnew, IDA_mem->ida_delnew,
811 IDA_mem->ida_user_data);
812 IDA_mem->ida_nre++;
813 if(retval < 0) return(IDA_RES_FAIL);
814 if(retval > 0) return(IC_FAIL_RECOV);
815
816 N_VScale(ONE, IDA_mem->ida_delnew, IDA_mem->ida_savres);
817
818 /* Call the linear solve function to get J-inverse F; return if failed. */
819 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delnew, IDA_mem->ida_ewt,
820 IDA_mem->ida_ynew, IDA_mem->ida_ypnew,
821 IDA_mem->ida_savres);
822 if(retval < 0) return(IDA_LSOLVE_FAIL);
823 if(retval > 0) return(IC_FAIL_RECOV);
824
825 /* Compute the WRMS-norm. */
826 *fnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delnew, IDA_mem->ida_ewt, SUNFALSE);
827
828
829 /* Are we computing SENSITIVITIES with the IDA_SIMULTANEOUS approach? */
830
831 if(IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS)) {
832
833 /* Evaluate the residual for sensitivities. */
834 retval = IDA_mem->ida_resS(IDA_mem->ida_Ns, IDA_mem->ida_t0,
835 IDA_mem->ida_ynew, IDA_mem->ida_ypnew,
836 IDA_mem->ida_savres,
837 IDA_mem->ida_yyS0new,
838 IDA_mem->ida_ypS0new,
839 IDA_mem->ida_delnewS,
840 IDA_mem->ida_user_dataS,
841 IDA_mem->ida_tmpS1, IDA_mem->ida_tmpS2,
842 IDA_mem->ida_tmpS3);
843 IDA_mem->ida_nrSe++;
844 if(retval < 0) return(IDA_RES_FAIL);
845 if(retval > 0) return(IC_FAIL_RECOV);
846
847 /* Save delnewS in savresS. */
848 for(is=0; is<IDA_mem->ida_Ns; is++)
849 N_VScale(ONE, IDA_mem->ida_delnewS[is], IDA_mem->ida_savresS[is]);
850
851 /* Call the linear solve function to get J-inverse deltaS. */
852 for(is=0; is<IDA_mem->ida_Ns; is++) {
853
854 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delnewS[is],
855 IDA_mem->ida_ewtS[is],
856 IDA_mem->ida_ynew,
857 IDA_mem->ida_ypnew,
858 IDA_mem->ida_savres);
859 if(retval < 0) return(IDA_LSOLVE_FAIL);
860 if(retval > 0) return(IC_FAIL_RECOV);
861 }
862
863 /* Include sensitivities in norm. */
864 *fnorm = IDASensWrmsNormUpdate(IDA_mem, *fnorm, IDA_mem->ida_delnewS,
865 IDA_mem->ida_ewtS, SUNFALSE);
866 }
867
868 /* Rescale norm if index = 0. */
869 if(IDA_mem->ida_sysindex == 0)
870 (*fnorm) *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
871
872 return(IDA_SUCCESS);
873
874 }
875
876 /*
877 * -----------------------------------------------------------------
878 * IDANewyyp
879 * -----------------------------------------------------------------
880 * IDANewyyp updates the vectors ynew and ypnew from yy0 and yp0,
881 * using the current step vector lambda*delta, in a manner
882 * depending on icopt and the input id vector.
883 *
884 * The return value is always IDA_SUCCESS = 0.
885 * -----------------------------------------------------------------
886 */
887
IDANewyyp(IDAMem IDA_mem,realtype lambda)888 static int IDANewyyp(IDAMem IDA_mem, realtype lambda)
889 {
890 int retval;
891
892 retval = IDA_SUCCESS;
893
894 /* IDA_YA_YDP_INIT case: ynew = yy0 - lambda*delta where id_i = 0
895 ypnew = yp0 - cj*lambda*delta where id_i = 1. */
896 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
897
898 N_VProd(IDA_mem->ida_id, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
899 N_VLinearSum(ONE, IDA_mem->ida_yp0, -IDA_mem->ida_cj*lambda,
900 IDA_mem->ida_dtemp, IDA_mem->ida_ypnew);
901 N_VLinearSum(ONE, IDA_mem->ida_delta, -ONE,
902 IDA_mem->ida_dtemp, IDA_mem->ida_dtemp);
903 N_VLinearSum(ONE, IDA_mem->ida_yy0, -lambda,
904 IDA_mem->ida_dtemp, IDA_mem->ida_ynew);
905
906 }else if(IDA_mem->ida_icopt == IDA_Y_INIT) {
907
908 /* IDA_Y_INIT case: ynew = yy0 - lambda*delta. (ypnew = yp0 preset.) */
909 N_VLinearSum(ONE, IDA_mem->ida_yy0, -lambda, IDA_mem->ida_delta,
910 IDA_mem->ida_ynew);
911 }
912
913 if(IDA_mem->ida_sensi && (IDA_mem->ida_ism==IDA_SIMULTANEOUS))
914 retval = IDASensNewyyp(IDA_mem, lambda);
915
916 return(retval);
917
918 }
919
920 /*
921 * -----------------------------------------------------------------
922 * IDANewy
923 * -----------------------------------------------------------------
924 * IDANewy updates the vector ynew from yy0,
925 * using the current step vector delta, in a manner
926 * depending on icopt and the input id vector.
927 *
928 * The return value is always IDA_SUCCESS = 0.
929 * -----------------------------------------------------------------
930 */
931
IDANewy(IDAMem IDA_mem)932 static int IDANewy(IDAMem IDA_mem)
933 {
934
935 /* IDA_YA_YDP_INIT case: ynew = yy0 - delta where id_i = 0. */
936 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
937 N_VProd(IDA_mem->ida_id, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
938 N_VLinearSum(ONE, IDA_mem->ida_delta, -ONE,
939 IDA_mem->ida_dtemp, IDA_mem->ida_dtemp);
940 N_VLinearSum(ONE, IDA_mem->ida_yy0, -ONE,
941 IDA_mem->ida_dtemp, IDA_mem->ida_ynew);
942 return(IDA_SUCCESS);
943 }
944
945 /* IDA_Y_INIT case: ynew = yy0 - delta. */
946 N_VLinearSum(ONE, IDA_mem->ida_yy0, -ONE, IDA_mem->ida_delta,
947 IDA_mem->ida_ynew);
948 return(IDA_SUCCESS);
949
950 }
951 /*
952 * -----------------------------------------------------------------
953 * Sensitivity I.C. functions
954 * -----------------------------------------------------------------
955 */
956
957 /*
958 * -----------------------------------------------------------------
959 * IDASensNlsIC
960 * -----------------------------------------------------------------
961 * IDASensNlsIC solves nonlinear systems for sensitivities consistent
962 * initial conditions. It mainly relies on IDASensNewtonIC.
963 *
964 * The return value is IDA_SUCCESS = 0 if no error occurred.
965 * The error return values (positive) considered recoverable are:
966 * IC_FAIL_RECOV if res, lsetup, or lsolve failed recoverably
967 * IC_CONSTR_FAILED if the constraints could not be met
968 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
969 * or on maxbacks test)
970 * IC_CONV_FAIL if the Newton iterations failed to converge
971 * IC_SLOW_CONVRG if the iterations are converging slowly
972 * (failed the convergence test, but showed
973 * norm reduction or convergence rate < 1)
974 * The error return values (negative) considered non-recoverable are:
975 * IDA_RES_FAIL if res had a non-recoverable error
976 * IDA_FIRST_RES_FAIL if res failed recoverably on the first call
977 * IDA_LSETUP_FAIL if lsetup had a non-recoverable error
978 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
979 * -----------------------------------------------------------------
980 */
IDASensNlsIC(IDAMem IDA_mem)981 static int IDASensNlsIC(IDAMem IDA_mem)
982 {
983 int retval;
984 int is, nj;
985
986 retval = IDA_mem->ida_resS(IDA_mem->ida_Ns, IDA_mem->ida_t0,
987 IDA_mem->ida_yy0, IDA_mem->ida_yp0,
988 IDA_mem->ida_delta, IDA_mem->ida_yyS0,
989 IDA_mem->ida_ypS0,
990 IDA_mem->ida_deltaS,
991 IDA_mem->ida_user_dataS,
992 IDA_mem->ida_tmpS1, IDA_mem->ida_tmpS2,
993 IDA_mem->ida_tmpS3);
994 IDA_mem->ida_nrSe++;
995 if(retval < 0) return(IDA_RES_FAIL);
996 if(retval > 0) return(IDA_FIRST_RES_FAIL);
997
998 /* Save deltaS */
999 for(is=0; is<IDA_mem->ida_Ns; is++)
1000 N_VScale(ONE, IDA_mem->ida_deltaS[is], IDA_mem->ida_savresS[is]);
1001
1002 /* Loop over nj = number of linear solve Jacobian setups. */
1003
1004 for(nj = 1; nj <= 2; nj++) {
1005
1006 /* Call the Newton iteration routine */
1007 retval = IDASensNewtonIC(IDA_mem);
1008 if(retval == IDA_SUCCESS) return(IDA_SUCCESS);
1009
1010 /* If converging slowly and lsetup is nontrivial and this is the first pass,
1011 update Jacobian and retry. */
1012 if(retval == IC_SLOW_CONVRG && IDA_mem->ida_lsetup && nj==1) {
1013
1014 /* Restore deltaS. */
1015 for(is=0; is<IDA_mem->ida_Ns; is++)
1016 N_VScale(ONE, IDA_mem->ida_savresS[is], IDA_mem->ida_deltaS[is]);
1017
1018 IDA_mem->ida_nsetupsS++;
1019 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy0, IDA_mem->ida_yp0,
1020 IDA_mem->ida_delta, IDA_mem->ida_tmpS1,
1021 IDA_mem->ida_tmpS2, IDA_mem->ida_tmpS3);
1022 if(retval < 0) return(IDA_LSETUP_FAIL);
1023 if(retval > 0) return(IC_FAIL_RECOV);
1024
1025 continue;
1026 } else {
1027 return(retval);
1028 }
1029 }
1030
1031 return(IDA_SUCCESS);
1032 }
1033
1034 /*
1035 * -----------------------------------------------------------------
1036 * IDASensNewtonIC
1037 * -----------------------------------------------------------------
1038 * IDANewtonIC performs the Newton iteration to solve for
1039 * sensitivities consistent initial conditions. It calls
1040 * IDASensLineSrch within each iteration.
1041 * On return, savresS contains the current residual vectors.
1042 *
1043 * The return value is IDA_SUCCESS = 0 if no error occurred.
1044 * The error return values (positive) considered recoverable are:
1045 * IC_FAIL_RECOV if res or lsolve failed recoverably
1046 * IC_CONSTR_FAILED if the constraints could not be met
1047 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
1048 * or on maxbacks test)
1049 * IC_CONV_FAIL if the Newton iterations failed to converge
1050 * IC_SLOW_CONVRG if the iterations appear to be converging slowly.
1051 * They failed the convergence test, but showed
1052 * an overall norm reduction (by a factor of < 0.1)
1053 * or a convergence rate <= ICRATEMAX).
1054 * The error return values (negative) considered non-recoverable are:
1055 * IDA_RES_FAIL if res had a non-recoverable error
1056 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
1057 * -----------------------------------------------------------------
1058 */
IDASensNewtonIC(IDAMem IDA_mem)1059 static int IDASensNewtonIC(IDAMem IDA_mem)
1060 {
1061 int retval, is, mnewt;
1062 realtype delnorm, fnorm, fnorm0, oldfnrm, rate;
1063
1064 for(is=0;is<IDA_mem->ida_Ns;is++) {
1065
1066 /* Call the linear solve function to get the Newton step, delta. */
1067 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_deltaS[is],
1068 IDA_mem->ida_ewtS[is], IDA_mem->ida_yy0,
1069 IDA_mem->ida_yp0, IDA_mem->ida_delta);
1070 if(retval < 0) return(IDA_LSOLVE_FAIL);
1071 if(retval > 0) return(IC_FAIL_RECOV);
1072
1073 }
1074 /* Compute the norm of the step and return if it is small enough */
1075 fnorm = IDASensWrmsNorm(IDA_mem, IDA_mem->ida_deltaS,
1076 IDA_mem->ida_ewtS, SUNFALSE);
1077 if(IDA_mem->ida_sysindex == 0)
1078 fnorm *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
1079 if(fnorm <= IDA_mem->ida_epsNewt) return(IDA_SUCCESS);
1080 fnorm0 = fnorm;
1081
1082 rate = ZERO;
1083
1084 /* Newton iteration loop */
1085 for(mnewt = 0; mnewt < IDA_mem->ida_maxnit; mnewt++) {
1086
1087 IDA_mem->ida_nniS++;
1088 delnorm = fnorm;
1089 oldfnrm = fnorm;
1090
1091 /* Call the Linesearch function and return if it failed. */
1092 retval = IDASensLineSrch(IDA_mem, &delnorm, &fnorm);
1093 if(retval != IDA_SUCCESS) return(retval);
1094
1095 /* Set the observed convergence rate and test for convergence. */
1096 rate = fnorm/oldfnrm;
1097 if(fnorm <= IDA_mem->ida_epsNewt) return(IDA_SUCCESS);
1098
1099 /* If not converged, copy new step vectors, and loop. */
1100 for(is=0; is<IDA_mem->ida_Ns; is++)
1101 N_VScale(ONE, IDA_mem->ida_delnewS[is], IDA_mem->ida_deltaS[is]);
1102
1103 } /* End of Newton iteration loop */
1104
1105 /* Return either IC_SLOW_CONVRG or recoverable fail flag. */
1106 if(rate <= ICRATEMAX || fnorm < PT1*fnorm0) return(IC_SLOW_CONVRG);
1107 return(IC_CONV_FAIL);
1108 }
1109
1110 /*
1111 * -----------------------------------------------------------------
1112 * IDASensLineSrch
1113 * -----------------------------------------------------------------
1114 * IDASensLineSrch performs the Linesearch algorithm with the
1115 * calculation of consistent initial conditions for sensitivities
1116 * systems.
1117 *
1118 * On entry, yyS0 and ypS0 contain the current values, the Newton
1119 * steps are contained in deltaS, the current residual vectors FS are
1120 * savresS, delnorm is sens-WRMS-norm(deltaS), and fnorm is
1121 * max { WRMS-norm( J-inverse FS[is] ) : is=1,2,...,Ns }
1122 *
1123 * On a successful return, yy0, yp0, and savres have been updated,
1124 * delnew contains the current values of J-inverse FS, and fnorm is
1125 * max { WRMS-norm(delnewS[is]) : is = 1,2,...Ns }
1126 *
1127 * The return value is IDA_SUCCESS = 0 if no error occurred.
1128 * The error return values (positive) considered recoverable are:
1129 * IC_FAIL_RECOV if res or lsolve failed recoverably
1130 * IC_CONSTR_FAILED if the constraints could not be met
1131 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
1132 * or on maxbacks test)
1133 * The error return values (negative) considered non-recoverable are:
1134 * IDA_RES_FAIL if res had a non-recoverable error
1135 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
1136 * -----------------------------------------------------------------
1137 */
1138
IDASensLineSrch(IDAMem IDA_mem,realtype * delnorm,realtype * fnorm)1139 static int IDASensLineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm)
1140 {
1141 int is, retval, nbacks;
1142 realtype f1norm, fnormp, f1normp, slpi, minlam;
1143 realtype lambda, ratio;
1144
1145 /* Set work space pointer. */
1146 IDA_mem->ida_dtemp = IDA_mem->ida_phi[3];
1147
1148 f1norm = (*fnorm)*(*fnorm)*HALF;
1149
1150 /* Initialize local variables. */
1151 ratio = ONE;
1152 slpi = -TWO*f1norm*ratio;
1153 minlam = IDA_mem->ida_steptol / (*delnorm);
1154 lambda = ONE;
1155 nbacks = 0;
1156
1157 for(;;) {
1158
1159 if (nbacks == IDA_mem->ida_maxbacks)
1160 return(IC_LINESRCH_FAILED);
1161 /* Get new iteration in (ySnew, ypSnew). */
1162 IDASensNewyyp(IDA_mem, lambda);
1163
1164 /* Get the norm of new function value. */
1165 retval = IDASensfnorm(IDA_mem, &fnormp);
1166 if (retval!=IDA_SUCCESS) return retval;
1167
1168 /* If lsoff option is on, break out. */
1169 if(IDA_mem->ida_lsoff) break;
1170
1171 /* Do alpha-condition test. */
1172 f1normp = fnormp*fnormp*HALF;
1173 if(f1normp <= f1norm + ALPHALS*slpi*lambda) break;
1174 if(lambda < minlam) return(IC_LINESRCH_FAILED);
1175 lambda /= TWO;
1176 IDA_mem->ida_nbacktr++; nbacks++;
1177 }
1178
1179 /* Update yyS0, ypS0 and fnorm and return. */
1180 for(is=0; is<IDA_mem->ida_Ns; is++) {
1181 N_VScale(ONE, IDA_mem->ida_yyS0new[is], IDA_mem->ida_yyS0[is]);
1182 }
1183
1184 if (IDA_mem->ida_icopt == IDA_YA_YDP_INIT)
1185 for(is=0; is<IDA_mem->ida_Ns; is++)
1186 N_VScale(ONE, IDA_mem->ida_ypS0new[is], IDA_mem->ida_ypS0[is]);
1187
1188 *fnorm = fnormp;
1189 return(IDA_SUCCESS);
1190 }
1191
1192 /*
1193 * -----------------------------------------------------------------
1194 * IDASensfnorm
1195 * -----------------------------------------------------------------
1196 * IDASensfnorm computes the norm of the current function value, by
1197 * evaluating the sensitivity residual function, calling the linear
1198 * system solver, and computing a WRMS-norm.
1199 *
1200 * On return, savresS contains the current residual vectors FS, and
1201 * delnewS contains J-inverse FS.
1202 *
1203 * The return value is IDA_SUCCESS = 0 if no error occurred, or
1204 * IC_FAIL_RECOV if res or lsolve failed recoverably, or
1205 * IDA_RES_FAIL if res had a non-recoverable error, or
1206 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error.
1207 * -----------------------------------------------------------------
1208 */
1209
IDASensfnorm(IDAMem IDA_mem,realtype * fnorm)1210 static int IDASensfnorm(IDAMem IDA_mem, realtype *fnorm)
1211 {
1212 int is, retval;
1213
1214 /* Get sensitivity residual */
1215 retval = IDA_mem->ida_resS(IDA_mem->ida_Ns, IDA_mem->ida_t0,
1216 IDA_mem->ida_yy0, IDA_mem->ida_yp0,
1217 IDA_mem->ida_delta,
1218 IDA_mem->ida_yyS0new,
1219 IDA_mem->ida_ypS0new,
1220 IDA_mem->ida_delnewS,
1221 IDA_mem->ida_user_dataS,
1222 IDA_mem->ida_tmpS1, IDA_mem->ida_tmpS2,
1223 IDA_mem->ida_tmpS3);
1224 IDA_mem->ida_nrSe++;
1225 if(retval < 0) return(IDA_RES_FAIL);
1226 if(retval > 0) return(IC_FAIL_RECOV);
1227
1228 for(is=0; is<IDA_mem->ida_Ns; is++)
1229 N_VScale(ONE, IDA_mem->ida_delnewS[is], IDA_mem->ida_savresS[is]);
1230
1231 /* Call linear solve function */
1232 for(is=0; is<IDA_mem->ida_Ns; is++) {
1233
1234 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delnewS[is],
1235 IDA_mem->ida_ewtS[is],
1236 IDA_mem->ida_yy0,
1237 IDA_mem->ida_yp0,
1238 IDA_mem->ida_delta);
1239 if(retval < 0) return(IDA_LSOLVE_FAIL);
1240 if(retval > 0) return(IC_FAIL_RECOV);
1241 }
1242
1243 /* Compute the WRMS-norm; rescale if index = 0. */
1244 *fnorm = IDASensWrmsNorm(IDA_mem, IDA_mem->ida_delnewS, IDA_mem->ida_ewtS, SUNFALSE);
1245 if(IDA_mem->ida_sysindex == 0)
1246 (*fnorm) *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
1247
1248 return(IDA_SUCCESS);
1249 }
1250
1251 /*
1252 * -----------------------------------------------------------------
1253 * IDASensNewyyp
1254 * -----------------------------------------------------------------
1255 * IDASensNewyyp computes the Newton updates for each of the
1256 * sensitivities systems using the current step vector lambda*delta,
1257 * in a manner depending on icopt and the input id vector.
1258 *
1259 * The return value is always IDA_SUCCESS = 0.
1260 * -----------------------------------------------------------------
1261 */
1262
IDASensNewyyp(IDAMem IDA_mem,realtype lambda)1263 static int IDASensNewyyp(IDAMem IDA_mem, realtype lambda)
1264 {
1265 int is;
1266
1267 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
1268
1269 /* IDA_YA_YDP_INIT case:
1270 - ySnew = yS0 - lambda*deltaS where id_i = 0
1271 - ypSnew = ypS0 - cj*lambda*delta where id_i = 1. */
1272
1273 for(is=0; is<IDA_mem->ida_Ns; is++) {
1274
1275 /* It is ok to use dtemp as temporary vector here. */
1276 N_VProd(IDA_mem->ida_id, IDA_mem->ida_deltaS[is], IDA_mem->ida_dtemp);
1277 N_VLinearSum(ONE, IDA_mem->ida_ypS0[is], -IDA_mem->ida_cj*lambda,
1278 IDA_mem->ida_dtemp, IDA_mem->ida_ypS0new[is]);
1279 N_VLinearSum(ONE, IDA_mem->ida_deltaS[is], -ONE,
1280 IDA_mem->ida_dtemp, IDA_mem->ida_dtemp);
1281 N_VLinearSum(ONE, IDA_mem->ida_yyS0[is], -lambda,
1282 IDA_mem->ida_dtemp, IDA_mem->ida_yyS0new[is]);
1283 } /* end loop is */
1284 }else {
1285
1286 /* IDA_Y_INIT case:
1287 - ySnew = yS0 - lambda*deltaS. (ypnew = yp0 preset.) */
1288
1289 for(is=0; is<IDA_mem->ida_Ns; is++)
1290 N_VLinearSum(ONE, IDA_mem->ida_yyS0[is], -lambda,
1291 IDA_mem->ida_deltaS[is], IDA_mem->ida_yyS0new[is]);
1292 } /* end loop is */
1293 return(IDA_SUCCESS);
1294 }
1295
1296 /*
1297 * -----------------------------------------------------------------
1298 * IDAICFailFlag
1299 * -----------------------------------------------------------------
1300 * IDAICFailFlag prints a message and sets the IDACalcIC return
1301 * value appropriate to the flag retval returned by IDANlsIC.
1302 * -----------------------------------------------------------------
1303 */
1304
IDAICFailFlag(IDAMem IDA_mem,int retval)1305 static int IDAICFailFlag(IDAMem IDA_mem, int retval)
1306 {
1307
1308 /* Depending on retval, print error message and return error flag. */
1309 switch(retval) {
1310
1311 case IDA_RES_FAIL:
1312 IDAProcessError(IDA_mem, IDA_RES_FAIL, "IDAS", "IDACalcIC", MSG_IC_RES_NONREC);
1313 return(IDA_RES_FAIL);
1314
1315 case IDA_FIRST_RES_FAIL:
1316 IDAProcessError(IDA_mem, IDA_FIRST_RES_FAIL, "IDAS", "IDACalcIC", MSG_IC_RES_FAIL);
1317 return(IDA_FIRST_RES_FAIL);
1318
1319 case IDA_LSETUP_FAIL:
1320 IDAProcessError(IDA_mem, IDA_LSETUP_FAIL, "IDAS", "IDACalcIC", MSG_IC_SETUP_FAIL);
1321 return(IDA_LSETUP_FAIL);
1322
1323 case IDA_LSOLVE_FAIL:
1324 IDAProcessError(IDA_mem, IDA_LSOLVE_FAIL, "IDAS", "IDACalcIC", MSG_IC_SOLVE_FAIL);
1325 return(IDA_LSOLVE_FAIL);
1326
1327 case IC_FAIL_RECOV:
1328 IDAProcessError(IDA_mem, IDA_NO_RECOVERY, "IDAS", "IDACalcIC", MSG_IC_NO_RECOVERY);
1329 return(IDA_NO_RECOVERY);
1330
1331 case IC_CONSTR_FAILED:
1332 IDAProcessError(IDA_mem, IDA_CONSTR_FAIL, "IDAS", "IDACalcIC", MSG_IC_FAIL_CONSTR);
1333 return(IDA_CONSTR_FAIL);
1334
1335 case IC_LINESRCH_FAILED:
1336 IDAProcessError(IDA_mem, IDA_LINESEARCH_FAIL, "IDAS", "IDACalcIC", MSG_IC_FAILED_LINS);
1337 return(IDA_LINESEARCH_FAIL);
1338
1339 case IC_CONV_FAIL:
1340 IDAProcessError(IDA_mem, IDA_CONV_FAIL, "IDAS", "IDACalcIC", MSG_IC_CONV_FAILED);
1341 return(IDA_CONV_FAIL);
1342
1343 case IC_SLOW_CONVRG:
1344 IDAProcessError(IDA_mem, IDA_CONV_FAIL, "IDAS", "IDACalcIC", MSG_IC_CONV_FAILED);
1345 return(IDA_CONV_FAIL);
1346
1347 case IDA_BAD_EWT:
1348 IDAProcessError(IDA_mem, IDA_BAD_EWT, "IDAS", "IDACalcIC", MSG_IC_BAD_EWT);
1349 return(IDA_BAD_EWT);
1350
1351 }
1352 return -99;
1353 }
1354