How to add (energy) constraint when using NDSolve to Equation of Motion

From the docs,

solConstraint2[x0_, p0_, m_, ω_, time_] := 
 NDSolve[{x'[t] == p[t]/m, p'[t] == -m ω^2 x[t], x[0] == x0, p[0] == p0},
  {x, p}, {t, 0, time}, 
  Method -> {"TimeIntegration" -> {"Projection", "Invariants" -> energy[x[t], p[t], m, ω]}}]

This keeps the energy to within about 0.12 of its starting value. Not great, but it does not drift.

Note: Using the option WorkingPrecision -> 50 on the original sol also reigns in the drift and keeps it around 6 * 10^-10.

I do not know offhand what is a good way to recast as a DAE. One way to enforce the algebraic constraint, without getting an overdetermined (albeit consistent) system, is to use the projection method. I cribbed some of the submethod settings from advanced documentation in tutorial/NDSolveProjection.

sol2[x0_, p0_, m_, ω_, time_] := 
 NDSolve[{x'[t] == p[t]/m, p'[t] == -m ω^2 x[t], x[0] == x0, 
   p[0] == p0}, {x[t], p[t]}, {t, 0, time}, MaxSteps -> ∞, 
  Method -> {"Projection", 
    Method -> {"SymplecticPartitionedRungeKutta", 
      "DifferenceOrder" -> 4, "PositionVariables" -> {x[t]}}, 
    "Invariants" -> energy[x[t], p[t], m, ω]}, PrecisionGoal -> 12]

This seems to remove that drift off the constraint manifold and also to keep fluctiations to around 10^(-9). Could perhaps do better with more tuning of the various options and suboptions.