Skip to content

Commit

Permalink
switch to using Newton polynomial
Browse files Browse the repository at this point in the history
  • Loading branch information
@gardner48
gardner48 committed Jul 6, 2023
1 parent 205f885 commit 35951d3
Show file tree
Hide file tree
Showing 3 changed files with 90 additions and 34 deletions.
119 changes: 86 additions & 33 deletions src/cvode/cvode.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -959,9 +959,41 @@ int CVodeResizeHistory(void *cvode_mem, sunrealtype* t_hist, N_Vector* y_hist,
N_VScale(ONE, cv_mem->cv_tempv, cv_mem->cv_zn[cv_mem->cv_qmax]);
}

/* ------------------------- *
* Construct Nordsieck Array *
* ------------------------- */
/* ------------------------------- *
* Construct Nordsieck Array (new) *
* ------------------------------- */

/* Compute interpolation coefficients */
retval = NewtonPolyCoef(t_hist, y_hist, n_hist, cv_mem->resize_wrk);
if (retval)
{
cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeResizeHistory",
"NewtonPolyCoef failed");
return CV_ILL_INPUT;
}

retval = NewtonPolyMultiDerEval(t_hist, cv_mem->resize_wrk, n_hist, t_hist[0],
cv_mem->cv_qprime + 1, cv_mem->cv_zn);
if (retval)
{
cvProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeResizeHistory",
"NewtonPolyMultiDerEval failed");
return CV_ILL_INPUT;
}

sunrealtype scale[L_MAX];
scale[0] = ONE;
for (int i = 1; i < L_MAX; i++)
{
scale[i] *= cv_mem->cv_hscale / ((sunrealtype) i);
}

/* N_VScale((cv_mem->cv_hscale * cv_mem->cv_hscale) / TWO, */
/* cv_mem->cv_zn[2], cv_mem->cv_zn[2]); */

/* ------------------------------- *
* Construct Nordsieck Array (old) *
* ------------------------------- */

/* Test 0: Input y_hist is actually a copy of zn. Output should match that
from the problem without resizing except for the duplicated values. */
Expand Down Expand Up @@ -989,20 +1021,20 @@ int CVodeResizeHistory(void *cvode_mem, sunrealtype* t_hist, N_Vector* y_hist,

/* Test 3: Copy yn, evaluate fn, and limit to first order (externally) */

/* zn[0] = y_n-1 */
N_VScale(ONE, y_hist[0], cv_mem->cv_zn[0]);
/* /\* zn[0] = y_n-1 *\/ */
/* N_VScale(ONE, y_hist[0], cv_mem->cv_zn[0]); */

/* zn[1] = h_n-1 * y'(t_n-1, y_n-1) */
retval = cv_mem->cv_f(t_hist[0], y_hist[0],
cv_mem->cv_zn[1], cv_mem->cv_user_data);
cv_mem->cv_nfe++;
if (retval)
{
cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODE", "CVode",
MSGCV_RHSFUNC_FAILED, cv_mem->cv_tn);
return CV_RHSFUNC_FAIL;
}
N_VScale(cv_mem->cv_hscale, cv_mem->cv_zn[1], cv_mem->cv_zn[1]);
/* /\* zn[1] = h_n-1 * y'(t_n-1, y_n-1) *\/ */
/* retval = cv_mem->cv_f(t_hist[0], y_hist[0], */
/* cv_mem->cv_zn[1], cv_mem->cv_user_data); */
/* cv_mem->cv_nfe++; */
/* if (retval) */
/* { */
/* cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODE", "CVode", */
/* MSGCV_RHSFUNC_FAILED, cv_mem->cv_tn); */
/* return CV_RHSFUNC_FAIL; */
/* } */
/* N_VScale(cv_mem->cv_hscale, cv_mem->cv_zn[1], cv_mem->cv_zn[1]); */

/* for (int j = 2; j < maxord; j++) */
/* { */
Expand All @@ -1011,25 +1043,25 @@ int CVodeResizeHistory(void *cvode_mem, sunrealtype* t_hist, N_Vector* y_hist,

/* Test 3: Copy yn and evaluate fn (like above). Compute higher order
derivatives from polynomial interpolant */
sunrealtype a[4];
/* sunrealtype a[4]; */

/* (>= 2nd order) zn[2] = ((h_n-1)^2 / 2) * y''(t_n-1, y_n-1) */
if (cv_mem->cv_qprime >= 2)
{
for (int j = 0; j < 3; j++) a[j] = LBasisD2(j, t_hist[0], t_hist);
N_VLinearCombination(3, a, y_hist, cv_mem->cv_zn[2]);
N_VScale((cv_mem->cv_hscale * cv_mem->cv_hscale) / TWO,
cv_mem->cv_zn[2], cv_mem->cv_zn[2]);
}
/* /\* (>= 2nd order) zn[2] = ((h_n-1)^2 / 2) * y''(t_n-1, y_n-1) *\/ */
/* if (cv_mem->cv_qprime >= 2) */
/* { */
/* for (int j = 0; j < 3; j++) a[j] = LBasisD2(j, t_hist[0], t_hist); */
/* N_VLinearCombination(3, a, y_hist, cv_mem->cv_zn[2]); */
/* N_VScale((cv_mem->cv_hscale * cv_mem->cv_hscale) / TWO, */
/* cv_mem->cv_zn[2], cv_mem->cv_zn[2]); */
/* } */

/* (>= 3rd order) zn[3] = ((h_n-1)^3 / 6) * y'''(t_n-1, y_n-1) */
if (cv_mem->cv_qprime >= 3)
{
for (int j = 0; j < 4; j++) a[j] = LBasisD3(j, t_hist[0], t_hist);
N_VLinearCombination(4, a, y_hist, cv_mem->cv_zn[3]);
N_VScale((cv_mem->cv_hscale * cv_mem->cv_hscale * cv_mem->cv_hscale) / SUN_RCONST(6.0),
cv_mem->cv_zn[3], cv_mem->cv_zn[3]);
}
/* /\* (>= 3rd order) zn[3] = ((h_n-1)^3 / 6) * y'''(t_n-1, y_n-1) *\/ */
/* if (cv_mem->cv_qprime >= 3) */
/* { */
/* for (int j = 0; j < 4; j++) a[j] = LBasisD3(j, t_hist[0], t_hist); */
/* N_VLinearCombination(4, a, y_hist, cv_mem->cv_zn[3]); */
/* N_VScale((cv_mem->cv_hscale * cv_mem->cv_hscale * cv_mem->cv_hscale) / SUN_RCONST(6.0), */
/* cv_mem->cv_zn[3], cv_mem->cv_zn[3]); */
/* } */

/* >>>
Do we need to adjust q and qprime and not apply order change update?
Expand Down Expand Up @@ -2171,6 +2203,26 @@ static booleantype cvAllocVectors(CVodeMem cv_mem, N_Vector tmpl)
}
}

for (j=0; j <= cv_mem->cv_qmax; j++)
{
cv_mem->resize_wrk[j] = N_VClone(tmpl);
N_VConst(NAN, cv_mem->resize_wrk[j]);
if (cv_mem->resize_wrk[j] == NULL)
{
N_VDestroy(cv_mem->cv_ewt);
N_VDestroy(cv_mem->cv_acor);
N_VDestroy(cv_mem->cv_tempv);
N_VDestroy(cv_mem->cv_ftemp);
N_VDestroy(cv_mem->cv_vtemp1);
N_VDestroy(cv_mem->cv_vtemp2);
N_VDestroy(cv_mem->cv_vtemp3);
for (i=0; i <= cv_mem->cv_qmax; i++) N_VDestroy(cv_mem->cv_zn[i]);
for (i=0; i < j; i++) N_VDestroy(cv_mem->resize_wrk[i]);
return(SUNFALSE);
}
}


/* Update solver workspace lengths */
cv_mem->cv_lrw += (cv_mem->cv_qmax + 8)*cv_mem->cv_lrw1;
cv_mem->cv_liw += (cv_mem->cv_qmax + 8)*cv_mem->cv_liw1;
Expand Down Expand Up @@ -2201,6 +2253,7 @@ static void cvFreeVectors(CVodeMem cv_mem)
N_VDestroy(cv_mem->cv_vtemp2);
N_VDestroy(cv_mem->cv_vtemp3);
for (j=0; j <= maxord; j++) N_VDestroy(cv_mem->cv_zn[j]);
for (j=0; j <= maxord; j++) N_VDestroy(cv_mem->resize_wrk[j]);

cv_mem->cv_lrw -= (maxord + 8)*cv_mem->cv_lrw1;
cv_mem->cv_liw -= (maxord + 8)*cv_mem->cv_liw1;
Expand Down
3 changes: 3 additions & 0 deletions src/cvode/cvode_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -474,6 +474,9 @@ typedef struct CVodeMemRec {

booleantype cv_usefused; /* flag indicating if CVODE specific fused kernels should be used */

/* Resize history workspace */
N_Vector resize_wrk[L_MAX];

} *CVodeMem;

/*
Expand Down
2 changes: 1 addition & 1 deletion test/unit_tests/cvode/CXX_serial/cv_test_resize_history.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -255,7 +255,7 @@ int main(int argc, char* argv[])
// 14 steps - reach 3rd order
// 22 steps - reach 4th order
// 27 steps - reach 5th order
for (int i = 1; i <= 22; i++)
for (int i = 1; i <= 30; i++)
{
flag = CVode(cvode_mem, tf, y, &(tret[i]), CV_ONE_STEP);
if (check_flag(flag, "CVode")) return 1;
Expand Down

0 comments on commit 35951d3

Please sign in to comment.