Runge-Kutta Method

[conorhef...@gmail.com]

Hello,

I have recently implemented the non-Stiffly Accurate 4th/5th Order Runge-Kutta Methods from the recent Carpenter & Kennedy paper:

https://doi.org/10.1016/j.apnum.2018.10.007

I have found an error with implementing the final stage with the b coefficients (equation (4) in the paper). The methods appear to work if this final step is omitted but causes problems when implemented. I suspect there is some problem with how I have done the final step but can't see where it is going wrong. My edits to the current RK methods in Dedalus to implement the higher order methods is shown below  in green. Changing the line for i in range(1, self.stages+2): to for i in range(1, self.stages+1): omits the final step mentioned above.

I have also attached plots of the Kinetic Energy from a sample problem (Doubly Periodic NSE in a box) which shows that the RK443 method appears to agree with the RK4 method but causes an error when the final step is included (it jumps then rapidly decays).

If anyone can see where I may have made a mistake with my edit of the RK method source code I would greatly appreciate it.

Cheers,

Conor

Code:

def __init__(self, pencil_length, domain):

        self.RHS = CoeffSystem(pencil_length, domain)

        # Create coefficient systems for multistep history
        self.MX0 = CoeffSystem(pencil_length, domain)
        self.LX = LX = [CoeffSystem(pencil_length, domain)
                        for i in range(self.stages+2)]
        self.F = F = [CoeffSystem(pencil_length, domain)
                      for i in range(self.stages+2)]

        self._LHS_params = None

    def step(self, solver, dt):
        """Advance solver by one timestep."""

        # Solver references
        pencils = solver.pencils
        evaluator = solver.evaluator
        state = solver.state

        evaluator_kw = {}
        evaluator_kw['world_time'] = world_time = solver.get_world_time()
        evaluator_kw['wall_time'] = world_time - solver.start_time
        evaluator_kw['sim_time'] = sim_time_0 = solver.sim_time
        evaluator_kw['timestep'] = dt
        evaluator_kw['iteration'] = solver.iteration

        # Other references
        RHS = self.RHS
        MX0 = self.MX0
        LX = self.LX
        F = self.F
        A = self.A
        H = self.H
        c = self.c
        b = self.b
        k = dt

        # Check on updating LHS
        update_LHS = (k != self._LHS_params)
        self._LHS_params = k
        final_stage = False

        # Compute M.X(n,0)
        MX0.data.fill(0)
        for p in pencils:
            fast_csr_matvec(p.M, state.get_pencil(p), MX0.get_pencil(p))
            if update_LHS:
                p.LHS_solvers = [None] * (self.stages+1)
        # Compute stages
        # (M + k Hii L).X(n,i) = M.X(n,0) + k Aij F(n,j) - k Hij L.X(n,j)
        # The extra loop sets final_stage = True and computes eq 4 from Carp/Kenn
        for i in range(1, self.stages+2):
            # Compute F(n,i-1), L.X(n,i-1)
            if i == self.stages+1 :
                final_stage = True

            state.scatter()
            evaluator_kw['sim_time'] = solver.sim_time
            if i == 1:
                evaluator.evaluate_scheduled(**evaluator_kw)
            else:
                evaluator.evaluate_group('F', **evaluator_kw)
            LX[i-1].data.fill(0)
            F[i-1].data.fill(0)
            for p in pencils:
                fast_csr_matvec(p.L, state.get_pencil(p),
                                LX[i-1].get_pencil(p))
                fast_csr_matvec(p.pre_left, solver.F.get_pencil(
                    p), F[i-1].get_pencil(p))

            # Construct RHS(n,i)
            np.copyto(RHS.data, MX0.data)
            for j in range(0, i):
                if final_stage == False:
                    RHS.data += (k * A[i, j]) * F[j].data
                    RHS.data -= (k * H[i, j]) * LX[j].data
                else :
                    # Compute M.X(n+1,0) = M.X(n,0) + k bj F(n,j) - k bj L.X(n,j)
                    RHS.data += (k * b[j]) * F[j].data
                    RHS.data -= (k * b[j]) * LX[j].data

            state.data.fill(0)
            for p in pencils:
                pRHS = RHS.get_pencil(p)
                # Construct LHS(n,i)
                if update_LHS and final_stage == False:
                    np.copyto(p.LHS.data, p.M_exp.data +
                              (k*H[i, i])*p.L_exp.data)
                    # Remove old solver reference before building new solver
                    p.LHS_solvers[i] = None
                    p.LHS_solvers[i] = solver.matsolver(p.LHS, solver)
                if final_stage == False :
                    pX = p.LHS_solvers[i].solve(pRHS)
                else :
                    # All stages have been computed
                    pX = pRHS

                if p.pre_right is None :
                    state.set_pencil(p, pX)
                else:
                    fast_csr_matvec(p.pre_right, pX, state.get_pencil(p))
            if final_stage == False :
                solver.sim_time = sim_time_0 + k*c[i]