[petsc-users] Convergence rate for spatially varying Helmholtz system

[Barry Smith]

> On Sep 30, 2021, at 6:16 PM, Amneet Bhalla <mail2amneet at gmail.com> wrote:
> 
>  
> >> If you want to solve systems accurately, you should non-dimensionalize the system prior to discretization. This would mean that
> your C and b have elements in the [1, D] range, where D is the dynamic range of your problem, say 1e4, rather than these huge
> numbers you have now.
> 
> @Matt: We have done non-dimensionalization and the diagonal matrix ranges from 1 to 1e4 now. Still it takes 4-5 iterations to converge for the non-dimensional diagonal matrix. The convergence trend is looking much better now, though:
> 
> Residual norms for temperature_ solve.
>   0 KSP preconditioned resid norm 4.724547545716e-04 true resid norm 2.529423250889e+00 ||r(i)||/||b|| 4.397759655853e-05
>   1 KSP preconditioned resid norm 6.504853596318e-06 true resid norm 2.197130494439e-02 ||r(i)||/||b|| 3.820021755431e-07
>   2 KSP preconditioned resid norm 7.733420341215e-08 true resid norm 3.539290481432e-04 ||r(i)||/||b|| 6.153556501117e-09
>   3 KSP preconditioned resid norm 6.419092250844e-10 true resid norm 5.220398494466e-06 ||r(i)||/||b|| 9.076400273607e-11
>   4 KSP preconditioned resid norm 5.095955157158e-12 true resid norm 2.484163999489e-08 ||r(i)||/||b|| 4.319070053474e-13
>   5 KSP preconditioned resid norm 6.828200916501e-14 true resid norm 2.499229854610e-10 ||r(i)||/||b|| 4.345264170970e-15
> Linear temperature_ solve converged due to CONVERGED_RTOL iterations 5
> 
> 
> Only when all the equations are scaled individually the convergence is achieved in a single iteration. In the above, all equations are scaled using the same non-dimensional parameter. Do you think this is reasonable or do you expect the diagonal system to converge in a single iteration irrespective of the range of diagonal entries? 

   For a diagonal system with this modest range of values Jacobi should converge in a single iteration. 

   The output below is confusing, it is a system with 1 variable and should definitely converge in one iterations. 

   I am concerned we may be talking apples and oranges here and your test may not be as simple as you think it is (with regard to the diagonal).

> 
> @Barry: 
>> 
> 
>    What is the result of -ksp_view on the solve? 
> 
> KSP Object: (temperature_) 1 MPI processes
>   type: gmres
>     restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed
>     happy breakdown tolerance 1e-30
>   maximum iterations=1000, nonzero initial guess
>   tolerances:  relative=1e-12, absolute=1e-50, divergence=10000.
>   left preconditioning
>   using PRECONDITIONED norm type for convergence test
> PC Object: (temperature_) 1 MPI processes
>   type: shell
>     IEPSemiImplicitHierarchyIntegrator::helmholtz_precond::Temperature
>   linear system matrix = precond matrix:
>   Mat Object: 1 MPI processes
>     type: shell
>     rows=1, cols=1 
>  
> 
>    The way you describe your implementation it does not sound like standard PETSc practice. 
> 
> Yes, we do it differently in IBAMR. Succinctly, the main solver is a matrix-free one, whereas the preconditioner is a FAC multigrid solver with its bottom solver formed on the coarsest level of AMR grid using PETSc (matrix-based KSP).  
> 
> In the above -ksp_view temperature_ is the matrix-free KSP solver and IEPSemiImplicitHierarchyIntegrator::helmholtz_precond is the FAC preconditioner.
> 
> With PETSc using a matrix-free operation mA and a matrix from which KSP will build the preconditioner  A one uses  KSPSetOperator(ksp,mA,A); and then just selects the preconditioner with -pc_type xxx  For example to use Jacobi preconditioning one uses -pc_type jacobi (note that this only uses the diagonal of A, the rest of A is never used).
> 
> We run -pc_type jacobi on the bottom solver of the FAC preconditioner.  
> 
> If you wish to precondition mA by fully solving with the matrix A one can use -ksp_monitor_true_residual -pc_type ksp -ksp_ksp_type yyy -ksp_pc_type xxx  -ksp_ksp_monitor_true_residual with, for example, yyy of richardson and xxx of jacobi
> 
> Yes, this is what we do.  
> 
>   Barry
> 
> 
> 
> 
>> To verify that I am indeed solving a diagonal system I printed the PETSc matrix from the preconditioner and viewed it in Matlab. It indeed shows it to be a diagonal system. Attached is the plot of the spy command on the printed matrix. The matrix in binary form is also attached. 
>> 
>> My understanding is that because the C coefficient is varying in 4 orders of magnitude, i.e., Max(C)/Min(C) ~ 10^4, the matrix is poorly scaled. When I rescale my matrix by 1/C then the system converges in 1 iteration as expected. Is my understanding correct, and that scaling 1/C should be done even for a diagonal system?
>> 
>> When D is non-zero, then scaling by 1/C seems to be very inconvenient as D is stored as side-centered data for the matrix free solver. 
>> 
>> In the case that I do not scale my equations by 1/C, is there some solver setting that improves the convergence rate? (With D as non-zero, I have also tried gmres as the ksp solver in the matrix-based preconditioner to get better performance, but it didn't matter much.)
>> 
>> 
>> Thanks,
>> Ramakrishnan Thirumalaisamy
>> San Diego State University.
>> <Temperature_fill.pdf><matrix_temperature>
> 
> 
> 
> -- 
> --Amneet 
> 
> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210930/30e8302f/attachment-0001.html>