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

[Matthew Knepley]

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

> 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?
>>
> Yes. ||b|| is 1e14.  We have verified that.
>
>>
>> 2) Your preconditioned residual is 11 orders of magnitude less than the
>> true residual. This usually indicates that the system is near singular.
>>
> The system is diagonal, as shown in Temperature_fill.pdf. So we don't
> think it is singular unless we are missing something very obvious.
> Diagonal elements range from 8.8e6 to 5.62896e+10, while off-diagonal terms
> are 0, as shown in the spy plot.
>
>>
>> 3) The disparity above does not seem possible if C only has elements ~
>> 1e4. The preconditioner consistently has norm around 1e-11.
>>
> The value of C in the Helmholtz system is computed as : *C =
> rho*specific_heat/dt* in which dt = 5e-5, specific_heat ~10^3 and rho
> ranges from 0.4 to 2700. Hence, C ranges from 8.8e6 to 5.62896e10.
> Max_diagonal(C)/Min_diagonal(C) ~ 10^4.
>
>>
>> 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.
>>
> Indeed the numbers are large because C =  rho*specific_heat/dt.
>

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.

  Thanks,

     Matt


> On Wed, Sep 29, 2021 at 3:58 PM Matthew Knepley <knepley at gmail.com> wrote:
>
>> 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/>
>>
>

-- 
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/20210930/af907835/attachment-0001.html>