Instability, Courant Condition and Robustness about solving 2D+1 PDE

Question

After several discussions, I would like to focus on the robustness of solving 2D+1 PDE by considering all suggested methods from @xzczd (see here) I found that the Ratio between the convection term and diffusion is crucial. Here is the code.

Clear["Global`*"]
Clear[tosameorder, fix]
tosameorder[state_NDSolve`StateData, order_] := 
 state /. a_NDSolve`FiniteDifferenceDerivativeFunction :> 
   NDSolve`FiniteDifferenceDerivative[a@"DerivativeOrder", 
    a@"Coordinates", "DifferenceOrder" -> order, 
    PeriodicInterpolation -> a@"PeriodicInterpolation"]

fix[endtime_, order_] := 
 Function[{ndsolve}, 
  Module[{state = 
     First[NDSolve`ProcessEquations @@ Unevaluated@ndsolve], 
    newstate}, newstate = tosameorder[state, order]; 
   NDSolve`Iterate[newstate, endtime];
   Unevaluated[ndsolve][[2]] /. NDSolve`ProcessSolutions@newstate], 
  HoldAll]
a = 1;
T = 1;
ωb = -15; ωt = 15;
A = 6.5;
γ = .1;
kT = 0.1;
φ = 0;

With[{u = u[t,θ, ω]}, 
eq = D[u, t] == -D[ω u,θ] - D[-A Sin[3θ] u, ω] + γ (1 + Sin[3θ])  kT  D[
   u, {ω, 2}] + γ  (1 + Sin[3θ]) D[ω u, ω];
ic = u == E^(-((ω^2 +(θ+Pi/4)^2)/(2 a^2))) 1/(2 π a) /. t -> 0];
ufun = NDSolveValue[{eq, ic, u[t, -π, ω] == u[t, π, ω], 
 u[t,θ, ωb] == 0, u[t,θ, ωt] == 0}, u, {t, 0, T}, {θ, -π, π}, {ω, ωb, ωt}, 
      Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MaxPoints" -> {81, 51}, 
         "MinPoints" -> {41, 31}}}]; // AbsoluteTiming
plots = Table[
Plot3D[Abs[ufun[t,θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, AxesLabel -> Automatic, 
 PlotPoints -> 30, BoxRatios -> {Pi, ωb, 1}, 
 ColorFunction -> "LakeColors", PlotRange -> All], {t, 0, T, 
 T/10}]; // AbsoluteTiming
ListAnimate[plots]

One can see that the coefficient of diffusion term (2nd-order term) is really small ($0.1*0.1*sin(3\theta)$) especially when $sin(3\theta)=-1$, how could it possible satisfy the Courant condition with such vanishing diffusion term?

The below is the result:

What I expect is the something similar to following (obtain by adding artificial diffusion)

The main question is the robustness about solving this partially convection-dominated problem in a efficient and stable way. Many thanks.

Note about artificial diffusion:

Max[γ[θ], 0.3] D[ u, {ω, 2}]

This is how I put the artificial diffusion. But for the problem, the angle-dependent diffusion and convection are different.

Full code is here, it's a little bit messy:

Clear["Global`*"]
(*//////////////////////////////////////         parameters         \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

n = 3; 
ϕ = π/2; vg0 = 5;
vg[t_] := vg[t] = vg0;
(*vg[t_]:=vg[t]=2vg0*1/(\[ExponentialE]^(k(t-Tp1))+1)-vg0;*)
(*vg[t_]:=vg[t]=2vg0*1/(\[ExponentialE]^(k(t-Tp1))+1)-2vg0*1/(\
\[ExponentialE]^(k(t-Tp2))+1)+vg0;*)


τi = 3; Γ = 10; 

Vb = 0; μ[α_] := (-1)^(α + 1) Vb/2  ;


XTicks1 = Table[2 π*j, {j, -10, 10}];
XTicks2 = Table[π/6*j, {j, -10, 10}];
YTicks = Table[2 π*j, {j, -10, 10}];



(*//////////////////////////////////////         Karrasch poles and \
coefficients          \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

Np = 25; M = 2 Np;
B = Normal[
   SparseArray[{Band[{2, 1}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}], 
     Band[{1, 2}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}]}, M]];
{bvals, bvecs} = Eigensystem[B];
Zp = Table[Abs[N[1/bvals[[2 p]]]], {p, 1, Np}];
Rp = Table[
   N[(Normalize[bvecs[[2 p]]][[1]]/(2 bvals[[2 p]]))^2], {p, 1, Np}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

σ0[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := σ0[θ, V, τ0, Γ, 
    Vg, ϕ] = 
   1/2 - I/(
     4 π) (PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]);
τ[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := τ[θ, V, τ0, Γ, 
    Vg, ϕ] = -τ0 n (Cos[n θ + ϕ] - 
      Cos[n  θ]) σ0[θ, 
     V, τ0, Γ, Vg, ϕ];
σ1[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := σ1[θ, V, τ0, Γ, 
    Vg, ϕ] = 
   1/(8 π^2 Γ)*(PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]) - \!\(
\*UnderoverscriptBox[\(∑\), \(p = 
       1\), \(Np\)]\(Rp[\([\)\(p\)\(]\)] \((
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ 
                 Sin[n\ θ])\) - V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ 
                  Sin[n\ θ])\) - V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]] + 
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ 
                 Sin[n\ θ])\) + V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ 
                  Sin[n\ θ])\) + V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]])\)\)\);
Uprime[θ_, τ0_, ϕ_] := τ0 n Cos[n θ];

mol[m_Integer, n_Integer] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> m, 
    "MinPoints" -> n}}
τeff = 
  FunctionInterpolation[-τi n Cos[n θ] + τ[θ, 
     Vb, τi, Γ, 
     vg0, ϕ], {θ, -π, π}];
γ = 
  FunctionInterpolation[(τi n (Cos[n θ + ϕ] - 
        Cos[n  θ]))^2 σ1[θ, 
     Vb, τi, Γ, 
     vg0, ϕ], {θ, -π, π}, AccuracyGoal -> 5, 
   InterpolationPrecision -> MachinePrecision];
Plot[{τeff[θ], γ[θ]}, {θ, -π, \
π}, PlotRange -> All]
a = 0.7;
T = 20;
ωb = -8;
ωt = 8;
θ0 = -π/4;
L = 1;
With[{u = u[t, θ, ω]}, 
  eq = D[u, t] == -ω D[ u, θ] - 
     1/L*τeff[θ] D[u, ω] + 
     Re[γ[θ]] D[ ω u, ω] + 
     Max[γ[θ], 0.3] D[ u, {ω, 2}];
  ic = u == 
     E^(-(ω^2 + (θ + π/
                4)^2)/(2 a^2))/(2  π a^2) /. t -> 0];

ufun = NDSolveValue[{eq, ic, 
     u[t, -π, ω] == u[t, π, ω], 
     u[t, θ, ωb] == 0, u[t, θ, ωt] == 0}, 
    u, {t, 0, 
     T}, {θ, -π, π}, {ω, ωb, ωt}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> {41, 61}, "MinPoints" -> {41, 61}, 
        "DifferenceOrder" -> 4}}]; // AbsoluteTiming
(*Plot3D[ufun[T,θ,ω],{θ,-π,π},{ω,\
ωb,ωt},PlotRange\[Rule]All,AxesLabel\[Rule]Automatic,\
PlotPoints\[Rule]30,ColorFunction\[Rule]"LakeColors"]*)
plots = Table[
    Plot3D[Abs[
      ufun[t, θ, ω]], {θ, -π, π}, {\
ω, ωb, ωt}, PlotRange -> All, 
     AxesLabel -> Automatic, PlotPoints -> 30, 
     ColorFunction -> "LakeColors"], {t, 0, T, 
     T/50}]; // AbsoluteTiming
ListAnimate[plots] // AbsoluteTiming

Stable Julia code for the same problem: Matrix of 2D+1 PDE

function F_eff(x, Gamma, Delta, Q, EpOm, A)
    return -A*sin(x + Delta)
end

# effective friction
function gamma_eff(x, Gamma, Delta, Q, EpOm, A)
    return Gamma
end

# effective diffusion
function D_eff(x, Gamma, Delta, Q, EpOm, A, kBT)
    return 2.0*kBT*Gamma
end

function make_FPE_matrix(xi,vj,Gamma, Delta, Q, EpOm, A, kBT)

    Nx = length(xi)
    Nv = length(vj)
    hx = xi[2]-xi[1]
    hv = vj[2]-vj[1]

    mat = zeros(Float64, Nx*Nv+1, Nx*Nv)

    for i=0:(Nx-1)
        for j=0:(Nv-1)

            mat[i*(Nv)+j+1,i*(Nv)+j+1] = -D_eff( xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)

            # -d/dx (v P)
            if i == 0
                mat[i*(Nv)+j+1, (Nx-1)*(Nv)+j+1] = vj[j+1]/(2*hx) # PBC
            end
            if i > 0
                mat[i*(Nv)+j+1, (i-1)*(Nv)+j+1] = vj[j+1]/(2*hx)
            end
            if i < Nx-1
                mat[i*(Nv)+j+1, (i+1)*(Nv)+j+1] = -vj[j+1]/(2*hx)
            end
            if i == Nx-1
                mat[i*(Nv)+j+1, (0)*(Nv)+j+1] = -vj[j+1]/(2*hx)   # PBC
            end


            if j > 0
                mat[i*(Nv)+j+1, (i)*(Nv)+(j-1)+1] = F_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)/(2*hv) -
                    gamma_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)*vj[j+1-1]/(2*hv)
                if j> 1
                    mat[i*(Nv)+j+1, (i)*(Nv)+(j-2) + 1] = 0.5*D_eff(xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)
                end
            end
            if j < Nv-1
                mat[i*(Nv)+j+1, (i)*(Nv)+(j+1)+1] = -F_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)/(2*hv) +
                    gamma_eff(xi[i+1], Gamma, Delta, Q, EpOm, A)*vj[j+1+1]/(2*hv)
                if j < Nv-2
                    mat[i*(Nv)+j+1, (i)*(Nv)+(j+2)+1] = 0.5*D_eff(xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)
                end
            end

        end
    end

    for i = 1:Nx*Nv
         mat[end,i] = hx*hv
    end

    return mat
end

Update (9/4) I try to run the compiled code (with VS2017 community). Here is my code:

Clear["Global`*"]
(*//////////////////////////////////////         parameters         \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

n = 3; 
ϕ = π/2; vg0 = 5;
vg[t_] := vg[t] = vg0;


τi = 3; Γ = 10; 

Vb = 0; μ[α_] := (-1)^(α + 1) Vb/2  ;


XTicks1 = Table[2 π*j, {j, -10, 10}];
XTicks2 = Table[π/6*j, {j, -10, 10}];
YTicks = Table[2 π*j, {j, -10, 10}];



(*//////////////////////////////////////         Karrasch poles and \
coefficients          \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

Np = 25; M = 2 Np;
B = Normal[
   SparseArray[{Band[{2, 1}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}], 
     Band[{1, 2}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}]}, M]];
{bvals, bvecs} = Eigensystem[B];
Zp = Table[Abs[N[1/bvals[[2 p]]]], {p, 1, Np}];
Rp = Table[
   N[(Normalize[bvecs[[2 p]]][[1]]/(2 bvals[[2 p]]))^2], {p, 1, Np}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

σ0[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := σ0[θ, V, τ0, Γ, 
    Vg, ϕ] = 
   1/2 - I/(
     4 π) (PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]);
τ[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := τ[θ, V, τ0, Γ, 
    Vg, ϕ] = -τ0 n (Cos[n θ + ϕ] - 
      Cos[n  θ]) σ0[θ, 
     V, τ0, Γ, Vg, ϕ];
σ1[θ_, V_, τ0_, Γ_, 
   Vg_, ϕ_] := σ1[θ, V, τ0, Γ, 
    Vg, ϕ] = 
   1/(8 π^2 Γ)*(PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + ϕ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]) - \!\(
\*UnderoverscriptBox[\(∑\), \(p = 
       1\), \(Np\)]\(Rp[\([\)\(p\)\(]\)] \((
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ 
                 Sin[n\ θ])\) - V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ 
                  Sin[n\ θ])\) - V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]] + 
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ 
                 Sin[n\ θ])\) + V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ 
                  Sin[n\ θ])\) + V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]])\)\)\);
Uprime[θ_, τ0_, ϕ_] := τ0 n Cos[n θ];

mol[m_Integer, n_Integer] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> m, 
    "MinPoints" -> n}}
τeff = 
  FunctionInterpolation[-τi n Cos[n θ] + τ[θ, 
     Vb, τi, Γ, 
     vg0, ϕ], {θ, -π, π}];
γ = 
  FunctionInterpolation[(τi n (Cos[n θ + ϕ] - 
        Cos[n  θ]))^2 σ1[θ, 
     Vb, τi, Γ, 
     vg0, ϕ], {θ, -π, π}, AccuracyGoal -> 5, 
   InterpolationPrecision -> MachinePrecision];
Plot[{τeff[θ], γ[θ]}, {θ, -π, \
π}, PlotRange -> All]
a = 0.7;
T = 20;
ωb = -8;
ωt = 8;
θ0 = -π/4;
L = 1;
Clear[τeff, γ]
points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
τefflst = 
  Chop@N@Array[
     Function[θ, -τi n Cos[n θ] + τ[θ, Vb, τi, Γ, vg0, ϕ]], points@θ, {-π, π}];
γlst = 
  Chop@Array[
     Function[θ, (τi n (Cos[n θ + ϕ] - Cos[n θ]))^2 σ1[θ, Vb, τi, Γ, vg0, ϕ]], 
     points@θ, {-π, π}];

With[{u = u[θ, ω]}, 
  rhs2 = -ω ct@D[u, θ] - 1/L τeff[θ] ct@D[u, ω] + 
         γ[θ] ct@D[ω u, ω] + γ[θ] fw@D[bw@D[u, ω], ω]; 
  iclst2 = Table[E^(-((ω^2 + (θ + π/4)^2)/(2 a^2)))/(2 π a^2), 
                 {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];
rttheta = RescalingTransform[{{-Pi, Pi}}, {{1, points@θ}}];

With[{rc = RuleCondition, cg = Compile`GetElement}, 
  rhsfunc2 = Hold@
            Compile[{{u, _Real, 2}, {τeff, _Real, 1}, {γ, _Real, 1}}, 
             Table[rhs2, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}], 
             RuntimeOptions -> EvaluateSymbolically -> False, CompilationTarget -> C] /. 
           OwnValues@rhs2 /. 
          u[theta_, omega_] :> 
           rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
               Round@rt@{theta, omega}) /. (coef : τeff | γ)[theta_] :> 
          rc@(cg[coef, First@Round@rttheta@{theta}]) /. DownValues@delta /. 
       DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] // 
     ReleaseHold // Last];
T = 30;
ulstfunc2 = 
   NDSolveValue[{u'[t] == rhsfunc2[u[t], τefflst, γlst], u[0] == iclst2}, 
    u, {t, 0, T}, MaxSteps -> Infinity]; // AbsoluteTiming
(* {36.177260, Null} *)

lst = Table[
   ListPlot3D[ulstfunc2[t]\[Transpose], PlotRange -> All, 
    DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, 5, 1/20}];    
ListAnimate@lst

There is an error message:

CompiledFunction::rterr: -- Message text not found -- (compiledFunction5) (8)

Answer 1

OK, since neither of the default difference scheme nor the fix function works properly on this PDE, let's discretize the spatial derivative all by ourselves in a way different from that in NDSolve or in fix[…, …]@NDSolve.

The spatial discretization will be done in the following way: \begin{aligned} \frac{\partial^2 f}{\partial x^2}\Biggl|_{x=x_i} \approx & \frac{f(x_i+h)-2f(x_i)+f(x_i-h)}{h^2} \\ \frac{\partial f}{\partial x}\Biggl|_{x=x_i} \approx & \frac{f(x_i+h)-f(x_i-h)}{2h} \\ \frac{\partial (x f)}{\partial x}\Biggl|_{x=x_i} \approx & \frac{(x_i+h) f(x_i+h)-(x_i-h)f(x_i-h)}{2h} \end{aligned}

Notice the difference formula for $\frac{\partial (x f)}{\partial x}\Bigl|_{x=x_i}$ is critical here, because it produces a result different from

$$\frac{\partial (x f)}{\partial x}\Biggl|_{x=x_i}=\left(x\frac{\partial f}{\partial x}+f\frac{\partial x}{\partial x}\right)\Biggl|_{x=x_i}\approx x_i \frac{f(x_i+h)-f(x_i-h)}{2h}+f(x_i)$$

Periodic b.c. will be set in both directions because it's relatively easier to implement. Since the Dirichlet b.c. in $\omega$ direction is just an approximation for b.c. at infinity, this should not have any significant influence on the result.

The remaining work is just coding:

ωb = -5; ωt = 5;
a = 1; A = 6.5; γ = .1; kT = 0.1; φ = 0;

ClearAll[fw, bw, ct]
SetAttributes[#, HoldAll] & /@ {fw, bw, ct};
fw@D[expr_, x_] := Subtract @@ (expr /. {{x -> x + delta@x}, {x -> x}})/delta@x
bw@D[expr_, x_] := Subtract @@ (expr /. {{x -> x}, {x -> x - delta@x}})/delta@x
ct@D[expr_, x_] := 
 Subtract @@ (expr /. {{x -> x + delta@x}, {x -> x - delta@x}})/(2 delta@x)

Clear[delta]
delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0

points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
With[{u = u[θ, ω]}, 
  rhs =-ct@D[ω u, θ] - ct@D[-A Sin[3 θ] u, ω] + 
       γ (1 + Sin[3 θ]) kT bw@D[fw@D[u, ω], ω] + γ (1 + Sin[3 θ]) ct@D[ω u, ω];
  iclst = Table[
                E^(-((ω^2 + (θ + Pi/4)^2)/(2. a^2))) 1/(2 π a), 
                {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]
  ];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];

With[{rc = RuleCondition, cg = Compile`GetElement}, 
  rhsfunc = Hold@Compile[{{u, _Real, 2}},            
             Table[rhs, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}], 
             RuntimeOptions -> EvaluateSymbolically -> False, 
             CompilationTarget -> C] /. OwnValues@rhs /. 
          u[theta_, omega_] :> 
           rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
               Round /@ rt@{theta, omega}) /. DownValues@delta /. 
       DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] // 
     ReleaseHold // Last];

T = 30;
ulstfunc = NDSolveValue[{u'[t] == rhsfunc[u[t]], u[0] == iclst}, u, {t, 0, T}, 
    MaxSteps -> Infinity]; // AbsoluteTiming
(* {33.1812583, Null} *)

lst = Table[
   ListPlot3D[ulstfunc[t]\[Transpose], PlotRange -> All, 
    DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, T, 1}];   
ListAnimate@lst

The result seems to be what you're expecting.

If you don't have a C compiler installed, take the CompilationTarget option and // Last away. I do recommend you to install one though.

Some advanced technique has been used in this code piece to make the discretization less tedious. To understand it better, you may want to read the following posts:

When should I, and when should I not, set the HoldAll attribute on a function I define?

Is there a way to simplify this replacement rule

How to make the code inside Compile conciser without hurting performance?

Why is CompilationTarget -> C slower than directly writing with C?

Replacement inside held expression

---

Update: Uncompilable coefficient function treatment

The coefficients of equation solved above is compilable. When they're not, such as in your 2nd code sample, some further modification is needed. The key idea is calculating coefficient values at grid points first and pass the value list to Compile:

Clear[τeff, γ]
points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
τefflst = 
  Chop@N@Array[
     Function[θ, -τi n Cos[n θ] + τ[θ, Vb, τi, Γ, vg0, ϕ]], points@θ, {-π, π}];
γlst = 
  Chop@Array[
     Function[θ, (τi n (Cos[n θ + ϕ] - Cos[n θ]))^2 σ1[θ, Vb, τi, Γ, vg0, ϕ]], 
     points@θ, {-π, π}];

With[{u = u[θ, ω]}, 
  rhs2 = -ω ct@D[u, θ] - 1/L τeff[θ] ct@D[u, ω] + 
         γ[θ] ct@D[ω u, ω] + γ[θ] fw@D[bw@D[u, ω], ω]; 
  iclst2 = Table[E^(-((ω^2 + (θ + π/4)^2)/(2 a^2)))/(2 π a^2), 
                 {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];
rttheta = RescalingTransform[{{-Pi, Pi}}, {{1, points@θ}}];

With[{rc = RuleCondition, cg = Compile`GetElement}, 
  rhsfunc2 = Hold@
            Compile[{{u, _Real, 2}, {τeff, _Real, 1}, {γ, _Real, 1}}, 
             Table[rhs2, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}], 
             RuntimeOptions -> EvaluateSymbolically -> False, CompilationTarget -> C] /. 
           OwnValues@rhs2 /. 
          u[theta_, omega_] :> 
           rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
               Round@rt@{theta, omega}) /. (coef : τeff | γ)[theta_] :> 
          rc@(cg[coef, First@Round@rttheta@{theta}]) /. DownValues@delta /. 
       DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] // 
     ReleaseHold // Last];
T = 30;
ulstfunc2 = 
   NDSolveValue[{u'[t] == rhsfunc2[u[t], τefflst, γlst], u[0] == iclst2}, 
    u, {t, 0, T}, MaxSteps -> Infinity]; // AbsoluteTiming
(* {36.177260, Null} *)

lst = Table[
   ListPlot3D[ulstfunc2[t]\[Transpose], PlotRange -> All, 
    DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, 5, 1/20}];    
ListAnimate@lst

Notice definitions of σ1 etc. are not included here, please find them in the body of the question.