Python chemreac.ReactionDiffusion类代码示例

def test_n_jac_diags(n_jac_diags):
    N, n, nstencil = 10, 1, 7
    rd = ReactionDiffusion(n, [], [], [], N=N, nstencil=nstencil,
                           n_jac_diags=n_jac_diags, D=[9])
    assert np.allclose(rd.xcenters,
                       [.05, .15, .25, .35, .45, .55, .65, .75, .85, .95])
    y0 = np.ones(N)

    # Dense
    jref_cdns = np.zeros((n*N, n*N), order='F')
    jout_cdns = np.zeros((n*N, n*N), order='F')
    sm = SymRD.from_rd(rd)
    sm.dense_jac(0.0, y0.flatten(), jref_cdns)
    rd.dense_jac_cmaj(0.0, y0.flatten(), jout_cdns)
    assert np.allclose(jout_cdns, jref_cdns)

    # Banded
    jref_cbnd = rd.alloc_jout(order='F', pad=0)
    jout_cbnd = rd.alloc_jout(order='F')
    sm.banded_jac(0.0, y0.flatten(), jref_cbnd)
    rd.banded_jac_cmaj(0.0, y0.flatten(), jout_cbnd)
    assert np.allclose(jout_cbnd[rd.n*rd.n_jac_diags:, :], jref_cbnd)

    # Compressed
    jref_cmprs = rd.alloc_jout_compressed()
    jout_cmprs = rd.alloc_jout_compressed()
    sm.compressed_jac(0.0, y0.flatten(), jref_cmprs)
    rd.compressed_jac_cmaj(0.0, y0.flatten(), jout_cmprs)
    assert np.allclose(jout_cmprs, jref_cmprs)

def test_chained_parameter_variation():
    # A -> B
    names = ['A', 'B']
    rd = ReactionDiffusion(len(names), [], [], k=[],
                           substance_names=names, g_value_parents=[0], g_values=[[0, 1]],
                           param_names=['doserate'])
    durations = [1., 3., 2.]
    y0 = [13., 7.]
    ic = dict(zip(names, y0))
    doserates = [.3, .11, .7]
    npoints = 3
    odesys = rd._as_odesys(variables_from_params=dict(
        density=lambda self, params: 1.0
    ))
    res = odesys.chained_parameter_variation(
        durations, ic, {'doserate': doserates}, npoints=npoints,
        integrate_kwargs=dict(atol={k: 1e-8 for k in odesys.names}))
    assert res.xout.size == npoints*len(durations) + 1
    assert res.xout[0] == 0
    assert np.all(res.yout[0, :] == y0)
    expected = [.3]*npoints + [.11]*npoints + [.7]*(npoints+1)
    assert np.all(res.params[:, odesys.param_names.index('doserate')] == expected)
    cumulative = 0.0
    for dr, dur in zip(doserates, durations):
        mask = (cumulative <= res.xout) & (res.xout <= cumulative + dur)
        cumulative += dur
        t, y = res.xout[mask], res.yout[mask, :]
        a, b = y[:, 0], y[:, 1]
        refa = a[0]
        refb = b[0] + (t - t[0])*dr*a[0]
        assert np.allclose(refa, a)
        assert np.allclose(refb, b)
    res.extend_by_integration(np.sum(durations)+1, dict(doserate=doserates[-1]), integrator='cvode')
    assert abs(res.yout[-1, 1] - (refb[-1] + doserates[-1]*a[0])) < 1e-8

def test_integrate_nondimensionalisation__g_values(from_rsys):
    from chempy import Reaction, ReactionSystem
    from chempy.units import allclose, default_units as u
    rstr = "-> H + OH; Radiolytic({'radiolytic_yield': 2.1e-7*mol/J})"
    if from_rsys:
        rxn = Reaction.from_string(rstr, None)
        rsys = ReactionSystem([rxn], 'H OH')
        rd = ReactionDiffusion.from_ReactionSystem(
            rsys, unit_registry=SI_base_registry,
            variables=dict(doserate=0.15*u.Gy/u.s, density=0.998*u.kg/u.dm3))
        assert rd.g_value_parents == [-1]
        assert rd.g_values == [[2.1e-7]*2]
        assert abs(rd.fields[0][0] - 0.15*998) < 1e-14
    else:
        rd = ReactionDiffusion.nondimensionalisation(
            2, [[]], [[0, 1]], [2.1e-7*0.15*0.998*u.molar/u.second],
            unit_registry=SI_base_registry)
    C0 = [3*molar, 4*molar]
    tout = np.linspace(0, 1)*day
    integr = Integration(rd, C0, tout, integrator='scipy')
    k_m3_p_mol_p_sec = 0.15*998*2.1e-7
    t_sec = np.linspace(0, 24*3600)
    C0_mol_p_m3 = [3000, 4000]
    Cref_mol_p_m3 = np.empty(integr.Cout.squeeze().shape)
    Cref_mol_p_m3[:, 0] = C0_mol_p_m3[0] + k_m3_p_mol_p_sec*t_sec
    Cref_mol_p_m3[:, 1] = C0_mol_p_m3[1] + k_m3_p_mol_p_sec*t_sec
    print(integr.with_units('Cout').squeeze())
    print(integr.with_units('Cout').squeeze() - Cref_mol_p_m3*u.mole/u.metre**3)
    assert allclose(integr.with_units('tout'), t_sec*u.s)
    assert allclose(integr.with_units('Cout').squeeze(),
                    Cref_mol_p_m3*u.mole/u.metre**3)

def test_integrate_nondimensionalisation(from_rsys):
    from chempy import Reaction, ReactionSystem
    from chempy.units import allclose, default_units as u

    # 2A -> B
    if from_rsys:
        rxn = Reaction.from_string('2 A -> B; 2e-3*metre**3/mol/hour', None)
        rsys = ReactionSystem([rxn], 'A B')
        rd = ReactionDiffusion.from_ReactionSystem(rsys, unit_registry=SI_base_registry)
    else:
        rd = ReactionDiffusion.nondimensionalisation(
            2, [[0, 0]], [[1]], [2e-9/(umol/metre**3)/hour],
            unit_registry=SI_base_registry)
    C0 = [3*molar, 4*molar]
    tout = np.linspace(0, 1)*day
    integr = Integration(rd, C0, tout, integrator='scipy')

    k_m3_p_mol_p_sec = 2e-3/3600
    t_sec = np.linspace(0, 24*3600)
    C0_mol_p_m3 = [3000, 4000]
    Cref_mol_p_m3 = np.empty(integr.Cout.squeeze().shape)
    Cref_mol_p_m3[:, 0] = 1/(C0_mol_p_m3[0]**-1 + 2*k_m3_p_mol_p_sec*t_sec)
    missing_A = (C0_mol_p_m3[0] - Cref_mol_p_m3[:, 0])
    Cref_mol_p_m3[:, 1] = C0_mol_p_m3[1] + missing_A/2
    assert allclose(integr.with_units('tout'), t_sec*u.s)
    assert allclose(integr.with_units('Cout').squeeze(),
                    Cref_mol_p_m3*u.mol/u.metre**3, rtol=1e-6)

def test_ReactionDiffusion__f__wrong_fout_dimension():
    y0 = np.array([2.0, 3.0])
    k = 5.0
    # A -> B
    rd = ReactionDiffusion(2, [[0]], [[1]], [k])
    fout = np.ones((1,))*99  # fout too small
    with pytest.raises(AssertionError):
        rd.f(0.0, y0, fout)

def test_ReactionDiffusion__only_1_reaction__logy__logt(N):
    # See <test_ReactionDiffusion__only_1_reaction__logy_logt.png>
    t0 = 3.0
    y0 = np.array([2.0, 3.0]*N)
    k = 5.0
    # A -> B
    rd = ReactionDiffusion(2, [[0]], [[1]], [k], N, D=[0.0, 0.0],
                           logy=True, logt=True)
    fref = np.array([(-k*t0, t0*k*y0[i*2]/y0[i*2+1])
                     for i in range(N)]).flatten()
    _test_f_and_dense_jac_rmaj(rd, rd.logb(t0), np.log(y0), fref)

def test_ReactionDiffusion__only_1_species_diffusion_3bins(log):
    # Diffusion without reaction
    # 3 bins
    t0 = 3.0
    logy, logt = log
    D = 17.0
    y0 = np.array([23.0, 27.0, 37.0])
    N = 3
    x = [5.0, 7.0, 13.0, 15.0]
    xc = [4.0, 6.0, 10.0, 14.0, 16.0]
    nstencil = 3
    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x, N=N, logy=logy,
                           logt=logt, lrefl=False, rrefl=False,
                           nstencil=nstencil)
    assert np.allclose(rd.xc, xc)

    w = [1/16, -1/8, 1/16]  # finite diff. weights for 2nd order deriv
    for i in range(N):
        assert np.allclose(rd.D_weight[i*nstencil:(i+1)*nstencil], w)
    J = D*(w[0]*y0[0] + w[1]*y0[1] + w[2]*y0[2])
    fref = np.array([J, J, J])

    if logy:
        fref /= y0
    if logt:
        fref *= t0

    if logy:
        jref = D*np.array([  # jref[i, k] = ...
            [w[k]*y0[k]/y0[i] if k != i else -1/y0[k]*sum([
                w[j]*y0[j] if j != k else 0 for j in range(3)
            ]) for k in range(3)] for i in range(3)
        ])
        jref[0, 2] = 0.0  # dense_jac_rmaj only computes banded approx.
        jref[2, 0] = 0.0  # same as above.
    else:
        jref = D*np.array([[w[k] if abs(k-i) < 2 else 0.0 for
                          k in range(3)] for i in range(3)])

    if logt:
        jref *= t0

    y = rd.logb(y0) if logy else y0
    t = rd.logb(t0) if logt else t0
    _test_f_and_dense_jac_rmaj(rd, t, y, fref, jref)

    jout_bnd = np.zeros((4, 3), order='F')
    rd.banded_jac_cmaj(t, y, jout_bnd)
    jref_bnd = get_banded(jref, 1, 3)
    assert np.allclose(jout_bnd[1:, :], jref_bnd)

def test_integrators(log):
    logy, logt, use_log2 = log
    t0, tend, nt = 5.0, 17.0, 42
    tout = np.linspace(t0, tend, nt+1)

    # Update the dict if more integrators are added:
    solver_kwargs = {
        'scipy1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'tout': tout
        },
        'scipy2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'tout': (t0, tend),
            'dense_output': True
        },
        'cvode1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': tout
        },
        'cvode2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'method': 'adams',
            'tout': tout,
            'C0_is_log': True
        }
    }

    # A -> B
    n = 2
    k0 = 0.13
    rd = ReactionDiffusion(n, [[0]], [[1]], k=[k0], logy=logy, logt=logt, use_log2=use_log2)
    y0 = [3.0, 1.0]

    results = []
    for solver, kwargs in solver_kwargs.items():
        _y0 = rd.logb(y0) if kwargs.get('C0_is_log', False) else y0
        integr = Integration(rd, _y0, integrator=solver[:-1], **kwargs)
        if not kwargs.get('dense_output', False):
            results.append(integr.Cout)

    for result in results[1:]:
        assert np.allclose(results[0][0], result[0])

def test_from_ReactionSystem__g_values__multiple_types():
    from chempy import Reaction, ReactionSystem as RS
    from chempy.kinetics.rates import mk_Radiolytic
    RABG = mk_Radiolytic('alpha', 'beta', 'gamma')
    dens, yields_k, yields_v = .7, 'ya yb yg'.split(), [3, 5, 7]
    rxn = Reaction({}, {'H': 2}, RABG(yields_k))
    doserates = {'doserate_alpha': 11, 'doserate_beta': 13, 'doserate_gamma': 17}
    yields = dict(zip(yields_k, yields_v))
    params = dict(doserates)
    params.update(yields)
    params['density'] = dens
    ref = .7*2*(3*11 + 5*13 + 7*17)
    rat = rxn.rate(params)
    assert abs(rat['H'] - ref) < 1e-13
    assert RABG.parameter_keys == ('density', 'doserate_alpha', 'doserate_beta', 'doserate_gamma')
    assert RABG.argument_names == tuple('radiolytic_yield_%s' % k for k in 'alpha beta gamma'.split())

    rs = RS([rxn], checks=())
    rd = ReactionDiffusion.from_ReactionSystem(rs, variables=params)
    gv = rd.g_values
    assert len(gv) == 3
    assert np.allclose(sorted(gv), [[v*2] for v in sorted(yields_v)])
    assert len(rd.fields) == 3
    assert len(rd.fields[0]) == 1
    assert np.allclose(sorted(np.array(rd.fields).squeeze()), sorted([drat*dens for drat in doserates.values()]))
    fout = rd.alloc_fout()
    rd.f(0, np.array([0.0]), fout)
    assert np.allclose(fout, ref)

def test_ReactionDiffusion__only_1_field_dep_reaction_logy_logt(N):
    t0 = 3.0
    y0 = np.concatenate([np.array([2.0, 3.0])/(x+1) for x in range(N)])
    k = 5.0
    # A -> B

    rd = ReactionDiffusion(2, [], [], [], N, D=[0.0, 0.0],
                           fields=[[x+1 for x in range(N)]],
                           g_values=[[-k, k]], g_value_parents=[0],
                           logy=True, logt=True)

    def k_(bi):
        return k*(bi+1)

    fref = np.array([(-k_(i)*t0, k_(i)*t0*y0[i*2]/y0[i*2+1])
                     for i in range(N)]).flatten()
    _test_f_and_dense_jac_rmaj(rd, rd.logb(t0), np.log(y0), fref)

def test_ReactionDiffusion__only_1_reaction__logy(N):
    # See <test_ReactionDiffusion__only_1_reaction__logy.png>
    t0 = 3.0
    y0 = np.array([2.0, 3.0]*N)
    k = 5.0
    # A -> B
    rd = ReactionDiffusion(2, [[0]], [[1]], [k], N, D=[0.0, 0.0], logy=True)
    fref = np.array([(-k, k*y0[i*2]/y0[i*2+1]) for i in range(N)]).flatten()
    _test_f(rd, t0, rd.logb(y0), fref)

    jref = np.zeros((2*N, 2*N))
    for i in range(N):
        A = y0[i*2]
        B = y0[i*2+1]
        jref[i*2+1, i*2] = k/B*A
        jref[i*2+1, i*2+1] = -k/B*A
    _test_dense_jac_rmaj(rd, t0, rd.logb(y0), jref)

def test_ReactionSystem__to_ReactionDiffusion():
    sbstncs = mk_sn_dict_from_names('AB')
    r1 = Reaction({'A': 2}, {'B': 1}, 3.0)
    rsys = ReactionSystem([r1], sbstncs)
    rd = ReactionDiffusion.from_ReactionSystem(rsys)
    assert rd.stoich_active == [[0, 0]]
    assert rd.stoich_prod == [[1]]
    assert rd.k == [3.0]

def test_radyields2pdf_table():
    rsys = _get_rsys()
    rd = ReactionDiffusion.from_ReactionSystem(rsys)
    tempdir = tempfile.mkdtemp()
    try:
        radyields2pdf_table(rd, tempdir)
    finally:
        shutil.rmtree(tempdir)

def test_dc():
    n = 7
    k = [1]*(n-1)
    rd = ReactionDiffusion(n, [[i] for i in range(n-1)],
                           [[i] for i in range(1, n)], k=k)
    fout = rd.alloc_fout()
    y0 = [0]*n
    y0[0] = 1
    y_ = [0] + y0
    y0 = np.asarray(y0, dtype=np.float64)
    rd.f(0, y0, fout)
    pos, neg = [0]+k, k+[0]
    _test_f(rd, 0, y0)
    fref = [pos[i]*y_[i] - neg[i]*y_[i+1] for i in range(n)]
    assert np.allclose(fref, fout)

    jout = rd.alloc_jout(order='C')
    rd.dense_jac_rmaj(0, y0, jout)
    jref = np.zeros_like(jout)
    for i in range(n):
        if i < n - 1:
            jref[i, i] = -k[i]
        if i > 0:
            jref[i, i-1] = k[i-1]
    assert np.allclose(jref, jout)
    _test_dense_jac_rmaj(rd, 0, y0)

def test_chained_parameter_variation_from_ReactionSystem():
    g_E_mol_J = 2.1e-7
    rsys = ReactionSystem.from_string(
        """
        (H2O) -> e-(aq) + H+ + OH; Radiolytic(%.2e*mol/J)
        2 OH -> H2O2; 3.6e9/M/s
        H+ + OH- -> H2O; 1.4e11/M/s
        H2O -> H+ + OH-; 1.4e-3/s
        N2O + e-(aq) -> N2 + O-; 9.6e9/M/s
        O- + H+ -> OH; 1e11/M/s
        """ % g_E_mol_J  # neglecting a large body of reactions (just a test-case after all)
    )
    ureg = SI_base_registry
    field_u = get_derived_unit(ureg, 'doserate') * get_derived_unit(ureg, 'density')
    rd = ReactionDiffusion.from_ReactionSystem(rsys, fields=[[0*field_u]], unit_registry=ureg,
                                               param_names=['doserate'])
    dens_kg_dm3 = 0.998
    odesys = rd._as_odesys(
        variables_from_params=dict(
            density=lambda self, params: dens_kg_dm3*1e3*u.kg/u.m**3
        )
    )
    npoints = 5
    durations = [59*u.second, 42*u.minute, 2*u.hour]
    doserates = [135*u.Gy/u.s, 11*u.Gy/u.s, 180*u.Gy/u.minute]
    M = u.molar
    ic = defaultdict(lambda: 0*M, {'H2O': 55.4*M, 'H+': 1e-7*M, 'OH-': 1e-7*M, 'N2O': 20e-3*M})

    result = odesys.chained_parameter_variation(durations, ic, {'doserate': doserates}, npoints=npoints)
    ref_xout_s = [0]
    for dur in map(lambda dur: to_unitless(dur, u.s), durations):
        ref_xout_s += list(np.linspace(ref_xout_s[-1], ref_xout_s[-1] + dur, npoints+1)[1:])
    assert allclose(result.xout, ref_xout_s*u.s)

    N2_M = to_unitless(result.named_dep('N2'), u.M)
    H2O2_M = to_unitless(result.named_dep('H2O2'), u.M)

    e_accum_molar = 0
    for i, (dur, dr) in enumerate(zip(durations, doserates)):
        dur_s = to_unitless(dur, u.s)
        dr_Gy_s = to_unitless(dr, u.Gy/u.s)
        local_ts = np.linspace(0, dur_s, npoints+1)
        # local_ic = {k: result.named_dep(k)[i*npoints] for k in odesys.names}
        for j, (lt, ld) in enumerate(zip(local_ts[1:], np.diff(local_ts))):
            e_accum_molar += ld*g_E_mol_J*dr_Gy_s*dens_kg_dm3
            assert abs(N2_M[i*npoints + j + 1] - e_accum_molar)/e_accum_molar < 1e-3
            assert abs(H2O2_M[i*npoints + j + 1] - e_accum_molar)/e_accum_molar < 1e-3

    res2 = odesys.integrate(durations[0], ic, {'doserate': doserates[0]}, integrator='cvode')
    dr2 = res2.params[res2.odesys.param_names.index('doserate')]
    assert np.asarray(res2.params).shape[-1] == len(odesys.param_names)
    assert allclose(dr2, doserates[0])
    assert allclose(res2.xout[-1], durations[0])
    assert allclose(res2.named_dep('N2')[-1], durations[0]*doserates[0]*g_E_mol_J*u.mol/u.J*dens_kg_dm3*u.kg/u.dm3)
    to_unitless(res2.xout, u.s)
    to_unitless(res2.yout, u.molar)
    to_unitless(dr2, u.Gy/u.s)