1 /*
2 * -----------------------------------------------------------------
3 * $Revision$
4 * $Date$
5 * -----------------------------------------------------------------
6 * Programmers: Alan C. Hindmarsh, and Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * SUNDIALS Copyright Start
9 * Copyright (c) 2002-2020, 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 IDA.
19 * It is independent of the linear solver in use.
20 * -----------------------------------------------------------------
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "ida_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, N_Vector w,
66 booleantype mask);
67
68 static int IDAnlsIC(IDAMem IDA_mem);
69 static int IDANewtonIC(IDAMem IDA_mem);
70 static int IDALineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm);
71 static int IDAfnorm(IDAMem IDA_mem, realtype *fnorm);
72 static int IDANewyyp(IDAMem IDA_mem, realtype lambda);
73 static int IDANewy(IDAMem IDA_mem);
74 static int IDAICFailFlag(IDAMem IDA_mem, int retval);
75
76 /*
77 * =================================================================
78 * EXPORTED FUNCTIONS IMPLEMENTATION
79 * =================================================================
80 */
81
82 /*
83 * -----------------------------------------------------------------
84 * IDACalcIC
85 * -----------------------------------------------------------------
86 * IDACalcIC computes consistent initial conditions, given the
87 * user's initial guess for unknown components of yy0 and/or yp0.
88 *
89 * The return value is IDA_SUCCESS = 0 if no error occurred.
90 *
91 * The error return values (fully described in ida.h) are:
92 * IDA_MEM_NULL ida_mem is NULL
93 * IDA_NO_MALLOC ida_mem was not allocated
94 * IDA_ILL_INPUT bad value for icopt, tout1, or id
95 * IDA_LINIT_FAIL the linear solver linit routine failed
96 * IDA_BAD_EWT zero value of some component of ewt
97 * IDA_RES_FAIL res had a non-recoverable error
98 * IDA_FIRST_RES_FAIL res failed recoverably on the first call
99 * IDA_LSETUP_FAIL lsetup had a non-recoverable error
100 * IDA_LSOLVE_FAIL lsolve had a non-recoverable error
101 * IDA_NO_RECOVERY res, lsetup, or lsolve had a recoverable
102 * error, but IDACalcIC could not recover
103 * IDA_CONSTR_FAIL the inequality constraints could not be met
104 * IDA_LINESEARCH_FAIL the linesearch failed (either on steptol test
105 * or on the maxbacks test)
106 * IDA_CONV_FAIL the Newton iterations failed to converge
107 * -----------------------------------------------------------------
108 */
109
IDACalcIC(void * ida_mem,int icopt,realtype tout1)110 int IDACalcIC(void *ida_mem, int icopt, realtype tout1)
111 {
112 int ewtsetOK;
113 int ier, nwt, nh, mxnh, icret, retval=0;
114 realtype tdist, troundoff, minid, hic, ypnorm;
115 IDAMem IDA_mem;
116
117 /* Check if IDA memory exists */
118
119 if(ida_mem == NULL) {
120 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDACalcIC", MSG_NO_MEM);
121 return(IDA_MEM_NULL);
122 }
123 IDA_mem = (IDAMem) ida_mem;
124
125 /* Check if problem was malloc'ed */
126
127 if(IDA_mem->ida_MallocDone == SUNFALSE) {
128 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDACalcIC", MSG_NO_MALLOC);
129 return(IDA_NO_MALLOC);
130 }
131
132 /* Check inputs to IDA for correctness and consistency */
133
134 ier = IDAInitialSetup(IDA_mem);
135 if(ier != IDA_SUCCESS) return(IDA_ILL_INPUT);
136 IDA_mem->ida_SetupDone = SUNTRUE;
137
138 /* Check legality of input arguments, and set IDA memory copies. */
139
140 if(icopt != IDA_YA_YDP_INIT && icopt != IDA_Y_INIT) {
141 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDACalcIC", MSG_IC_BAD_ICOPT);
142 return(IDA_ILL_INPUT);
143 }
144 IDA_mem->ida_icopt = icopt;
145
146 if(icopt == IDA_YA_YDP_INIT && (IDA_mem->ida_id == NULL)) {
147 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDACalcIC", MSG_IC_MISSING_ID);
148 return(IDA_ILL_INPUT);
149 }
150
151 tdist = SUNRabs(tout1 - IDA_mem->ida_tn);
152 troundoff = TWO * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(tout1));
153 if(tdist < troundoff) {
154 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDACalcIC", MSG_IC_TOO_CLOSE);
155 return(IDA_ILL_INPUT);
156 }
157
158 /* Allocate space and initialize temporary vectors */
159
160 IDA_mem->ida_yy0 = N_VClone(IDA_mem->ida_ee);
161 IDA_mem->ida_yp0 = N_VClone(IDA_mem->ida_ee);
162 IDA_mem->ida_t0 = IDA_mem->ida_tn;
163 N_VScale(ONE, IDA_mem->ida_phi[0], IDA_mem->ida_yy0);
164 N_VScale(ONE, IDA_mem->ida_phi[1], IDA_mem->ida_yp0);
165
166 /* For use in the IDA_YA_YP_INIT case, set sysindex and tscale. */
167
168 IDA_mem->ida_sysindex = 1;
169 IDA_mem->ida_tscale = tdist;
170 if(icopt == IDA_YA_YDP_INIT) {
171 minid = N_VMin(IDA_mem->ida_id);
172 if(minid < ZERO) {
173 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDACalcIC", MSG_IC_BAD_ID);
174 return(IDA_ILL_INPUT);
175 }
176 if(minid > HALF) IDA_mem->ida_sysindex = 0;
177 }
178
179 /* Set the test constant in the Newton convergence test */
180
181 IDA_mem->ida_epsNewt = IDA_mem->ida_epiccon;
182
183 /* Initializations:
184 cjratio = 1 (for use in direct linear solvers);
185 set nbacktr = 0; */
186
187 IDA_mem->ida_cjratio = ONE;
188 IDA_mem->ida_nbacktr = 0;
189
190 /* Set hic, hh, cj, and mxnh. */
191
192 hic = PT001*tdist;
193 ypnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_yp0,
194 IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
195 if(ypnorm > HALF/hic) hic = HALF/ypnorm;
196 if(tout1 < IDA_mem->ida_tn) hic = -hic;
197 IDA_mem->ida_hh = hic;
198 if(icopt == IDA_YA_YDP_INIT) {
199 IDA_mem->ida_cj = ONE/hic;
200 mxnh = IDA_mem->ida_maxnh;
201 }
202 else {
203 IDA_mem->ida_cj = ZERO;
204 mxnh = 1;
205 }
206
207 /* Loop over nwt = number of evaluations of ewt vector. */
208
209 for(nwt = 1; nwt <= 2; nwt++) {
210
211 /* Loop over nh = number of h values. */
212 for(nh = 1; nh <= mxnh; nh++) {
213
214 /* Call the IC nonlinear solver function. */
215 retval = IDAnlsIC(IDA_mem);
216
217 /* Cut h and loop on recoverable IDA_YA_YDP_INIT failure; else break. */
218 if(retval == IDA_SUCCESS) break;
219 IDA_mem->ida_ncfn++;
220 if(retval < 0) break;
221 if(nh == mxnh) break;
222 /* If looping to try again, reset yy0 and yp0 if not converging. */
223 if(retval != IC_SLOW_CONVRG) {
224 N_VScale(ONE, IDA_mem->ida_phi[0], IDA_mem->ida_yy0);
225 N_VScale(ONE, IDA_mem->ida_phi[1], IDA_mem->ida_yp0);
226 }
227 hic *= PT1;
228 IDA_mem->ida_cj = ONE/hic;
229 IDA_mem->ida_hh = hic;
230 } /* End of nh loop */
231
232 /* Break on failure; else reset ewt, save yy0, yp0 in phi, and loop. */
233 if(retval != IDA_SUCCESS) break;
234 ewtsetOK = IDA_mem->ida_efun(IDA_mem->ida_yy0, IDA_mem->ida_ewt,
235 IDA_mem->ida_edata);
236 if(ewtsetOK != 0) {
237 retval = IDA_BAD_EWT;
238 break;
239 }
240 N_VScale(ONE, IDA_mem->ida_yy0, IDA_mem->ida_phi[0]);
241 N_VScale(ONE, IDA_mem->ida_yp0, IDA_mem->ida_phi[1]);
242
243 } /* End of nwt loop */
244
245 /* Free temporary space */
246
247 N_VDestroy(IDA_mem->ida_yy0);
248 N_VDestroy(IDA_mem->ida_yp0);
249
250 /* Load the optional outputs. */
251
252 if(icopt == IDA_YA_YDP_INIT) IDA_mem->ida_hused = hic;
253
254 /* On any failure, print message and return proper flag. */
255
256 if(retval != IDA_SUCCESS) {
257 icret = IDAICFailFlag(IDA_mem, retval);
258 return(icret);
259 }
260
261 /* Otherwise return success flag. */
262
263 return(IDA_SUCCESS);
264
265 }
266
267 /*
268 * =================================================================
269 * PRIVATE FUNCTIONS IMPLEMENTATION
270 * =================================================================
271 */
272
273 /*
274 * -----------------------------------------------------------------
275 * IDAnlsIC
276 * -----------------------------------------------------------------
277 * IDAnlsIC solves a nonlinear system for consistent initial
278 * conditions. It calls IDANewtonIC to do most of the work.
279 *
280 * The return value is IDA_SUCCESS = 0 if no error occurred.
281 * The error return values (positive) considered recoverable are:
282 * IC_FAIL_RECOV if res, lsetup, or lsolve failed recoverably
283 * IC_CONSTR_FAILED if the constraints could not be met
284 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
285 * or on maxbacks test)
286 * IC_CONV_FAIL if the Newton iterations failed to converge
287 * IC_SLOW_CONVRG if the iterations are converging slowly
288 * (failed the convergence test, but showed
289 * norm reduction or convergence rate < 1)
290 * The error return values (negative) considered non-recoverable are:
291 * IDA_RES_FAIL if res had a non-recoverable error
292 * IDA_FIRST_RES_FAIL if res failed recoverably on the first call
293 * IDA_LSETUP_FAIL if lsetup had a non-recoverable error
294 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
295 * -----------------------------------------------------------------
296 */
297
IDAnlsIC(IDAMem IDA_mem)298 static int IDAnlsIC (IDAMem IDA_mem)
299 {
300 int retval, nj;
301 N_Vector tv1, tv2, tv3;
302
303 tv1 = IDA_mem->ida_ee;
304 tv2 = IDA_mem->ida_tempv2;
305 tv3 = IDA_mem->ida_phi[2];
306
307 retval = IDA_mem->ida_res(IDA_mem->ida_t0, IDA_mem->ida_yy0,
308 IDA_mem->ida_yp0, IDA_mem->ida_delta,
309 IDA_mem->ida_user_data);
310 IDA_mem->ida_nre++;
311 if(retval < 0) return(IDA_RES_FAIL);
312 if(retval > 0) return(IDA_FIRST_RES_FAIL);
313
314 N_VScale(ONE, IDA_mem->ida_delta, IDA_mem->ida_savres);
315
316 /* Loop over nj = number of linear solve Jacobian setups. */
317
318 for(nj = 1; nj <= IDA_mem->ida_maxnj; nj++) {
319
320 /* If there is a setup routine, call it. */
321 if(IDA_mem->ida_lsetup) {
322 IDA_mem->ida_nsetups++;
323 retval = IDA_mem->ida_lsetup(IDA_mem, IDA_mem->ida_yy0,
324 IDA_mem->ida_yp0, IDA_mem->ida_delta,
325 tv1, tv2, tv3);
326 if(retval < 0) return(IDA_LSETUP_FAIL);
327 if(retval > 0) return(IC_FAIL_RECOV);
328 }
329
330 /* Call the Newton iteration routine, and return if successful. */
331 retval = IDANewtonIC(IDA_mem);
332 if(retval == IDA_SUCCESS) return(IDA_SUCCESS);
333
334 /* If converging slowly and lsetup is nontrivial, retry. */
335 if(retval == IC_SLOW_CONVRG && IDA_mem->ida_lsetup) {
336 N_VScale(ONE, IDA_mem->ida_savres, IDA_mem->ida_delta);
337 continue;
338 } else {
339 return(retval);
340 }
341
342 } /* End of nj loop */
343
344 /* No convergence after maxnj tries; return with retval=IC_SLOW_CONVRG */
345 return(retval);
346
347 }
348
349 /*
350 * -----------------------------------------------------------------
351 * IDANewtonIC
352 * -----------------------------------------------------------------
353 * IDANewtonIC performs the Newton iteration to solve for consistent
354 * initial conditions. It calls IDALineSrch within each iteration.
355 * On return, savres contains the current residual vector.
356 *
357 * The return value is IDA_SUCCESS = 0 if no error occurred.
358 * The error return values (positive) considered recoverable are:
359 * IC_FAIL_RECOV if res or lsolve failed recoverably
360 * IC_CONSTR_FAILED if the constraints could not be met
361 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
362 * or on maxbacks test)
363 * IC_CONV_FAIL if the Newton iterations failed to converge
364 * IC_SLOW_CONVRG if the iterations appear to be converging slowly.
365 * They failed the convergence test, but showed
366 * an overall norm reduction (by a factor of < 0.1)
367 * or a convergence rate <= ICRATEMAX).
368 * The error return values (negative) considered non-recoverable are:
369 * IDA_RES_FAIL if res had a non-recoverable error
370 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
371 * -----------------------------------------------------------------
372 */
373
IDANewtonIC(IDAMem IDA_mem)374 static int IDANewtonIC(IDAMem IDA_mem)
375 {
376 int retval, mnewt;
377 realtype delnorm, fnorm, fnorm0, oldfnrm, rate;
378
379 /* Set pointer for vector delnew */
380 IDA_mem->ida_delnew = IDA_mem->ida_phi[2];
381
382 /* Call the linear solve function to get the Newton step, delta. */
383 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delta,
384 IDA_mem->ida_ewt, IDA_mem->ida_yy0,
385 IDA_mem->ida_yp0, IDA_mem->ida_savres);
386 if(retval < 0) return(IDA_LSOLVE_FAIL);
387 if(retval > 0) return(IC_FAIL_RECOV);
388
389 /* Compute the norm of the step; return now if this is small. */
390 fnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delta, IDA_mem->ida_ewt, SUNFALSE);
391 if(IDA_mem->ida_sysindex == 0)
392 fnorm *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
393 if(fnorm <= IDA_mem->ida_epsNewt)
394 return(IDA_SUCCESS);
395 fnorm0 = fnorm;
396
397 /* Initialize rate to avoid compiler warning message */
398 rate = ZERO;
399
400 /* Newton iteration loop */
401
402 for(mnewt = 0; mnewt < IDA_mem->ida_maxnit; mnewt++) {
403
404 IDA_mem->ida_nni++;
405 delnorm = fnorm;
406 oldfnrm = fnorm;
407
408 /* Call the Linesearch function and return if it failed. */
409 retval = IDALineSrch(IDA_mem, &delnorm, &fnorm);
410 if(retval != IDA_SUCCESS) return(retval);
411
412 /* Set the observed convergence rate and test for convergence. */
413 rate = fnorm/oldfnrm;
414 if(fnorm <= IDA_mem->ida_epsNewt) return(IDA_SUCCESS);
415
416 /* If not converged, copy new step vector, and loop. */
417 N_VScale(ONE, IDA_mem->ida_delnew, IDA_mem->ida_delta);
418
419 } /* End of Newton iteration loop */
420
421 /* Return either IC_SLOW_CONVRG or recoverable fail flag. */
422 if(rate <= ICRATEMAX || fnorm < PT1*fnorm0) return(IC_SLOW_CONVRG);
423 return(IC_CONV_FAIL);
424
425 }
426
427
428 /*
429 * -----------------------------------------------------------------
430 * IDALineSrch
431 * -----------------------------------------------------------------
432 * IDALineSrch performs the Linesearch algorithm with the
433 * calculation of consistent initial conditions.
434 *
435 * On entry, yy0 and yp0 are the current values of y and y', the
436 * Newton step is delta, the current residual vector F is savres,
437 * delnorm is WRMS-norm(delta), and fnorm is the norm of the vector
438 * J-inverse F.
439 *
440 * On a successful return, yy0, yp0, and savres have been updated,
441 * delnew contains the current value of J-inverse F, and fnorm is
442 * WRMS-norm(delnew).
443 *
444 * The return value is IDA_SUCCESS = 0 if no error occurred.
445 * The error return values (positive) considered recoverable are:
446 * IC_FAIL_RECOV if res or lsolve failed recoverably
447 * IC_CONSTR_FAILED if the constraints could not be met
448 * IC_LINESRCH_FAILED if the linesearch failed (either on steptol test
449 * or on maxbacks test)
450 * The error return values (negative) considered non-recoverable are:
451 * IDA_RES_FAIL if res had a non-recoverable error
452 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error
453 * -----------------------------------------------------------------
454 */
455
IDALineSrch(IDAMem IDA_mem,realtype * delnorm,realtype * fnorm)456 static int IDALineSrch(IDAMem IDA_mem, realtype *delnorm, realtype *fnorm)
457 {
458 booleantype conOK;
459 int retval, nbacks;
460 realtype f1norm, fnormp, f1normp, ratio, lambda, minlam, slpi;
461 N_Vector mc;
462
463 /* Initialize work space pointers, f1norm, ratio.
464 (Use of mc in constraint check does not conflict with ypnew.) */
465 mc = IDA_mem->ida_ee;
466 IDA_mem->ida_dtemp = IDA_mem->ida_phi[3];
467 IDA_mem->ida_ynew = IDA_mem->ida_tempv2;
468 IDA_mem->ida_ypnew = IDA_mem->ida_ee;
469 f1norm = (*fnorm)*(*fnorm)*HALF;
470 ratio = ONE;
471
472 /* If there are constraints, check and reduce step if necessary. */
473 if(IDA_mem->ida_constraintsSet) {
474
475 /* Update y and check constraints. */
476 IDANewy(IDA_mem);
477 conOK = N_VConstrMask(IDA_mem->ida_constraints, IDA_mem->ida_ynew, mc);
478
479 if(!conOK) {
480 /* Not satisfied. Compute scaled step to satisfy constraints. */
481 N_VProd(mc, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
482 ratio = PT99*N_VMinQuotient(IDA_mem->ida_yy0, IDA_mem->ida_dtemp);
483 (*delnorm) *= ratio;
484 if((*delnorm) <= IDA_mem->ida_steptol) return(IC_CONSTR_FAILED);
485 N_VScale(ratio, IDA_mem->ida_delta, IDA_mem->ida_delta);
486 }
487
488 } /* End of constraints check */
489
490 slpi = -TWO*f1norm*ratio;
491 minlam = IDA_mem->ida_steptol / (*delnorm);
492 lambda = ONE;
493 nbacks = 0;
494
495 /* In IDA_Y_INIT case, set ypnew = yp0 (fixed) for linesearch. */
496 if(IDA_mem->ida_icopt == IDA_Y_INIT)
497 N_VScale(ONE, IDA_mem->ida_yp0, IDA_mem->ida_ypnew);
498
499 /* Loop on linesearch variable lambda. */
500
501 for(;;) {
502
503 if (nbacks == IDA_mem->ida_maxbacks) return(IC_LINESRCH_FAILED);
504 /* Get new (y,y') = (ynew,ypnew) and norm of new function value. */
505 IDANewyyp(IDA_mem, lambda);
506 retval = IDAfnorm(IDA_mem, &fnormp);
507 if(retval != IDA_SUCCESS) return(retval);
508
509 /* If lsoff option is on, break out. */
510 if(IDA_mem->ida_lsoff) break;
511
512 /* Do alpha-condition test. */
513 f1normp = fnormp*fnormp*HALF;
514 if(f1normp <= f1norm + ALPHALS*slpi*lambda) break;
515 if(lambda < minlam) return(IC_LINESRCH_FAILED);
516 lambda /= TWO;
517 IDA_mem->ida_nbacktr++; nbacks++;
518
519 } /* End of breakout linesearch loop */
520
521 /* Update yy0, yp0, and fnorm, then return. */
522 N_VScale(ONE, IDA_mem->ida_ynew, IDA_mem->ida_yy0);
523 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT)
524 N_VScale(ONE, IDA_mem->ida_ypnew, IDA_mem->ida_yp0);
525 *fnorm = fnormp;
526 return(IDA_SUCCESS);
527
528 }
529
530 /*
531 * -----------------------------------------------------------------
532 * IDAfnorm
533 * -----------------------------------------------------------------
534 * IDAfnorm computes the norm of the current function value, by
535 * evaluating the DAE residual function, calling the linear
536 * system solver, and computing a WRMS-norm.
537 *
538 * On return, savres contains the current residual vector F, and
539 * delnew contains J-inverse F.
540 *
541 * The return value is IDA_SUCCESS = 0 if no error occurred, or
542 * IC_FAIL_RECOV if res or lsolve failed recoverably, or
543 * IDA_RES_FAIL if res had a non-recoverable error, or
544 * IDA_LSOLVE_FAIL if lsolve had a non-recoverable error.
545 * -----------------------------------------------------------------
546 */
547
IDAfnorm(IDAMem IDA_mem,realtype * fnorm)548 static int IDAfnorm(IDAMem IDA_mem, realtype *fnorm)
549 {
550
551 int retval;
552
553 /* Get residual vector F, return if failed, and save F in savres. */
554 retval = IDA_mem->ida_res(IDA_mem->ida_t0, IDA_mem->ida_ynew,
555 IDA_mem->ida_ypnew, IDA_mem->ida_delnew,
556 IDA_mem->ida_user_data);
557 IDA_mem->ida_nre++;
558 if(retval < 0) return(IDA_RES_FAIL);
559 if(retval > 0) return(IC_FAIL_RECOV);
560
561 N_VScale(ONE, IDA_mem->ida_delnew, IDA_mem->ida_savres);
562
563 /* Call the linear solve function to get J-inverse F; return if failed. */
564 retval = IDA_mem->ida_lsolve(IDA_mem, IDA_mem->ida_delnew,
565 IDA_mem->ida_ewt, IDA_mem->ida_ynew,
566 IDA_mem->ida_ypnew, IDA_mem->ida_savres);
567 if(retval < 0) return(IDA_LSOLVE_FAIL);
568 if(retval > 0) return(IC_FAIL_RECOV);
569
570 /* Compute the WRMS-norm; rescale if index = 0. */
571 *fnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delnew, IDA_mem->ida_ewt, SUNFALSE);
572 if(IDA_mem->ida_sysindex == 0)
573 (*fnorm) *= IDA_mem->ida_tscale * SUNRabs(IDA_mem->ida_cj);
574
575 return(IDA_SUCCESS);
576
577 }
578
579 /*
580 * -----------------------------------------------------------------
581 * IDANewyyp
582 * -----------------------------------------------------------------
583 * IDANewyyp updates the vectors ynew and ypnew from yy0 and yp0,
584 * using the current step vector lambda*delta, in a manner
585 * depending on icopt and the input id vector.
586 *
587 * The return value is always IDA_SUCCESS = 0.
588 * -----------------------------------------------------------------
589 */
590
IDANewyyp(IDAMem IDA_mem,realtype lambda)591 static int IDANewyyp(IDAMem IDA_mem, realtype lambda)
592 {
593
594 /* IDA_YA_YDP_INIT case: ynew = yy0 - lambda*delta where id_i = 0
595 ypnew = yp0 - cj*lambda*delta where id_i = 1. */
596 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
597 N_VProd(IDA_mem->ida_id, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
598 N_VLinearSum(ONE, IDA_mem->ida_yp0, -IDA_mem->ida_cj*lambda,
599 IDA_mem->ida_dtemp, IDA_mem->ida_ypnew);
600 N_VLinearSum(ONE, IDA_mem->ida_delta, -ONE,
601 IDA_mem->ida_dtemp, IDA_mem->ida_dtemp);
602 N_VLinearSum(ONE, IDA_mem->ida_yy0, -lambda,
603 IDA_mem->ida_dtemp, IDA_mem->ida_ynew);
604 return(IDA_SUCCESS);
605 }
606
607 /* IDA_Y_INIT case: ynew = yy0 - lambda*delta. (ypnew = yp0 preset.) */
608 N_VLinearSum(ONE, IDA_mem->ida_yy0, -lambda,
609 IDA_mem->ida_delta, IDA_mem->ida_ynew);
610 return(IDA_SUCCESS);
611
612 }
613
614 /*
615 * -----------------------------------------------------------------
616 * IDANewy
617 * -----------------------------------------------------------------
618 * IDANewy updates the vector ynew from yy0,
619 * using the current step vector delta, in a manner
620 * depending on icopt and the input id vector.
621 *
622 * The return value is always IDA_SUCCESS = 0.
623 * -----------------------------------------------------------------
624 */
625
IDANewy(IDAMem IDA_mem)626 static int IDANewy(IDAMem IDA_mem)
627 {
628
629 /* IDA_YA_YDP_INIT case: ynew = yy0 - delta where id_i = 0. */
630 if(IDA_mem->ida_icopt == IDA_YA_YDP_INIT) {
631 N_VProd(IDA_mem->ida_id, IDA_mem->ida_delta, IDA_mem->ida_dtemp);
632 N_VLinearSum(ONE, IDA_mem->ida_delta, -ONE,
633 IDA_mem->ida_dtemp, IDA_mem->ida_dtemp);
634 N_VLinearSum(ONE, IDA_mem->ida_yy0, -ONE,
635 IDA_mem->ida_dtemp, IDA_mem->ida_ynew);
636 return(IDA_SUCCESS);
637 }
638
639 /* IDA_Y_INIT case: ynew = yy0 - delta. */
640 N_VLinearSum(ONE, IDA_mem->ida_yy0, -ONE,
641 IDA_mem->ida_delta, IDA_mem->ida_ynew);
642 return(IDA_SUCCESS);
643
644 }
645
646 /*
647 * -----------------------------------------------------------------
648 * IDAICFailFlag
649 * -----------------------------------------------------------------
650 * IDAICFailFlag prints a message and sets the IDACalcIC return
651 * value appropriate to the flag retval returned by IDAnlsIC.
652 * -----------------------------------------------------------------
653 */
654
IDAICFailFlag(IDAMem IDA_mem,int retval)655 static int IDAICFailFlag(IDAMem IDA_mem, int retval)
656 {
657
658 /* Depending on retval, print error message and return error flag. */
659 switch(retval) {
660
661 case IDA_RES_FAIL:
662 IDAProcessError(IDA_mem, IDA_RES_FAIL, "IDA", "IDACalcIC", MSG_IC_RES_NONREC);
663 return(IDA_RES_FAIL);
664
665 case IDA_FIRST_RES_FAIL:
666 IDAProcessError(IDA_mem, IDA_FIRST_RES_FAIL, "IDA", "IDACalcIC", MSG_IC_RES_FAIL);
667 return(IDA_FIRST_RES_FAIL);
668
669 case IDA_LSETUP_FAIL:
670 IDAProcessError(IDA_mem, IDA_LSETUP_FAIL, "IDA", "IDACalcIC", MSG_IC_SETUP_FAIL);
671 return(IDA_LSETUP_FAIL);
672
673 case IDA_LSOLVE_FAIL:
674 IDAProcessError(IDA_mem, IDA_LSOLVE_FAIL, "IDA", "IDACalcIC", MSG_IC_SOLVE_FAIL);
675 return(IDA_LSOLVE_FAIL);
676
677 case IC_FAIL_RECOV:
678 IDAProcessError(IDA_mem, IDA_NO_RECOVERY, "IDA", "IDACalcIC", MSG_IC_NO_RECOVERY);
679 return(IDA_NO_RECOVERY);
680
681 case IC_CONSTR_FAILED:
682 IDAProcessError(IDA_mem, IDA_CONSTR_FAIL, "IDA", "IDACalcIC", MSG_IC_FAIL_CONSTR);
683 return(IDA_CONSTR_FAIL);
684
685 case IC_LINESRCH_FAILED:
686 IDAProcessError(IDA_mem, IDA_LINESEARCH_FAIL, "IDA", "IDACalcIC", MSG_IC_FAILED_LINS);
687 return(IDA_LINESEARCH_FAIL);
688
689 case IC_CONV_FAIL:
690 IDAProcessError(IDA_mem, IDA_CONV_FAIL, "IDA", "IDACalcIC", MSG_IC_CONV_FAILED);
691 return(IDA_CONV_FAIL);
692
693 case IC_SLOW_CONVRG:
694 IDAProcessError(IDA_mem, IDA_CONV_FAIL, "IDA", "IDACalcIC", MSG_IC_CONV_FAILED);
695 return(IDA_CONV_FAIL);
696
697 case IDA_BAD_EWT:
698 IDAProcessError(IDA_mem, IDA_BAD_EWT, "IDA", "IDACalcIC", MSG_IC_BAD_EWT);
699 return(IDA_BAD_EWT);
700
701 }
702 return -99;
703 }
704
705