Following the example set by xzczd's answer: an alternative to an external library would be to implement a simple shock capturing finite volume scheme like the one put forward in New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations, 2000, A. Kurganov and E. Tadmor -- KT FV scheme -- see also this wikipedia article. I have implemented this scheme in Mathematica a few years ago. The scheme is purpose build for hyperbolic conservation laws and includes extensions for the treatment of parabolic conservation laws (e.g. PDEs including diffusion terms). The scheme can be applied directly to non-linear PDEs (initial value problems) of the type
$$
\begin{align}
\partial_t u(t,x)+\partial_x F(t,x,u)&=\partial_x Q(t,x,u,\partial_x u)\\
u(0,x)&=u_0(x)
\end{align}
$$
with a broad range of boundary conditions (ones which can be implemented easily with ghost-cells). Below you find an implementation of the scheme using CompiledFunction for performance.
(** Periodic boundary conditions for ghost cells **)
BCperiodic=Function[{vin},Module[{v},
v=vin;
v[[1]]=v[[-5]];(* Subscript[u, -2 ]= Subscript[u, n-3] *)
v[[2]]=v[[-4]];(* Subscript[u, -1 ]= Subscript[u, n-2] *)
v[[-3]]=v[[3]];(* Subscript[u, n-1 ]= Subscript[u, 0] *)
v[[-2]]=v[[4]];(* Subscript[u, n] = Subscript[u, 1] *)
v[[-1]]=v[[5]];(* Subscript[u, n+1] = Subscript[u, 2] *)
v
]];
(** Selected MUSCL scheme Flux limiter **)
ϕMinMod=Function[{Δuj,Δujp1},If[Abs@Δujp1==0.,1.,Max[0.,Min[1.,Δuj/Δujp1]]]]; (* MinMod Limiter [https://en.wikipedia.org/wiki/Flux_limiter] *)
(** Main step method **)
ClearAll[KTstepper];
KTstepper[nx_][Δx_,x_,x12_][F_,dFdu_,Q_,BC_,pars___]:=
Compile[{{t,_Real},{v,_Real,1}},
Module[{
u,(* Conserved quanties, dim=nx+4 *)
Δu,(* Differences of conserved quanties, dim=nx+3 *)
ux,(* Reconstructed slopes, dim =nx+2 *)
up12T,(* Right intermediate values at the cell interface, dim =nx+1 *)
um12T,(* Left intermediate values at the cell interfaces, dim=nx+1 *)
Fp12T,(* Advection flux through the cell interfaces, dim=nx+1 *)
Fm12T,(* Advection flux through the cell interfaces, dim=nx+1 *)
λp12T, (* Jacobian eigenvalues (single value) the cell interfaces, dim=nx+1 *)
λm12T, (* Jacobian eigenvalues (single value)the cell interfaces, dim=nx+1 *)
a12T,(* Approxmiate velocities at the cell interfaces, dim=nx+1 *)
H12 ,(* Numerical advection fluxes at the cell interfaces, dim=nx+1 *)
P12(* Diffusion fluxe through the cell interfaces, dim =nx+1 *)
},
(* **** Boundary condition **** *)
u=BC[Join[{0.,0.},Take[v,{1,nx}],{0.,0.}]];
(* **** PDE Convection Flux **** *)
Δu=Differences[u];
ux=MapThread[0.5*#2(ϕMinMod[#1,#2])&,{Take[Δu,{1,-2}],Take[Δu,{2,-1}]}];(* [KTO2-0, eq. (2.4)*0.5*Δx]: modified to be compatible with generic flux limiters *)
up12T=Take[u,{3,-2}]-Take[ux,{2,-1}];(* [KTO2-0, eq. (4.5)]: modified to be compatible with generic flux limiters *)
um12T=Take[u,{2,-3}]+Take[ux,{1,-2}];(* [KTO2-0, eq. (4.5)]: modified to be compatible with generic flux limiters *)
λp12T=MapThread[dFdu[t,#1,#2,pars]&,{x12,up12T}];
λm12T=MapThread[dFdu[t,#1,#2,pars]&,{x12,um12T}];
a12T=MapThread[Max[Max@Abs@#1,Max@Abs@#2]&,{λp12T,λm12T}]; (* [KTO2-0, eq. (3.2)] [KTO2-0, footnote 2] *)
Fp12T=MapThread[F[t,#1,#2,pars]&,{x12,up12T}];
Fm12T=MapThread[F[t,#1,#2,pars]&,{x12,um12T}];
H12=0.5*(Fp12T+Fm12T-a12T*(up12T-um12T)); (* [KTO2-0, eq. (4.4)] *)
(* **** PDE Diffusion Flux **** *)
P12=0.5*MapThread[Q[t,#1,#3,#5,pars]+Q[t,#2,#4,#5,pars]&,{
Join[{x12[[1]]-Δx*0.5},x],
Join[x,{x12[[-1]]+Δx*0.5}],
Take[u,{2,-3}],
Take[u,{3,-2}],
Take[Δu,{2,-2}]/Δx
}]; (* [KTO2-0, eq. (4.14)] *)
(* **** Result **** *)
Return[Flatten[(-Differences[H12]+Differences[P12])/Δx,1]] (* [KTO2-0,eq.(4.13)] *)
], RuntimeOptions->"Speed", CompilationOptions->{"InlineExternalDefinitions"->True,"ExpressionOptimization"->True}
]
(** Solver/time-stepper using NDSolve to solve the ODE (MOL) system of KTstepper *)
ClearAll[KTsolver]
KTsolver/:Options[KTsolver]={AccuracyGoal->Automatic,PrecisionGoal->Automatic,WorkingPrecision->MachinePrecision,Method->Automatic,MaxSteps->Automatic};
KTsolver[OptionsPattern[]][x0_,x1_,n_Integer][F_,dFdu_,Q_,BC_,params___][u0_][t0_,t1_]:=Module[{
no,
Δx,xi,xi12,
v0,
PDEcompileFkt,PDEfkt,PDEsystem,PDEv,PDEt,PDEsolver,PDEsolution,
nt,tw,tm},
(* x0 and x1 are the first and last cell centers *)
Δx=(x1-x0)/(n-1);
xi=Table[x0+Δx*i,{i,0,n-1}];
xi12=Table[x0+Δx*(i-0.5),{i,0,n}];
v0=u0[#]&/@xi;(* approimate cell averages with mid point value *)
PDEcompileFkt=KTstepper[n][Δx,xi,xi12][F,dFdu,Q,BC,params];
PDEfkt[t_?NumericQ,v_]:=PDEcompileFkt[t,v];
PDEsystem={Equal[PDEv[t0],v0],Equal[Derivative[1][PDEv][PDEt],PDEfkt[PDEt,PDEv[PDEt]]]};
nt=0;
tw=-AbsoluteTime[];
PDEsolver=Inactive[NDSolveValue][PDEsystem,PDEv,{PDEt,t0,t1},
StepMonitor:>{nt++,tm=PDEt},
Method->OptionValue[Method],
AccuracyGoal->OptionValue[AccuracyGoal],
PrecisionGoal->OptionValue[PrecisionGoal],
WorkingPrecision->OptionValue[WorkingPrecision],
MaxSteps->OptionValue[MaxSteps]
];
PDEsolution=Activate@PDEsolver;
tw+=AbsoluteTime[];
tm=PDEsolution["Domain"]//Last//Last;
KTsolution[<|"solution"->PDEsolution,"n"->n,"no"->no,"Δx"->Δx,"xi"->xi,"nt"->nt,"tw"->tw,"t0"->t0,"t1"->t1,"tm"->tm|>]
]
(** Output wrapper **)
ClearAll[KTsolution];
KTsolution[asoc_][key_]:=asoc[key] /; MemberQ[Keys[asoc], key]
KTsolution[asoc_][key_Symbol]:=asoc[ToString@key] /; MemberQ[Keys[asoc], ToString@key]
KTsolution[asoc_][t_]:=With[{sol=asoc["solution"],xi=asoc["xi"]},{xi,sol[t]}\[Transpose]/;!MissingQ[sol]]
MakeBoxes[KTsolution[asoc_Association],StandardForm]/;BoxForm`UseIcons:=Module[{xi,n,no,Δx,nt},
(* [https://mathematica.stackexchange.com/a/99914] *)
(* [https://mathematica.stackexchange.com/a/79891] *)
xi=With[{a=asoc["xi"]},If[MissingQ[a],{},{a[[1]],a[[2]],Skeleton[Length[a]-3],Last@a}]];
n=With[{a=asoc["n"]},If[MissingQ[a],0,a]];
nt=With[{a=asoc["nt"]},If[MissingQ[a],0,a]];
Δx=With[{a=asoc["Δx"]},If[MissingQ[a],Missing,a]];
BoxForm`ArrangeSummaryBox[KTsolution,
KTsolution[asoc],
None,
{{"n = "<>ToString[n],"Δx = "<>ToString[Δx]},{Row[{"Subscript[x, i] = ",xi}],SpanFromLeft},{"nt="<>ToString[nt],SpanFromLeft}},
{},
StandardForm,
"Interpretable"->True
]
]
The only interesting part of the code above is the one for KTstepper as it implements the semi-discrete method-of-lines (MOL) finite volume discretion of Kurganov and Tadmore. KTsolver generates and solves the MOL ODE system using NDSolve and KTsolution is just an output wrapper.
In the above notation the Burger's equation is given by $F(t,x,u)=F(u)=u^2$ and $Q(t,x,u,\partial_x u)=Q(\partial_x u)=\nu \partial_x u$ in case of non-vanishing diffusion $\nu>0$. Apart from initial and boundary (in this example periodic ones) conditions the scheme requires the Eigenvalues of the Jacobian $(\partial F/\partial u)$, which in case of the Burger's equation is simply $\partial F/\partial u = u$. The following code computes three solutions at variable $\nu$:
burgersEq=Sequence[{t,x,u}|->0.5u^2,{t,x,u}|->u,{t,x,u,dudx,nu}|->nu*dudx,BCperiodic];
νi={0,0.01,0.05};
νiSols=Table[KTsolver[AccuracyGoal->6,PrecisionGoal->6][0.,2π,512][burgersEq,ν][x|->1/3+2 Sin[x]/3][0,2π],{ν,νi}]
Manipulate[ListLinePlot[#[t]&/@νiSols,GridLines->Automatic,Frame->True,FrameLabel->{"x","u["<>ToString[NumberForm[t,3]]<>",x]"},PlotLegends->Placed[Row[{"ν=",#}]&/@νi,{Left,Top}],PlotRange->{{0,2π},{-0.5,1.25}},ImagePadding->{{45,10},{35,10}}],{t,0,2π}]
The method works well even on small spatial grids:
ni={64,128,256,512};
niSols=Table[KTsolver[AccuracyGoal->6,PrecisionGoal->6][0.,2π,n][burgersEq,0][x|->1/3+2 Sin[x]/3][0,2π],{n,ni}]
Manipulate[ListLinePlot[#[t]&/@niSols,GridLines->Automatic,Frame->True,FrameLabel->{"x","u["<>ToString[NumberForm[t,3]]<>",x]"},PlotLegends->Placed[Row[{"Subscript[n, x]=",#}]&/@ni,{Left,Top}],PlotRange->{{0,2π},{-0.5,1.25}},ImagePadding->{{45,10},{35,10}}],{t,0,2π}]
Note that the solutions are completely free of spurious oscillations around the shock wave. The reason for this is the fact that the KT FV scheme is shock capturing and second order accurate in $\Delta x$. It employs a MUSCL reconstruction. In my experience with the numerical solution of this kind of conservation laws the KT FV is a very good compromise between simplicity and performance. It works shockingly well and can be implemented in only a few lines of code (most of the code above is just setup, NDSolve options and the output wrapper) using mainly list manipulations. The performance in Mathematica (especially when using the CompiledFunction) is also not bad. If you want to see the method flex its "MUSCLs" have a look into the examples studied in the paper or on wikipedia.
Obviously the exact solution cannot be like that, it seems the above code is no longer working for large $t$, but I don't know how to edit the code to make the numerical results reasonable. May I ask what should I add/change to the code to solve the results for $t=2\pi$? Thanks in advance!
and the result looks nice too :)
