slip output inconsistency

[ydluo]

Dear QDYN developers,

During the first SCEC-SEAS benchmark practice we discovered the QDYN slip output is updated with only first order accuracy by outputting v*dt instead of the high-order RK/BS method as used in computation.
It is only a output issue and does not affect the actual computation, yet it better be fixed.

Yingdi

[martijnende]

The RK45 (4th order) method that is currently implemented in QDYN can perform a single integration step from t0 to t1, and it will automatically determine how many intermediate steps are needed to remain within a certain tolerance. However, for this we need a continuous function that gives us velocity as a function of time (which we don't have). This is also a problem for the BS solver. Using the RK45 method, we can make a 4th order reconstruction of the slip at time t1 if the intermediate values of V are stored. I'm not sure how this works for the BS solver.

[jpampuero]

A possible fix is to add a new variable "slip" and the following equation to the ODE system:
d(slip)/dt = (slip velocity)
This would solve for "slip" with the same order as the other variables, whatever ODE solver we choose.

[martijnende]

A possible fix is to add a new variable "slip" and the following equation to the ODE system:
d(slip)/dt = (slip velocity)
This would solve for "slip" with the same order as the other variables, whatever ODE solver we choose.

That would work, though I'm not sure if this approach exerts any undesired controls on the step size. In the RK45 solver, the step size is determined by the maximum over all vector elements of the difference between the 4th and 5th order approximations, so we may be paying a computational price for resolving the slip within the specified tolerance.

Before diving into the code, I would suggest waiting for my pull requests for v. 2.0 to get merged into the master branch. I have generalised the procedure for packing and unpacking the input/output vector for the ODE solver, which makes it more convenient (and potentially less buggy) to add new variables to the ODE system. Once my computer has arrived I will do some final testing of the code changes and perform the merge.

[ydluo]

Isn't it true that with the full history, full accuracy V and DT, we can calculate high order D with proper RK/BS weights with existing nodes of V? It will be one order less acuurate than V, explicit rather than implicit, but still should be way more accurate than current 1st order so might be able to mitigate the D output artifact we found during the benchmark. There will also be (almost) no increase in computation.

[martijnende]

Yes, in the case of the RK solver, if you store the intermediate outputs of V, you can reconstruct the slip using the same RK integration scheme. I don't know how this works for the BS solver.

[jpampuero]

Coefficients depend on the method. The strategy I proposed has the advantage that it applies to all methods (RK, BS, etc) and does not require storing additional "intermediate outputs". I think it's more developer-friendly (any "blackbox" ODE solver can be used). I don't expect it to affect the time step, but that's hard to tell without trying it.

[martijnende]

All right. Since it is easy enough to simply append another vector to the ODE input vector, we can just try it and see how much it affects performance. Since the stiffest part of the ODE is in the computation of v (exponential dependence on stress), v will most likely control the time step. In that case there's just a bit of additional overhead of computing the error (N extra floating point divisions).

I will execute this fix once I have tested the changes to v2.0 and merged everything into the master branch.

[ydluo]

I am thinking of adding the functionality of outputting off-fault observation points at surface for 3D models so we can get slip / slip rate at surface that can e compared with observation.
It seems quite straightforward (just calculate unilateral Vx, Vy and Vz at Obs. point, run-time from Okada) but should have same issue of accuracy for slip. Any thoughts? We don't want to put the Obs. points in the ODE system for sure.

[jpampuero]

No need to add ground displacements/velocities in the ODE solver. Their computation is a post-processing step.
Fixing the accuracy of slip will automatically make the ground displacements accurate.
Ground velocities can be computed from slip rates, hence they should be accurate already in the current version.

[martijnende]

I've moved the slip integration into the main ODE routine (#38), which solves this issue