This phase is meant to be a more accurate representation of what you might actually do if the project problem came up in the course of your research. Rather than manually writing a rudimentary leapfrog algorithm, you’ll leverage existing infrastructure, taking advantage of the huge amount of time and research that others have put into fast, flexible differential equation solvers. You’ll modify WaveSim.jl to replace the leapfrog algorithm with a SecondOrderODEProblem (or whichever “problem” you’d like) from SciML’s DifferentialEquations package. Since this means a higher performance ceiling than a fixed-time-step leapfrog algorithm, it will be required to run on 2d-medium-in.wo in 30 seconds. I recommend using your phase 4 code as a starting point since you’ve already done some of the optimization work you’ll need.

Preparation

You’ll need to add DifferentialEquations to your package before using it. To do so, activate WaveSim.jl and run add DifferentialEquations in the package manager.

Using an ODE Problem

The self-contained example code uses the one-dimensional mountain range example problem; it’s similar to what you’ll do, but there are a couple of tricks to using a second order ODE problem.

The Derivative Function

Defining the function that you’ll pass to SecondOrderODEProblem itself isn’t too bad, and dhdt! in the example offers some guidance. Notice that dhdt!’s main computation is just the differential equation that defines the example problem–you’ll do something similar with the damped wave equation differential equation. Here’s the gist:

function dvdt!(a, v, u, c, t) # u: displacement; v: velocity; a: acceleration
    @. a = laplacian(u) - c*v # just pseudo-code, needs significant overhaul
    return nothing
end

ContinuousCallback

The function you provide to ContinuousCallback is a bit tricky since what’s passed to the function is an object that holds both the displacement and the velocity in its member x. To compensate for this, you can overload energy to take u and v, then do a callback along the lines of:

# When this returns 0, iteration stops
terminationcondition(vu, args...) = energy(vu.x[2], vu.x[1])-energythreshold
cb = ContinuousCallback(terminationcondition, terminate!)

vu in the example is an ArrayPartition which contains two arrays in its x member; the first element (vu.x[1]) is velocity, the second (vu.x[2]) displacement.

Since energy decays consistently for the wave plates you’ll be dealing with, you can pass interp_points=0 to ContinuousCallback to prevent energy from being evaluated so much–this often results in dramatic speedups.

SecondOrderODEProblem

A SecondOrderODEProblem is only slightly different from ODEProblem–in addition to specifying an initial value (displacement in your case), you’ll pass in its initial derivative (velocity):

ODEProblem{true}(           derivative!,     u0, timespan, c, callback=cb)
SecondOrderODEProblem{true}(derivative!, v0, u0, timespan, c, callback=cb)

You can pass save_on=false and save_start=false to SecondOrderODEProblem to save time and memory at the expense of diagnostic information.

Updating the Wave Orthotope

Once you get a solution using SecondOrderODEProblem, you’ll need to update your wave orthotope with the results. That will look something like:

solution = solve(problem)
w.t[] += last(solution.t)
w.u   .= last(solution.u).x[2] # get displacement from ArrayPartition
w.v   .= last(solution.u).x[1] # get velocity from ArrayPartition

Results

Different solving algorithms will yield slightly different results, but they shouldn’t differ from previous results by more than a few percent. As an example, using the default solver should result in a final simulation time of just over 104 for 2d-medium-in.wo.

Submission

You’ll submit a tarball named WaveSim.jl.tar.gz that I’ll grade exactly as I did in phase 4.

Grading

This phase is worth 20 points. The following deductions, up to 20 points total, will apply for a program that doesn’t work as laid out by the spec:

Defect Deduction
Failure to run correctly on 2d-medium-in.wo in 30 seconds on m9 5 points
Failure to run correctly on 2d-medium-in.wo in 60 seconds on m9 5 points
Failure of your WaveSim.jl to work correctly on each of 3 2-dimensional test files 4 points each
WaveSim.jl doesn’t use an external differential equation solver from DifferentialEquations 1-20 points