I think I have I found the flaw in my methods. And it stems from my failure to do a proper dimensional analysis.
The units of my source terms were in volumetric heat generation (W/m^3) whereas all differential terms have units of (K/s). To correct this, I had to divide the source terms (including the boundary value constant flux) by (ro * Cp). This corrected my originally strange results, and the time evolution of the problem now seems to follow an appropriate scale. My apologies for bombarding the Fipy list with this avoidable problem, but I always like to post the resolution for future googlers. Perhaps this example might make for a nice tutorial addition in the Fipy manual. I would be willing to write it up if you all are willing to accept it. Also, if I still have glaring holes in my logic or method, please let me know! Cheers, Chris
