Python ReactionDiffusion.efield方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import efield [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 efield [as 别名]
def integrate_rd(D=2e-3, t0=3., tend=7., x0=0.0, xend=1.0, mu=None, N=32,
                 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,
                 atol=1e-6, rtol=1e-6, efield=False, random_seed=42,
                 verbose=False, use_log2=False):
    if random_seed:
        np.random.seed(random_seed)
    n = 1
    mu = float(mu or x0)
    tout = np.linspace(t0, tend, nt)

    assert geom in 'fcs'

    # 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)

    mob = 0.3
    # Initial conditions
    y0 = {
        'f': y0_flat_cb,
        'c': y0_cylindrical_cb,
        's': y0_spherical_cb
    }[geom](x, logx)

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

    assert not lrefl
    assert not rrefl

    rd = ReactionDiffusion(
        n, stoich_active, stoich_prod, k, N,
        D=[D],
        z_chg=[1],
        mobility=[mob],
        x=x,
        geom=geom,
        logy=logy,
        logt=logt,
        logx=logx,
        nstencil=nstencil,
        lrefl=lrefl,
        rrefl=rrefl,
        use_log2=use_log2
    )

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

    # Analytic reference values
    t = tout.copy().reshape((nt, 1))
    Cref = np.repeat(y0[np.newaxis, :, np.newaxis], nt, axis=0)
    if efield:
        Cref += t.reshape((nt, 1, 1))*mob

    # Run the integration
    integr = run(rd, y0, tout, atol=atol, rtol=rtol,
                 with_jacobian=(not num_jacobian), method=method)
    Cout, info = integr.Cout, integr.info

    if verbose:
        print(info)

    def lin_err(i=slice(None), j=slice(None)):
        return integr.Cout[i, :, j] - Cref[i, :, j]

    rmsd = np.sum(lin_err()**2 / N, axis=1)**0.5
    ave_rmsd_over_atol = np.average(rmsd, axis=0)/info['atol']

    # Plot results
    if plot:
        import matplotlib.pyplot as plt

        def _plot(y, c, ttl=None, apply_exp_on_y=False):
            plt.plot(rd.xcenters, rd.expb(y) if apply_exp_on_y else y, c=c)
            if N < 100:
                plt.vlines(rd.x, 0, np.ones_like(rd.x)*max(y), linewidth=.1,
                           colors='gray')
            plt.xlabel('x / m')
            plt.ylabel('C / M')
            if ttl:
                plt.title(ttl)

        for i in range(nt):
            c = 1-tout[i]/tend
            c = (1.0-c, .5-c/2, .5-c/2)  # over time: dark red -> light red

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

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import efield [as 别名]
def integrate_rd(D=2e-3, t0=3., tend=7., x0=0.0, xend=1.0, mu=None, N=64,
                 nt=42, geom='f', logt=False, logy=False, logx=False,
                 random=False, k=0.0, nstencil=3, linterpol=False,
                 rinterpol=False, num_jacobian=False, method='bdf',
                 scale_x=False, atol=1e-6, rtol=1e-6,
                 efield=False, random_seed=42):
    if t0 == 0.0:
        raise ValueError("t0==0 => Dirac delta function C0 profile.")
    if random_seed:
        np.random.seed(random_seed)
    decay = (k != 0.0)
    n = 2 if decay else 1
    mu = float(mu or x0)
    tout = np.linspace(t0, tend, nt)

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

    # Setup the grid
    _x0 = log(x0) if logx else x0
    _xend = log(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)

    rd = ReactionDiffusion(
        2 if decay else 1,
        [[0]] if decay else [],
        [[1]] if decay else [],
        [k] if decay else [],
        N,
        D=[D]*(2 if decay else 1),
        z_chg=[1]*(2 if decay else 1),
        mobility=[0.01]*(2 if decay else 1),
        x=x,
        geom=geom,
        logy=logy,
        logt=logt,
        logx=logx,
        nstencil=nstencil,
        lrefl=not linterpol,
        rrefl=not rinterpol,
        xscale=1/(x[1]-x[0]) if scale_x else 1.0
    )

    if efield:
        if geom != FLAT:
            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, mu, x0, xend,
                    0.01 if efield else 0, logy, logx).reshape(nt, N, 1)

    if decay:
        yref = np.concatenate((yref, yref), axis=2)
        if logy:
            yref[:, :, 0] += -k*t
            yref[:, :, 1] += np.log(1-np.exp(-k*t))
        else:
            yref[:, :, 0] *= np.exp(-k*t)
            yref[:, :, 1] *= 1-np.exp(-k*t)

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