Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
@davidshepherd7

This comment has been minimized.

Copy link
Owner Author

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().

@davidshepherd7

This comment has been minimized.

Copy link
Owner Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.