Figure 9.35:

Predictions of LC measurement for reduced model.

Code for Figure 9.35

Text of the GNU GPL.

flow.dat

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
   9.000000      9.7817838E-02   4.786000    
   19.00000       1.290601       14.21800    
   29.00000       1.262779       13.99800    
   39.00000       1.240016       13.81800    
   49.00000       1.215736       13.62600    
   59.00000       1.187156       13.40000    
   69.00000       1.160599       13.19000    
   79.00000      0.5545961       8.398001    
   89.00000       1.109508       12.78600    
   99.00000       1.361925       14.78200    
   109.0000       1.332333       14.54800    
   119.0000       1.294648       14.25000    
   129.0000       1.271885       14.07000    
   139.0000       1.244569       13.85400    
   149.0000       1.220794       13.66600    
   159.0000       1.197778       13.48400    
   169.0000       1.170463       13.26800    
   179.0000       1.150988       13.11400    
   189.0000       1.128983       12.94000    
   199.0000       1.105968       12.75800    
   209.0000       1.092563       12.65200    
   219.0000       1.360155       14.76800    
   229.0000       1.344473       14.64400    
   239.0000       1.329045       14.52200    
   249.0000       1.305524       14.33600    
   259.0000       1.275426       14.09800    
   269.0000       1.259239       13.97000    
   279.0000       1.244063       13.85000    
   289.0000       1.219783       13.65800    
   299.0000       1.191961       13.43800    
   309.0000       1.173751       13.29400    
   319.0000       1.158322       13.17200    
   329.0000      0.3163429       6.514000    
   339.0000      0.5495375       8.358001    
   349.0000       1.224841       13.69800    
   359.0000      0.3228684       6.565600    
   369.0000      9.4845297E-04   4.020000    
   379.0000      4.6794413E-04   4.016201    
   389.0000      5.4383278E-04   4.016800    
   399.0000      4.9326418E-04   4.016400    
   409.0000      3.6680698E-04   4.015400    
   419.0000      0.2701340       6.148601    
   429.0000       1.185891       13.39000    
   439.0000      0.9299334       11.36600    
   449.0000      0.5259147       8.171202    
   459.0000      8.9793204E-04   4.019600    
   469.0000      8.7256433E-04   4.019400    
   479.0000      7.9672335E-04   4.018800    
   489.0000      4.9319270E-04   4.016400    
   499.0000      4.6799183E-04   4.016200    
   509.0000      3.6678315E-04   4.015400    
   519.0000      0.2330809       5.855600    
   529.0000      0.7784327       10.16800    
   539.0000      7.5244429E-03   4.072000    
   549.0000      7.5244186E-03   4.072000    
   559.0000      5.7539940E-03   4.058000    
   569.0000      8.0302954E-03   4.076001    
   579.0000      0.5019882       7.982001    
   589.0000     -3.0983209E-03   3.988000    
   599.0000      9.4850064E-04   4.020000    
   609.0000     -6.3252446E-05   4.012000    
   619.0000     -1.3278484E-03   4.002000    
   629.0000     -4.6158792E-03   3.976000    
   639.0000     -6.3204767E-05   4.012000    
   649.0000     -2.0865917E-03   3.996000    
   659.0000      0.1752122       5.398000    
   669.0000      0.3024322       6.404000    
   679.0000      4.9952744E-03   4.052001    
   689.0000      3.9835931E-03   4.044001    
   699.0000      0.5114223       8.056601    
   709.0000      3.4151078E-04   4.015200    
   719.0000      3.1619071E-04   4.015000    
   729.0000      1.1386871E-04   4.013400    
   739.0000      6.3300133E-05   4.013000    
   749.0000      3.4151078E-04   4.015200    
   759.0000      0.1059619       4.850400    
   769.0000      6.7019463E-04   4.017800    
   779.0000      6.8099439E-02   4.551000    
   789.0000      8.7261200E-04   4.019400    
   799.0000      9.7372534E-04   4.020200    
   809.0000      8.7256433E-04   4.019400    
   819.0000      7.9669955E-04   4.018800    
   829.0000      7.2083471E-04   4.018200    
   839.0000      5.6912901E-04   4.017000    
   849.0000      4.6801567E-04   4.016200    
   859.0000      5.1856041E-04   4.016600    
   869.0000      4.4267176E-04   4.016000    

lc.dat

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
    414        0.17050
    424        0.16040
    434        0.13250
    444        0.10840
    493        0.10140
    503        0.10420
    513        0.10450
    523        0.09700
    533        0.08240
    582        0.06780
    592        0.06970
    602        0.07020
    612        0.06920
    621        0.06840
    631        0.06750
    641        0.06730
    651        0.06730
    661        0.06540
    671        0.05770
    681        0.05710
    690        0.05630
    700        0.04660
    710        0.04550
    720        0.04490
    730        0.04560
    740        0.04500
    750        0.04340
    759        0.04280
    769        0.04260
    779        0.04090
    789        0.03990
    799        0.03970
    848        0.03940
    858        0.03740
    868        0.03710

find_time_index.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# Converted from find_time_index.m - helper function for lc_reduced
import numpy as np

def find_time_index(t, now):
    """Return (index, is_switching_time) for the last t[index] <= now.
    t must be sorted with no duplicates; now must be >= t[0].
    """
    t = np.asarray(t)
    if now < t[0]:
        raise ValueError(f'find_time_index: now={now} out of range [{t[0]}, {t[-1]}]')
    if now > t[-1]:
        index = len(t) - 1
    else:
        index = int(np.max(np.where(t <= now)[0]))
    is_switching_time = bool(t[index] == now)
    return index, is_switching_time

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# Converted from lc_reduced.m - reactor LC prediction with reduced model
# [depends] flow.dat lc.dat find_time_index.py
# This py-file loads data file 'flow.dat'.
# This py-file loads data file 'lc.dat'.
import numpy as np
from scipy.integrate import solve_ivp
from misc import gnuplotsave
from find_time_index import find_time_index

vrdiclo   = 2.3497
k1k2ratio = 2.0
vro       = 2370.0
teaf      = 0.00721
teaden    = 0.728

flow = np.loadtxt('flow.dat')
lc   = np.loadtxt('lc.dat')
tflow    = np.concatenate([[0.0], flow[:, 0]])
flowrate = np.concatenate([[0.0], flow[:, 1] / teaden])
tlc    = lc[:, 0]
lcmeas = lc[:, 1]

ntimes = 80
tlin   = np.linspace(0, tflow[-1], ntimes)
timeout = np.unique(np.concatenate([tlin, tflow, tlc]))
timeout = np.sort(timeout)

idx = np.array([int(np.where(timeout == t)[0][0]) for t in tlc])
assert np.allclose(timeout[idx], tlc), 'index extraction failed'

def res_reduced(time, x, vrdiclo_, k1k2ratio_):
    vr, eps2 = x
    if vr <= 0.0:
        raise ValueError(f'vr is zero at time={time}')
    badded = (vr - vro) * teaf
    eps1   = badded - eps2
    index, _ = find_time_index(tflow, time)
    fteaf    = flowrate[index] * teaf
    dvr_dt   = flowrate[index]
    if badded == 0.0:
        deps2_dt = 0.0
    else:
        deps2_dt = fteaf / (1.0 + k1k2ratio_ * (vrdiclo_ - (badded - eps2)) /
                             (badded - 2.0*eps2))
    return [dvr_dt, deps2_dt]

x0   = [vro, 0.0]
ks   = [0.5, 1.0, 2.0, 8.0, 32.0]
nks  = len(ks)

lcpred = np.zeros((len(timeout), nks))

for i, k12 in enumerate(ks):
    sol = solve_ivp(res_reduced, [timeout[0], timeout[-1]], x0,
                    args=(vrdiclo, k12), t_eval=timeout,
                    method='Radau', rtol=1.5e-8, atol=1.5e-8,
                    dense_output=False)
    vr   = sol.y[0]
    eps2 = sol.y[1]
    badded = (vr - vro) * teaf
    eps1   = badded - eps2
    cc     = (eps1 - eps2) / vr
    ccd    = eps2 / vr
    with np.errstate(invalid='ignore', divide='ignore'):
        lcp = np.where((cc + 2*ccd) > 0, cc / (cc + 2*ccd), 1.0)
    # First index where flowrate is positive (0-indexed); before this,
    # no TEA has been added so cc=ccd=0 and the lc ratio is 0/0.
    firstposflow = int(np.min(np.where(flowrate > 0)[0]))
    lcp[:firstposflow+1] = 1.0
    lcpred[:, i] = lcp

table = np.column_stack([timeout, lcpred])
data  = {'lc': lc, 'table': table}
gnuplotsave('lc_reduced.dat', data)