what is the proper way to solve it efficiently?
What are you using currently? And why is not giving you appropriate performance?
will PETSc reuse the inversion of the matrix?
Yes.
Unfortunately, I got a blow up for my new code, every variable looks weird even in the first time step.
I am afraid that my boundary condition is imposed in a wrong way. It seems that constructing rhs by
varf rhs(u,v)=int2d(Th)(uold*v)+on(1,u=0); real[int] b=rhs(0,Vh,tgv=-1);
is wrong. 
Maybe a mass matrix multiply uold M*uold[] is correct, I just thought they were the same thing.
Actually I am not sure with the time derivative term discretized by finite element.
Anyway, I will figure out what’s going wrong as soon as possible.