changed param values in IRForcing.jl

examples/twodturb/IRForcing.jl

@@ -2,16 +2,16 @@ using PyPlot, FourierFlows
 import FourierFlows.TwoDTurb
 import FourierFlows.TwoDTurb: energy, enstrophy, dissipation, injection, drag
 
- n, L  = 128, 2π
- ν, nν = 1e-6, 2
+ n, L  = 256, 2π
+ ν, nν = 1e-7, 2
  μ, nμ = 1e-1, 0
-dt, tf = 0.01, 2/μ
+dt, tf = 0.01, 0.2/μ
 nt = round(Int, tf/dt)
 ns = 4
 
 # Forcing
 kf, dkf = 12.0, 2.0
-σ = 0.005
+σ = 0.1
 
 gr  = TwoDGrid(n, L)
 
@@ -51,13 +51,16 @@ function makeplot(prob, diags)
 
   E, D, I, R = diags
 
+
+  t = round(μ*prob.state.t, 2)
   sca(axs[1]); cla()
   pcolormesh(prob.grid.X, prob.grid.Y, prob.vars.q)
   xlabel(L"$x$")
   ylabel(L"$y$")
+  title("\$\\nabla^2\\psi(x,y,\\mu t= $t )\$")
   axis("square")
 
-  sca(axs[2]); cla()
+  sca(axs[3]); cla()
 
   i₀ = 1
   dEdt = (E[(i₀+1):E.count] - E[i₀:E.count-1])/prob.ts.dt
@@ -76,26 +79,26 @@ function makeplot(prob, diags)
 
   # If the Ito interpretation was used for the work
   # then we need to add the drift term
-  plot(E.time[ii], I[ii2] + σ, label="injection (\$I\$)")   # Ito
-  # plot(E.time[ii], I[ii2] , label="injection (\$I\$)")      # Stratonovich
-  plot(E.time[ii], σ+0*E.time[ii], "--")
-  plot(E.time[ii], -D[ii], label="dissipation (\$D\$)")
-  plot(E.time[ii], -R[ii], label="drag (\$R\$)")
+  plot(μ*E.time[ii], I[ii2] + σ, label=L"work ($P$)")   # Ito
+  # plot(μ*E.time[ii], I[ii2] , label=L"work ($P$)")      # Stratonovich
+  plot(μ*E.time[ii], σ+0*E.time[ii], "--", label=L"ensemble mean  work ($\langle P\rangle $)")
+  # plot(μ*E.time[ii], -D[ii], label="dissipation (\$D\$)")
+  plot(μ*E.time[ii], -R[ii], label=L"drag ($D=2\mu E$)")
+  plot(μ*E.time[ii], 0*E.time[ii], "k:", linewidth=0.5)
   ylabel("Energy sources and sinks")
-  xlabel(L"$t$")
+  xlabel(L"$\mu t$")
   legend(fontsize=10)
 
-  sca(axs[3]); cla()
-  plot(E.time[ii], total[ii], label=L"computed $dE/dt$")
-  plot(E.time[ii], dEdt, "--k", label=L"numerical $dE/dt$")
+  sca(axs[2]); cla()
+  plot(μ*E.time[ii], total[ii], label=L"computed $P-D$")
+  plot(μ*E.time[ii], dEdt, "--k", label=L"numerical $dE/dt$")
   ylabel(L"$dE/dt$")
-  xlabel(L"$t$")
+  xlabel(L"$\mu t$")
   legend(fontsize=10)
 
   sca(axs[4]); cla()
-  plot(E.time[ii], residual, "c-", label=L"residual = computed $-$ numerical")
-  xlabel(L"$t$")
-  ylabel(L"$dE/dt$")
+  plot(μ*E.time[ii], residual, "c-", label=L"residual $dE/dt$ = computed $-$ numerical")
+  xlabel(L"$\mu t$")
   legend(fontsize=10)
 
   residual
@@ -128,6 +131,8 @@ for i = 1:ns
   # savefig(savename, dpi=240)
 end
 
+# savename = @sprintf("./plots/stochastictest_kf%d_%06d.png", kf, prob.step)
+# savefig(savename, dpi=240)
 
 # savediagnostic(E, "energy", out.filename)
 # savediagnostic(D, "dissipation", out.filename)