7
$\begingroup$

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: enter image description here

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

enter image description here

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)
$\endgroup$
  • 1
    $\begingroup$ See section "Stabilization of Convection-Dominated Equations" in this guide. Note that this is also an ongoing research topic. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 9:17
  • $\begingroup$ @HenrikSchumacher Thanks for the suggestion. Actually, the second animation is obtained from adding artificial diffusion as suggested in the documentation. In question is "How far away is solution relative to the actual solution" and "Can we actually improve the stability by specifying a proper grids or boundary condition". Of course, the Courant condition is an important issue but hard to control since time-step is adaptive, which is not able to modify manually. $\endgroup$ – Bob Lin Aug 29 '18 at 10:11
  • $\begingroup$ "What I expect is the something similar to following (obtain by adding artificial diffusion)" How do you add the artificial diffusion? Can you include the specific setting? $\endgroup$ – xzczd Aug 29 '18 at 11:09
  • $\begingroup$ @xzczd I put the code in my question. The idea is to keep some minimal value of diffusion such that the solution is still stable. $\endgroup$ – Bob Lin Aug 29 '18 at 13:08
  • 1
    $\begingroup$ The Julia code is incomplete. Anyway, I managed to find a hint from it. $\endgroup$ – xzczd Aug 30 '18 at 11:26
8
$\begingroup$

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

enter image description here

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

enter image description here

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

$\endgroup$
  • $\begingroup$ I notice that you use the same difference scheme in the d(xf)/dx as in my friend's Julia code. Such a subtle detail can really affect the stability very much? Anyway, thanks a lot! $\endgroup$ – Bob Lin Aug 30 '18 at 14:14
  • $\begingroup$ @BobLin Yeah, surprising, but it's true. Notice this phenomenon is rather rare, before your question came out, this seems to be the only example in this site. $\endgroup$ – xzczd Aug 30 '18 at 14:28
  • $\begingroup$ @BobLin If you're going to be doing a lot of numerical solution of PDEs for your work, you should get more general background on the topic. One nice starting place is Numerical Recipes by Press et al. $\endgroup$ – Chris K Aug 30 '18 at 21:11
  • $\begingroup$ @xzczd Is it possible to apply your code to the PDE with f, which are interpolating function or user define function? In my second code, there are interpolating function like $τ_{eff}(θ)$ and $γ(θ)$. I want to apply your method to get the exact solution in the second example, but it seems to run quite too long. $\endgroup$ – Bob Lin Aug 31 '18 at 8:38
  • $\begingroup$ I notice that the functions are always outside the the derivative. Something like With[{u = u[t, θ, ω]}, eq = D[u, t] == -ω D[ u, θ] - 1/L*τeff[θ] D[u, ω] + γ[θ] D[ ω u, ω] + γ[θ]D[ u, {ω, 2}]; ic = u == E^(-(ω^2 + (θ + π/ 4)^2)/(2 a^2))/(2 π a^2) /. t -> 0]; But it still very slow, when I run NDSolve. $\endgroup$ – Bob Lin Aug 31 '18 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.