Initialization and write simulation information

Clear["Global`*"]
Needs["pc`"]
Get["saffman`"]
pcInitialize[]

dir="/Users/hzhou/Work/pencilCode/runs/hongzhe_pc_dardel/saffman/";
outputDir="/Users/hzhou/Dropbox/Apps/Overleaf/Saffman_helicity_invariant/";
runs=FileNames["k*",dir];

label=AssociationThread[runs->(runs//StringDrop[#,60]&//StringReplace[#,
{"k"->"K","del2"->"D1","del4"->"D2","del6"->"D3",
"ac"->"c","a"->"t"}]&)];
Print["Loaded runs:\n",TableForm[MapThread[StringDrop[#1,52]<>" --> "<>#2&,{runs,label/@runs}]]]
Do[group["D"<>n<>sp]=Cases[runs,x_String/;StringMatchQ[label[x],___~~"D"<>n<>sp]],
{n,ToString/@{2,4,6}},{sp,{"","c"}}]

Read and compute everything, except (Subscript[p, {H,UB}]) which needs (Subscript[r, 0])

Do[
Nx=getParam[sim,"Nx"];
Lx=(2π)/(2π);
dx=Lx/Nx//N;

(*read spectral+cube method*)
{t[sim,1],Ivspc[sim]}=read1D[sim,"Iv_spc_saffman_mag.dat"];
Ivspc[sim]=Map[Transpose[{dx*Range[Nx/2],#}]&,Ivspc[sim]];

(*read power spectrum of h*)
{t[sim,2],sph[sim]}=read1D[sim,"power_saffman_mag.dat"];
sph[sim]=Map[Transpose[{Range[0,Nx/2-1],#}]&,sph[sim]];

(*read helicity spectrum A.B*)
{t[sim,3],spAB[sim]}=read1D[sim,"powerhel_mag.dat"];

(*read power spectrum of B*)
{t[sim,4],spB[sim]}=read1D[sim,"power_mag.dat"];

(*make sure time is correct*)
n=Min[Length[t[sim,#]]&/@Range[4]];
{t[sim,1],Ivspc[sim],t[sim,2],sph[sim],t[sim,3],spAB[sim],t[sim,4],spB[sim]}=
Take[#,n]&/@{t[sim,1],Ivspc[sim],t[sim,2],sph[sim],t[sim,3],spAB[sim],t[sim,4],spB[sim]};
t[sim]=t[sim,1];

(*Lundquist number*)
Lu[sim]=LuNspec[sim,1,n,"l1/n"->False];
Lu2[sim]=LuNspec[sim,1,n,"l1/n"->True];

(*compute decay exponents*)
Em[sim]=Total/@spB[sim];
pE[sim]=pDecay[t[sim],Em[sim]];
pAB[sim]=pDecay[t[sim],Total/@spAB[sim]];
xim[sim]=lB[sim,t[sim],spB[sim]];
q[sim]=pDecay[t[sim],xim[sim],-1];

(*color scale with time*)
hue[sim]=Directive[AbsoluteThickness[1],Hue[0.9#]]&/@Rescale[t[sim]],

{sim,runs}];

ClearAll[round,sciForm,mathForm]
round[x_]/;x<1:=DecimalForm[x,1]
round[x_]/;x>=1:=Round[x]
Attributes[sciForm]={Listable};
sciForm[x_]:=x//pcNumberToTex[#,2]&//StringTake[#,{2,-2}]&//StringReplace[#,"10^1"->"10"]&//StringReplace[#,StartOfString~~xx_~~"\\times":>StringJoin[xx,".0\\times"]]&
mathForm[x_]:=StringJoin["$",ToString@x,"$"]

tb=Table[{
label[sim],"$"<>ToString@readDim[sim]["nx"]<>"^3$",
read1D[sim,"power_mag.dat"][[1,2]]//pcNumberToTex,
readParamNml[sim,"start.in","KPEAK_AA"]//Round//mathForm,
Max[readParamNml[sim,"run.in","ETA"],readParamNml[sim,"run.in","ETA_HYPER3"]]//pcNumberToTex,
readParamNml[sim,"run.in","IRESISTIVITY"][[1]]/.{"'hyper3'"->3,"'hyper3-tdep'"->3,"'hyper2-tdep'"->2,"'eta-const'"->1,"'eta-tdep'"->1}//mathForm,
readParamNml[sim,"run.in","ETA_TDEP_EXPONENT"]//Round[#,0.01]&//#/.{-0.43->"$-3/7$",0.->"$0$"}&,
Lu[sim][[;;,2]]//DeleteCases[#,Indeterminate]&//MinMax//sciForm//StringJoin["$",#[[1]],"\\to",#[[2]],"$"]&
},{sim,runs}]//Quiet;

head= {"Run","$N^3$","$\\EEM(0)$","$k_\\text{peak}$","$\\nu_n=\\eta_n$","$n$","$r$","$\\Lu_n$"};
tb//TableForm[#,TableHeadings->{None,head}]&
pcWriteTexTable[outputDir<>"runs.tex",head,tb,"lAlignment"->"l"]

Find the asymptotic range and check spectrum

Table[{
infoMatrix[sim,"folder"],
ListLogLogPlot[Ivspc[sim],PlotStyle->hue[sim],
GridLines->{{r0=0.115},None},FrameLabel->{"R/2π","(Subscript[I, H])"}
],
With[{spec=read1D[sim,"power_mag.dat"]},
ListLogLogPlot[spec[[2,;;,2;;]],PlotStyle->(Hue[0.9#]&/@Rescale[spec[[1]]]),FrameLabel->{"k","(Subscript[ℰ, M])"},
GridLines->{{{1/(2r0),Black}},None}
]]
},{sim,runs}]//Grid

sim="/Users/hzhou/Work/pencilCode/runs/hongzhe_pc_dardel/saffman/k60del2b"//Echo[#,"Run: ",label]&;
xticks=Table[{10^i,Superscript["10",ToString@i]},{i,-20,10}];
yticks=Table[{10^i,ToString@i},{i,-20,10}];
style={PlotStyle->hue[sim],GridLines->{{{Echo[1/(2r0),"2π/(2R)=",N],Black}},None},
FrameTicks->{{yticks,Automatic},{xticks,Automatic}},ImageSize->{imgSz=300},
ImagePadding->{{60,10},{50,30}}};
{
(*a*)
Show[ListLogLogPlot[Rest/@spB[sim],FrameLabel->{"k/(Subscript[k, 1])",Subscript["(Subscript[Log, 10]) E","M"]},PlotRange->{All,{All,10^(-1.7)}},style],
ListLogLogPlot[Table[{kk,1*^-10 *kk^4},{kk,3,12}],PlotStyle->Directive[Black,Dashed]],
ListLogLogPlot[Table[{kk,1*^1 *kk^(-5/3)},{kk,150,500,50}],PlotStyle->Directive[Black,Dashed]],
Epilog->{Inset[Style["(k^4)",pcLabelStyle],Scaled[{0.36,0.35}]],
Inset[Style["(k^(-5/3))",pcLabelStyle],Scaled[{0.87,0.91}]],
Inset[Style["(a)",pcLabelStyle],Scaled[{0.08,0.9}]]}
],

(*b*)
Show[ListLogLogPlot[Rest/@sph[sim],FrameLabel->{"k/(Subscript[k, 1])","(Subscript[Log, 10]) Sp(h)"},PlotRange->{All,{All,1*^-6}},style],
ListLogLogPlot[Table[{kk,2.9*^-10 *kk^2},{kk,1,5}],PlotStyle->Directive[Black,Dashed]],
ListLogLogPlot[Table[{kk,1.1*^-7 *(kk/150)^(-5/3)},{kk,150,480,50}],PlotStyle->Directive[Black,Dashed]],
Epilog->{Inset[Style["(k^2)",pcLabelStyle],Scaled[{0.18,0.72}]],
Inset[Style["(k^(-5/3))",pcLabelStyle],Scaled[{0.87,0.9}]],
Inset[Style["(b)",pcLabelStyle],Scaled[{0.08,0.9}]]}
],

(*legend*)
DensityPlot[y,{x,-1,1},{y,Sequence@@MinMax[t[sim]]},ColorFunction->(Hue[0.9#]&),
PlotLegends->False,ImageSize->{imgSz/3,imgSz}/1.5,AspectRatio->8,
PlotRangePadding->None,PlotLabel->Style["t",pcLabelStyle],
FrameTicks->{{None,Subdivide[Sequence@@MinMax[t[sim]],4]//Round},{None,None}},
ImagePadding->{{10,30},{30,10}},FrameStyle->None
],

(*c*)
norm=20480/969*Em[sim]^2*xim[sim]^5;
ListLogLogPlot[MapThread[
((Rest[#1])/(#2*Range[getParam[sim,"Nx"]/2-1]^2))&,
{sph[sim][[;;,;;,2]],norm}
],FrameLabel->{"k/(Subscript[k, 1])","Compensated Sp(h)"},PlotRange->{0.001,20},
Epilog->{Inset[Style["(c)",pcLabelStyle],Scaled[{0.08,0.1}]]},
FrameTicks->{{xticks,Automatic},{xticks,Automatic}},
style],

(*d*)
ListLogLogPlot[
pos=Nearest[Ivspc[sim][[1,;;,1]]->"Index",r0][[1]];
Sort@Transpose[{t[sim],Ivspc[sim][[;;,pos,2]]/(π*norm)}],
FrameLabel->{"t","(Subscript[(OverTilde[I]), H])"},
Epilog->{Inset[Style["(d)",pcLabelStyle],Scaled[{0.08,0.9}]]},
FrameTicks->{{Automatic,Automatic},{xticks,Automatic}},
GridLines->None,PlotStyle->Directive[Black,AbsoluteThickness[1]],
style
]
}//Partition[#,UpTo@3]&//Grid[#,Spacings->{{0,1,1,0},0},Alignment->Center]&//Echo//
Export[outputDir<>"/figures/spectra.eps",#]&

Time series, pq, and (pH,q)

Get["saffman`"]
Table[{
infoMatrix[sim,"folder"],
pos=Nearest[Ivspc[sim][[1,;;,1]]->"Index",r0][[1]];
ListLogLogPlot[tmp=Transpose[{t[sim],Ivspc[sim][[;;,pos,2]]}],FrameLabel->{"t","(Subscript[I, H])"},ImagePadding->{{100,10},{100,10}},ImageSize->400],
pH[sim]=pDecay[Sequence@@Transpose[tmp]];
pqDiag[sim,pE[sim],q[sim],Most@t[sim],{FrameLabel->{"q","p"},PlotLabel->label[sim]}],
pqDiag[sim,pH[sim],q[sim],Most@t[sim],{FrameLabel->{"q","(Subscript[p, H])"}},"lnotPQ"->True]
},{sim,runs}]//Grid

With[{runs=Flatten[(Position[label,#]&/@{"K200D3t","K60D1c","K60D1b","K60D1bc"})/.Key->List]},
Print["Plotting runs: ",label/@runs];
Get["saffman`"];
panelLabel=AssociationThread[runs->{"(a)","(b)","(c)","(d)"}];
Table[
pqDiag[sim,pE[sim],q[sim],Most@t[sim],{FrameLabel->{"q","p"}},
"lPanelLabel"->panelLabel[sim]
],{sim,runs}]//Partition[#,2]&//Grid
]//Echo//Export[outputDir<>"/figures/pq.eps",#]&

(p,pH) vs Lu

style={ImagePadding->{{60,10},{50,10}},ImageSize->350};
markers=Charting`CommonDump`GraphicsOpenPlotMarkers[][[1;;2]];

(*rainbow color*)
colors0=Hue[0.9#]&/@Subdivide[0,1,Length[runs]-1];

(*same color in group*)
(*colors0=With[{l=label[#]},
Apply[Which,Flatten[
Transpose[{
Flatten@Table[StringMatchQ[l,"*D"<>sp1<>sp2],{sp1,{"1","2","3"}},{sp2,{"t","c"}}],
Hue[0.8#]&/@Subdivide[0,1,5]
}]
]]]&/@runs;*)


colors=ConstantArray[colors0,2]//Transpose//Flatten;
xLabel=Table[{10^i,Superscript["10",ToString[i]]},{i,Range[-20,20]}];

dataP=Table[{Lu[sim][[1;;-2,2]],pE[sim]}//Transpose,{sim,runs}];
dataPAB=Table[{Lu[sim][[1;;-2,2]],pAB[sim]}//Transpose,{sim,runs}];
dataPH=Table[{Lu[sim][[1;;-2,2]],pH[sim]}//Transpose,{sim,runs}];

Clear[normalize]
normalize[l_]:=Module[{x,y},{x,y}=Transpose[l];Transpose[{x,y/y[[1]]}]]

Clear[panelLabel]
panelLabel[x_]:=Inset[Style["("<>x<>")",pcLabelStyle],Scaled@{0.05,0.9}]
{
(*a*)
ListLogLogPlot[
Table[pos=Nearest[Ivspc[sim][[1,;;,1]]->"Index",r0][[1]];normalize@Transpose[{t[sim],Ivspc[sim][[;;,pos,2]]}],{sim,runs}],
FrameLabel->{"t","(Subscript[I, H])/(Subscript[I, H])(0)"},PlotStyle->colors0,
FrameTicks->{{Automatic,Automatic},{xLabel,Automatic}},
PlotLegends->Placed[LineLegend[label/@runs,LegendLayout->{"Column",2}],{Left,Bottom}],
Epilog->{panelLabel["a"]},
style
],

(*b*)
Show[ListLogLogPlot[Riffle[dataP,dataPH],Joined->False,FrameLabel->{"(Subscript[Lu, n])","p and (Subscript[p, H])"},
PlotMarkers->markers,PlotStyle->colors,
PlotLegends->Placed[PointLegend[{Black,Black},{"p","(Subscript[p, H])"},LegendFunction->(Framed[#1, FrameMargins -> 0] & )],{Left,Bottom}],
PlotRange->{Automatic,2},
FrameTicks->{{Table[{10^i,Superscript["10",ToString[i]]},{i,Range[-10,0]}],Automatic},{xLabel,Automatic}},
style
],
Plot[0.8x^(-0.25),{x,15,10^8},PlotStyle->Directive[Black,Dashed],ScalingFunctions->{"Log","Log"}],
Epilog->{Inset[Style[Superscript["0.8(Subscript[Lu, n])","-1/4"],pcLabelStyle],Scaled[{0.68,0.48}]],
panelLabel["b"]}
],

(*c*)
Show[ListLogLogPlot[Riffle[dataPAB,dataPH],Joined->False,FrameLabel->{"(Subscript[Lu, n])","(Subscript[p, AB]) and (Subscript[p, H])"},
PlotMarkers->markers,PlotStyle->colors,
PlotLegends->Placed[PointLegend[{Black,Black},{"(Subscript[p, AB])","(Subscript[p, H])"},LegendFunction->(Framed[#1, FrameMargins -> 0] & )],{Left,Bottom}],
PlotRange->{Automatic,2},
FrameTicks->{{Table[{10^i,Superscript["10",ToString[i]]},{i,Range[-10,0]}],Automatic},{xLabel,Automatic}},
style
],
Plot[0.8x^(-0.25),{x,15,10^8},PlotStyle->Directive[Black,Dashed],ScalingFunctions->{"Log","Log"}],
Epilog->{Inset[Style[Superscript["0.8(Subscript[Lu, n])","-1/4"],pcLabelStyle],Scaled[{0.68,0.48}]],
panelLabel["c"]}
],

(*d*)
Show[
With[{sim="/Users/hzhou/Work/pencilCode/runs/hongzhe_pc_dardel/saffman/k60del6bc"//Echo[#,"(d): Run ",label]&},
ListLogLogPlot[Ivspc[sim],PlotStyle->hue[sim],
GridLines->{{{r0//Echo[#,"(Subscript[R, 0])/2π="]&,Directive[AbsoluteThickness[1],Black]}},None},FrameLabel->{"R/2π","(Subscript[I, H])"},
FrameTicks->{{Table[{10^i,Superscript["10",ToString[i]]},{i,Range[-30,0]}],Automatic},
{Table[{5*10^i,Superscript["5\[Cross]10",ToString[i]]},{i,Range[-30,0]}],Automatic}},
style
]],
Plot[5*10^(-13)*(x/0.0002)^3,{x,0.0005,0.0015},PlotStyle->Directive[Black,Dashed],ScalingFunctions->{"Log","Log"}],
Epilog->{(*Inset[
pos=Nearest[Ivspc[sim][[1,;;,1]]->"Index",r0][[1]];
ListLogLogPlot[
Transpose[{t[sim],10^9*Ivspc[sim][[;;,pos,2]]}],FrameLabel->{"t","(10^9)(Subscript[I, H])"},ImagePadding->{{50,10},{40,10}},ImageSize->200,
FrameTicks->{{Automatic,Automatic},{xLabel,Automatic}}
],
Scaled[{0.64,0.32}]
],*)panelLabel["d"],
Inset[Style["(R^3)",pcLabelStyle],Scaled@{0.12,0.68}]}
]
}//Partition[#,UpTo@2]&//Grid[#,Spacings->{3, 1}]&//Echo//Export[outputDir<>"/figures/Lu.eps",#]&

ListLogLogPlot[Lu2/@runs,ImageSize->Large,PlotLegends->(label/@runs),
FrameLabel->{"t","(Subscript[Lu, n])"},
PlotStyle->(Hue[0.9#]&/@Subdivide[0,1,Length[runs]-1])]

Do[
(*length scale*)
l[sim]=lB[sim,t[sim],spB[sim]];
(*total helicity*)
H[sim]=Total/@spAB[sim];
(*total energy*)
e[sim]=2Total/@spB[sim];
(*h^2*)
eh2[sim]=Total/@sph[sim][[;;,;;,2]];
,
{sim,runs}]

Clear[normPlot]
normPlot[norm_,yLabel_]:=ListLogLogPlot[
Table[pos=Nearest[Ivspc[sim][[1,;;,1]]->"Index",r0][[1]];
Transpose[{t[sim],(Ivspc[sim][[;;,pos,2]])/norm[sim]}],{sim,runs}],
FrameLabel->{"t",yLabel},PlotStyle->colors0,
PlotLegends->Placed[LineLegend[label/@runs,LegendLayout->{"Column",3}],{Left,Bottom}],
ImagePadding->{{100,5},{50,5}},
style
]
{V=(2π)^3;
normPlot[(eh2[#]*l[#]^3)&,"(Subscript[I, H])/(⟨(h^2)⟩*(L^3))"],
normPlot[(H[#]^2*l[#]^3)&,"(Subscript[I, H])/(⟨h⟩^2*(L^3))"],
normPlot[(e[#]^2*l[#]^5)&,"(Subscript[I, H])/(⟨(B^2)⟩^2*(L^5))"]
}//List//Grid[#,Spacings->{1, 1}]&

Gauge dependence

sim="/Users/hzhou/Work/pencilCode/runs/hongzhe_pc_dardel/saffman/k60del2ac"//Echo[#,"Run: ",label]&;

Nx=getParam[sim,"Nx"];
Lx=(2π)/(2π);
dx=Lx/Nx//N;
(*read spectral+cube method in Coulomb gauge*)
{t[sim],IvspcC[sim]}=read1D[sim,"Iv_spc_saffman_mag_c.dat"];
IvspcC[sim]=Map[Transpose[{dx*Range[Nx/2],#}]&,IvspcC[sim]/.Around[x_,y_]:>x];

n=1;;-1;;3;
t[sim][[n]]//Round[#,0.1]&//Echo[#,"Time: "]&;
ListLogPlot[Riffle[Ivspc[sim][[n]],IvspcC[sim][[n]]],
PlotStyle->{{Black,AbsoluteThickness[1]},{Red,Dashed,AbsoluteThickness[1]}},
PlotRange->{1*^-15,Automatic},
FrameLabel->{"R/L","(Subscript[Log, 10]) (Subscript[I, H])"},
PlotLegends->Placed[LineLegend@{"Resistive","Coulomb"},Scaled@{0.35,0.2}],
Epilog->{
Arrow[{Scaled@{0.8,0.98},Scaled@{0.8,0.4}}]
},
FrameTicks->{{Table[{10^i,ToString@i},{i,-20,-10}],Automatic},{Automatic,Automatic}}
]//Echo//Export[outputDir<>"/figures/gauge.eps",#]&

Compare four methods

runs

sim="/Users/hzhou/Work/pencilCode/runs/hongzhe_pc_dardel/saffman/k200del6a";

Nx=getParam[sim,"Nx"];
Lx=(2π)/(2π);
dx=Lx/Nx//N;
(*read all 4 methods*)
Clear[m]
{t[sim],m[3]}=read1D[sim,"Iv_spc_saffman_mag.dat"];
{t[sim],m[4]}=read1D[sim,"Iv_sps_saffman_mag.dat"];
{t[sim],m[1]}=read1D[sim,"Iv_bcc_saffman_mag.dat"];
{t[sim],m[2]}=read1D[sim,"Iv_bcs_saffman_mag.dat"];

Do[
m[i]=Map[Transpose[{dx*Range[Nx/2],#}]&,m[i]];
pos=Nearest[m[i][[1,;;,1]]->"Index",r0][[1]];
m[i]=Transpose[{t[sim],m[i][[;;,pos,2]]}];
ph[i]=pDecay@@Transpose[m[i]],
{i,4}];

{t[sim],m[5]}=read1D[sim,"power_saffman_mag.dat"];
m[5]=Transpose[{t[sim],Map[a[1]/.FindFit[#[[2;;5]],a[1]*k^2/(2π^2)+a[2]*k^4,{a[1],a[2]},k]&,m[5]]}];
ph[5]=pDecay@@Transpose[m[5]];

ticks1=Table[{10^i,Superscript["10",ToString@i]},{i,-20,20}];
ticks2=Join[Table[{10^i,ToString@i},{i,-20,-10}],Table[{10^i,PaddedForm[i,{2,1}]},{i,-5,1,0.5}]];
style={ImagePadding->{{80,10},{50,10}},FrameTicks->{{ticks2,Automatic},{ticks1,Automatic}},
PlotStyle->({Black,#}&/@{AbsoluteThickness[1.5],Dashed,Dotted,DotDashed,AbsoluteThickness[0.5]}),
ImageSize->400
};
Print["Plotting for run: ",label[sim]];
{
ListLogLogPlot[m/@Range[5],FrameLabel->{"t","(Subscript[Log, 10]) (Subscript[I, H])"},
Epilog->{Inset[Style["(a)",pcLabelStyle],Scaled@{0.93,0.9}]},
style],
ListLogLogPlot[Transpose[{Most@t[sim],ph[#]}]&/@Range[5],
FrameLabel->{"t","(Subscript[Log, 10]) (Subscript[p, H])"},
PlotRange->{1*^-3,All},
PlotLegends->Placed[
LineLegend[{"BC, cubic","BC, spherical","CI, cubic","CI, spherical","Fitting"},LegendLayout->{"Row",3}],
{Right,Bottom}],
Epilog->{Inset[Style["(b)",pcLabelStyle],Scaled@{0.07,0.9}]},
style]}//List//Grid//Echo//Export[outputDir<>"/figures/4methods.eps",#]&