def stress_util2d_for_joe_the_file_king(ml, nlay, nrow, ncol):
    dis = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol)
    hk = np.ones((nlay, nrow, ncol))
    vk = np.ones((nlay, nrow, ncol)) + 1.0
    # save hk up one dir from model_ws
    fnames = []
    for i, h in enumerate(hk):
        fname = os.path.join("test_{0}.ref".format(i))
        fnames.append(fname)
        np.savetxt(fname, h, fmt="%15.6e", delimiter='')
        vk[i] = i + 1.

    lpf = flopy.modflow.ModflowLpf(ml, hk=fnames, vka=vk)
    ml.lpf.vka[0].format.binary = True
    ml.lpf.vka.cnstnt = 2.0
    ml.write_input()

    assert np.array_equal(ml.lpf.hk.array,hk)
    assert np.array_equal(ml.lpf.vka.array,vk * 2.0)

    ml1 = flopy.modflow.Modflow.load(ml.namefile,
                                     model_ws=ml.model_ws,
                                     verbose=True, forgive=False)
    print("testing load")
    assert ml1.load_fail == False
    assert np.array_equal(ml1.lpf.vka.array, vk * 2.0)
    assert np.array_equal(ml1.lpf.hk.array, hk)
    assert np.array_equal(ml1.lpf.vka.array, ml.lpf.vka.array)
    assert np.array_equal(ml1.lpf.hk.array, ml.lpf.hk.array)



    # more binary testing
    ml.lpf.vka[0]._array[0,0] *= 3.0
    ml.write_input()
    ml1 = flopy.modflow.Modflow.load(ml.namefile,
                                     model_ws=ml.model_ws,
                                     verbose=True, forgive=False)
    assert np.array_equal(ml.lpf.vka.array,ml1.lpf.vka.array)
    assert np.array_equal(ml.lpf.hk.array,ml1.lpf.hk.array)
dataset/ETHPy150Open modflowpy/flopy/autotest/t004_test.py/stress_util2d_for_joe_the_file_king/original