Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function? I can try to use it to test. Thanks a lot, Jed.
Best,
Yujie
Can you send a reproducible test case? (Instructions for how to run a Libmesh example that would cause this problem is okay, producing it with a PETSc example would be better.)