I am experiencing poor residuals with my iterative solver in PETSc when solving the elasticity equation. I normally use -pc_gamg and -ksp_gmres, and I have not had issues with solving 3D meshes when my Young’s modulus was an arbitrary normalized value. Now, however, I cannot seem to achieve good residuals using realistic properties (MPa) and under an applied traction (kN) and steady state thermal load. I can only achieve good residuals when my load is much smaller (small traction and no thermal), but that does not satisfy the goal of my research. From my understanding, there is no apparent singularity in the [K] matrix, and finite element analysis should be feasible to solve given the units of the problem.
I have attached an example script with my formulation. I would appreciate any guidance on how I can solve the system of elasticity using an appropriate iterative solver in PETSc.