After I set a=pi/4<pi/2, the plots look like what I expected to see. Thanks again.

Thank you so much!

Do you need a boundary condition like $F(\infty)=0$?

Hi Alex: I noticed an interesting difference between Thomas numerical solution in Ref.1 and your solution in the answer. If you let $z1=4*10^{-3}} and do Plot[u[40][x,z1],{x,-0.02,0.02}], you will see that u[t(40)][x=0,z1]>1. But from inspection of Figure (3.1) in Ref.1, you will see that |u[t(40)],x=0,z1)|<1/10. This also true for w. As you now that (Thomas)u1==(your)u, (Thomas)omega1==(your)w.

In the answer to question "mathematica.stackexchange.com/questions/188101/…", @xzczd mentioned that "We can use r=0 as the left boundary because "FiniteElement" is able to handle the removable singularity there". This probably explained why FEM worked here as well.

Alex: How can I improve the quality of your NDSolve/FEM solution. (1) set "MaxCellMeasure" to 0.0001 or even smaller; reduce dt to 1/10 of its original value and increase the total time steps tp 400. Can I set "InterpolationOrder" to 4? Best regards-

I also have question about DirichletCondition[{w[x, z] == 0, f[x, z] == 0, u[x, z] == 0}, True]. Is this an initial condition? Why do we not need boundary conditions?

You are right. Since doubling the domain in $x$ will not automatically reduce the error near $x=0$, I go back to original domain in Hou's paper $D=\{(x,z):0\leq x\leq 1,0\leq z\leq 1/2\}$. The solution looks similarly nice. $1/x$ term will not cause trouble.

I see. These are the lower-left and up-right corners. Have you checked plots like "Plot[U[40][x, 0], {x, -0.05, 0.05}, PlotRange -> All]" and "Plot[U[40][x, 0], {x, -0.05, 0.05}, PlotRange -> All]"? Why are they anti-symmetric in x?

The thing people want to now is whether the singularity will form on the z axis (x=0) in finite time. Thus the thing that is less ideal in Hou's paper, in my opinion, is that x=0 is the boundary of their simulation. This may introduce large error. In my equations above, (x=0,z=0) is in the middle of the domain $D_2={(x,z): -1<=x<=1,-1/2<=z<=1/2}$, I was hoping to see that by putting the potential singular region in the middle of the simulation domain, we would see more accurate behavior of the solutions (finite time blowup singlarity formation)

Very impressive! I would suggest that you publish it. I do not quite understand that "Rectangle[{-1, -1/2}, {1, 1/2}]". The solutions are symmetric in $x$ and anti-symmetric in $z$. I can follow you if you solve the PDEs for $-1<=x<=0,0<=z<=1/2$ and then symmetrize it wrt $x$ and anti-symmetrized wrt $z$. But why "Rectangle[{-1, -1/2}, {1, 1/2}]"? What do is mean?

@AlexTrounev: Their adaptive mesh method looks too complicated to me. I hope to take short cut with NDSolve. :-)

@MichaelE2: I tried. It does not help. Thanks again!