Python ReactionDiffusion.from_ReactionSystem方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
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)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
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)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
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)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
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]

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_radyields2pdf_table():
    rsys = _get_rsys()
    rd = ReactionDiffusion.from_ReactionSystem(rsys)
    tempdir = tempfile.mkdtemp()
    try:
        radyields2pdf_table(rd, tempdir)
    finally:
        shutil.rmtree(tempdir)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
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)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_from_ReactionSystem__g_values():
    from chempy import ReactionSystem as RS
    rs = RS.from_string('-> H + OH; Radiolytic(2.1e-7)', checks=())
    rd = ReactionDiffusion.from_ReactionSystem(rs, variables={'density': 998, 'doserate': 0.15})
    gv = rd.g_values
    assert len(gv) == 1
    assert np.allclose(gv[0], rs.as_per_substance_array({'H': 2.1e-7, 'OH': 2.1e-7}))
    assert len(rd.fields) == 1
    assert len(rd.fields[0]) == 1
    assert np.allclose(rd.fields[0][0], 998*0.15)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_chemistry():
    sbstncs = mk_sn_dict_from_names('ABC', D=[.1, .2, .3])
    r1 = Reaction({'A': 1, 'B': 1}, {'C': 1}, 0.3)
    rsys = ReactionSystem([r1], sbstncs)
    rd = ReactionDiffusion.from_ReactionSystem(rsys)
    serialized_rd = load(JSON_PATH)
    assert rd.stoich_active == serialized_rd.stoich_active
    assert rd.stoich_prod == serialized_rd.stoich_prod
    assert rd.stoich_inact == serialized_rd.stoich_inact
    assert np.allclose(rd.k, serialized_rd.k)
    assert np.allclose(rd.D, serialized_rd.D)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_autobinary():
    from chemreac.chemistry import (
        Reaction, ReactionSystem, mk_sn_dict_from_names
    )
    sbstncs = mk_sn_dict_from_names('AB')
    k = 3.0
    r1 = Reaction({'A': 2}, {'B': 1}, k)
    rsys = ReactionSystem([r1], sbstncs)
    rd = ReactionDiffusion.from_ReactionSystem(rsys)

    _test_f_and_dense_jac_rmaj(rd, 0, [1, 37], [-2*3, 3])

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_from_ReactionSystem__g_values__units():
    from chempy import ReactionSystem as RS
    from chempy.units import SI_base_registry, default_units as u
    rs = RS.from_string('-> H + OH; Radiolytic(2.1*per100eV)', checks=())
    variables = {'density': .998 * u.kg/u.dm3, 'doserate': 0.15*u.Gy/u.s}
    rd = ReactionDiffusion.from_ReactionSystem(rs, variables=variables, unit_registry=SI_base_registry)
    gv = rd.g_values
    per100eV_as_mol_per_joule = 1.0364268556366418e-07
    ref = 2.1 * per100eV_as_mol_per_joule
    assert len(gv) == 1
    assert np.allclose(gv[0], rs.as_per_substance_array({'H': ref, 'OH': ref}))
    assert len(rd.fields) == 1
    assert len(rd.fields[0]) == 1
    assert np.allclose(rd.fields[0][0], 998*0.15)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_autodimerization(arrhenius):
    # A + A -> B
    from chemreac.chemistry import (
        Reaction, ReactionSystem, mk_sn_dict_from_names
    )
    sbstncs = mk_sn_dict_from_names('AB')
    if arrhenius:
        from chempy.kinetics.arrhenius import ArrheniusParam
        k = ArrheniusParam(7e11, -8.314472*298.15*(np.log(3) - np.log(7) - 11*np.log(10)))
        variables = {'temperature': 298.15}
        param = 3.0
    else:
        param = k = 3.0
        variables = None
    r1 = Reaction({'A': 2}, {'B': 1}, k)
    rsys = ReactionSystem([r1], sbstncs)
    rd = ReactionDiffusion.from_ReactionSystem(rsys, variables=variables)
    t = np.linspace(0, 5, 3)
    A0, B0 = 1.0, 0.0
    integr = run(rd, [A0, B0], t)
    Aref = 1/(1/A0+2*param*t)
    yref = np.vstack((Aref, (A0-Aref)/2)).transpose()
    assert np.allclose(integr.yout[:, 0, :], yref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def test_from_ReactionSystem__fields():
    from chempy import Reaction as R, ReactionSystem as RS
    from chempy.kinetics.rates import mk_Radiolytic
    gvals = [2, 3, 5]
    rs = RS([R({}, {'H', 'OH'}, Rad(v)) for v, Rad in zip(
        gvals, map(mk_Radiolytic, 'alpha beta gamma'.split()))])
    for s in rs.substances.values():
        s.data['D'] = 0.0  # no diffusion
    nbins = 73
    fields = OrderedDict([('alpha', np.linspace(7, 11, nbins)),
                          ('beta', np.linspace(13, 17, nbins)),
                          ('gamma', np.linspace(19, 23, nbins))])
    rd = ReactionDiffusion.from_ReactionSystem(rs, variables={'density': 998}, fields=fields, N=nbins)

    C0 = {'H': np.linspace(29, 31, nbins), 'OH': np.linspace(37, 41, nbins)}
    integr = Integration(rd, np.array([C0[k]*np.ones(nbins) for k in rd.substance_names]).T,
                         [0, 43], integrator='cvode')
    assert integr.tout[-1] == 43
    assert integr.tout.size > 2
    Cref = np.array([C0[sk] + integr.tout.reshape((-1, 1))*reduce(add, [
        fields[fk]*g for fk, g in zip(fields, gvals)
    ]).reshape((1, -1)) for sk in rd.substance_names]).transpose(1, 2, 0)
    assert np.allclose(integr.Cout, Cref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import from_ReactionSystem [as 别名]
def integrate_rd(tend=1e2, A0=1.0, B0=0.0, C0=0.0, k1=0.04, k2=1e4, k3=3e7,
                 t0=1e2, nt=100, N=1, nstencil=3, logt=False, logy=False,
                 plot=False, savefig='None', verbose=False, dump_expr='False',
                 use_chempy=False, D=2e-3):
    if N == 1:
        init_conc = (A0, B0, C0)
    else:
        init_conc = np.tile((A0, B0, C0), (N, 1))
        init_conc /= np.linspace(1, 2, N).reshape((N, 1))**.5

    rsys = ReactionSystem(get_reactions((k1, k2, k3)), 'ABC')
    if verbose:
        print([str(_) for _ in rsys.rxns])
    if use_chempy:
        from chempy.kinetics.ode import get_odesys
        odesys = get_odesys(rsys, include_params=True)
        if N != 1:
            raise ValueError("ChemPy does not support diffusion")
        odesys.integrate(np.logspace(log10(t0), log10(tend)), init_conc)
        if plot:
            odesys.plot_result(xscale='log', yscale='log')
        result = None
    else:
        rd = ReactionDiffusion.from_ReactionSystem(
            rsys, N=N, nstencil=1 if N == 1 else nstencil, logt=logt,
            logy=logy, D=[D/2, D/3, D/5])

        if dump_expr.lower() not in ('false', '0'):
            from chemreac.symbolic import SymRD
            import sympy as sp
            cb = {'latex': sp.latex,
                  'ccode': sp.ccode}.get(dump_expr.lower(), str)
            srd = SymRD.from_rd(rd, k=sp.symbols('k:3'))
            print('dydx:')
            print('\n'.join(map(cb, srd._f)))
            print('jac:')
            for ri, row in enumerate(srd.jacobian.tolist()):
                for ci, expr in enumerate(row):
                    if expr == 0:
                        continue
                    print(ri, ci, cb(expr))
            return None

        if t0 == 0 and logt:
            t0 = 1e-3*suggest_t0(rd, init_conc)
            if verbose:
                print("Using t0 = %12.5g" % t0)

        t = np.logspace(np.log10(t0), np.log10(tend), nt)
        print(t[0], t[-1])
        integr = run(rd, init_conc, t)

        if verbose:
            import pprint
            pprint.pprint(integr.info)

        if plot:
            if N == 1:
                plot_C_vs_t(integr, xscale='log', yscale='log')
            else:
                import matplotlib.pyplot as plt
                for idx, name in enumerate('ABC', 1):
                    plt.subplot(1, 3, idx)
                    rgb = [.5, .5, .5]
                    rgb[idx-1] = 1
                    plot_faded_time(integr, name, rgb=rgb, log_color=True)
        result = integr

    if plot:
        save_and_or_show_plot(savefig=savefig)

    return result