Python hysop.Field類代碼示例

def check_subset(discr, fileref, class_name):
    Cartesian.reset_counter()
    dim = len(discr.resolution)
    dom = Box(dimension=dim, length=ldom[:dim],
              origin=xdom[:dim])
    topo = dom.create_topology(discr)
    vref = Field(domain=topo.domain, name='vref')
    vd = vref.discretize(topo)
    iop = IOParams(working_dir + fileref)
    reader = HDF_Reader(variables={vref: topo}, io_params=iop, restart=0)
    reader.discretize()
    reader.setup()
    reader.apply()
    reader.finalize()
    #return topo, dom, vd
    rd = ldom[0] * 0.4
    subs = class_name(origin=xpos[:dim], radius=rd, parent=dom)
    assert subs.radius == rd
    assert (subs.origin[:dim] == xpos[:dim]).all()
    assert hasattr(subs.chi, '__call__')
    subs.discretize(topo)
    assert subs.ind.keys()[0] == topo
    assert len(subs.ind[topo][0]) == dom.dimension
    scal = Field(domain=topo.domain, is_vector=False, name='s0')
    sd = scal.discretize(topo)
    sd[0][subs.ind[topo][0]] = 2.
    icompute = topo.mesh.iCompute
    return np.allclose(sd[0][icompute], vd[0][icompute])

def test_analytical_field_9():
    box = Box()
    topo = box.create_topology(d3d)
    res2 = Discretization([65, 33, 65], [1, 1, 1])
    topo2 = box.create_topology(res2, dim=2)
    coords = topo.mesh.coords
    coords2 = topo2.mesh.coords
    caf = Field(box, name='f1', formula=func_vec_1, is_vector=True)
    ref = Field(box, name='f2', is_vector=True)
    refd = ref.discretize(topo)
    cafd = caf.discretize(topo)
    cafd2 = caf.discretize(topo2)
    refd2 = ref.discretize(topo2)
    ids = [0, ] * 3
    for i in xrange(3):
        ids[i] = id(cafd2.data[i])
        # init on topo2
    caf.initialize(topo=topo2)
    refd2.data = func_vec_1(refd2.data, *(coords2 + (0.,)))
    refd.data = func_vec_1(refd.data, *(coords + (0.,)))
    for i in xrange(caf.nb_components):
        assert allclose(cafd2[i], refd2.data[i])
        assert id(cafd2.data[i]) == ids[i]
        assert not allclose(cafd[i], refd.data[i])
    caf.initialize(topo=topo)
    for i in xrange(caf.nb_components):
        assert allclose(cafd[i], refd.data[i])

def test_analytical():
    """Test formula"""
    dom = Box()
    caf = Field(dom, name='f1', formula=formula, vectorize_formula=True)
    assert caf.value(0., 0., 0.) == (0., 1., 0.)
    assert caf.value(1., 1., 1.) == (1., 1., -1.)
    assert caf.value(1., 2., 3.) == (6., 1., -6.)

def test_cb_3d_fulldomain():
    """
    A cb defined on the whole domain.
    Integrals on upper surfaces must fail,
    because of periodic bc.
    """
    topo, cb = check_control_box(discr3D, orig, ldom)
    velo = Field(domain=topo.domain, is_vector=True, formula=v3d, name='velo')
    velo.discretize(topo)
    velo.initialize(topo=topo)
    i0 = cb.integrate_field_allc(velo, topo)
    vref = np.prod(cb.real_length[topo])
    assert (np.abs(i0 - vref) < 1e-6).all()
    nbc = velo.nb_components
    sref = npw.zeros(nbc)
    dirs = np.arange(nbc)
    for i in xrange(nbc):
        ilist = np.where(dirs != i)[0]
        sref = np.prod(cb.real_length[topo][ilist])
        isurf = cb.integrate_on_faces(velo, topo, [2 * i])
        assert np.abs(isurf - sref) < 1e-6
        res = False
        try:
            isurf = cb.integrate_on_faces(velo, topo, [2 * i + 1])
        except:
            res = True
        assert res

def test_randomize():
    box = Box()
    topo = box.create_topology(d3d)
    f2 = Field(box, name='target', is_vector=True)
    f2.randomize(topo)
    f2d = f2.discretize(topo)
    for d in xrange(f2.nb_components):
        assert not np.allclose(f2d[d], 0.)

def test_integ():
    topo, cub = check_subset(discr3D, xdef, ldef)
    velo = Field(domain=topo.domain, is_vector=True, formula=v3d, name='velo')
    vd = velo.discretize(topo)
    velo.initialize(topo=topo)
    i0 = cub.integrate_field_allc(velo, topo)
    i1 = cub.integrate_dfield_allc(vd)
    vref = np.prod(cub.real_length[topo])
    assert (np.abs(i0 - vref) < 1e-6).all()
    assert (np.abs(i1 - vref) < 1e-6).all()

def test_integ_2():
    topo, cub = check_subset(discr3D, xdef, ldef)
    velo = Field(domain=topo.domain, is_vector=True, formula=v3d, name='velo')
    vd = velo.discretize(topo)
    velo.initialize(topo=topo)
    vref = np.prod(cub.real_length[topo])
    for i in xrange(velo.nb_components):
        i0 = cub.integrate_field(velo, topo, component=i)
        assert np.abs(i0 - vref) < 1e-6
        i1 = cub.integrate_dfield(vd, component=i)
        assert np.abs(i1 - vref) < 1e-6

def test_div_rho_v_2d():
    # Reference field
    ref = Field(domain=topo2.domain, formula=analyticalDivWV2D,
                name='Analytical', is_vector=False)
    rd = ref.discretize(topo2)
    ref.initialize(topo=topo2)
    op, result = build_op(diffop.DivRhoV, 1, topo2, work2)
    tol = np.max(topo2.mesh.space_step ** 4)
    result = op(v2d.data, w2d.data[0:1], result)
    # Numerical VS analytical
    assert np.allclose(rd[0][ic2], result[0][ic2], atol=tol)

def test_laplacian():
    ref = Field(domain=topo3.domain, formula=laplacian_func,
                name='Analytical', is_vector=True)
    rd = ref.discretize(topo3)
    ref.initialize(topo=topo3)

    op, result = build_op(diffop.Laplacian, 3, topo3, work3)
    result = op(w3d.data, result)
    tol = np.max(topo3.mesh.space_step ** 2)
    for i in xrange(3):
        assert np.allclose(rd.data[i][ic3], result[i][ic3], rtol=tol)

def test_div_advection():
    # Reference scalar field
    ref = Field(domain=topo3.domain, formula=analyticalDivAdvection,
                name='Analytical')
    rd = ref.discretize(topo3)
    ref.initialize(topo=topo3)
    op, result = build_op(diffop.DivAdvection, 1, topo3, work3)
    result = op(v3d.data, result)

    # Numerical VS analytical
    errx = (topo3.domain.length[0] / (Nx - 1)) ** 4
    assert np.allclose(rd[0][ic3], result[0][ic3], rtol=errx)

def test_div_rho_v():
    # Reference field
    ref = Field(domain=topo3.domain, formula=analyticalDivWV,
                name='Analytical', is_vector=True)
    rd = ref.discretize(topo3)
    ref.initialize(topo=topo3)
    op, result = build_op(diffop.DivRhoV, 3, topo3, work3)
    tol = np.max(topo3.mesh.space_step ** 4)
    for i in xrange(3):
        result[i:i + 1] = op(v3d.data, w3d.data[i:i + 1], result[i:i + 1])
        # Numerical VS analytical
        assert np.allclose(rd[i][ic3], result[i][ic3], atol=tol)

def test_2d_subbox():
    """
    subset in 2D domain
    """
    topo, cub = check_subset(discr2D, xdef[:2], ldef[:2])
    velo = Field(domain=topo.domain, is_vector=True, formula=v2d, name='velo')
    vd = velo.discretize(topo)
    velo.initialize(topo=topo)
    i0 = cub.integrate_field_allc(velo, topo)
    vref = np.prod(cub.real_length[topo])
    assert (np.abs(i0 - vref) < 1e-6).all()
    i1 = cub.integrate_dfield_allc(vd)
    assert (np.abs(i1 - vref) < 1e-6).all()

def test_density():
    """
    Todo : write proper tests.
    Here we just check if discr/setup/apply process goes well.
    """
    box = pp.Box(length=[2., 1., 0.9], origin=[0.0, -1., -0.43])
    density = Field(domain=box, name='density')
    viscosity = Field(domain=box, formula=computeVisco, name='visco')
    op = DensityVisco(density=density, viscosity=viscosity, discretization=d3d)
    op.discretize()
    op.setup()
    topo = op.variables[viscosity]
    viscosity.initialize(topo=topo)

def test_2d_line():
    """
    subset == line in 2D domain
    Test integration on this subset.
    """
    topo, cub = check_subset(discr2D, xdef[:2], lr=[0., 1.0])
    velo = Field(domain=topo.domain, is_vector=True, formula=v2d, name='velo')
    vd = velo.discretize(topo)
    velo.initialize(topo=topo)
    i0 = cub.integrate_field_allc(velo, topo)
    vref = cub.real_length[topo][1]
    assert (np.abs(i0 - vref) < 1e-6).all()
    i1 = cub.integrate_dfield_allc(vd)
    assert (np.abs(i1 - vref) < 1e-6).all()

def test_subbox_2d():
    """
    subset == plane in 3D domain
    Test integration on this subset.
    """
    topo, cub = check_subset(discr3D, xdef, [0., 0.2, 0.8])
    velo = Field(domain=topo.domain, is_vector=True, formula=v3d, name='velo')
    vd = velo.discretize(topo)
    velo.initialize(topo=topo)
    i0 = cub.integrate_field_allc(velo, topo)
    vref = np.prod(cub.real_length[topo][1:])
    assert (np.abs(i0 - vref) < 1e-6).all()
    i1 = cub.integrate_dfield_allc(vd)
    assert (np.abs(i1 - vref) < 1e-6).all()