1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Radu Serban @ LLNL
4 * -----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, Lawrence Livermore National Security
7 * and Southern Methodist University.
8 * All rights reserved.
9 *
10 * See the top-level LICENSE and NOTICE files for details.
11 *
12 * SPDX-License-Identifier: BSD-3-Clause
13 * SUNDIALS Copyright End
14 * -----------------------------------------------------------------
15 * This is the implementation file for the CVDIAG linear solver.
16 * -----------------------------------------------------------------
17 */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21
22 #include "cvodes_diag_impl.h"
23 #include "cvodes_impl.h"
24
25 /* Other Constants */
26
27 #define FRACT RCONST(0.1)
28 #define ONE RCONST(1.0)
29
30 /* CVDIAG linit, lsetup, lsolve, and lfree routines */
31
32 static int CVDiagInit(CVodeMem cv_mem);
33
34 static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
35 N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
36 N_Vector vtemp2, N_Vector vtemp3);
37
38 static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
39 N_Vector ycur, N_Vector fcur);
40
41 static int CVDiagFree(CVodeMem cv_mem);
42
43
44 /*
45 * ================================================================
46 *
47 * PART I - forward problems
48 *
49 * ================================================================
50 */
51
52
53 /* Readability Replacements */
54
55 #define lrw1 (cv_mem->cv_lrw1)
56 #define liw1 (cv_mem->cv_liw1)
57 #define f (cv_mem->cv_f)
58 #define uround (cv_mem->cv_uround)
59 #define tn (cv_mem->cv_tn)
60 #define h (cv_mem->cv_h)
61 #define rl1 (cv_mem->cv_rl1)
62 #define gamma (cv_mem->cv_gamma)
63 #define ewt (cv_mem->cv_ewt)
64 #define nfe (cv_mem->cv_nfe)
65 #define zn (cv_mem->cv_zn)
66 #define linit (cv_mem->cv_linit)
67 #define lsetup (cv_mem->cv_lsetup)
68 #define lsolve (cv_mem->cv_lsolve)
69 #define lfree (cv_mem->cv_lfree)
70 #define lmem (cv_mem->cv_lmem)
71 #define vec_tmpl (cv_mem->cv_tempv)
72 #define setupNonNull (cv_mem->cv_setupNonNull)
73
74 #define gammasv (cvdiag_mem->di_gammasv)
75 #define M (cvdiag_mem->di_M)
76 #define bit (cvdiag_mem->di_bit)
77 #define bitcomp (cvdiag_mem->di_bitcomp)
78 #define nfeDI (cvdiag_mem->di_nfeDI)
79 #define last_flag (cvdiag_mem->di_last_flag)
80
81
82 /*
83 * -----------------------------------------------------------------
84 * CVDiag
85 * -----------------------------------------------------------------
86 * This routine initializes the memory record and sets various function
87 * fields specific to the diagonal linear solver module. CVDense first
88 * calls the existing lfree routine if this is not NULL. Then it sets
89 * the cv_linit, cv_lsetup, cv_lsolve, cv_lfree fields in (*cvode_mem)
90 * to be CVDiagInit, CVDiagSetup, CVDiagSolve, and CVDiagFree,
91 * respectively. It allocates memory for a structure of type
92 * CVDiagMemRec and sets the cv_lmem field in (*cvode_mem) to the
93 * address of this structure. It sets setupNonNull in (*cvode_mem) to
94 * SUNTRUE. Finally, it allocates memory for M, bit, and bitcomp.
95 * The CVDiag return value is SUCCESS = 0, LMEM_FAIL = -1, or
96 * LIN_ILL_INPUT=-2.
97 * -----------------------------------------------------------------
98 */
99
CVDiag(void * cvode_mem)100 int CVDiag(void *cvode_mem)
101 {
102 CVodeMem cv_mem;
103 CVDiagMem cvdiag_mem;
104
105 /* Return immediately if cvode_mem is NULL */
106 if (cvode_mem == NULL) {
107 cvProcessError(NULL, CVDIAG_MEM_NULL, "CVDIAG", "CVDiag", MSGDG_CVMEM_NULL);
108 return(CVDIAG_MEM_NULL);
109 }
110 cv_mem = (CVodeMem) cvode_mem;
111
112 /* Check if N_VCompare and N_VInvTest are present */
113 if(vec_tmpl->ops->nvcompare == NULL ||
114 vec_tmpl->ops->nvinvtest == NULL) {
115 cvProcessError(cv_mem, CVDIAG_ILL_INPUT, "CVDIAG", "CVDiag", MSGDG_BAD_NVECTOR);
116 return(CVDIAG_ILL_INPUT);
117 }
118
119 if (lfree != NULL) lfree(cv_mem);
120
121 /* Set four main function fields in cv_mem */
122 linit = CVDiagInit;
123 lsetup = CVDiagSetup;
124 lsolve = CVDiagSolve;
125 lfree = CVDiagFree;
126
127 /* Get memory for CVDiagMemRec */
128 cvdiag_mem = NULL;
129 cvdiag_mem = (CVDiagMem) malloc(sizeof(CVDiagMemRec));
130 if (cvdiag_mem == NULL) {
131 cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
132 return(CVDIAG_MEM_FAIL);
133 }
134
135 last_flag = CVDIAG_SUCCESS;
136
137 /* Allocate memory for M, bit, and bitcomp */
138
139 M = N_VClone(vec_tmpl);
140 if (M == NULL) {
141 cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
142 free(cvdiag_mem); cvdiag_mem = NULL;
143 return(CVDIAG_MEM_FAIL);
144 }
145 bit = N_VClone(vec_tmpl);
146 if (bit == NULL) {
147 cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
148 N_VDestroy(M);
149 free(cvdiag_mem); cvdiag_mem = NULL;
150 return(CVDIAG_MEM_FAIL);
151 }
152 bitcomp = N_VClone(vec_tmpl);
153 if (bitcomp == NULL) {
154 cvProcessError(cv_mem, CVDIAG_MEM_FAIL, "CVDIAG", "CVDiag", MSGDG_MEM_FAIL);
155 N_VDestroy(M);
156 N_VDestroy(bit);
157 free(cvdiag_mem); cvdiag_mem = NULL;
158 return(CVDIAG_MEM_FAIL);
159 }
160
161 /* Attach linear solver memory to integrator memory */
162 lmem = cvdiag_mem;
163
164 return(CVDIAG_SUCCESS);
165 }
166
167 /*
168 * -----------------------------------------------------------------
169 * CVDiagGetWorkSpace
170 * -----------------------------------------------------------------
171 */
172
CVDiagGetWorkSpace(void * cvode_mem,long int * lenrwLS,long int * leniwLS)173 int CVDiagGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS)
174 {
175 CVodeMem cv_mem;
176
177 /* Return immediately if cvode_mem is NULL */
178 if (cvode_mem == NULL) {
179 cvProcessError(NULL, CVDIAG_MEM_NULL, "CVDIAG", "CVDiagGetWorkSpace", MSGDG_CVMEM_NULL);
180 return(CVDIAG_MEM_NULL);
181 }
182 cv_mem = (CVodeMem) cvode_mem;
183
184 *lenrwLS = 3*lrw1;
185 *leniwLS = 3*liw1;
186
187 return(CVDIAG_SUCCESS);
188 }
189
190 /*
191 * -----------------------------------------------------------------
192 * CVDiagGetNumRhsEvals
193 * -----------------------------------------------------------------
194 */
195
CVDiagGetNumRhsEvals(void * cvode_mem,long int * nfevalsLS)196 int CVDiagGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS)
197 {
198 CVodeMem cv_mem;
199 CVDiagMem cvdiag_mem;
200
201 /* Return immediately if cvode_mem is NULL */
202 if (cvode_mem == NULL) {
203 cvProcessError(NULL, CVDIAG_MEM_NULL, "CVDIAG", "CVDiagGetNumRhsEvals", MSGDG_CVMEM_NULL);
204 return(CVDIAG_MEM_NULL);
205 }
206 cv_mem = (CVodeMem) cvode_mem;
207
208 if (lmem == NULL) {
209 cvProcessError(cv_mem, CVDIAG_LMEM_NULL, "CVDIAG", "CVDiagGetNumRhsEvals", MSGDG_LMEM_NULL);
210 return(CVDIAG_LMEM_NULL);
211 }
212 cvdiag_mem = (CVDiagMem) lmem;
213
214 *nfevalsLS = nfeDI;
215
216 return(CVDIAG_SUCCESS);
217 }
218
219 /*
220 * -----------------------------------------------------------------
221 * CVDiagGetLastFlag
222 * -----------------------------------------------------------------
223 */
224
CVDiagGetLastFlag(void * cvode_mem,long int * flag)225 int CVDiagGetLastFlag(void *cvode_mem, long int *flag)
226 {
227 CVodeMem cv_mem;
228 CVDiagMem cvdiag_mem;
229
230 /* Return immediately if cvode_mem is NULL */
231 if (cvode_mem == NULL) {
232 cvProcessError(NULL, CVDIAG_MEM_NULL, "CVDIAG", "CVDiagGetLastFlag", MSGDG_CVMEM_NULL);
233 return(CVDIAG_MEM_NULL);
234 }
235 cv_mem = (CVodeMem) cvode_mem;
236
237 if (lmem == NULL) {
238 cvProcessError(cv_mem, CVDIAG_LMEM_NULL, "CVDIAG", "CVDiagGetLastFlag", MSGDG_LMEM_NULL);
239 return(CVDIAG_LMEM_NULL);
240 }
241 cvdiag_mem = (CVDiagMem) lmem;
242
243 *flag = last_flag;
244
245 return(CVDIAG_SUCCESS);
246 }
247
248 /*
249 * -----------------------------------------------------------------
250 * CVDiagGetReturnFlagName
251 * -----------------------------------------------------------------
252 */
253
CVDiagGetReturnFlagName(long int flag)254 char *CVDiagGetReturnFlagName(long int flag)
255 {
256 char *name;
257
258 name = (char *)malloc(30*sizeof(char));
259
260 switch(flag) {
261 case CVDIAG_SUCCESS:
262 sprintf(name,"CVDIAG_SUCCESS");
263 break;
264 case CVDIAG_MEM_NULL:
265 sprintf(name,"CVDIAG_MEM_NULL");
266 break;
267 case CVDIAG_LMEM_NULL:
268 sprintf(name,"CVDIAG_LMEM_NULL");
269 break;
270 case CVDIAG_ILL_INPUT:
271 sprintf(name,"CVDIAG_ILL_INPUT");
272 break;
273 case CVDIAG_MEM_FAIL:
274 sprintf(name,"CVDIAG_MEM_FAIL");
275 break;
276 case CVDIAG_INV_FAIL:
277 sprintf(name,"CVDIAG_INV_FAIL");
278 break;
279 case CVDIAG_RHSFUNC_UNRECVR:
280 sprintf(name,"CVDIAG_RHSFUNC_UNRECVR");
281 break;
282 case CVDIAG_RHSFUNC_RECVR:
283 sprintf(name,"CVDIAG_RHSFUNC_RECVR");
284 break;
285 case CVDIAG_NO_ADJ:
286 sprintf(name,"CVDIAG_NO_ADJ");
287 break;
288 default:
289 sprintf(name,"NONE");
290 }
291
292 return(name);
293 }
294
295 /*
296 * -----------------------------------------------------------------
297 * CVDiagInit
298 * -----------------------------------------------------------------
299 * This routine does remaining initializations specific to the diagonal
300 * linear solver.
301 * -----------------------------------------------------------------
302 */
303
CVDiagInit(CVodeMem cv_mem)304 static int CVDiagInit(CVodeMem cv_mem)
305 {
306 CVDiagMem cvdiag_mem;
307
308 cvdiag_mem = (CVDiagMem) lmem;
309
310 nfeDI = 0;
311
312 last_flag = CVDIAG_SUCCESS;
313 return(0);
314 }
315
316 /*
317 * -----------------------------------------------------------------
318 * CVDiagSetup
319 * -----------------------------------------------------------------
320 * This routine does the setup operations for the diagonal linear
321 * solver. It constructs a diagonal approximation to the Newton matrix
322 * M = I - gamma*J, updates counters, and inverts M.
323 * -----------------------------------------------------------------
324 */
325
CVDiagSetup(CVodeMem cv_mem,int convfail,N_Vector ypred,N_Vector fpred,booleantype * jcurPtr,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)326 static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
327 N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
328 N_Vector vtemp2, N_Vector vtemp3)
329 {
330 realtype r;
331 N_Vector ftemp, y;
332 booleantype invOK;
333 CVDiagMem cvdiag_mem;
334 int retval;
335
336 cvdiag_mem = (CVDiagMem) lmem;
337
338 /* Rename work vectors for use as temporary values of y and f */
339 ftemp = vtemp1;
340 y = vtemp2;
341
342 /* Form y with perturbation = FRACT*(func. iter. correction) */
343 r = FRACT * rl1;
344 N_VLinearSum(h, fpred, -ONE, zn[1], ftemp);
345 N_VLinearSum(r, ftemp, ONE, ypred, y);
346
347 /* Evaluate f at perturbed y */
348 retval = f(tn, y, M, cv_mem->cv_user_data);
349 nfeDI++;
350 if (retval < 0) {
351 cvProcessError(cv_mem, CVDIAG_RHSFUNC_UNRECVR, "CVDIAG", "CVDiagSetup", MSGDG_RHSFUNC_FAILED);
352 last_flag = CVDIAG_RHSFUNC_UNRECVR;
353 return(-1);
354 }
355 if (retval > 0) {
356 last_flag = CVDIAG_RHSFUNC_RECVR;
357 return(1);
358 }
359
360 /* Construct M = I - gamma*J with J = diag(deltaf_i/deltay_i) */
361 N_VLinearSum(ONE, M, -ONE, fpred, M);
362 N_VLinearSum(FRACT, ftemp, -h, M, M);
363 N_VProd(ftemp, ewt, y);
364 /* Protect against deltay_i being at roundoff level */
365 N_VCompare(uround, y, bit);
366 N_VAddConst(bit, -ONE, bitcomp);
367 N_VProd(ftemp, bit, y);
368 N_VLinearSum(FRACT, y, -ONE, bitcomp, y);
369 N_VDiv(M, y, M);
370 N_VProd(M, bit, M);
371 N_VLinearSum(ONE, M, -ONE, bitcomp, M);
372
373 /* Invert M with test for zero components */
374 invOK = N_VInvTest(M, M);
375 if (!invOK) {
376 last_flag = CVDIAG_INV_FAIL;
377 return(1);
378 }
379
380 /* Set jcur = SUNTRUE, save gamma in gammasv, and return */
381 *jcurPtr = SUNTRUE;
382 gammasv = gamma;
383 last_flag = CVDIAG_SUCCESS;
384 return(0);
385 }
386
387 /*
388 * -----------------------------------------------------------------
389 * CVDiagSolve
390 * -----------------------------------------------------------------
391 * This routine performs the solve operation for the diagonal linear
392 * solver. If necessary it first updates gamma in M = I - gamma*J.
393 * -----------------------------------------------------------------
394 */
395
CVDiagSolve(CVodeMem cv_mem,N_Vector b,N_Vector weight,N_Vector ycur,N_Vector fcur)396 static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
397 N_Vector ycur, N_Vector fcur)
398 {
399 booleantype invOK;
400 realtype r;
401 CVDiagMem cvdiag_mem;
402
403 cvdiag_mem = (CVDiagMem) lmem;
404
405 /* If gamma has changed, update factor in M, and save gamma value */
406
407 if (gammasv != gamma) {
408 r = gamma / gammasv;
409 N_VInv(M, M);
410 N_VAddConst(M, -ONE, M);
411 N_VScale(r, M, M);
412 N_VAddConst(M, ONE, M);
413 invOK = N_VInvTest(M, M);
414 if (!invOK) {
415 last_flag = CVDIAG_INV_FAIL;
416 return (1);
417 }
418 gammasv = gamma;
419 }
420
421 /* Apply M-inverse to b */
422 N_VProd(b, M, b);
423
424 last_flag = CVDIAG_SUCCESS;
425 return(0);
426 }
427
428 /*
429 * -----------------------------------------------------------------
430 * CVDiagFree
431 * -----------------------------------------------------------------
432 * This routine frees memory specific to the diagonal linear solver.
433 * -----------------------------------------------------------------
434 */
435
CVDiagFree(CVodeMem cv_mem)436 static int CVDiagFree(CVodeMem cv_mem)
437 {
438 CVDiagMem cvdiag_mem;
439
440 cvdiag_mem = (CVDiagMem) lmem;
441
442 N_VDestroy(M);
443 N_VDestroy(bit);
444 N_VDestroy(bitcomp);
445 free(cvdiag_mem);
446 cv_mem->cv_lmem = NULL;
447
448 return(0);
449 }
450
451
452 /*
453 * ================================================================
454 *
455 * PART II - backward problems
456 *
457 * ================================================================
458 */
459
460
461 /*
462 * CVDiagB
463 *
464 * Wrappers for the backward phase around the corresponding
465 * CVODES functions
466 */
467
CVDiagB(void * cvode_mem,int which)468 int CVDiagB(void *cvode_mem, int which)
469 {
470 CVodeMem cv_mem;
471 CVadjMem ca_mem;
472 CVodeBMem cvB_mem;
473 void *cvodeB_mem;
474 int flag;
475
476 /* Check if cvode_mem exists */
477 if (cvode_mem == NULL) {
478 cvProcessError(NULL, CVDIAG_MEM_NULL, "CVSDIAG", "CVDiagB", MSGDG_CVMEM_NULL);
479 return(CVDIAG_MEM_NULL);
480 }
481 cv_mem = (CVodeMem) cvode_mem;
482
483 /* Was ASA initialized? */
484 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
485 cvProcessError(cv_mem, CVDIAG_NO_ADJ, "CVSDIAG", "CVDiagB", MSGDG_NO_ADJ);
486 return(CVDIAG_NO_ADJ);
487 }
488 ca_mem = cv_mem->cv_adj_mem;
489
490 /* Check which */
491 if ( which >= ca_mem->ca_nbckpbs ) {
492 cvProcessError(cv_mem, CVDIAG_ILL_INPUT, "CVSDIAG", "CVDiagB", MSGDG_BAD_WHICH);
493 return(CVDIAG_ILL_INPUT);
494 }
495
496 /* Find the CVodeBMem entry in the linked list corresponding to which */
497 cvB_mem = ca_mem->cvB_mem;
498 while (cvB_mem != NULL) {
499 if ( which == cvB_mem->cv_index ) break;
500 cvB_mem = cvB_mem->cv_next;
501 }
502
503 cvodeB_mem = (void *) (cvB_mem->cv_mem);
504
505 flag = CVDiag(cvodeB_mem);
506
507 return(flag);
508 }
509
510