Python ReactionDiffusion.nondimensionalisation方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import nondimensionalisation [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 nondimensionalisation [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 nondimensionalisation [as 别名]
def main(logy=False, logt=False, unit_registry=None):
    # A -> B
    n = 2
    k0 = 0.13 * second**-1
    if unit_registry is None:
        unit_registry = SI_base_registry
    rd = ReactionDiffusion.nondimensionalisation(
        n, [[0]], [[1]], k=[k0], logy=logy, logt=logt,
        unit_registry=unit_registry)
    y0 = [3.0e3 * mole/metre**3, 1.0*molar]
    t0, tend, nt = 5.0*second, 17.0*second, 42
    tout = linspace(t0, tend, nt)

    Cref = molar*np.array([3.0*np.exp(-0.13*second**-1*(tout-t0)),
                           1.0 + 3.0*(1 - np.exp(
                               -0.13*second**-1*(tout-t0)))]).transpose()

    # scipy
    integr1 = Integration(
        rd, y0, tout, integrator='scipy')
    return integr1, Cref, rd

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import nondimensionalisation [as 别名]
def integrate_rd(
        tend=1.9, A0=4.2, B0=3.1, C0=1.4, nt=100, t0=0.0, kf=0.9, kb=0.23,
        atol='1e-7,1e-6,1e-5', rtol='1e-6', integrator='scipy', method='bdf',
        logy=False, logt=False, num_jac=False, plot=False, savefig='None',
        splitplots=False, plotlogy=False, plotsymlogy=False, plotlogt=False,
        scale_err=1.0, scaling=1.0, verbose=False):
    """
    Runs the integration and (optionally) plots:

    - Individual concentrations as function of time
    - Reaction Quotient vs. time (with equilibrium constant as reference)
    - Numerical error commited (with tolerance span plotted)
    - Excess error committed (deviation outside tolerance span)

    Concentrations (A0, B0, C0) are taken to be in "M" (molar),
    kf in "M**-1 s**-1" and kb in "s**-1", t0 and tend in "s"
    """

    rtol = float(rtol)
    atol = list(map(float, atol.split(',')))
    if len(atol) == 1:
        atol = atol[0]
    registry = SI_base_registry.copy()
    registry['amount'] = 1.0/scaling*registry['amount']
    registry['length'] = registry['length']/10  # decimetre

    kf = kf/molar/second
    kb = kb/second

    rd = ReactionDiffusion.nondimensionalisation(
        3, [[0, 1], [2]], [[2], [0, 1]], [kf, kb], logy=logy, logt=logt,
        unit_registry=registry)

    C0 = np.array([A0, B0, C0])*molar
    if plotlogt:
        eps = 1e-16
        tout = np.logspace(np.log10(t0+eps), np.log10(tend+eps), nt)*second
    else:
        tout = np.linspace(t0, tend, nt)*second

    integr = Integration(
        rd, C0, tout, integrator=integrator, atol=atol, rtol=rtol,
        with_jacobian=not num_jac, method=method)
    Cout = integr.with_units('Cout')
    yout, info = integr.yout, integr.info
    try:
        import mpmath
        assert mpmath  # silence pyflakes
    except ImportError:
        use_mpmath = False
    else:
        use_mpmath = True
    time_unit = get_derived_unit(registry, 'time')
    conc_unit = get_derived_unit(registry, 'concentration')
    Cref = _get_Cref(
        to_unitless(tout - tout[0], time_unit),
        to_unitless(C0, conc_unit),
        [to_unitless(kf, 1/time_unit/conc_unit),
         to_unitless(kb, 1/time_unit)],
        use_mpmath
    ).reshape((nt, 1, 3))*conc_unit
    if verbose:
        print(info)

    if plot:
        npltcols = 3 if splitplots else 1
        import matplotlib.pyplot as plt
        plt.figure(figsize=(18 if splitplots else 6, 10))

        def subplot(row=0, idx=0, adapt_yscale=True, adapt_xscale=True,
                    span_all_x=False):
            offset = idx if splitplots else 0
            ax = plt.subplot(4, 1 if span_all_x else npltcols,
                             1 + row*npltcols + offset)
            if adapt_yscale:
                if plotlogy:
                    ax.set_yscale('log')
                elif plotsymlogy:
                    ax.set_yscale('symlog')
            if adapt_xscale and plotlogt:
                ax.set_xscale('log')
            return ax

        tout_unitless = to_unitless(tout, second)
        c = 'rgb'
        for i, l in enumerate('ABC'):
            # Plot solution trajectory for i:th species
            ax_sol = subplot(0, i)
            ax_sol.plot(tout_unitless, to_unitless(Cout[:, 0, i], molar),
                        label=l, color=c[i])

            if splitplots:
                # Plot relative error
                ax_relerr = subplot(1, 1)
                ax_relerr.plot(
                    tout_unitless, Cout[:, 0, i]/Cref[:, 0, i] - 1.0,
                    label=l, color=c[i])
                ax_relerr.set_title("Relative error")
                ax_relerr.legend(loc='best', prop={'size': 11})

#.........这里部分代码省略.........

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import nondimensionalisation [as 别名]
def test_with_units():
    rd = ReactionDiffusion.nondimensionalisation(
        2, [[0, 0]], [[1]], [2/molar/second], unit_registry=SI_base_registry)
    assert allclose(rd.with_units('k'), [2e-3 * metre**3/mole/second])

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import nondimensionalisation [as 别名]
def test_nondimensionalisation():
    rd = ReactionDiffusion.nondimensionalisation(
        2, [[0, 0]], [[1]], [2/molar/second], unit_registry=SI_base_registry)
    assert rd.k == [2e-3]

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import nondimensionalisation [as 别名]
def integrate_rd(D=2e-3, t0=1., tend=13., x0=1e-10, xend=1.0, N=256,
                 nt=42, logt=False, logy=False, logx=False,
                 random=False, k=1.0, nstencil=3, linterpol=False,
                 rinterpol=False, num_jacobian=False, method='bdf',
                 integrator='scipy', iter_type='default',
                 linear_solver='default', atol=1e-8, rtol=1e-10, factor=1e5,
                 random_seed=42, plot=False, savefig='None', verbose=False,
                 scaling=1.0, ilu_limit=1000.0, first_step=0.0, n_jac_diags=0, use_log2=False):
    """
    Solves the time evolution of diffusion from a constant (undepletable)
    source term. Optionally plots the results. In the plots time is represented
    by color scaling from black (:math:`t_0`) to red (:math:`t_{end}`)
    """
    if t0 == 0.0:
        raise ValueError("t0==0 => Dirac delta function C0 profile.")
    if random_seed:
        np.random.seed(random_seed)
    tout = np.linspace(t0, tend, nt)

    units = defaultdict(lambda: 1)
    units['amount'] = 1.0/scaling

    # Setup the grid
    x = generate_grid(x0, xend, N, logx, random=random)
    modulation = [1 if (i == 0) else 0 for i in range(N)]

    rd = ReactionDiffusion.nondimensionalisation(
        2,
        [[0], [1]],
        [[1], [0]],
        [k, factor*k],
        N=N,
        D=[0, D],
        x=x,
        logy=logy,
        logt=logt,
        logx=logx,
        nstencil=nstencil,
        lrefl=not linterpol,
        rrefl=not rinterpol,
        modulated_rxns=[0, 1],
        modulation=[modulation, modulation],
        unit_registry=units,
        ilu_limit=ilu_limit,
        n_jac_diags=n_jac_diags,
        faraday_const=1,
        vacuum_permittivity=1,
        use_log2=use_log2
    )

    # Calc initial conditions / analytic reference values
    Cref = analytic(rd.xcenters, tout, D, x0, xend, logx, use_log2=use_log2).reshape(
        nt, N, 1)
    source = np.zeros_like(Cref[0, ...])
    source[0, 0] = factor
    y0 = np.concatenate((source, Cref[0, ...]), axis=1)

    # Run the integration
    integr = Integration(
        rd, y0, tout, integrator=integrator, atol=atol, rtol=rtol,
        with_jacobian=(not num_jacobian), method=method,
        iter_type=iter_type,
        linear_solver=linear_solver, first_step=first_step)
    Cout = integr.with_units('Cout')
    if verbose:
        import pprint
        pprint.pprint(integr.info)
    spat_ave_rmsd_over_atol = spat_ave_rmsd_vs_time(
        Cout[:, :, 1], Cref[:, :, 0]) / atol
    tot_ave_rmsd_over_atol = np.average(spat_ave_rmsd_over_atol)

    if plot:
        # Plot results
        import matplotlib.pyplot as plt
        plt.figure(figsize=(6, 10))

        def _plot(C, c, ttl=None, vlines=False,
                  smooth=True):
            if vlines:
                plt.vlines(rd.x, 0, np.ones_like(rd.x)*max(C),
                           linewidth=1, colors='gray')
            if smooth:
                plt.plot(rd.xcenters, C, c=c)
            else:
                for i, _C in enumerate(C):
                    plt.plot([rd.x[i], rd.x[i+1]], [_C, _C], c=c)

            plt.xlabel('x / m')
            plt.ylabel('C / M')
            if ttl:
                plt.title(ttl)

        for i in range(nt):
            kwargs = dict(smooth=(N >= 20),
                          vlines=(i == 0 and N < 20))

            c = 1-tout[i]/tend
            c = (1.0-c, .5-c/2, .5-c/2)
            plt.subplot(4, 1, 1)
            _plot(Cout[i, :, 1], c, 'Simulation (N={})'.format(rd.N),
#.........这里部分代码省略.........