From 7cd161f995e6c5fe72386c309009dde360e38900 Mon Sep 17 00:00:00 2001 From: maggul Date: Thu, 19 Sep 2024 18:35:38 -0500 Subject: [PATCH] revision --- doc/arkode/guide/source/Usage/index.rst | 2 +- .../arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp | 10 +++-- examples/arkode/C_serial/lsrk_analytic.c | 2 +- .../arkode/C_serial/lsrk_analytic_varjac.c | 2 +- src/arkode/arkode_lsrkstep.c | 39 +++++++------------ 5 files changed, 24 insertions(+), 31 deletions(-) diff --git a/doc/arkode/guide/source/Usage/index.rst b/doc/arkode/guide/source/Usage/index.rst index e57475ced2..88cbfa7d3d 100644 --- a/doc/arkode/guide/source/Usage/index.rst +++ b/doc/arkode/guide/source/Usage/index.rst @@ -77,6 +77,6 @@ ARKBBDPRE can only be used with NVECTOR_PARALLEL. Preconditioners ARKStep/index.rst ERKStep/index.rst + LSRKStep/index.rst SPRKStep/index.rst MRIStep/index.rst - LSRKStep/index.rst diff --git a/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp b/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp index d35ee3940a..dd9c592fe2 100644 --- a/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp +++ b/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp @@ -196,8 +196,9 @@ struct UserData static int f(sunrealtype t, N_Vector u, N_Vector f, void* user_data); // Spectral radius estimation routine -static int eig(sunrealtype t, N_Vector y, sunrealtype* lambdaR, - sunrealtype* lambdaI, void* user_data); +static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3); // ----------------------------------------------------------------------------- // Helper functions @@ -1120,8 +1121,9 @@ static int WaitRecv(UserData* udata) } // Spectral radius estimation routine -static int eig(sunrealtype t, N_Vector y, sunrealtype* lambdaR, - sunrealtype* lambdaI, void* user_data) +static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3) { // Access problem data UserData* udata = (UserData*)user_data; diff --git a/examples/arkode/C_serial/lsrk_analytic.c b/examples/arkode/C_serial/lsrk_analytic.c index 98bfadfc4b..32c47df970 100644 --- a/examples/arkode/C_serial/lsrk_analytic.c +++ b/examples/arkode/C_serial/lsrk_analytic.c @@ -127,7 +127,7 @@ int main(void) if (check_flag(&flag, "LSRKStepSetDomEigFn", 1)) { return 1; } /* Specify after how many successful steps dom_eig is recomputed - Note that nsteps = 0 refers to Constant Jacobian */ + Note that nsteps = 0 refers to constant dominant eigenvalue */ flag = LSRKStepSetDomEigFrequency(arkode_mem, 0); if (check_flag(&flag, "LSRKStepSetDomEigFrequency", 1)) { return 1; } diff --git a/examples/arkode/C_serial/lsrk_analytic_varjac.c b/examples/arkode/C_serial/lsrk_analytic_varjac.c index 712de34943..6ba695ad47 100644 --- a/examples/arkode/C_serial/lsrk_analytic_varjac.c +++ b/examples/arkode/C_serial/lsrk_analytic_varjac.c @@ -157,7 +157,7 @@ int main(void) if (check_flag(&flag, "LSRKStepSetDomEigSafetyFactor", 1)) { return 1; } /* Specify the Runge--Kutta--Chebyshev LSRK method */ - flag = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKC_2); + flag = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKL_2); if (check_flag(&flag, "LSRKStepSetMethod", 1)) { return 1; } /* Open output stream for results, output comment line */ diff --git a/src/arkode/arkode_lsrkstep.c b/src/arkode/arkode_lsrkstep.c index f1628a91bd..577613b798 100644 --- a/src/arkode/arkode_lsrkstep.c +++ b/src/arkode/arkode_lsrkstep.c @@ -408,13 +408,11 @@ int lsrkStep_Init(ARKodeMem ark_mem, int init_type) ark_mem->e_data = ark_mem; } - /* Set default parameters for the default stepper */ + /* Set the default stepper */ if (ark_mem->step == lsrkStep_TakeStepRKC) { - step_mem->LSRKmethod = ARKODE_LSRK_RKC_2; - step_mem->nfusedopvecs = 5; - step_mem->q = ark_mem->hadapt_mem->q = 2; - step_mem->p = ark_mem->hadapt_mem->p = 2; + retval = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKC_2); + if (retval != ARK_SUCCESS) { return (retval); } } /* Allocate ARK RHS vector memory, update storage requirements */ @@ -595,21 +593,14 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } /* determine the number of required stages */ - int ss; - for (ss = 1; ss <= step_mem->stagemaxlimit; ss++) + int ss = SUNRceil(SUNRsqrt(onep54 * SUNRabs(ark_mem->h) * step_mem->sprad)); + step_mem->reqstages = SUNMAX(ss, 2); + if (step_mem->reqstages == step_mem->stagemaxlimit) { - if (SUNSQR(ss) >= (onep54 * SUNRabs(ark_mem->h) * step_mem->sprad)) - { - step_mem->reqstages = SUNMAX(ss, 2); - break; - } - if (ss == step_mem->stagemaxlimit) - { - hmax = SUN_RCONST(0.95) * SUNSQR(ss) / (onep54 * step_mem->sprad); - ark_mem->eta = hmax / ark_mem->h; - *nflagPtr = ARK_RETRY_STEP; - return (ARK_RETRY_STEP); - } + hmax = SUN_RCONST(0.95) * SUNSQR(ss) / (onep54 * step_mem->sprad); + ark_mem->eta = hmax / ark_mem->h; + *nflagPtr = ARK_RETRY_STEP; + return (ARK_RETRY_STEP); } step_mem->stagemax = SUNMAX(step_mem->reqstages, step_mem->stagemax); @@ -798,17 +789,17 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } /* determine the number of required stages */ - int ss; - for (ss = 1; ss <= step_mem->stagemaxlimit; ss++) + int ss = floor(SUNRsqrt(TWO * SUNRabs(ark_mem->h) * step_mem->sprad + THREE) - ONE); + for (ss; ss <= step_mem->stagemaxlimit; ss++) { - if ((SUNSQR(ss) + ss - 2) >= 2 * (SUNRabs(ark_mem->h) * step_mem->sprad)) + if ((SUNSQR(ss) + ss - TWO) >= TWO * (SUNRabs(ark_mem->h) * step_mem->sprad)) { step_mem->reqstages = SUNMAX(ss, 2); break; } if (ss == step_mem->stagemaxlimit) { - hmax = SUN_RCONST(0.95) * (SUNSQR(ss) + ss - 2) / (2.0 * step_mem->sprad); + hmax = SUN_RCONST(0.95) * (SUNSQR(ss) + ss - TWO) / (TWO * step_mem->sprad); ark_mem->eta = hmax / ark_mem->h; *nflagPtr = ARK_RETRY_STEP; return (ARK_RETRY_STEP); @@ -1295,7 +1286,7 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt ark_mem->fn, ark_mem->tempv1); } } - N_VLinearSum(SUN_RCONST(1.0 / 25.0), ark_mem->tempv2, SUN_RCONST(9.0 / 25.0), + N_VLinearSum(SUN_RCONST(1.0) / SUN_RCONST(25.0), ark_mem->tempv2, SUN_RCONST(9.0) / SUN_RCONST(25.0), ark_mem->ycur, ark_mem->tempv2); N_VLinearSum(SUN_RCONST(15.0), ark_mem->tempv2, SUN_RCONST(-5.0), ark_mem->ycur, ark_mem->ycur);