> On Sep 30, 2021, at 6:16 PM, Amneet Bhalla <mail2amn...@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 > > >
