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

[Matthew Knepley]

On Wed, Sep 29, 2021 at 6:03 PM Ramakrishnan Thirumalaisamy <
rthirumalaisam1857 at sdsu.edu> wrote:

> Hi all,
>
> I am trying to solve the Helmholtz equation for temperature T:
>
> (C I  + Div D grad) T = f
>
> in IBAMR, in which C is the spatially varying diagonal entries, and D is
> the spatially varying diffusion coefficient.   I use a matrix-free solver
> with matrix-based PETSc preconditioner. For the matrix-free solver, I use
> gmres solver and for the matrix based preconditioner, I use Richardson ksp
> + Jacobi as a preconditioner. As the simulation progresses, the iterations
> start to increase. To understand the cause, I set D to be zero, which
> results in a diagonal system:
>
> C T = f.
>
> This should result in convergence within a single iteration, but I get
> convergence in 3 iterations.
>
> Residual norms for temperature_ solve.
>
>   0 KSP preconditioned resid norm 4.590811647875e-02 true resid norm
> 2.406067589273e+09 ||r(i)||/||b|| 4.455533946945e-05
>
>   1 KSP preconditioned resid norm 2.347767895880e-06 true resid norm
> 1.210763896685e+05 ||r(i)||/||b|| 2.242081505717e-09
>
>   2 KSP preconditioned resid norm 1.245406571896e-10 true resid norm
> 6.328828824310e+00 ||r(i)||/||b|| 1.171966730978e-13
>
> Linear temperature_ solve converged due to CONVERGED_RTOL iterations 2
>

Several things look off here:

1) Your true residual norm is 2.4e9, but r_0/b is 4.4e-5. That seems to
indicate that ||b|| is 1e14. Is this true?

2) Your preconditioned residual is 11 orders of magnitude less than the
true residual. This usually indicates that the system is near singular.

3) The disparity above does not seem possible if C only has elements ~ 1e4.
The preconditioner consistently has norm around 1e-11.

4) Using numbers that large can be a problem. You lose precision, so that
you really only have 3-4 correct digits each time, as you see above. It
appears to be
     doing iterative refinement with a very low precision solve.

  Thanks,

     Matt


> 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.
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210929/25223c1d/attachment-0001.html>