def main(opt):
    from aero_kernels import dirichlet, dotPV, dotR, init_cg, res_calc, spMV, \
        update, updateP, updateUR
    try:
        with h5py.File(opt['mesh'], 'r') as f:
            # sets
            nodes = op2.Set.fromhdf5(f, 'nodes')
            bnodes = op2.Set.fromhdf5(f, 'bedges')
            cells = op2.Set.fromhdf5(f, 'cells')

            # maps
            pbnodes = op2.Map.fromhdf5(bnodes, nodes, f, 'pbedge')
            pcell = op2.Map.fromhdf5(cells, nodes, f, 'pcell')
            pvcell = op2.Map.fromhdf5(cells, nodes, f, 'pcell')

            # dats
            p_xm = op2.Dat.fromhdf5(nodes ** 2, f, 'p_x')
            p_phim = op2.Dat.fromhdf5(nodes, f, 'p_phim')
            p_resm = op2.Dat.fromhdf5(nodes, f, 'p_resm')
            p_K = op2.Dat.fromhdf5(cells ** 16, f, 'p_K')
            p_V = op2.Dat.fromhdf5(nodes, f, 'p_V')
            p_P = op2.Dat.fromhdf5(nodes, f, 'p_P')
            p_U = op2.Dat.fromhdf5(nodes, f, 'p_U')
    except __HOLE__:
        import sys
        print "Failed reading mesh: Could not read from %s\n" % opt['mesh']
        sys.exit(1)

    # Constants

    gam = 1.4
    gm1 = op2.Const(1, gam - 1.0, 'gm1', dtype=np.double)
    op2.Const(1, 1.0 / gm1.data, 'gm1i', dtype=np.double)
    op2.Const(2, [0.5, 0.5], 'wtg1', dtype=np.double)
    op2.Const(2, [0.211324865405187, 0.788675134594813], 'xi1',
              dtype=np.double)
    op2.Const(4, [0.788675134594813, 0.211324865405187,
                  0.211324865405187, 0.788675134594813],
              'Ng1', dtype=np.double)
    op2.Const(4, [-1, -1, 1, 1], 'Ng1_xi', dtype=np.double)
    op2.Const(4, [0.25] * 4, 'wtg2', dtype=np.double)
    op2.Const(16, [0.622008467928146, 0.166666666666667,
                   0.166666666666667, 0.044658198738520,
                   0.166666666666667, 0.622008467928146,
                   0.044658198738520, 0.166666666666667,
                   0.166666666666667, 0.044658198738520,
                   0.622008467928146, 0.166666666666667,
                   0.044658198738520, 0.166666666666667,
                   0.166666666666667, 0.622008467928146],
              'Ng2', dtype=np.double)
    op2.Const(32, [-0.788675134594813, 0.788675134594813,
                   -0.211324865405187, 0.211324865405187,
                   -0.788675134594813, 0.788675134594813,
                   -0.211324865405187, 0.211324865405187,
                   -0.211324865405187, 0.211324865405187,
                   -0.788675134594813, 0.788675134594813,
                   -0.211324865405187, 0.211324865405187,
                   -0.788675134594813, 0.788675134594813,
                   -0.788675134594813, -0.211324865405187,
                   0.788675134594813, 0.211324865405187,
                   -0.211324865405187, -0.788675134594813,
                   0.211324865405187, 0.788675134594813,
                   -0.788675134594813, -0.211324865405187,
                   0.788675134594813, 0.211324865405187,
                   -0.211324865405187, -0.788675134594813,
                   0.211324865405187, 0.788675134594813],
              'Ng2_xi', dtype=np.double)
    minf = op2.Const(1, 0.1, 'minf', dtype=np.double)
    op2.Const(1, minf.data ** 2, 'm2', dtype=np.double)
    op2.Const(1, 1, 'freq', dtype=np.double)
    op2.Const(1, 1, 'kappa', dtype=np.double)
    op2.Const(1, 0, 'nmode', dtype=np.double)
    op2.Const(1, 1.0, 'mfan', dtype=np.double)

    niter = 20

    for i in xrange(1, niter + 1):

        op2.par_loop(res_calc, cells,
                     p_xm(op2.READ, pvcell),
                     p_phim(op2.READ, pcell),
                     p_K(op2.WRITE),
                     p_resm(op2.INC, pcell))

        op2.par_loop(dirichlet, bnodes,
                     p_resm(op2.WRITE, pbnodes[0]))

        c1 = op2.Global(1, data=0.0, name='c1')
        c2 = op2.Global(1, data=0.0, name='c2')
        c3 = op2.Global(1, data=0.0, name='c3')
        # c1 = R' * R
        op2.par_loop(init_cg, nodes,
                     p_resm(op2.READ),
                     c1(op2.INC),
                     p_U(op2.WRITE),
                     p_V(op2.WRITE),
                     p_P(op2.WRITE))

        # Set stopping criteria
        res0 = sqrt(c1.data)
        res = res0
        res0 *= 0.1
        it = 0
        maxiter = 200

        while res > res0 and it < maxiter:

            # V = Stiffness * P
            op2.par_loop(spMV, cells,
                         p_V(op2.INC, pcell),
                         p_K(op2.READ),
                         p_P(op2.READ, pcell))

            op2.par_loop(dirichlet, bnodes,
                         p_V(op2.WRITE, pbnodes[0]))

            c2.data = 0.0

            # c2 = P' * V
            op2.par_loop(dotPV, nodes,
                         p_P(op2.READ),
                         p_V(op2.READ),
                         c2(op2.INC))

            alpha = op2.Global(1, data=c1.data / c2.data, name='alpha')

            # U = U + alpha * P
            # resm = resm - alpha * V
            op2.par_loop(updateUR, nodes,
                         p_U(op2.INC),
                         p_resm(op2.INC),
                         p_P(op2.READ),
                         p_V(op2.RW),
                         alpha(op2.READ))

            c3.data = 0.0
            # c3 = resm' * resm
            op2.par_loop(dotR, nodes,
                         p_resm(op2.READ),
                         c3(op2.INC))

            beta = op2.Global(1, data=c3.data / c1.data, name="beta")
            # P = beta * P + resm
            op2.par_loop(updateP, nodes,
                         p_resm(op2.READ),
                         p_P(op2.RW),
                         beta(op2.READ))

            c1.data = c3.data
            res = sqrt(c1.data)
            it += 1

        rms = op2.Global(1, data=0.0, name='rms')

        # phim = phim - Stiffness \ Load
        op2.par_loop(update, nodes,
                     p_phim(op2.RW),
                     p_resm(op2.WRITE),
                     p_U(op2.READ),
                     rms(op2.INC))

        print "rms = %10.5e iter: %d" % (sqrt(rms.data) / sqrt(nodes.size), it)
dataset/ETHPy150Open OP2/PyOP2/demo/aero.py/main