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