Fixing the difference order for space derivatives in a PDE when using NDSolve

Question

I am trying to implement the fix to make the difference order uniform for all space derivatives in a PDE, as described in this post: NDSolve uses different difference order for different spatial derivative when solving PDE

It does not work, and I am out of my wits as to why. (Actually, it is the system of 2 coupled PDEs that I am solving.) After the system is set up I included the standard call to NDSolve to show that the problem was set up correctly, and then in the end there is a minimal implementation of the fix. As I wrote it, it does not work, although the the output indicates that time stepping was performed. Can anyone please help?

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]

RHSCeq = D[H[t, x], x] D[Subscript[C, Bb][t, x], x, x]; 
RHSHeq = Subscript[C, Bb][t, x] D[H[t, x], x, x]; 
Nkmax = 6.482250811894023; 
Nλmax = 2 Pi/Nkmax;

domainlength = 20 Nλmax; center1 = domainlength/8; w1 = 0.05; 
IniShapeH = 1 - 0.02 Exp[-(x - center1)^2/w1^2];
IniShapeC = 1/2;

Δt = 0.1

solution = 
NDSolve[{D[Subscript[C, Bb][t, x], t] == RHSCeq, 
D[H[t, x], t] == RHSHeq, 
Subscript[C, Bb][t, 0] == Subscript[C, Bb][t, domainlength], 
H[t, 0] == H[t, domainlength], Subscript[C, Bb][0, x] == IniShapeC,
H[0, x] == IniShapeH, 
WhenEvent[Mod[t, Δt] == 0, Print[t]]}, {Subscript[C, 
Bb], H}, {x, 0, domainlength}, {t, 0, 1}, 
AccuracyGoal -> MachinePrecision/2, 
PrecisionGoal -> MachinePrecision/2, InterpolationOrder -> All, 
Method -> {"MethodOfLines", 
"SpatialDiscretization" -> {"TensorProductGrid", 
  "DifferenceOrder" -> 2, MinPoints -> 1000, MaxPoints -> 1000, 
  AccuracyGoal -> MachinePrecision/2, 
  PrecisionGoal -> MachinePrecision/2}, 
Method -> {"ImplicitRungeKutta", 
  "Coefficients" -> "ImplicitRungeKuttaRadauIIACoefficients", 
  "DifferenceOrder" -> 2}}]

Plot[{IniShapeH, First[H[1, x] /. solution]}, {x, 0, domainlength}, 
PlotRange -> All, 
PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}]

Plot[{IniShapeC, First[Subscript[C, Bb][1, x] /. solution]}, {x, 0, 
domainlength}, PlotRange -> All, 
PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}]

endtime = 1; nspacepoints = 1000; xdifforder = 4; 

mol[n_Integer, o_] := {"MethodOfLines", 
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
"MinPoints" -> n, "DifferenceOrder" -> o}}

nd := NDSolveValue[{D[Subscript[C, Bb][t, x], t] == RHSCeq, 
D[H[t, x], t] == RHSHeq, 
Subscript[C, Bb][t, 0] == Subscript[C, Bb][t, domainlength], 
H[t, 0] == H[t, domainlength], Subscript[C, Bb][0, x] == IniShapeC,
H[0, x] == IniShapeH, 
WhenEvent[Mod[t, Δt] == 0, Print[t]]}, {Subscript[C, 
Bb], H}, {x, 0, domainlength}, {t, 0, endtime}, 
Method -> mol[nspacepoints, xdifforder]]

sold = fix[endtime, xdifforder]@nd

Answer 1

It's actually a matter of evaluation order. Notice that in my answer, I've defined nd with a With[{nd:= ……}, ……] structure. Why? Because the last "step" inside fix is

Unevaluated[ndsolve][[2]] /. NDSolve`ProcessSolutions@newstate

The Unevaluated[ndsolve][[2]] part is for extracting the dependent variable from ndsolve, which requires the dependent variables to explictly exist. If you simply define nd:=…… without With, nd will still be nd inside Unevaluated. In other words, what happened here can be reproduced with the following example:

lst := {a, b}
Unevaluated[lst][[2]]

Part::partd

In this example, though the OwnValues of lst is {a, b}, when it's not evaluated, it's no more than a Symbol in Part's eyes, so the code doesn't work as expected.

Then why it'll work if I define nd with a With? Because the… Er… scoping mechanism of With is quite similar to rule substitution. When e.g.

With[{lst := {a, b}}, Unevaluated[lst][[2]]]

is executed, lst in Unevaluated[lst][[2]] is simply replaced by {a, b}.

So, to fix your code, we just need to replace nd with the code inside nd e.g. modify the last line to:

sold = Unevaluated@fix[endtime, xdifforder][nd] /. OwnValues@nd