Instantly share code, notes, and snippets.

Created Aug 26, 2013
 def emr_step(dt_n, y_n, dy_n, dt_nm1, y_nm1): """Take a single step of the explicit midpoint rule. From G&S pg. 715 and Prinja's thesis pg.45. """ dtr = dt_n / dt_nm1 y_np1 = (1 - dtr**2)*y_n + (1 + dtr)*dt_n*dy_n + (dtr**2)*(y_nm1) return y_np1 def bdf2_dydt(ts, ys): """Get dy/dt at time ts[-1] (allowing for varying dt). Gresho & Sani, pg. 715""" dt_n = ts[-1] - ts[-2] dt_nm1 = ts[-2] - ts[-3] y_np1 = ys[-1] y_n = ys[-2] y_nm1 = ys[-3] # Copied from oomph-lib (algebraic rearrangement of G&S forumla). weights = [1.0/dt_n + 1.0/(dt_n + dt_nm1), - (dt_n + dt_nm1)/(dt_n * dt_nm1), dt_n / ((dt_n + dt_nm1) * dt_nm1)] dydt = sp.dot(weights, [y_np1, y_n, y_nm1]) return dydt def bdf2_mp_gs_lte_estimate(ts, ys): """Estimate LTE using combination of bdf2 and explicit midpoint. From oomph-lib and G&S. """ # Get local values (makes maths more readable) dt_n = ts[-1] - ts[-2] dt_nm1 = ts[-2] - ts[-3] dtr = dt_n / dt_nm1 dtrinv = 1.0 / dtr y_np1 = ys[-1] y_n = ys[-2] y_nm1 = ys[-3] # Invert bdf2 to get predictor data (using the exact same function as # was used in the residual calculation). dy_n = bdf2_dydt(ts, ys) # Calculate predictor value (variable dt explicit mid point rule) y_np1_EMR = emr_step(dt_n, y_n, dy_n, dt_nm1, y_nm1) error_weight = ((1.0 + dtrinv)**2) / \ (1.0 + 3.0*dtrinv + 4.0 * dtrinv**2 + 2.0 * dtrinv**3) # Calculate truncation error -- oomph-lib error = (y_np1 - y_np1_EMR) * error_weight return error
Owner Author

davidshepherd7 commented Aug 26, 2013

 emr_step corresponds to the calculations done by BDF<2>::calculate_predicted_values with the weights set by BDF<2>::set_predictor_weights(). bdf2_dydt corresponds to TimeStepper::time_derivative(1,...) with the weights set by BDF<2>::set_weights(). bdf2_mp_gs_lte_estimate corresponds to BDF<2>::temporal_error_in_value with the error weight set by BDF<2>::set_error_weights().
Owner Author

davidshepherd7 commented Aug 26, 2013

 In python -1 gets the last element of the list, -2 the one before last, etc. I store time and y values in time order, i.e. ts[0] = t_0, ts[-2] = t_n and ts[-1] = t_{n+1}. ** is the power operator