Does phi = project( ((c3**2)/((c3**2)+(thr**2))) , Vphi) provide a good approximation in your application?
(If there are sharp layers in the quantity to be projected, the L^2 projection onto CG1 may have some oscillations, but you could instead try a lumped-mass projection, as demonstrated here, which also avoids doing a linear solve.)