In[]:=
SetAttributes[compile,HoldAll];​​compile[argu__]:=​​With[{cg=Compile`GetElement},​​Hold@Compile[argu,RuntimeOptions"Speed"(*,CompilationTarget"C"*)]/.​​Partcg/.HoldPattern[cg[a__]=rhs_](Part[a]=rhs)]//ReleaseHold​​​​vortex=compile[{{v,_Real,3}},​​Module[{nx,ny},​​{nx,ny}=Dimensions@v//Most;​​Table[((v[[i+1,j,2]]-v[[i-1,j,2]])-(v[[i,j+1,1]]-​​v[[i,j-1,1]]))/2,​​{j,2,ny-1},{i,2,nx-1}]]];
In[]:=
force…and…advection=​​With[{inter=Function[{valueL,valueR,scale},(1-scale)valueL+scalevalueR]},​​compile[{{arg,_Real,3}(*,{nx,_Integer},{ny,_Integer}*),{dt,_Real}},​​Module[{v=arg,nx,ny,inew,jnew,iL,jL,iR,jR},​​{nx,ny}=Dimensions@v//Most;​​Do[v[[1,j]]={.1,0.},{j,ny}];​​Do[(*If[i<1+nx/16,v[[i,j]]={.1,0.}]*)​​If[(i-nx/4)^2+(j-ny/2)^2<(nx/16)^2,v[[i,j]]={0.,0.}],{i,nx},{j,​​ny}];​​Table[​​(*{inew,jnew}=Mod[{i,j}-nxdtv[[i,j]],{nx,ny}-1,1];*)​​inew=Mod[i-nxdtv[[i,j,1]],nx-1,1];​​jnew=Mod[j-nxdtv[[i,j,2]],ny-1,1];​​(*{iL,jL}=Floor@{inew,jnew};​​{iR,jR}=Ceiling@{inew,jnew};*)​​iL=Floor@inew;jL=Floor@jnew;​​iR=iL+1(*Ceiling@inew*);jR=jL+1(*Ceiling@jnew*);​​Table[​​inter[​​inter[v[[iL,jL,index]],v[[iL,jR,index]],jnew-jL],​​inter[v[[iR,jL,index]],v[[iR,jR,index]],jnew-jL],​​inew-iL],​​{index,2}]​​,{i,nx},{j,ny}]​​]]];​​​​viscosity…and…conservation=​​compile[{{v,_Complex,3}(*,{nx,_Integer},{ny,_Integer}*),dt,mu},​​Module[{x,y,nx,ny,k},​​{nx,ny}=Dimensions@v//Most;​​Table[(*{x,y}=Mod[{i,j}-1,{nx,ny},-{nx,ny}/2];*)​​(*x=Mod[i-1+nx/2,nx]-nx/2;​​y=Mod[j-1+ny/2,ny]-ny/2;*)​​x=Mod[i-1,nx,-nx/2];​​y=Mod[j-1,ny,-ny/2];​​k=x^2+y^2;​​Exp[-kdtmu]If[k>0,​​With[{proj=(-yv[[i,j,1]]+xv[[i,j,2]])/k(*v[[i,j]].{-y,x}/​​k*)},{-yproj,xproj}],v[[i,j]]](*If[k>0,(v[[i,j]].{-y,x}/k){-y,x},v[[​​i,j]]]*),{i,nx},{j,ny}]]];
In[]:=
(*nx=128;ny=nx/2;dt=1.;mu=0.0005;nt=250;*)​​nx=64(*128*);ny=nx/2;dt=0.3;mu=0.001;nt=800;​​(*cut=nx/8;*)​​v=Table[{0.,0.},{nx},{ny}];​​separate=Transpose[#,{2,3,1}]&;​​combine=Transpose[#,{3,1,2}]&;​​​​vlst=Table[​​(*v=add…force[v];​​v=advection[v,dt];*)​​v=force…and…advection[v(*,nx,ny*),dt];​​v=Fourier/@separate@v//combine;​​v=viscosity…and…conservation[v(*,nx,ny*),dt,mu];​​v=InverseFourier/@separate@v//combine//Re,{t,nt}];//AbsoluteTiming​​(*{24.4562,Null}*)​​disk=Graphics[{Opacity[1/2],Disk[{nx/4(*-cut*),ny/2},nx/16]}];​​arrayplot=​​Show[ArrayPlot[#,DataReversedTrue,ColorFunction"Rainbow",​​ImageSizenx,PlotRange{0,1},AxesTrue],disk]&;​​​​piclst=arrayplot/@​​Rescale[vortex(*@#[[cut+1;;]]&*)/@vlst[[;;;;10]]];//AbsoluteTiming
Out[]=
{2.94563,Null}
Out[]=
{1.08019,Null}
In[]:=
piclst//ListAnimate
In[]:=
MaxMemoryUsed[]/1024^2."MB"
Out[]=
407.852MB