Python ReactionDiffusion.banded_jac_cmaj方法代码示例

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import banded_jac_cmaj [as 别名]
def test_n_jac_diags(n_jac_diags):
    N, n, nstencil = 10, 1, 7
    rd = ReactionDiffusion(n, [], [], [], N=N, nstencil=nstencil,
                           n_jac_diags=n_jac_diags, D=[9])
    assert np.allclose(rd.xcenters,
                       [.05, .15, .25, .35, .45, .55, .65, .75, .85, .95])
    y0 = np.ones(N)

    # Dense
    jref_cdns = np.zeros((n*N, n*N), order='F')
    jout_cdns = np.zeros((n*N, n*N), order='F')
    sm = SymRD.from_rd(rd)
    sm.dense_jac(0.0, y0.flatten(), jref_cdns)
    rd.dense_jac_cmaj(0.0, y0.flatten(), jout_cdns)
    assert np.allclose(jout_cdns, jref_cdns)

    # Banded
    jref_cbnd = rd.alloc_jout(order='F', pad=0)
    jout_cbnd = rd.alloc_jout(order='F')
    sm.banded_jac(0.0, y0.flatten(), jref_cbnd)
    rd.banded_jac_cmaj(0.0, y0.flatten(), jout_cbnd)
    assert np.allclose(jout_cbnd[rd.n*rd.n_jac_diags:, :], jref_cbnd)

    # Compressed
    jref_cmprs = rd.alloc_jout_compressed()
    jout_cmprs = rd.alloc_jout_compressed()
    sm.compressed_jac(0.0, y0.flatten(), jref_cmprs)
    rd.compressed_jac_cmaj(0.0, y0.flatten(), jout_cmprs)
    assert np.allclose(jout_cmprs, jref_cmprs)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import banded_jac_cmaj [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 banded_jac_cmaj [as 别名]

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

    def cflux(si, bi):
        f = 0.0
        for k in range(nstencil):
            f += rd.D_weight[bi*nstencil+k]*y0[pxci2bi[lb[bi]+k], si]
        return D[si]*f

    r = [
        [k[0]*modulation[0][bi]*y0[bi, 0]*y0[bi, 1] for
         bi in range(N)],
        [k[1]*modulation[1][bi]*y0[bi, 3]*y0[bi, 2] for
         bi in range(N)],
        [k[2]*y0[bi, 1]**2 for bi in range(N)],
    ]

    fref = np.array([[
        -r[0][bi] + r[1][bi] + cflux(0, bi),
        -r[0][bi] + r[1][bi] - 2*r[2][bi] + cflux(1, bi),
        r[0][bi] - r[1][bi] + cflux(2, bi),
        r[2][bi] + cflux(3, bi)
    ] for bi in range(N)]).flatten()

    # Now let's check that the Jacobian is correctly computed.
    def dfdC(bi, lri, lci):
        v = 0.0
        for ri in range(len(stoich_active)):
            totl = (stoich_prod[ri].count(lri) -
                    stoich_active[ri].count(lri))
            if totl == 0:
                continue
            actv = stoich_active[ri].count(lci)
            if actv == 0:
                continue
            v += actv*totl*r[ri][bi]/y0[bi, lci]
        return v

    def jac_elem(ri, ci):
        bri, bci = ri // n, ci // n
        lri, lci = ri % n,  ci % n
        elem = 0.0

        def _diffusion():
            _elem = 0.0
            for k in range(nstencil):
                if pxci2bi[lb[bri]+k] == bci:
                    _elem += D[lri]*rd.D_weight[bri*nstencil+k]
            return _elem

        if bri == bci:
            # on block diagonal
            elem += dfdC(bri, lri, lci)
            if lri == lci:
                elem += _diffusion()
        elif bri == bci - 1:
            if lri == lci:
                elem = _diffusion()
        elif bri == bci + 1:
            if lri == lci:
                elem = _diffusion()
        return elem

    jref = np.zeros((n*N, n*N), order='C')
    for ri, ci in np.ndindex(n*N, n*N):
        jref[ri, ci] = jac_elem(ri, ci)

    # Compare to what is calculated using our C++ callback
    _test_f_and_dense_jac_rmaj(rd, 0, y0.flatten(), fref, jref)

    jout_cmaj = np.zeros((n*N, n*N), order='F')
    rd.dense_jac_cmaj(0.0, y0.flatten(), jout_cmaj)
    assert np.allclose(jout_cmaj, jref)

    ref_banded_j = get_banded(jref, n, N)

    ref_banded_j_symbolic = rd.alloc_jout(order='F', pad=0)
    symrd = SymRD.from_rd(rd)
    symrd.banded_jac(0.0, y0.flatten(), ref_banded_j_symbolic)
    assert np.allclose(ref_banded_j_symbolic, ref_banded_j)

    jout_bnd_packed_cmaj = np.zeros((3*n+1, n*N), order='F')
    rd.banded_jac_cmaj(0.0, y0.flatten(), jout_bnd_packed_cmaj)

    if os.environ.get('plot_tests', False):
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        from chemreac.util.plotting import coloured_spy
        fig = plt.figure()
        ax = fig.add_subplot(3, 1, 1)
        coloured_spy(ref_banded_j, ax=ax)
        plt.title('ref_banded_j')
        ax = fig.add_subplot(3, 1, 2)
        coloured_spy(jout_bnd_packed_cmaj[n:, :], ax=ax)
        plt.title('jout_bnd_packed_cmaj')
        ax = fig.add_subplot(3, 1, 3)
        coloured_spy(ref_banded_j-jout_bnd_packed_cmaj[n:, :], ax=ax)
        plt.title('diff')
        plt.savefig(__file__+'.png')

    assert np.allclose(jout_bnd_packed_cmaj[n:, :], ref_banded_j)

# 需要导入模块: from chemreac import ReactionDiffusion [as 别名]
# 或者: from chemreac.ReactionDiffusion import banded_jac_cmaj [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)