Python ReactionDiffusion.logb方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_reaction__logy(N):
    # See <test_ReactionDiffusion__only_1_reaction__logy.png>
    t0 = 3.0
    y0 = np.array([2.0, 3.0]*N)
    k = 5.0
    # A -> B
    rd = ReactionDiffusion(2, [[0]], [[1]], [k], N, D=[0.0, 0.0], logy=True)
    fref = np.array([(-k, k*y0[i*2]/y0[i*2+1]) for i in range(N)]).flatten()
    _test_f(rd, t0, rd.logb(y0), fref)

    jref = np.zeros((2*N, 2*N))
    for i in range(N):
        A = y0[i*2]
        B = y0[i*2+1]
        jref[i*2+1, i*2] = k/B*A
        jref[i*2+1, i*2+1] = -k/B*A
    _test_dense_jac_rmaj(rd, t0, rd.logb(y0), jref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_reaction__logy__logt(N):
    # See <test_ReactionDiffusion__only_1_reaction__logy_logt.png>
    t0 = 3.0
    y0 = np.array([2.0, 3.0]*N)
    k = 5.0
    # A -> B
    rd = ReactionDiffusion(2, [[0]], [[1]], [k], N, D=[0.0, 0.0],
                           logy=True, logt=True)
    fref = np.array([(-k*t0, t0*k*y0[i*2]/y0[i*2+1])
                     for i in range(N)]).flatten()
    _test_f_and_dense_jac_rmaj(rd, rd.logb(t0), np.log(y0), fref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_species_diffusion_3bins(log):
    # Diffusion without reaction
    # 3 bins
    t0 = 3.0
    logy, logt = log
    D = 17.0
    y0 = np.array([23.0, 27.0, 37.0])
    N = 3
    x = [5.0, 7.0, 13.0, 15.0]
    xc = [4.0, 6.0, 10.0, 14.0, 16.0]
    nstencil = 3
    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x, N=N, logy=logy,
                           logt=logt, lrefl=False, rrefl=False,
                           nstencil=nstencil)
    assert np.allclose(rd.xc, xc)

    w = [1/16, -1/8, 1/16]  # finite diff. weights for 2nd order deriv
    for i in range(N):
        assert np.allclose(rd.D_weight[i*nstencil:(i+1)*nstencil], w)
    J = D*(w[0]*y0[0] + w[1]*y0[1] + w[2]*y0[2])
    fref = np.array([J, J, J])

    if logy:
        fref /= y0
    if logt:
        fref *= t0

    if logy:
        jref = D*np.array([  # jref[i, k] = ...
            [w[k]*y0[k]/y0[i] if k != i else -1/y0[k]*sum([
                w[j]*y0[j] if j != k else 0 for j in range(3)
            ]) for k in range(3)] for i in range(3)
        ])
        jref[0, 2] = 0.0  # dense_jac_rmaj only computes banded approx.
        jref[2, 0] = 0.0  # same as above.
    else:
        jref = D*np.array([[w[k] if abs(k-i) < 2 else 0.0 for
                          k in range(3)] for i in range(3)])

    if logt:
        jref *= t0

    y = rd.logb(y0) if logy else y0
    t = rd.logb(t0) if logt else t0
    _test_f_and_dense_jac_rmaj(rd, t, y, fref, jref)

    jout_bnd = np.zeros((4, 3), order='F')
    rd.banded_jac_cmaj(t, y, jout_bnd)
    jref_bnd = get_banded(jref, 1, 3)
    assert np.allclose(jout_bnd[1:, :], jref_bnd)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_integrators(log):
    logy, logt, use_log2 = log
    t0, tend, nt = 5.0, 17.0, 42
    tout = np.linspace(t0, tend, nt+1)

    # Update the dict if more integrators are added:
    solver_kwargs = {
        'scipy1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'tout': tout
        },
        'scipy2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'tout': (t0, tend),
            'dense_output': True
        },
        'cvode1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': tout
        },
        'cvode2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'method': 'adams',
            'tout': tout,
            'C0_is_log': True
        }
    }

    # A -> B
    n = 2
    k0 = 0.13
    rd = ReactionDiffusion(n, [[0]], [[1]], k=[k0], logy=logy, logt=logt, use_log2=use_log2)
    y0 = [3.0, 1.0]

    results = []
    for solver, kwargs in solver_kwargs.items():
        _y0 = rd.logb(y0) if kwargs.get('C0_is_log', False) else y0
        integr = Integration(rd, _y0, integrator=solver[:-1], **kwargs)
        if not kwargs.get('dense_output', False):
            results.append(integr.Cout)

    for result in results[1:]:
        assert np.allclose(results[0][0], result[0])

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_field_dep_reaction_logy_logt(N):
    t0 = 3.0
    y0 = np.concatenate([np.array([2.0, 3.0])/(x+1) for x in range(N)])
    k = 5.0
    # A -> B

    rd = ReactionDiffusion(2, [], [], [], N, D=[0.0, 0.0],
                           fields=[[x+1 for x in range(N)]],
                           g_values=[[-k, k]], g_value_parents=[0],
                           logy=True, logt=True)

    def k_(bi):
        return k*(bi+1)

    fref = np.array([(-k_(i)*t0, k_(i)*t0*y0[i*2]/y0[i*2+1])
                     for i in range(N)]).flatten()
    _test_f_and_dense_jac_rmaj(rd, rd.logb(t0), np.log(y0), fref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_field_dep_reaction_logy(N):
    y0 = np.concatenate([np.array([2.0, 3.0])/(x+1) for x in range(N)])
    k = 5.0
    # A -> B

    rd = ReactionDiffusion(2, [], [], [], N, D=[0.0, 0.0],
                           fields=[[x+1 for x in range(N)]],
                           g_values=[[-k, k]],
                           g_value_parents=[0], logy=True)

    def k_(bi):
        return k*(bi+1)

    fref = np.array([(-k_(i), k_(i)*y0[i*2]/y0[i*2+1])
                     for i in range(N)]).flatten()
    if N == 1:
        jref = np.array([[0, 0],
                         [k*y0[0]/y0[1], -k*y0[0]/y0[1]]])
    else:
        jref = None
    _test_f_and_dense_jac_rmaj(rd, 0, rd.logb(y0), fref, jref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [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 logb [as 别名]
def test_integrators(log):
    logy, logt, use_log2 = log
    t0, tend, nt = 5.0, 17.0, 42
    tout = np.linspace(t0, tend, nt+1)

    # Update the dict if more integrators are added:
    solver_kwargs = {
        'scipy1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'tout': tout
        },
        'scipy2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'tout': (t0, tend),
            'dense_output': True
        },
        'cvode1': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': tout
        },
        'cvode2': {
            'atol': 1e-8,
            'rtol': 1e-8,
            'method': 'adams',
            'tout': tout,
            'C0_is_log': True,
            'ew_ele': True
        },
        'cvode3': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': (t0, tend),
            'ew_ele': True
        },
        'cvode4': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': tend-t0,
        },
        'cvode5': {
            'atol': [1e-8, 1e-8],
            'rtol': 1e-8,
            'method': 'bdf',
            'tout': (t0, tend),
            'constraints': [1.0, 1.0]
        }
    }
    import pycvodes
    if logy or pycvodes.sundials_version < (3, 2, 0):
        solver_kwargs.pop('cvode5')  # sundials >=3.2.0 required for constraints

    # A -> B
    n = 2
    k0 = 0.13
    rd = ReactionDiffusion(n, [[0]], [[1]], k=[k0], logy=logy, logt=logt, use_log2=use_log2)
    y0 = [3.0, 1.0]

    results = []
    for solver, kwargs in solver_kwargs.items():
        _y0 = rd.logb(y0) if kwargs.get('C0_is_log', False) else y0
        integr = Integration(rd, _y0, integrator=solver[:-1], **kwargs)
        if not kwargs.get('dense_output', False):
            results.append(integr.Cout)
        ew_ele = integr.info.get('ew_ele', None)
        if ew_ele is not None:
            assert np.all(np.abs(np.prod(ew_ele, axis=1)) < 2)

    for result in results[1:]:
        assert np.allclose(results[0][0], result[0])

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__only_1_species_diffusion_7bins(log):
    # Diffusion without reaction
    N = 7
    nstencil = 5
    nsidep = (nstencil-1)//2
    t0 = 3.0
    logy, logt = log
    D = 2.0
    y0 = np.array([12, 8, 11, 5, 7, 4, 9], dtype=np.float64)
    x = np.array([3, 5, 13, 17, 23, 25, 35, 37], dtype=np.float64)
    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x,
                           logy=logy, logt=logt, nstencil=nstencil,
                           lrefl=False, rrefl=False)
    weights = [
        [951/8800, -716/2475, 100/297, -75/352, 311/5400],
        [321/8800, -161/2475, 7/297, 3/352, -19/5400],
        [-39/8800, 109/2475, -127/1485, 87/1760, -19/5400],
        [-2/693, 38/675, -129/1100, 7/108, -1/1050],
        [0, 9/160, -7/72, 2/45, -1/288],
        [-8/1575, 9/400, 0, -19/450, 25/1008],
        [16/315, -9/32, 31/72, -13/45, 179/2016]
    ]
    assert np.allclose(rd.D_weight, np.array(weights).flatten())

    lb = stencil_pxci_lbounds(nstencil, N)
    yi = pxci_to_bi(nstencil, N)
    fref = np.array([sum([D*weights[i][j]*y0[yi[j+lb[i]]] for j
                          in range(nstencil)]) for i in range(N)])

    if logy:
        fref /= y0
    if logt:
        fref *= t0

    jref = np.zeros((N, N))
    for i in range(N):
        for j in range(max(0, i-1), min(N, i+2)):
            if logy:
                if j == i+1 or j == i-1:
                    for k in range(nstencil):
                        if yi[k+lb[i]] == j:
                            jref[i, j] += D*weights[i][k]*y0[j]/y0[i]
                else:  # j == i
                    assert i == j
                    for k in range(nstencil):
                        cyi = yi[k+lb[i]]
                        if i == cyi:
                            continue
                        jref[i, i] -= D*weights[i][k]*y0[cyi]/y0[i]
            else:
                if i-1 <= j and j <= i+1:
                    jref[i, j] = D*weights[i][j-lb[i]+nsidep]
    if logt:
        jref *= t0
    t = rd.logb(t0) if logt else t0
    y = rd.logb(y0) if logy else y0
    _test_f_and_dense_jac_rmaj(rd, t, y, fref, jref)

    jout_bnd = np.zeros((4, N), order='F')
    rd.banded_jac_cmaj(t, y, jout_bnd)
    jref_bnd = get_banded(jref, 1, N)
    assert np.allclose(jout_bnd[1:, :], jref_bnd)

    # compressed_jac_cmaj actually use all diagonals
    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x,
                           logy=logy, logt=logt, nstencil=nstencil,
                           lrefl=False, rrefl=False, n_jac_diags=2)
    jout_cmprs = rd.alloc_jout_compressed()
    rd.compressed_jac_cmaj(t, y, jout_cmprs)
    from block_diag_ilu import Compressed_from_dense

    jref2 = np.zeros((N, N))
    for i in range(N):
        for j in range(max(0, i-2), min(N, i+3)):
            if logy:
                if i-2 <= j <= i+2:
                    if i == j:
                        for k in range(nstencil):
                            cyi = yi[k+lb[i]]
                            if i == cyi:
                                continue
                            jref2[i, i] -= D*weights[i][k]*y0[cyi]/y0[i]
                    else:
                        for k in range(nstencil):
                            if yi[k+lb[i]] == j:
                                jref2[i, j] += D*weights[i][k]*y0[j]/y0[i]

            else:
                if i-2 <= j <= i+2:
                    jref2[i, j] = D*weights[i][j-lb[i]+nsidep]
    if logt:
        jref2 *= t0
    jref_cmprs = Compressed_from_dense(jref2, N, 1, nsidep).data
    assert np.allclose(jout_cmprs, jref_cmprs)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__lrefl_7(log):
    # Diffusion without reaction (7 bins)
    from sympy import finite_diff_weights
    lrefl, rrefl = True, False
    N = 7
    t0 = 3.0
    logy, logt = log
    D = 17.0
    nstencil = 5
    nsidep = (nstencil-1)//2
    x = np.array([3, 5, 13, 17, 23, 25, 35, 37], dtype=np.float64)
    y0 = np.array([12, 8, 11, 5, 7, 4, 9], dtype=np.float64)
    xc_ = padded_centers(x, nsidep)
    lb = stencil_pxci_lbounds(nstencil, N, lrefl=lrefl, rrefl=rrefl)

    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x, logy=logy,
                           logt=logt, N=N, nstencil=nstencil,
                           lrefl=lrefl, rrefl=rrefl)

    y = rd.logb(y0) if logy else y0
    t = rd.logb(t0) if logt else t0
    le = nsidep if lrefl else 0
    D_weight_ref = np.array([
        finite_diff_weights(
            2, xc_[lb[i]:lb[i]+nstencil], x0=xc_[le+i])[-1][-1]
        for i in range(N)], dtype=np.float64)
    assert np.allclose(rd.D_weight, D_weight_ref.flatten())

    yi = pxci_to_bi(nstencil, N)

    fref = D*np.array([
        sum([rd.D_weight[i*nstencil+j]*y0[yi[lb[i]+j]]
             for j in range(nstencil)])
        for i in range(N)
    ])

    if logy:
        fref /= y0
    if logt:
        fref *= t0

    _test_f(rd, t, y, fref)

    if logy:
        def cb(i, j):
            if abs(i-j) > 1:
                return 0  # imperfect Jacobian
            elif i == j:
                res = 0
                for k in range(nstencil):
                    cyi = yi[lb[i] + k]  # current y index
                    if cyi != i:
                        res -= y0[cyi]/y0[i]*rd.D_weight[i*nstencil + k]
                return res
            else:
                res = 0
                for k in range(nstencil):
                    cyi = yi[lb[i] + k]  # current y index
                    if cyi == j:
                        res += y0[j]/y0[i]*rd.D_weight[i*nstencil + k]
                return res
    else:
        def cb(i, j):
            if abs(i-j) > 1:
                return 0  # imperfect Jacobian
            res = 0
            for k in range(nstencil):
                if yi[lb[i]+k] == j:
                    res += rd.D_weight[i*nstencil + k]
            return res
    jref = D*np.array([cb(i, j) for i, j in product(
        range(N), range(N))]).reshape(N, N)
    if logt:
        jref *= t0

    _test_dense_jac_rmaj(rd, t, y, jref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [as 别名]
def test_ReactionDiffusion__rrefl_3(log):
    # Diffusion without reaction
    # 3 bins
    t0 = 3.0
    logy, logt, use_log2 = log
    D = 17.0
    y0 = np.array([23.0, 27.0, 37.0])
    x = [5.0, 9.0, 13.0, 15.0]
    nstencil = 3
    xc = [3.0, 7.0, 11.0, 14.0, 16.0]

    rd = ReactionDiffusion(1, [], [], [], D=[D], x=x, logy=logy,
                           nstencil=nstencil, logt=logt,
                           lrefl=False, rrefl=True, use_log2=use_log2)
    assert np.allclose(rd.xc, xc)

    # In [7]: xlst=[3, 7, 11, 14, 16]
    # In [8]: print(finite_diff_weights(2, xlst[1:4], x0=xlst[1])[-1][-1])
    # [1/14, -1/6, 2/21]
    # In [9]: print(finite_diff_weights(2, xlst[1:4], x0=xlst[2])[-1][-1])
    # [1/14, -1/6, 2/21]
    # In [10]: print(finite_diff_weights(2, xlst[2:5], x0=xlst[3])[-1][-1])
    # [2/15, -1/3, 1/5]
    D_weight_ref = np.array([1/14, -1/6, 2/21, 1/14,
                             -1/6, 2/21, 2/15, -1/3, 1/5])
    assert np.allclose(rd.D_weight, D_weight_ref)

    fref = np.array([
        1/14*y0[0] - 1/6*y0[1] + 2/21*y0[2],  # lrefl=False
        1/14*y0[0] - 1/6*y0[1] + 2/21*y0[2],
        2/15*y0[1] - 1/3*y0[2] + 1/5*y0[2],  # rrefl=True
    ])*D

    if logy:
        fref /= y0
        if not logt and use_log2:
            fref /= np.log(2)
    if logt:
        fref *= t0
        if not logy and use_log2:
            fref *= np.log(2)

    y = rd.logb(y0) if logy else y0
    t = rd.logb(t0) if logt else t0
    _test_f(rd, t, y, fref)

    if logy:
        jref = D*np.array([
            [
                1/6*y0[1]/y0[0] - 2/21*y0[2]/y0[0],
                -1/6*y0[1]/y0[0],
                2/21*y0[2]/y0[0]
            ],
            [
                1/14*y0[0]/y0[1],
                -1/14*y0[0]/y0[1] - 2/21*y0[2]/y0[1],
                2/21*y0[2]/y0[1]
            ],
            [
                0,
                2/15*y0[1]/y0[2],
                -2/15*y0[1]/y0[2]
            ]
        ])
    else:
        jref = D*np.array([
            [1/14, -1/6, 2/21],
            [1/14, -1/6, 2/21],
            [0,    2/15, 1/5-1/3]
        ])

    if logt:
        jref *= t0
        if use_log2:
            jref *= np.log(2)

    _test_dense_jac_rmaj(rd, t, y, jref)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import logb [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):
#.........这里部分代码省略.........