Python ReactionDiffusion.integrated_conc方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import integrated_conc [as 别名]
def test_integrated_conc(params):
    geom, logx = params
    N = 8192
    x0, xend = 0.11, 1.37
    x = np.linspace(x0, xend, N+1)
    rd = ReactionDiffusion(1, [], [], [], D=[0], N=N,
                           x=np.log(x) if logx else x, geom=geom, logx=logx)
    xc = rd.expb(rd.xcenters) if logx else rd.xcenters
    y = xc*np.exp(-xc)

    def primitive(t):
        if geom == 'f':
            return -(t+1)*np.exp(-t)
        elif geom == 'c':
            return 2*np.exp(-t)*np.pi*(-2 - 2*t - t**2)
        elif geom == 's':
            return 4*np.exp(-t)*np.pi*(-6 - 6*t - 3*t**2 - t**3)
        else:
            raise NotImplementedError
    res = rd.integrated_conc(y)
    ref = (primitive(xend) - primitive(x0))
    assert abs(res - ref) < 1e-8

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import integrated_conc [as 别名]

#.........这里部分代码省略.........
    info = integr.info

    if logy:
        def lin_err(i, j):
            linref = rd.expb(yref[i, :, j])
            linerr = rd.expb(integr.yout[i, :, j])-linref
            linatol = np.average(yref[i, :, j])
            linrtol = linatol
            return linerr/(linrtol*np.abs(linref)+linatol)

    if logy:
        rmsd = np.sum(lin_err(slice(None), slice(None))**2 / N, axis=1)**0.5
    else:
        rmsd = np.sum((yref-integr.yout)**2 / N, axis=1)**0.5
    ave_rmsd_over_atol = np.average(rmsd, axis=0)/atol

    if verbose:
        # Print statistics
        from pprint import pprint
        pprint(info)
        pprint(ave_rmsd_over_atol)

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

        # colors: (0.5, 0.5, 0.5), (0.5, 0.5, 1), ...
        base_colors = list(product([.5, 1], repeat=3))[1:-1]

        def color(ci, ti):
            return np.array(base_colors[ci % len(base_colors)])*tout[ti]/tend

        for ti in range(nt):
            plt.subplot(4, 1, 1)
            for si in range(nspecies):
                plt.plot(rd.xcenters, integr.Cout[ti, :, si], c=color(si, ti),
                         label=None if ti < nt - 1 else rd.substance_names[si])

            plt.subplot(4, 1, 2)
            for si in range(nspecies):
                plt.plot(rd.xcenters, rd.expb(yref[ti, :, si]) if logy
                         else yref[ti, :, si], c=color(si, ti))

            plt.subplot(4, 1, 3)
            if logy:
                for si in range(nspecies):
                    plt.plot(rd.xcenters, lin_err(ti, si)/atol,
                             c=color(si, ti))
            else:
                for si in range(nspecies):
                    plt.plot(
                        rd.xcenters,
                        (yref[ti, :, si] - integr.yout[ti, :, si])/atol,
                        c=color(si, ti))

        if N < vline_limit:
            for idx in range(1, 4):
                plt.subplot(4, 1, idx)
                for bi in range(N):
                    plt.axvline(rd.x[bi], color='gray')

        plt.subplot(4, 1, 1)
        plt.title('Simulation (N={})'.format(rd.N))
        plt.xlabel('x / m')
        plt.ylabel('C / M')
        plt.gca().set_yscale(yscale)
        plt.legend()

        plt.subplot(4, 1, 2)
        plt.title('Analytic solution')
        plt.gca().set_yscale(yscale)

        plt.subplot(4, 1, 3)
        plt.title('Linear rel. error / Abs. tol. (={})'.format(atol))

        plt.subplot(4, 1, 4)
        plt.title('RMS error vs. time'.format(atol))
        tspan = [tout[0], tout[-1]]
        for si in range(nspecies):
            plt.plot(tout, rmsd[:, si] / atol, c=color(si, -1))
            plt.plot(tspan, [ave_rmsd_over_atol[si]]*2,
                     c=color(si, -1), ls='--')

        plt.xlabel('Time / s')
        plt.ylabel(r'$\sqrt{\langle E^2 \rangle} / atol$')
        plt.tight_layout()
        save_and_or_show_plot(savefig=savefig)

    if check_conserv:
        tot_amount = np.zeros(tout.size)
        for ti in range(tout.size):
            for si in range(nspecies):
                tot_amount[ti] += rd.integrated_conc(integr.yout[ti, :, si])
        if plot:
            plt.plot(tout, tot_amount)
            plt.show()
        assert np.allclose(tot_amount[0], tot_amount[1:])

    return tout, integr.yout, info, ave_rmsd_over_atol, rd, rmsd

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import integrated_conc [as 别名]

#.........这里部分代码省略.........
        z_chg=[1, -1],
        mobility=[mobility, -mobility],
        x=x,
        geom=geom,
        logy=logy,
        logt=logt,
        logx=logx,
        nstencil=nstencil,
        lrefl=lrefl,
        rrefl=rrefl,
        auto_efield=True,
        surf_chg=surf_chg,
        eps_rel=eps_rel,  # water at 20 deg C
        faraday_const=1,
        vacuum_permittivity=1,
        use_log2=use_log2
    )

    # Initial conditions
    sigma = (xend-x0)/sigma_q
    sigma = [(1-sigma_skew)*sigma, sigma_skew*sigma]
    y0 = np.vstack(pair_of_gaussians(
        rd.xcenters, [base+offset, base-offset], sigma, logy, logx, geom, use_log2)).transpose()
    if logy:
        y0 = sigm(y0)

    if plot:
        # Plot initial E-field
        import matplotlib.pyplot as plt
        plt.figure(figsize=(6, 10))
        rd.calc_efield((rd.expb(y0) if logy else y0).flatten())
        plt.subplot(4, 1, 3)
        plt.plot(rd.xcenters, rd.efield, label="E at t=t0")
        plt.plot(rd.xcenters, rd.xcenters*0, label="0")

    # Run the integration
    tout = np.linspace(t0, tend, nt)
    integr = run(rd, y0, tout,
                 atol=atol, rtol=rtol, sigm_damp=True,
                 C0_is_log=logy,
                 with_jacobian=(not num_jacobian), method=method)
    Cout = integr.Cout

    if verbose:
        print(integr.info)
    # Plot results
    if plot:
        def _plot(y, ttl=None,  **kwargs):
            plt.plot(rd.xcenters, y, **kwargs)
            plt.xlabel((('log_%s({})' % ('2' if use_log2 else 'e')) if logx else '{}').format('x / m'))
            plt.ylabel('C / M')
            if ttl:
                plt.title(ttl)

        for i in range(nt):
            plt.subplot(4, 1, 1)
            c = 1-tout[i]/tend
            c = (1.0-c, .5-c/2, .5-c/2)
            _plot(Cout[i, :, 0], 'Simulation (N={})'.format(rd.N),
                  c=c, label='$z_A=1$' if i == nt-1 else None)
            _plot(Cout[i, :, 1], c=c[::-1],
                  label='$z_B=-1$' if i == nt-1 else None)
            plt.legend()

            plt.subplot(4, 1, 2)
            delta_y = Cout[i, :, 0] - Cout[i, :, 1]
            _plot(delta_y, 'Diff'.format(rd.N),
                  c=[c[2], c[0], c[1]],
                  label='A-B (positive excess)' if i == nt-1 else None)
            plt.legend(loc='best')
            plt.xlabel("$x~/~m$")
            plt.ylabel(r'Concentration / M')
        ylim = plt.gca().get_ylim()
        if N < 100:
            plt.vlines(rd.x, ylim[0], ylim[1],
                       linewidth=1.0, alpha=0.2, colors='gray')

        plt.subplot(4, 1, 3)
        plt.plot(rd.xcenters, rd.efield, label="E at t=tend")
        plt.xlabel("$x~/~m$")
        plt.ylabel("$E~/~V\cdot m^{-1}$")
        plt.legend()

        for i in range(3):
            plt.subplot(4, 1, i+1)
            ylim = plt.gca().get_ylim()
            for d in (-1, 1):
                center_loc = [x0+(base+d*offset)*(xend-x0)]*2
                plt.plot(rd.logb(center_loc) if logx else center_loc,
                         ylim, '--k')
        plt.subplot(4, 1, 4)
        for i in range(n):
            amount = [rd.integrated_conc(Cout[j, :, i]) for j in range(nt)]
            plt.plot(tout, amount, c=c[::(1, -1)[i]], label=chr(ord('A')+i))
        plt.xlabel('Time / s')
        plt.ylabel('Amount / mol')
        plt.legend(loc='best')
        plt.tight_layout()
        save_and_or_show_plot(savefig=savefig)
    return tout, Cout, integr.info, rd