<html><head><meta http-equiv="content-type" content="text/html; charset=utf-8"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br></div> I always like a challenge.<div><br></div><div> The MR <a href="https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/8675__;!!G_uCfscf7eWS!eEJyFKOE1fIr1zG138qFEBHnEAPDDHz6-nhpA3l8TnLtt8VK3i98Cc13eRy7zJGgmWBtGQ8IB0yIDNjR7M6rtACjJrE$">https://gitlab.com/petsc/petsc/-/merge_requests/8675</a> should eliminate both redundant TSComputeRHSFunction() calls in your code and thus decrease its run time. </div><div><br></div><div> Thanks for bringing the performance bug to our attention </div><div><br></div><div> Barry</div><div><br id="lineBreakAtBeginningOfMessage"><div><br><blockquote type="cite"><div>On Aug 30, 2025, at 5:42 PM, Aldo Bonfiglioli <aldo.bonfiglioli@unibas.it> wrote:</div><br class="Apple-interchange-newline"><div><div dir="auto"><div>Barry,</div><div dir="auto">thank you for the update.</div><div dir="auto">Aldo </div><div><br></div><div data-smartmail="gmail_signature"><div dir="ltr"><img width="200" height="116" src="https://ci3.googleusercontent.com/mail-sig/AIorK4xnXXTqJRS0ILQ8nMbbv0sxqjdpcbPqBphPJIRPxrTxJsSffBROVMincHK3lqUgs8EgQg_8rNdfT0Ea"><br><div><font face="garamond, times new roman, serif"><b>Aldo Bonfiglioli</b></font></div><div><font face="garamond, times new roman, serif">Professore associato di </font><span style="font-family:garamond,"times new roman",serif">Fluidodinamica </span><span style="font-family:garamond,"times new roman",serif">settore scientifico disciplinare IIND-01/F </span></div><div><span style="font-family:garamond,"times new roman",serif">Dipartimento di Ingegneria</span></div><div><font face="garamond, times new roman, serif">Viale dell'Ateneo Lucano, 10 85100 Potenza</font></div><div><font face="garamond, times new roman, serif">Tel: 0971.205203</font></div><div><font face="garamond, times new roman, serif"><br></font></div></div></div></div><br><div class="gmail_quote gmail_quote_container"><div dir="ltr" class="gmail_attr">Il ven 29 ago 2025, 20:52 Barry Smith <<a href="mailto:barryfsmith@icloud.com">barryfsmith@icloud.com</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="line-break:after-white-space"><div><br></div> Aldo,<div><br></div><div> I understand the situation you are facing. </div><div><br></div><div> With pseudo-continuation, when you provide an RHS function, SNES is computing udot -</div><div>F(u). Inside the SNESSolve(), specifically in TSComputeIFunction(), with the lines</div><div><br></div><div><div> if (!ts->ops->ifunction) {</div><div> ierr = TSComputeRHSFunction(ts,t,X,Y);CHKERRQ(ierr);</div></div><div> ierr = VecAYPX(Y,-1,Xdot);CHKERRQ(ierr);</div><div><br></div><div>But for the original system, we care about the norm for F(u). So in TSStep_Pseudo(), after the nonlinear solve</div><div>there is the line</div><div><br></div><div> pseudo->fnorm = -1; /* The current norm is no longer valid, monitor must recompute it. */</div><div><br></div><div>and further down</div><div><br></div><div><div> if (pseudo->fnorm < 0) {</div><div> PetscCall(VecZeroEntries(pseudo->xdot));</div><div> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));</div><div> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));</div></div><div><br></div><div>this norm is needed to determine of the nonlinear system has converged satisfactory (or if more pseudo time steps are needed).</div><div><br></div><div>Hence, your FormRHSFunction gets called twice with identical input. </div><div><br></div><div>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 </div><div><br></div><div><div> PetscCall(VecZeroEntries(pseudo->xdot));</div><div> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));</div><div> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));</div></div><div><br></div><div>And the same norm is also recomputed.</div><div><br></div><div>Now the diff </div><div><br></div><div><div>diff --git a/src/ts/impls/pseudo/posindep.c b/src/ts/impls/pseudo/posindep.c</div><div>index 7c013435f98..813a28311e7 100644</div><div>--- a/src/ts/impls/pseudo/posindep.c</div><div>+++ b/src/ts/impls/pseudo/posindep.c</div><div>@@ -660,9 +660,11 @@ PetscErrorCode TSPseudoTimeStepDefault(TS ts, PetscReal *newdt, void *dtctx)</div><div> PetscReal inc = pseudo->dt_increment;</div><div> </div><div> PetscFunctionBegin;</div><div>- PetscCall(VecZeroEntries(pseudo->xdot));</div><div>- PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));</div><div>- PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));</div><div>+ if (pseudo->fnorm < 0.0) {</div><div>+ PetscCall(VecZeroEntries(pseudo->xdot));</div><div>+ PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));</div><div>+ PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));</div><div>+ }</div><div> if (pseudo->fnorm_initial < 0) {</div><div> /* first time through so compute initial function norm */</div><div> pseudo->fnorm_initial = pseudo->fnorm</div></div><div><br></div><div>will prevent the THIRD computation of the same function evaluation. So this is a good optimization</div><div><br></div><div>But eliminating the second computation is not so trivial. To review, deep inside the TS code there is</div><div><br></div><div><div>TSComputeRHSFunction(ts,t,X,Y);</div><div>VecAYPX(Y,-1,Xdot);</div></div><div><br></div><div>This code knows nothing about TSPPseudo etc.</div><div><br></div><div>Then later in the Pseudo code</div><div><br></div><div><div> PetscCall(VecZeroEntries(pseudo->xdot));</div><div> PetscCall(TSComputeIFunction(ts, ts->ptime, ts->vec_sol, pseudo->xdot, pseudo->func, PETSC_FALSE));</div><div> PetscCall(VecNorm(pseudo->func, NORM_2, &pseudo->fnorm));</div></div><div><br></div><div> again calls TSComputeRHSFunction(ts,t,X,Y); with the same X. (since Xdot is not used directly by TSComputeRHSFunction()). </div><div><br></div><div>Now if all the code was custom for Pseudo we could simply put in a line </div><div><br></div><div><div><div>TSComputeRHSFunction(ts,t,X,Y);</div><div>VecNorm(Y,NORM_2, &pseudo->fnorm)</div><div>VecAYPX(Y,-1,Xdot);</div></div></div><div><br></div><div>into TSComputeIFunction() and have the norm available later for pseudo. </div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>Thanks for reporting this challenging problem. </div><div><br></div><div> Barry</div><div><br></div><div>Note: I don't think providing the Jacobian improves the efficiency, it is just better at hiding the redundant computation.</div><div> </div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br id="m_462023396166240265lineBreakAtBeginningOfMessage"><div><br><blockquote type="cite"><div>On Aug 26, 2025, at 4:21 AM, Aldo Bonfiglioli <<a href="mailto:aldo.bonfiglioli@unibas.it" target="_blank" rel="noreferrer">aldo.bonfiglioli@unibas.it</a>> wrote:</div><br><div>
<div><p>Hi there,</p><p>I am using -ts_type pseudo to find steady solutions of the
Euler/NS PDEs.</p><p>It looks like FormRHSFunction is called twice with the exactly
same value of the dependent variable (u)</p><p>before invoking TSRHSJacobianFn and then solving the linear
system with KSP (u is then updated)<br>
</p><div>
<br></div><blockquote type="cite"><span style="font-family:monospace"><span style="background-color:rgb(255,255,255)">1 TS dt
0.141197 time 0.141197</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
FormRHSFunction; ||u|| = 60.489405479528003
463.15635675129079 62.813532026019146
5.9789502719351155 </span><span style="background-color:rgb(255,255,255)"> </span><br>
<span style="background-color:rgb(255,255,255)"> Calling
FormRHSFunction; ||u|| = 60.489405479528003
463.15635675129079</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">62.813532026019146</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">5.9789502719351155</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
TSRHSJacobianFn with ||u|| = 60.489405479528003
463.15635675129079 62.813532026019146
5.9789502719351155</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)">
Linear solve converged due to CONVERGED_RTOL iterations
9</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
FormRHSFunction; ||u|| = 60.489405476878737</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">463.15635675602186</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">62.813532031616965</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">5.9789502710519811</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)">2 TS dt
0.155917 time 0.297115</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
FormRHSFunction; ||u|| = 60.489405476878737</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">463.15635675602186</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">62.813532031616965</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">5.9789502710519811</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
FormRHSFunction; ||u|| = 60.489405476878737</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">463.15635675602186</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">62.813532031616965</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">5.9789502710519811</span><span style="background-color:rgb(255,255,255)">
</span><br>
<span style="background-color:rgb(255,255,255)"> Calling
TSRHSJacobianFn with ||u|| = 60.489405476878737</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">463.15635675602186</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">62.813532031616965</span><span style="background-color:rgb(255,255,255)"> </span><span style="background-color:rgb(255,255,255)">5.9789502710519811</span><br>
<span style="background-color:rgb(255,255,255)">
</span><br>
</span></blockquote>
Why is it so? Am I possibly misusing (or missing) some options? <br><div><br></div><p>My .petscrc file is attached.</p><p>Thank you,</p><p>Aldo<br>
</p>
<pre cols="72">--
Dr. Aldo Bonfiglioli
Associate professor of Fluid Mechanics
Dipartimento di Ingegneria
Universita' della Basilicata
V.le <a href="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$">dell'Ateneo Lucano, 10 85100 Potenza ITALY</a>
tel:+39.0971.205203 fax:+39.0971.205215
web: <a href="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$" target="_blank" rel="noreferrer">http://docenti.unibas.it/site/home/docente.html?m=002423</a></pre>
</div>
</div></blockquote></div></div></div><div style="line-break:after-white-space"><div><div><blockquote type="cite"><div></div></blockquote></div><br></div></div></blockquote></div>
</div></blockquote></div><br></div></body></html>