<div dir="ltr"><div dir="ltr">On Wed, Sep 29, 2021 at 8:31 PM Ramakrishnan Thirumalaisamy <<a href="mailto:rthirumalaisam1857@sdsu.edu">rthirumalaisam1857@sdsu.edu</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div>Several things look off here:</div><div><br></div><div>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? </div></div></div></blockquote><div>Yes. ||b|| is 1e14.  We have verified that.</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div><br></div><div>2) Your preconditioned residual is 11 orders of magnitude less than the true residual. This usually indicates that the system is near singular.</div></div></div></blockquote><div>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. </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div><br></div><div>3) The disparity above does not seem possible if C only has elements ~ 1e4. The preconditioner consistently has norm around 1e-11.</div></div></div></blockquote><div>The value of C in the Helmholtz system is computed as : <b>C = rho*specific_heat/dt</b> 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.</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div><br></div><div>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</div><div>     doing iterative refinement with a very low precision solve.</div></div></div></blockquote><div>Indeed the numbers are large because C =  rho*specific_heat/dt.</div></div></div></blockquote><div><br></div><div>If you want to solve systems accurately, you should non-dimensionalize the system prior to discretization. This would mean that</div><div>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</div><div>numbers you have now.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Sep 29, 2021 at 3:58 PM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Wed, Sep 29, 2021 at 6:03 PM Ramakrishnan Thirumalaisamy <<a href="mailto:rthirumalaisam1857@sdsu.edu" target="_blank">rthirumalaisam1857@sdsu.edu</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hi all,<div><br></div><div>I am trying to solve the Helmholtz equation for temperature T:</div><div><br></div><div>(C I  + Div D grad) T = f</div><div><br></div><div>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:</div><div>  </div><div>C T = f.</div><div><br></div><div>This should result in convergence within a single iteration, but I get convergence in 3 iterations.</div><div><br></div><div><p style="margin:0px;font-stretch:normal;font-size:15px;line-height:normal;font-family:Monaco;color:rgb(242,242,242);background-color:rgba(0,0,0,0.85)"><span style="font-variant-ligatures:no-common-ligatures">Residual norms for temperature_ solve.</span></p>
<p style="margin:0px;font-stretch:normal;font-size:15px;line-height:normal;font-family:Monaco;color:rgb(242,242,242);background-color:rgba(0,0,0,0.85)"><span style="font-variant-ligatures:no-common-ligatures"><span>  </span>0 KSP preconditioned resid norm 4.590811647875e-02 true resid norm 2.406067589273e+09 ||r(i)||/||b|| 4.455533946945e-05</span></p>
<p style="margin:0px;font-stretch:normal;font-size:15px;line-height:normal;font-family:Monaco;color:rgb(242,242,242);background-color:rgba(0,0,0,0.85)"><span style="font-variant-ligatures:no-common-ligatures"><span>  </span>1 KSP preconditioned resid norm 2.347767895880e-06 true resid norm 1.210763896685e+05 ||r(i)||/||b|| 2.242081505717e-09</span></p>
<p style="margin:0px;font-stretch:normal;font-size:15px;line-height:normal;font-family:Monaco;color:rgb(242,242,242);background-color:rgba(0,0,0,0.85)"><span style="font-variant-ligatures:no-common-ligatures"><span>  </span>2 KSP preconditioned resid norm 1.245406571896e-10 true resid norm 6.328828824310e+00 ||r(i)||/||b|| 1.171966730978e-13</span></p>
<p style="margin:0px;font-stretch:normal;font-size:15px;line-height:normal;font-family:Monaco;color:rgb(242,242,242);background-color:rgba(0,0,0,0.85)"><span style="font-variant-ligatures:no-common-ligatures">Linear temperature_ solve converged due to CONVERGED_RTOL iterations 2</span></p></div></div></blockquote><div><br></div><div>Several things look off here:</div><div><br></div><div>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? </div></div></div></blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div><br></div><div>2) Your preconditioned residual is 11 orders of magnitude less than the true residual. This usually indicates that the system is near singular.</div><div><br></div><div>3) The disparity above does not seem possible if C only has elements ~ 1e4. The preconditioner consistently has norm around 1e-11.</div><div><br></div><div>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</div><div>     doing iterative refinement with a very low precision solve.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>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. </div><div><br></div><div>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?</div><div><br></div><div>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. </div><div><br></div><div>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.)</div><div><br></div><div><br></div><div>Thanks,</div><div>Ramakrishnan Thirumalaisamy</div><div>San Diego State University.</div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>