Python ReactionDiffusion.expb方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import expb [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 expb [as 别名]
def integrate_rd(N=64, geom='f', nspecies=1, nstencil=3,
                 D=2e-3, t0=3.0, tend=7., x0=0.0, xend=1.0, center=None,
                 nt=42, logt=False, logy=False, logx=False,
                 random=False, p=0, a=0.2,
                 linterpol=False, rinterpol=False, ilu_limit=5.0,
                 n_jac_diags=-1, num_jacobian=False,
                 method='bdf', integrator='cvode', iter_type='undecided',
                 linear_solver='default',
                 atol=1e-8, rtol=1e-10,
                 efield=False, random_seed=42, mobility=0.01,
                 plot=False, savefig='None', verbose=False, yscale='linear',
                 vline_limit=100, use_log2=False, Dexpr='[D]*nspecies', check_conserv=False
                 ):  # remember: anayltic_N_scaling.main kwargs
    # Example:
    # python3 analytic_diffusion.py --plot --Dexpr "D*np.exp(10*(x[:-1]+np.diff(x)/2))"
    if t0 == 0.0:
        raise ValueError("t0==0 => Dirac delta function C0 profile.")
    if random_seed:
        np.random.seed(random_seed)
    # decay = (nspecies > 1)
    # n = 2 if decay else 1
    center = float(center or x0)
    tout = np.linspace(t0, tend, nt)

    assert geom in 'fcs'
    analytic = {
        'f': flat_analytic,
        'c': cylindrical_analytic,
        's': spherical_analytic
    }[geom]

    # Setup the grid
    logx0 = math.log(x0) if logx else None
    logxend = math.log(xend) if logx else None
    if logx and use_log2:
        logx0 /= math.log(2)
        logxend /= math.log(2)
    _x0 = logx0 if logx else x0
    _xend = logxend if logx else xend
    x = np.linspace(_x0, _xend, N+1)
    if random:
        x += (np.random.random(N+1)-0.5)*(_xend-_x0)/(N+2)

    def _k(si):
        return (si+p)*math.log(a+1)
    k = [_k(i+1) for i in range(nspecies-1)]
    rd = ReactionDiffusion(
        nspecies,
        [[i] for i in range(nspecies-1)],
        [[i+1] for i in range(nspecies-1)],
        k,
        N,
        D=eval(Dexpr),
        z_chg=[1]*nspecies,
        mobility=[mobility]*nspecies,
        x=x,
        geom=geom,
        logy=logy,
        logt=logt,
        logx=logx,
        nstencil=nstencil,
        lrefl=not linterpol,
        rrefl=not rinterpol,
        ilu_limit=ilu_limit,
        n_jac_diags=n_jac_diags,
        use_log2=use_log2
    )

    if efield:
        if geom != 'f':
            raise ValueError("Only analytic sol. for flat drift implemented.")
        rd.efield = _efield_cb(rd.xcenters)

    # Calc initial conditions / analytic reference values
    t = tout.copy().reshape((nt, 1))
    yref = analytic(rd.xcenters, t, D, center, x0, xend,
                    -mobility if efield else 0, logy, logx, use_log2).reshape(nt, N, 1)

    if nspecies > 1:
        from batemaneq import bateman_parent
        bateman_out = np.array(bateman_parent(k, tout)).T
        terminal = (1 - np.sum(bateman_out, axis=1)).reshape((nt, 1))
        bateman_out = np.concatenate((bateman_out, terminal), axis=1).reshape(
            (nt, 1, nspecies))
        if logy:
            yref = yref + rd.logb(bateman_out)
        else:
            yref = yref * bateman_out

    # Run the integration
    integr = run(rd, yref[0, ...], tout, atol=atol, rtol=rtol,
                 with_jacobian=(not num_jacobian), method=method,
                 iter_type=iter_type, linear_solver=linear_solver,
                 C0_is_log=logy, integrator=integrator)
    info = integr.info

    if logy:
        def lin_err(i, j):
            linref = rd.expb(yref[i, :, j])
            linerr = rd.expb(integr.yout[i, :, j])-linref
#.........这里部分代码省略.........

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import expb [as 别名]
def integrate_rd(D=-3e-1, t0=0.0, tend=7., x0=0.1, xend=1.0, N=1024,
                 base=0.5, offset=0.25, nt=25, geom='f',
                 logt=False, logy=False, logx=False, random=False,
                 nstencil=3, lrefl=False, rrefl=False,
                 num_jacobian=False, method='bdf', plot=False,
                 savefig='None', atol=1e-6, rtol=1e-6, random_seed=42,
                 surf_chg=(0.0, 0.0), sigma_q=101, sigma_skew=0.5,
                 verbose=False, eps_rel=80.10, use_log2=False):
    """
    A negative D (diffusion coefficent) denotes:
        mobility := -D
        D := 0
    A positive D calculates mobility from Einstein-Smoluchowski relation

    """
    assert 0 <= base and base <= 1
    assert 0 <= offset and offset <= 1
    if random_seed:
        np.random.seed(random_seed)
    n = 2

    if D < 0:
        mobility = -D
        D = 0
    else:
        mobility = electrical_mobility_from_D(D, 1, 298.15)
        print(D, mobility)

    # Setup the grid
    logb = (lambda arg: log(arg)/log(2)) if use_log2 else log

    _x0 = logb(x0) if logx else x0
    _xend = logb(xend) if logx else xend
    x = np.linspace(_x0, _xend, N+1)
    if random:
        x += (np.random.random(N+1)-0.5)*(_xend-_x0)/(N+2)

    # Setup the system
    stoich_active = []
    stoich_prod = []
    k = []

    rd = ReactionDiffusion(
        n, stoich_active, stoich_prod, k, N,
        D=[D, D],
        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):
#.........这里部分代码省略.........