[petsc-users] Question about Pseudo Time Stepping
Barry Smith
barryfsmith at icloud.com
Sat Aug 30 22:18:24 CDT 2025
I always like a challenge.
The MR https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/8675__;!!G_uCfscf7eWS!eEJyFKOE1fIr1zG138qFEBHnEAPDDHz6-nhpA3l8TnLtt8VK3i98Cc13eRy7zJGgmWBtGQ8IB0yIDNjR7M6rtACjJrE$ should eliminate both redundant TSComputeRHSFunction() calls in your code and thus decrease its run time.
Thanks for bringing the performance bug to our attention
Barry
> On Aug 30, 2025, at 5:42 PM, Aldo Bonfiglioli <aldo.bonfiglioli at unibas.it> wrote:
>
> Barry,
> thank you for the update.
> Aldo
>
>
> Aldo Bonfiglioli
> Professore associato di Fluidodinamica settore scientifico disciplinare IIND-01/F
> Dipartimento di Ingegneria
> Viale dell'Ateneo Lucano, 10 85100 Potenza
> Tel: 0971.205203
>
>
> Il ven 29 ago 2025, 20:52 Barry Smith <barryfsmith at icloud.com <mailto:barryfsmith at icloud.com>> ha scritto:
>>
>> Aldo,
>>
>> I understand the situation you are facing.
>>
>> With pseudo-continuation, when you provide an RHS function, SNES is computing udot -
>> F(u). Inside the SNESSolve(), specifically in TSComputeIFunction(), with the lines
>>
>> if (!ts->ops->ifunction) {
>> ierr = TSComputeRHSFunction(ts,t,X,Y);CHKERRQ(ierr);
>> ierr = VecAYPX(Y,-1,Xdot);CHKERRQ(ierr);
>>
>> But for the original system, we care about the norm for F(u). So in TSStep_Pseudo(), after the nonlinear solve
>> there is the line
>>
>> pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */
>>
>> and further down
>>
>> if (pseudo->fnorm < 0) {
>> PetscCall(VecZeroEntries(pseudo->xdot));
>> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
>> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
>>
>> this norm is needed to determine of the nonlinear system has converged satisfactory (or if more pseudo time steps are needed).
>>
>> Hence, your FormRHSFunction gets called twice with identical input.
>>
>> To make matters even worse your FormRHSFunction is actually being called THREE times with identical input. Because at the start of TSStep_Pseudo() it calls TSPseudoComputeTimeStep(ts, &next_time_step) which by default uses TSPseudoTimeStepDefault() which has the same chunk of code
>>
>> PetscCall(VecZeroEntries(pseudo->xdot));
>> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
>> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
>>
>> And the same norm is also recomputed.
>>
>> Now the diff
>>
>> diff --git a/src/ts/impls/pseudo/posindep.c b/src/ts/impls/pseudo/posindep.c
>> index 7c013435f98..813a28311e7 100644
>> --- a/src/ts/impls/pseudo/posindep.c
>> +++ b/src/ts/impls/pseudo/posindep.c
>> @@ -660,9 +660,11 @@ PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)
>> PetscReal inc = pseudo->dt_increment;
>>
>> PetscFunctionBegin;
>> - PetscCall(VecZeroEntries(pseudo->xdot));
>> - PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
>> - PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
>> + if (pseudo->fnorm < 0.0) {
>> + PetscCall(VecZeroEntries(pseudo->xdot));
>> + PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
>> + PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
>> + }
>> if (pseudo->fnorm_initial < 0) {
>> /* first time through so compute initial function norm */
>> pseudo->fnorm_initial = pseudo->fnorm
>>
>> will prevent the THIRD computation of the same function evaluation. So this is a good optimization
>>
>> But eliminating the second computation is not so trivial. To review, deep inside the TS code there is
>>
>> TSComputeRHSFunction(ts,t,X,Y);
>> VecAYPX(Y,-1,Xdot);
>>
>> This code knows nothing about TSPPseudo etc.
>>
>> Then later in the Pseudo code
>>
>> PetscCall(VecZeroEntries(pseudo->xdot));
>> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));
>> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));
>>
>> again calls TSComputeRHSFunction(ts,t,X,Y); with the same X. (since Xdot is not used directly by TSComputeRHSFunction()).
>>
>> Now if all the code was custom for Pseudo we could simply put in a line
>>
>> TSComputeRHSFunction(ts,t,X,Y);
>> VecNorm(Y,NORM_2, &pseudo->fnorm)
>> VecAYPX(Y,-1,Xdot);
>>
>> into TSComputeIFunction() and have the norm available later for pseudo.
>>
>> It is not clear that there is any definitive benefit for pseudo to utilize the TS framework. It possibly gets some code reuse from TS but at the cost of extra code that won't be needed if written directly. Plus this performance hit.
>>
>> Anyways I will make a merge request to remove the third redundant computation of FormRHSFunction but will await input from others on if anything else could be done to remove the inefficiency.
>>
>> Thanks for reporting this challenging problem.
>>
>> Barry
>>
>> Note: I don't think providing the Jacobian improves the efficiency, it is just better at hiding the redundant computation.
>>
>>
>>
>>
>>
>>
>>
>>
>>> On Aug 26, 2025, at 4:21 AM, Aldo Bonfiglioli <aldo.bonfiglioli at unibas.it <mailto:aldo.bonfiglioli at unibas.it>> wrote:
>>>
>>> Hi there,
>>>
>>> I am using -ts_type pseudo to find steady solutions of the Euler/NS PDEs.
>>>
>>> It looks like FormRHSFunction is called twice with the exactly same value of the dependent variable (u)
>>>
>>> before invoking TSRHSJacobianFn and then solving the linear system with KSP (u is then updated)
>>>
>>>
>>>> 1 TS dt 0.141197 time 0.141197
>>>> Calling FormRHSFunction; ||u|| = 60.489405479528003 463.15635675129079 62.813532026019146 5.9789502719351155
>>>> Calling FormRHSFunction; ||u|| = 60.489405479528003 463.15635675129079 62.813532026019146 5.9789502719351155
>>>> Calling TSRHSJacobianFn with ||u|| = 60.489405479528003 463.15635675129079 62.813532026019146 5.9789502719351155
>>>> Linear solve converged due to CONVERGED_RTOL iterations 9
>>>> Calling FormRHSFunction; ||u|| = 60.489405476878737 463.15635675602186 62.813532031616965 5.9789502710519811
>>>> 2 TS dt 0.155917 time 0.297115
>>>> Calling FormRHSFunction; ||u|| = 60.489405476878737 463.15635675602186 62.813532031616965 5.9789502710519811
>>>> Calling FormRHSFunction; ||u|| = 60.489405476878737 463.15635675602186 62.813532031616965 5.9789502710519811
>>>> Calling TSRHSJacobianFn with ||u|| = 60.489405476878737 463.15635675602186 62.813532031616965 5.9789502710519811
>>>>
>>> Why is it so? Am I possibly misusing (or missing) some options?
>>>
>>> My .petscrc file is attached.
>>>
>>> Thank you,
>>>
>>> Aldo
>>>
>>> --
>>> Dr. Aldo Bonfiglioli
>>> Associate professor of Fluid Mechanics
>>> Dipartimento di Ingegneria
>>> Universita' della Basilicata
>>> V.le dell'Ateneo Lucano, 10 85100 Potenza ITALY <https://urldefense.us/v3/__https://www.google.com/maps/search/dell'Ateneo*Lucano,*10*85100*Potenza*ITALY?entry=gmail&source=g__;KysrKys!!G_uCfscf7eWS!eEJyFKOE1fIr1zG138qFEBHnEAPDDHz6-nhpA3l8TnLtt8VK3i98Cc13eRy7zJGgmWBtGQ8IB0yIDNjR7M6rNuzDtLU$ >
>>> tel:+39.0971.205203 fax:+39.0971.205215
>>> web: https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!eEJyFKOE1fIr1zG138qFEBHnEAPDDHz6-nhpA3l8TnLtt8VK3i98Cc13eRy7zJGgmWBtGQ8IB0yIDNjR7M6rchJhC0A$ <https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!ci_aySHdF2Qx_B-37IOzELFDtcQRthJ2NItK-JNbfXB0M2_21B_z97WIfjxZUyKCVRRhd1mzL6b-rnCed_i8LLA56JvFdz7JBrw$>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250830/10ee7c4f/attachment.html>
More information about the petsc-users
mailing list