Spectral Model for
Thermal Convection with
Free-Slip, Constant Temperature
Boundaries

Solves the equations for the spectral representation of thermal convection in a layer.  Certain simplifying symmetry restrictions are made on the convection.
Plots the results.

This note book should be executed before the subsequent notebooks are used.

Constants

(*The critical Rayleigh number*)
Rc=27*Pi^4/4

Subroutines

(*Dirac delta function*)
	id[m_,n_]:=If[m==n,1,0]

(*nonlinear interaction coefficient for buoyancy*)
		ct[mi_,ni_,mj_,nj_,mk_,nk_]:=
	Module[
		{norm},
				norm=If[mi==0,.125,.25];
	norm*(
		mj*nk*
		(id[mj+mk,mi]+id[mj-mk,mi]+id[mk-mj,mi])
*(id[nj+nk,ni]+id[nj-nk,ni]-id[nk-nj,ni])
+mk*nj*
		(-id[mj+mk,mi]+id[mj-mk,mi]+id[mk-mj,mi])
 *(id[nj+nk,ni]-id[nj-nk,ni]+id[nk-nj,ni])
		)
	]

(*nonlinear interaction coefficient for stream function*)cp[mi_,ni_,mj_,nj_,mk_,nk_]:=Module[{},
		.25*(
	mj*nk*
(id[mj+mk,mi]-id[mj-mk,mi]+id[mk-mj,mi])
*(id[nj+nk,ni]+id[nj-nk,ni]-id[nk-nj,ni])
-mk*nj*
(id[mj+mk,mi]+id[mj-mk,mi]-id[mk-mj,mi])
*(id[nj+nk,ni]-id[nj-nk,ni]+id[nk-nj,ni])
)
]

(*individual spectral adection term for buoyancy*)
advt[i_,j_,k_]:=p[j][t]*T[k][t]*
	ct[
			tmodes[[i,1]],tmodes[[i,2]],
			pmodes[[j,1]],pmodes[[j,2]],
			tmodes[[k,1]],tmodes[[k,2]]
		]

(*individual spectral adection term for vorticity*)advom[i_,j_,k_]:=-p[j][t]*p[k][t]*
		(a^2*pmodes[[k,1]]^2+b^2*pmodes[[k,2]]^2)*
	cp[
			pmodes[[i,1]],pmodes[[i,2]],
			pmodes[[j,1]],pmodes[[j,2]],
			pmodes[[k,1]],pmodes[[k,2]]
		]

(*buoyancy production term , (the one with the Rayleigh number)*)ray[i_,j_]:=p[j][t]*a*pmodes[[j,1]]*
		id[tmodes[[i,1]],pmodes[[j,1]]]*
		id[tmodes[[i,2]],pmodes[[j,2]]]

(*baroclinic production of vorticity*)baro[i_,j_]:=T[j][t]*a*tmodes[[j,1]]*
		id[pmodes[[i,1]],tmodes[[j,1]]]*
		id[pmodes[[i,2]],tmodes[[j,2]]]

(*sum of buoyancy advection terms*)
advectt[i_]:=Sum[
		Sum[advt[i,j,k],
			{j,1,nump}],
			{k,1,numt}]

(*sum of vorticity advection terms*)
advectom[i_]:=Sum[
		Sum[advom[i,j,k],
			{j,1,nump}],
			{k,1,nump}]

(*sum of buoyancy production terms*)
rayt[i_]:=Sum[
		ray[i,j],
			{j,1,nump}]

(*sum of baroclinc generation terms*)
barop[i_]:=Sum[
	baro[i,j],
			{j,1,numt}]

(*defines some global variables based on set of modes tm and pm*)
setup[tm_,pm_]:=Module[{},
		tmodes=tm;
		pmodes=pm;
			numt=Dimensions[tmodes][[1]];
		nump=Dimensions[pmodes][[1]];
		]

whatis[q_,t_,sol_]:=q[t]/.sol[[1,1]]

Programs (distinction between “programs” and “subroutines” is somewhat arbitrary)

(*A model is a list of things, see below...*)
makemodel[tm_,pm_]:=Module[{},
		setup[tm,pm];
	{makeode,makevar,tmodes,pmodes}
	]

(*Makes a list of ODE's compatible with NDSolve*)
makeode:=
			Join[Table[
		T[i]'[t]
				+(a^2*tmodes[[i,1]]^2+
				b^2*tmodes[[i,2]]^2)T[i][t]
						==
			-a*b*advectt[i]
		+rayt[i]
		,
				{i,1,numt}],
			Table[
		p[i]'[t]
				+P(a^2*pmodes[[i,1]]^2+
				b^2*pmodes[[i,2]]^2)p[i][t]
						==
			(-a*b*advectom[i]
		+P*R*barop[i])
		/(a^2*pmodes[[i,1]]^2+
				b^2*pmodes[[i,2]]^2)
		,
				{i,1,nump}]
		]

(*Makes a list of initialization equations compatible with NDSolve*)makeinit[themodel_,m_,n_,amp_]:=Module[{},
		setup[themodel[[3]],themodel[[4]]];
		Join[
				Table[
		T[i][0]==amp*id[tmodes[[i,1]],m]*id[tmodes[[i,2]],n],
					{i,1,numt}],
			Table[
		p[i][0]==0,
				{i,1,nump}]
	]
		]

(*Makes a list of the variables that will be solved for in NDSolve*)
makevar:=
		Join[
				Table[
		T[i],
					{i,1,numt}],
			Table[
		p[i],
				{i,1,nump}]
	]

(*integrates a model forward in time using  NDSolve*)runmodel[themodel_,init_,parms_,tstop_]:=
	Module[{},
		setup[themodel[[3]],themodel[[4]]];
		{NDSolve[Join[themodel[[1]]/.parms,init],themodel[[2]],{t,0,tstop},MaxSteps->100000],
			parms,themodel[[3]],themodel[[4]],maketfun,makepfun}
		]

(*makes a function b(x,z) that can be plotted*)
maketfun:=
	1-z*b/Pi+
		Sum[
			T[i][t]*Cos[a*tmodes[[i,1]]*x]
							*Sin[b*tmodes[[i,2]]*z],
							{i,1,numt}
						]

(*makes a function psi(x,z) that can bi plotted *)
makepfun:=
		Sum[
			p[i][t]*Sin[a*tmodes[[i,1]]*x]
							*Sin[b*tmodes[[i,2]]*z],
							{i,1,nump}
						]

Plotting Programs

What do these do?  They are best explained by seeing them in action in subsequent notebooks.

modeplot[tlist_,plist_,mmax_,nmax_]:=Module[{backpoints,rpoint},backpoints={};rpoint[q_]:=Text[Style["\[Cross]",FontSize->24],q];Do[Do[backpoints=Append[backpoints,{i,j}],{i,0,mmax,1}],{j,1,nmax,1}];theplot=ListPlot[{{0,0}},PlotRange->{{-1,mmax},{0,nmax}},AspectRatio->(nmax/(mmax+1)),AxesOrigin->{-1,0},Epilog->{{RGBColor[1,0,0],Point/@backpoints},{RGBColor[0.7`,0.7`,1],AbsolutePointSize[9],Point/@tlist},{RGBColor[0,0,0],rpoint/@plist}}]]

specplot[sol_,time_,what_,size_,color_,mmax_,nmax_]:=Module[{hack,splot,alltheblobs},splot[q_]:={AbsolutePointSize[size (Sqrt[Abs[q[[3]]]])],Point[{q[[1]],q[[2]]}]};If[what=="p",hack=Table[Flatten[{sol[[4,i]],whatis[p[i],time,sol]}],{i,1,nump}];,hack=Table[Flatten[{sol[[3,i]],whatis[T[i],time,sol]}],{i,1,numt}];];alltheblobs=splot/@hack;theplot=ListPlot[{{0,0}},PlotRange->{{-1,mmax},{0,nmax}},AspectRatio->(nmax/(mmax+1)),AxesOrigin->{-1,0},Epilog->{{color,alltheblobs}}]]

BlueGrayRed[z_]:=RGBColor[z,2*z*(1-z),1-z]

Grayy[z_]:=RGBColor[.3+.6*z,.3+.6*z,.3+.6*z]

MyTable=Grayy;

conp[sol_,time_]:=Module[{thefunc,theratio},thefunc=N[sol[[6]]/. sol[[1,1]]/. t->time/. sol[[2]]];theratio=(a/b)/. sol[[2]];theplot=ContourPlot[thefunc,{x,0,(π/a)/. sol[[2]]},{z,0,(π/b)/. sol[[2]]},ColorFunction->MyTable,ContourStyle->Automatic,AspectRatio->theratio,Frame->False]]

cont[sol_,time_]:=Module[{thefunc,theratio},thefunc=N[sol[[5]]/. sol[[1,1]]/. t->time/. sol[[2]]];theratio=(a/b)/. sol[[2]];theplot=ContourPlot[thefunc,{x,0,(π/a)/. sol[[2]]},{z,0,(π/b)/. sol[[2]]},ColorFunction->MyTable,ContourStyle->Automatic,AspectRatio->theratio,Frame->False,Contours->{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}]]

?Plot

hist[sol_,time1_,time2_,color_]:=Module[{thefunc},thefunc=T[1][t]/. sol[[1,1]];theplot=Plot[thefunc,{t,time1,time2},PlotRange->{-0.5,0.5},PlotStyle->{{color}},PerformanceGoal->"Quality",PlotPoints->(time2-time1) 100]]

butterfly[sol_,time1_,time2_,color_]:=Module[{tfunc,t2func},tfunc=T[1][t]/. sol[[1,1]];t2func=T[2][t]/. sol[[1,1]];theplot=ParametricPlot[{tfunc,t2func},{t,time1,time2},PlotPoints->(time2-time1) 100,PlotStyle->{{color}},PlotRange->{{-0.5,0.5},{-0.7,0.1}}]]