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