Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Sorry this code has not been changed. Barry > On Sep 30, 2019, at 4:24 PM, Sajid Ali > wrote: > > Hi PETSc-developers, > > Has this bug been fixed in the new 3.12 release ? > > Thank You, > Sajid Ali > Applied Physics > Northwestern University > s-sajid-ali.github.io
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi PETSc-developers, Has this bug been fixed in the new 3.12 release ? Thank You, Sajid Ali Applied Physics Northwestern University s-sajid-ali.github.io
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
There is no harm in having the GMRES there even if you use a direct solver (for testing) so just leave the GMRES. Changing to preonly every time you try LU is prone to error if you forget to change back. Barry > On May 22, 2019, at 2:45 PM, Sajid Ali via petsc-users > wrote: > > Hi Matt, > > Thanks for the explanation. That makes sense since I'd get reasonably close > convergence with preonly sometimes and not at other times which was confusing. > > Anyway, since there's no pc_tol (analogous to ksp_rtol/ksp_atol, etc), I'd > have to more carefully set the gamg preconditioner options to ensure that it > converges in one run, but since there's no guarantee that what works for one > problem might not work for another (or the same problem at a different grid > size), I'll stick with GMRES+gamg for now. > > Thank You, > Sajid Ali > Applied Physics > Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
> On May 22, 2019, at 2:26 PM, Sajid Ali via petsc-users > wrote: > > Hi Hong, > > Looks like this is my fault since I'm using -ksp_type preonly -pc_type gamg. > If I use the default ksp (GMRES) then everything works fine on a smaller > problem. > > Just to confirm, -ksp_type preonly is to be used only with direct-solve > preconditioners like LU,Cholesky, right ? You can use it any time you like but it only applies the preconditioner; thus unless your preconditioner is a really good approximation to the operator it won't give you much information. > > Thank You, > Sajid Ali > Applied Physics > Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi Matt, Thanks for the explanation. That makes sense since I'd get reasonably close convergence with preonly sometimes and not at other times which was confusing. Anyway, since there's no pc_tol (analogous to ksp_rtol/ksp_atol, etc), I'd have to more carefully set the gamg preconditioner options to ensure that it converges in one run, but since there's no guarantee that what works for one problem might not work for another (or the same problem at a different grid size), I'll stick with GMRES+gamg for now. Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi Hong, Looks like this is my fault since I'm using -ksp_type preonly -pc_type gamg. If I use the default ksp (GMRES) then everything works fine on a smaller problem. Just to confirm, -ksp_type preonly is to be used only with direct-solve preconditioners like LU,Cholesky, right ? Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Sajid, I have also rested the simpler problem you provided. The branch hongzh/fix-computejacobian gives exactly the same numerical results as the master branch does, but runs much faster. So the solver seems to work correctly. To rule out the possible compiler issues, you might want to try a different compiler or different optimization flags. Also you might want to try smaller stepsizes. Hong On May 20, 2019, at 4:47 PM, Sajid Ali mailto:sajidsyed2...@u.northwestern.edu>> wrote: Hi Hong, I tried running a simpler problem that solves the equation ` u_t = A*u_xx + A*u_yy; ` and the fix-computejacobian branch works for this on a coarse grid. The code for the same is here : https://github.com/s-sajid-ali/xwp_petsc/blob/master/2d/FD/free_space/ex_dmda.c and it requires no input file. It writes at all times steps and comparing the output at the last time step, everything looks fine. I want to eliminate a possible source of error which could be the fact that I installed both versions (3.11.2 and fix-computejacobian) with intel compilers and O3. Could floating point errors occur due to this ? I didn't specify -fp-model strict but since the results I got were reasonable I never bothered to run a test suite. Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi Hong, The solution has the right characteristics but it's off by many orders of magnitude. It is ~3.5x faster as before. Am I supposed to keep the TSRHSJacobianSetReuse function or not? Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
> On May 16, 2019, at 8:04 PM, Sajid Ali > wrote: > > While there is a ~3.5X speedup, deleting the aforementioned 20 lines also > leads the new version of petsc to give the wrong solution (off by orders of > magnitude for the same program). Ok, sorry about this. Unfortunately this stuff has been giving us headaches for years and we are struggling to get it right. > > I tried switching over the the IFunction/IJacobian interface as per the > manual (page 146) which the following lines : It is probably better to not switch to the IFunction/IJacobian, we are more likely to get the TS version working properly. > ``` > TSSetProblemType(ts,TSLINEAR); > TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL); > TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,NULL); > ``` > are equivalent to : > ``` > TSSetProblemType(ts,TSLINEAR); > TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,NULL); > TSSetIJacobian(ts,A,A,TSComputeIJacobianConstant,NULL); > ``` > But the example at src/ts/examples/tutorials/ex3.c employs a strategy of > setting a shift flag to prevent re-computation for time-independent problems. > Moreover, the docs say "using this function (TSComputeIFunctionLinear) is NOT > equivalent to using TSComputeRHSFunctionLinear()" and now I'm even more > confused. > > PS : Doing the simple switch is as slow as the original code and the answer > is wrong as well. > > Thank You, > Sajid Ali > Applied Physics > Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi Sajid, Can you please try this branch hongzh/fix-computejacobian quickly and see if it makes a difference? Thanks, Hong (Mr.) On May 16, 2019, at 8:04 PM, Sajid Ali via petsc-users mailto:petsc-users@mcs.anl.gov>> wrote: While there is a ~3.5X speedup, deleting the aforementioned 20 lines also leads the new version of petsc to give the wrong solution (off by orders of magnitude for the same program). I tried switching over the the IFunction/IJacobian interface as per the manual (page 146) which the following lines : ``` TSSetProblemType(ts,TSLINEAR); TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL); TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,NULL); ``` are equivalent to : ``` TSSetProblemType(ts,TSLINEAR); TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,NULL); TSSetIJacobian(ts,A,A,TSComputeIJacobianConstant,NULL); ``` But the example at src/ts/examples/tutorials/ex3.c employs a strategy of setting a shift flag to prevent re-computation for time-independent problems. Moreover, the docs say "using this function (TSComputeIFunctionLinear) is NOT equivalent to using TSComputeRHSFunctionLinear()" and now I'm even more confused. PS : Doing the simple switch is as slow as the original code and the answer is wrong as well. Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
While there is a ~3.5X speedup, deleting the aforementioned 20 lines also leads the new version of petsc to give the wrong solution (off by orders of magnitude for the same program). I tried switching over the the IFunction/IJacobian interface as per the manual (page 146) which the following lines : ``` TSSetProblemType(ts,TSLINEAR); TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL); TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,NULL); ``` are equivalent to : ``` TSSetProblemType(ts,TSLINEAR); TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,NULL); TSSetIJacobian(ts,A,A,TSComputeIJacobianConstant,NULL); ``` But the example at src/ts/examples/tutorials/ex3.c employs a strategy of setting a shift flag to prevent re-computation for time-independent problems. Moreover, the docs say "using this function (TSComputeIFunctionLinear) is NOT equivalent to using TSComputeRHSFunctionLinear()" and now I'm even more confused. PS : Doing the simple switch is as slow as the original code and the answer is wrong as well. Thank You, Sajid Ali Applied Physics Northwestern University
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Hi Barry, Thanks a lot for pointing this out. I'm seeing ~3X speedup in time ! Attached are the new log files. Does everything look right ? Thank You, Sajid Ali Applied Physics Northwestern University out_50 Description: Binary data out_100 Description: Binary data
Re: [petsc-users] Question about TSComputeRHSJacobianConstant
Sajid,
This is a huge embarrassing performance bug in PETSc
https://bitbucket.org/petsc/petsc/issues/293/refactoring-of-ts-handling-of-reuse-of
It is using 74 percent of the time to perform MatAXPY() on two large sparse
matrices, not knowing they have identical nonzero patterns and one of which has
all zeros off of the diagonal. This despite the fact that a few lines higher in
the code is special purpose code for exactly the case you have that only
stores one matrix and only ever shifts the diagonal of the matrix.
Please edit TSSetUp() and remove the lines
if (ts->rhsjacobian.reuse && rhsjac == TSComputeRHSJacobianConstant) {
Mat Amat,Pmat;
SNES snes;
ierr = TSGetSNES(ts,);CHKERRQ(ierr);
ierr = SNESGetJacobian(snes,,,NULL,NULL);CHKERRQ(ierr);
/* Matching matrices implies that an IJacobian is NOT set, because if it
had been set, the IJacobian's matrix would
* have displaced the RHS matrix */
if (Amat && Amat == ts->Arhs) {
/* we need to copy the values of the matrix because for the constant
Jacobian case the user will never set the numerical values in this new location
*/
ierr = MatDuplicate(ts->Arhs,MAT_COPY_VALUES,);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr);
ierr = MatDestroy();CHKERRQ(ierr);
}
if (Pmat && Pmat == ts->Brhs) {
ierr = MatDuplicate(ts->Brhs,MAT_COPY_VALUES,);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr);
ierr = MatDestroy();CHKERRQ(ierr);
}
}
You will be stunned by the improvement in time.
> On May 16, 2019, at 3:06 PM, Sajid Ali via petsc-users
> wrote:
>
> Hi PETSc developers,
>
> I have a question about TSComputeRHSJacobianConstant. If I create a TS (of
> type linear) for a problem where the jacobian does not change with time (set
> with the aforementioned option) and run it for different number of time
> steps, why does the time it takes to evaluate the jacobian change (as
> indicated by TSJacobianEval) ?
>
> To clarify, I run with the example with different TSSetTimeStep, but the same
> jacobian matrix. I see that the time spent in KSPSolve increases with
> increasing number of steps (which is as expected as this is a KSPOnly SNES
> solver). But surprisingly, the time spent in TSJacobianEval also increases
> with decreasing time-step (or increasing number of steps).
>
> For reference, I attach the log files for two cases which were run with
> different time steps and the source code.
>
> Thank You,
> Sajid Ali
> Applied Physics
> Northwestern University
>
