C++ scalargpuField::begin方法代码示例

void Foam::cyclicACMIFvPatch::makeWeights(scalargpuField& w) const
{
    if (coupled())
    {
        const cyclicACMIFvPatch& nbrPatch = neighbFvPatch();
        const fvPatch& nbrPatchNonOverlap = nonOverlapPatch();

        const scalargpuField deltas(nf() & coupledFvPatch::delta());

        const scalargpuField nbrDeltas
        (
            interpolate
            (
                nbrPatch.nf() & nbrPatch.coupledFvPatch::delta(),
                nbrPatchNonOverlap.nf() & nbrPatchNonOverlap.delta()
            )
        );

        thrust::transform
        (
            deltas.begin(),
            deltas.end(),
            nbrDeltas.begin(),
            w.begin(),
            cyclicACMIFvPatchMakeWeightsFunctor()
        );
    }
    else
    {
        // Behave as uncoupled patch
        fvPatch::makeWeights(w);
    }
}

void Foam::porosityModels::powerLaw::apply
(
    scalargpuField& Udiag,
    const scalargpuField& V,
    const RhoFieldType& rho,
    const vectorgpuField& U
) const
{
    const scalar C0 = C0_;
    const scalar C1m1b2 = (C1_ - 1.0)/2.0;

    forAll(cellZoneIDs_, zoneI)
    {
        const labelgpuList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]].getList();
  
        thrust::transform
        (
            thrust::make_permutation_iterator
            (
                Udiag.begin(),
                cells.begin()
            ),
            thrust::make_permutation_iterator
            (
                Udiag.begin(),
                cells.end()
            ),
            thrust::make_zip_iterator(thrust::make_tuple
            (
                thrust::make_permutation_iterator
                (
                V.begin(),cells.begin()
                ),
                thrust::make_permutation_iterator
                (
                rho.begin(),cells.begin()
                ),
                thrust::make_permutation_iterator
                (
                U.begin(),cells.begin()
                )
            )),
            thrust::make_permutation_iterator
            (
                Udiag.begin(),
                cells.begin()
            ),
            powerLawUdiagFunctor(C0,C1m1b2)
        );
    }
}

void Foam::diagonalPreconditioner::precondition
(
    scalargpuField& wA,
    const scalargpuField& rA,
    const direction
) const
{
    thrust::transform
    (
        rD.begin(),
        rD.end(),
        rA.begin(),
        wA.begin(),
        multiplyOperatorFunctor<scalar,scalar,scalar>()
    );
}

tmp<scalargpuField> nutURoughWallFunctionFvPatchScalarField::calcNut() const
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
            dimensionedInternalField().group()
        )
    );
    const scalargpuField& y = turbModel.y()[patchi];
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
    const tmp<scalargpuField> tnuw = turbModel.nu(patchi);
    const scalargpuField& nuw = tnuw();

    // The flow velocity at the adjacent cell centre
    const scalargpuField magUp(mag(Uw.patchInternalField() - Uw));

    tmp<scalargpuField> tyPlus = calcYPlus(magUp);
    scalargpuField& yPlus = tyPlus();

    tmp<scalargpuField> tnutw(new scalargpuField(patch().size(), 0.0));
    scalargpuField& nutw = tnutw();

    thrust::transform
    (
        y.begin(),
        y.end(),
        thrust::make_zip_iterator(thrust::make_tuple
        (
            yPlus.begin(),
            nuw.begin(),
            magUp.begin()
        )),
        nutw.begin(),
        nutURoughCalcNutFunctor(yPlusLam_)
    );

/*
    forAll(yPlus, facei)
    {
        if (yPlus[facei] > yPlusLam_)
        {
            const scalar Re = magUp[facei]*y[facei]/nuw[facei] + ROOTVSMALL;
            nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
        }
    }
*/
    return tnutw;
}

tmp<scalargpuField> nutURoughWallFunctionFvPatchScalarField::calcYPlus
(
    const scalargpuField& magUp
) const
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
            dimensionedInternalField().group()
        )
    );
    const scalargpuField& y = turbModel.y()[patchi];
    const tmp<scalargpuField> tnuw = turbModel.nu(patchi);
    const scalargpuField& nuw = tnuw();

    tmp<scalargpuField> tyPlus(new scalargpuField(patch().size(), 0.0));
    scalargpuField& yPlus = tyPlus();

    if (roughnessHeight_ > 0.0)
    {
        thrust::transform
        (
            y.begin(),
            y.end(),
            thrust::make_zip_iterator(thrust::make_tuple
            (
                nuw.begin(),
                magUp.begin()
            )),
            yPlus.begin(),
            nutURoughCalcYPlusRoughFunctor
            (
                yPlusLam_,
                kappa_,
                E_,
                roughnessHeight_,
                roughnessConstant_,
                roughnessFactor_
            )
        );
/*
        // Rough Walls
        const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
        static const scalar c_2 = 2.25/(90 - 2.25);
        static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
        static const scalar c_4 = c_3*log(2.25);

        //if (KsPlusBasedOnYPlus_)
        {
            // If KsPlus is based on YPlus the extra term added to the law
            // of the wall will depend on yPlus
            forAll(yPlus, facei)
            {
                const scalar magUpara = magUp[facei];
                const scalar Re = magUpara*y[facei]/nuw[facei];
                const scalar kappaRe = kappa_*Re;

                scalar yp = yPlusLam_;
                const scalar ryPlusLam = 1.0/yp;

                int iter = 0;
                scalar yPlusLast = 0.0;
                scalar dKsPlusdYPlus = roughnessHeight_/y[facei];

                // Additional tuning parameter - nominally = 1
                dKsPlusdYPlus *= roughnessFactor_;

                do
                {
                    yPlusLast = yp;

                    // The non-dimensional roughness height
                    scalar KsPlus = yp*dKsPlusdYPlus;

                    // The extra term in the law-of-the-wall
                    scalar G = 0.0;

                    scalar yPlusGPrime = 0.0;

                    if (KsPlus >= 90)
                    {
                        const scalar t_1 = 1 + roughnessConstant_*KsPlus;
                        G = log(t_1);
                        yPlusGPrime = roughnessConstant_*KsPlus/t_1;
                    }
                    else if (KsPlus > 2.25)
                    {
                        const scalar t_1 = c_1*KsPlus - c_2;
                        const scalar t_2 = c_3*log(KsPlus) - c_4;
                        const scalar sint_2 = sin(t_2);
                        const scalar logt_1 = log(t_1);
                        G = logt_1*sint_2;
                        yPlusGPrime =
                            (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
                    }

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

//.........这里部分代码省略.........
            }

            // Scale coarse-grid correction field
            // but not on the coarsest level because it evaluates to 1
            if
            (
                scaleCorrection_
                && (interpolateCorrection_ || leveli < coarsestLevel - 1)
            )
            {
                scale
                (
                    coarseCorrFields[leveli],
                    ACfRef,
                    matrixLevels_[leveli],
                    interfaceLevelsBouCoeffs_[leveli],
                    interfaceLevels_[leveli],
                    coarseSources[leveli],
                    cmpt
                );
            }

            // Only add the preSmoothedCoarseCorrField if pre-smoothing is
            // used
            if (nPreSweeps_)
            {
                coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
            }

            smoothers[leveli + 1].smooth
            (
                coarseCorrFields[leveli],
                coarseSources[leveli],
                cmpt,
                min
                (
                    nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
                    maxPostSweeps_
                )
            );
        }
    }

    // Prolong the finest level correction
    agglomeration_.prolongField
    (
        finestCorrection,
        coarseCorrFields[0],
        0
    );

    if (interpolateCorrection_)
    {
        interpolate
        (
            finestCorrection,
            Apsi,
            matrix_,
            interfaceBouCoeffs_,
            interfaces_,
            agglomeration_.restrictSortAddressing(0),
            agglomeration_.restrictTargetAddressing(0),
            agglomeration_.restrictTargetStartAddressing(0),
            coarseCorrFields[0],
            cmpt
        );
    }

    if (scaleCorrection_)
    {
        // Scale the finest level correction
        scale
        (
            finestCorrection,
            Apsi,
            matrix_,
            interfaceBouCoeffs_,
            interfaces_,
            finestResidual,
            cmpt
        );
    }

    thrust::transform
    (
        psi.begin(),
        psi.end(),
        finestCorrection.begin(),
        psi.begin(),
        thrust::plus<scalar>()
    );

    smoothers[0].smooth
    (
        psi,
        source,
        cmpt,
        nFinestSweeps_
    );
}

void Foam::GAMGSolver::scale
(
    scalargpuField& field,
    scalargpuField& Acf,
    const lduMatrix& A,
    const FieldField<gpuField, scalar>& interfaceLevelBouCoeffs,
    const lduInterfaceFieldPtrsList& interfaceLevel,
    const scalargpuField& source,
    const direction cmpt
) const
{
    A.Amul
    (
        Acf,
        field,
        interfaceLevelBouCoeffs,
        interfaceLevel,
        cmpt
    );

    scalar scalingFactorNum = 0.0;
    scalar scalingFactorDenom = 0.0;

    scalingFactorNum  = 
        thrust::reduce
        (
            thrust::make_transform_iterator
            (
                thrust::make_zip_iterator(thrust::make_tuple
                (
                    source.begin(),
                    field.begin()
                )),
                multiplyTupleFunctor()
            ),
            thrust::make_transform_iterator
            (
                thrust::make_zip_iterator(thrust::make_tuple
                (
                    source.end(),
                    field.end()
                )),
                multiplyTupleFunctor()
            ),
            0.0,
            thrust::plus<scalar>()
        );

    scalingFactorDenom  = 
        thrust::reduce
        (
            thrust::make_transform_iterator
            (
                thrust::make_zip_iterator(thrust::make_tuple
                (
                    Acf.begin(),
                    field.begin()
                )),
                multiplyTupleFunctor()
            ),
            thrust::make_transform_iterator
            (
                thrust::make_zip_iterator(thrust::make_tuple
                (
                    Acf.end(),
                    field.end()
                )),
                multiplyTupleFunctor()
            ),
            0.0,
            thrust::plus<scalar>()
        );

/*
    forAll(field, i)
    {
        scalingFactorNum += source[i]*field[i];
        scalingFactorDenom += Acf[i]*field[i];
    }
*/
    vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
    A.mesh().reduce(scalingVector, sumOp<vector2D>());

    scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);

    if (debug >= 2)
    {
        Pout<< sf << " ";
    }

    const scalargpuField& D = A.diag();

/*
    forAll(field, i)
    {
        field[i] = sf*field[i] + (source[i] - sf*Acf[i])/D[i];
    }
*/

    thrust::transform
//.........这里部分代码省略.........