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;
}

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)
        );
    }
}

Foam::tmp<Foam::scalargpuField> Foam::GAMGInterface::agglomerateCoeffs
(
    const scalargpuField& fineCoeffs
) const
{
    tmp<scalargpuField> tcoarseCoeffs(new scalargpuField(size(), 0.0));
    scalargpuField& coarseCoeffs = tcoarseCoeffs();

    if (fineCoeffs.size() != faceRestrictAddressing_.size())
    {
        FatalErrorIn
        (
            "GAMGInterface::agglomerateCoeffs(const scalarField&) const"
        )   << "Size of coefficients " << fineCoeffs.size()
            << " does not correspond to the size of the restriction "
            << faceRestrictAddressing_.size()
            << abort(FatalError);
    }
    if (debug && max(faceRestrictAddressing_) > size())
    {
        FatalErrorIn
        (
            "GAMGInterface::agglomerateCoeffs(const scalargpuField&) const"
        )   << "Face restrict addressing addresses outside of coarse interface"
            << " size. Max addressing:" << max(faceRestrictAddressing_)
            << " coarse size:" << size()
            << abort(FatalError);
    }

    thrust::transform
    (
        faceRestrictTargetStartAddressing_.begin(),
        faceRestrictTargetStartAddressing_.end()-1,
        faceRestrictTargetStartAddressing_.begin()+1,
        thrust::make_permutation_iterator
        (
            coarseCoeffs.begin(),
            faceRestrictTargetAddressing_.begin()
        ),
        GAMGInterfaceAgglomerateCoeffs
        (
            fineCoeffs.data(),
            faceRestrictSortAddressing_.data()
        )
    );

    return tcoarseCoeffs;
}

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>()
    );
}

void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
    gpuField<Type>& result,
    const gpuField<Type>& psiInternal,
    const scalargpuField& coeffs,
    const Pstream::commsTypes
) const
{
    // note: only applying coupled contribution

    const labelgpuList& nbrFaceCellsCoupled =
        cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().getFaceCells();

    gpuField<Type> pnf(psiInternal, nbrFaceCellsCoupled);

    // Transform according to the transformation tensors
    transformCoupleField(pnf);

    pnf = cyclicACMIPatch_.interpolate(pnf);

    matrixPatchOperation
    (
        this->patch().index(),
        result,
        this->patch().boundaryMesh().mesh().lduAddr(),
        matrixInterfaceFunctor<Type>
        (
            coeffs.data(),
            pnf.data()
        )
    );
}

void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
(
    gpuField<Type>& result,
    const gpuField<Type>& psiInternal,
    const scalargpuField& coeffs,
    const Pstream::commsTypes
) const
{
    const labelgpuList& nbrFaceCells =
        this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().getFaceCells();

    gpuField<Type> pnf(psiInternal, nbrFaceCells);

    if (this->cyclicAMIPatch().applyLowWeightCorrection())
    {
        pnf =
            this->cyclicAMIPatch().interpolate
            (
                pnf,
                this->patchInternalField()()
            );

    }
    else
    {
        pnf = this->cyclicAMIPatch().interpolate(pnf);
    }

    // only apply jump to original field
    if (&psiInternal == &this->internalField())
    {
        gpuField<Type> jf(this->jump());
        if (!this->cyclicAMIPatch().owner())
        {
            jf *= -1.0;
        }

        pnf -= jf;
    }

    // Transform according to the transformation tensors
    this->transformCoupleField(pnf);

    // Multiply the field by coefficients and add into the result
    matrixPatchOperation
    (
        this->patch().index(),
        result,
        this->patch().boundaryMesh().mesh().lduAddr(),
        matrixInterfaceFunctor<Type>
        (
            coeffs.data(),
            pnf.data()
        )
    );
}

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);
    }
}

Foam::solverPerformance Foam::smoothSolver::solve
(
    scalargpuField& psi,
    const scalargpuField& source,
    const direction cmpt
) const
{
    // Setup class containing solver performance data
    solverPerformance solverPerf(typeName, fieldName_);

    // If the nSweeps_ is negative do a fixed number of sweeps
    if (nSweeps_ < 0)
    {
        autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
        (
            fieldName_,
            matrix_,
            interfaceBouCoeffs_,
            interfaceIntCoeffs_,
            interfaces_,
            controlDict_
        );

        smootherPtr->smooth
        (
            psi,
            source,
            cmpt,
            -nSweeps_
        );

        solverPerf.nIterations() -= nSweeps_;
    }
    else
    {
        scalar normFactor = 0;

        {
            scalargpuField Apsi(psi.size());
            scalargpuField temp(psi.size());

            // Calculate A.psi
            matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);

            // Calculate normalisation factor
            normFactor = this->normFactor(psi, source, Apsi, temp);

            // Calculate residual magnitude
            solverPerf.initialResidual() = gSumMag
            (
                (source - Apsi)(),
                matrix().mesh().comm()
            )/normFactor;
            solverPerf.finalResidual() = solverPerf.initialResidual();
        }

        if (lduMatrix::debug >= 2)
        {
            Info.masterStream(matrix().mesh().comm())
                << "   Normalisation factor = " << normFactor << endl;
        }


        // Check convergence, solve if not converged
        if
        (
            minIter_ > 0
         || !solverPerf.checkConvergence(tolerance_, relTol_)
        )
        {
            autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
            (
                fieldName_,
                matrix_,
                interfaceBouCoeffs_,
                interfaceIntCoeffs_,
                interfaces_,
                controlDict_
            );

            // Smoothing loop
            do
            {
                smootherPtr->smooth
                (
                    psi,
                    source,
                    cmpt,
                    nSweeps_
                );

                // Calculate the residual to check convergence
                solverPerf.finalResidual() = gSumMag
                (
                    matrix_.residual
                    (
                        psi,
                        source,
                        interfaceBouCoeffs_,
                        interfaces_,
//.........这里部分代码省略.........

void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
(
    gpuField<Type>& result,
    const gpuField<Type>&,
    const scalargpuField& coeffs,
    const Pstream::commsTypes commsType
) const
{
    if (this->updatedMatrix())
    {
        return;
    }

    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
    {
        // Fast path.
        if
        (
            outstandingRecvRequest_ >= 0
         && outstandingRecvRequest_ < Pstream::nRequests()
        )
        {
            UPstream::waitRequest(outstandingRecvRequest_);
        }
        // Recv finished so assume sending finished as well.
        outstandingSendRequest_ = -1;
        outstandingRecvRequest_ = -1;

        // Consume straight from receiveBuf_

        // Transform according to the transformation tensor
        
        gpuReceiveBuf_ = receiveBuf_;
        transformCoupleField(gpuReceiveBuf_);

        // Multiply the field by coefficients and add into the result
        matrixPatchOperation
        (
            this->patch().index(),
            result,
            this->patch().boundaryMesh().mesh().lduAddr(),
            processorFvPatchFunctor<Type>
            (
                coeffs.data(),
                gpuReceiveBuf_.data()
            )
        );
    }
    else
    {
        gpuField<Type> pnf
        (
            procPatch_.compressedReceive<Type>(commsType, this->size())()
        );

        // Transform according to the transformation tensor
        transformCoupleField(pnf);

        // Multiply the field by coefficients and add into the result
        matrixPatchOperation
        (
            this->patch().index(),
            result,
            this->patch().boundaryMesh().mesh().lduAddr(),
            processorFvPatchFunctor<Type>
            (
                coeffs.data(),
                pnf.data()
            )
        );
    }

    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
}

Foam::solverPerformance Foam::GAMGSolver::solve
(
    scalargpuField& psi,
    const scalargpuField& source,
    const direction cmpt
) const
{
    // Setup class containing solver performance data
    solverPerformance solverPerf(typeName, fieldName_);

    // Calculate A.psi used to calculate the initial residual
    scalargpuField Apsi(psi.size());
    matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);

    // Create the storage for the finestCorrection which may be used as a
    // temporary in normFactor
    scalargpuField finestCorrection(psi.size());

    // Calculate normalisation factor
    scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);

    if (debug >= 2)
    {
        Pout<< "   Normalisation factor = " << normFactor << endl;
    }

    // Calculate initial finest-grid residual field
    scalargpuField finestResidual(source - Apsi);

    // Calculate normalised residual for convergence test
    solverPerf.initialResidual() = gSumMag
                                   (
                                       finestResidual,
                                       matrix().mesh().comm()
                                   )/normFactor;
    solverPerf.finalResidual() = solverPerf.initialResidual();


    // Check convergence, solve if not converged
    if
    (
        minIter_ > 0
        || !solverPerf.checkConvergence(tolerance_, relTol_)
    )
    {
        // Create coarse grid correction fields
        PtrList<scalargpuField> coarseCorrFields;

        // Create coarse grid sources
        PtrList<scalargpuField> coarseSources;

        // Create the smoothers for all levels
        PtrList<lduMatrix::smoother> smoothers;

        // Scratch fields if processor-agglomerated coarse level meshes
        // are bigger than original. Usually not needed
        scalargpuField scratch1;
        scalargpuField scratch2;

        // Initialise the above data structures
        initVcycle
        (
            coarseCorrFields,
            coarseSources,
            smoothers,
            scratch1,
            scratch2
        );

        do
        {
            Vcycle
            (
                smoothers,
                psi,
                source,
                Apsi,
                finestCorrection,
                finestResidual,

                (scratch1.size() ? scratch1 : Apsi),
                (scratch2.size() ? scratch2 : finestCorrection),

                coarseCorrFields,
                coarseSources,
                cmpt
            );

            // Calculate finest level residual field
            matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
            finestResidual = source;
            finestResidual -= Apsi;

            solverPerf.finalResidual() = gSumMag
                                         (
                                             finestResidual,
                                             matrix().mesh().comm()
                                         )/normFactor;

            if (debug >= 2)
//.........这里部分代码省略.........

void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
(
    scalargpuField& result,
    const scalargpuField&,
    const scalargpuField& coeffs,
    const direction cmpt,
    const Pstream::commsTypes commsType
) const
{
    if (updatedMatrix())
    {
        return;
    }

    label oldWarn = UPstream::warnComm;
    UPstream::warnComm = comm();

    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
    {
        // Fast path.
        if
        (
            outstandingRecvRequest_ >= 0
         && outstandingRecvRequest_ < Pstream::nRequests()
        )
        {
            UPstream::waitRequest(outstandingRecvRequest_);
        }
        // Recv finished so assume sending finished as well.
        outstandingSendRequest_ = -1;
        outstandingRecvRequest_ = -1;

        // Consume straight from scalarReceiveBuf_

        if( ! Pstream::gpuDirectTransfer)
        {
            scalargpuReceiveBuf_ = scalarReceiveBuf_;
        }

        // Transform according to the transformation tensor
        transformCoupleField(scalargpuReceiveBuf_, cmpt);

        // Multiply the field by coefficients and add into the result  
        GAMGUpdateInterfaceMatrix
        (
            result,
            coeffs,
            scalargpuReceiveBuf_,
            procInterface_
        );      
    }
    else
    {
        scalargpuReceiveBuf_.setSize(coeffs.size());
        procInterface_.compressedReceive<scalar>(commsType, scalargpuReceiveBuf_);

        transformCoupleField(scalargpuReceiveBuf_, cmpt);

        GAMGUpdateInterfaceMatrix
        (
            result,
            coeffs,
            scalargpuReceiveBuf_,
            procInterface_
        ); 
    }

    const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = true;

    UPstream::warnComm = oldWarn;
}

void Foam::GAMGSolver::Vcycle
(
    const PtrList<lduMatrix::smoother>& smoothers,
    scalargpuField& psi,
    const scalargpuField& source,
    scalargpuField& Apsi,
    scalargpuField& finestCorrection,
    scalargpuField& finestResidual,

    scalargpuField& scratch1,
    scalargpuField& scratch2,

    PtrList<scalargpuField>& coarseCorrFields,
    PtrList<scalargpuField>& coarseSources,
    const direction cmpt
) const
{
    //debug = 2;

    const label coarsestLevel = matrixLevels_.size() - 1;

    // Restrict finest grid residual for the next level up.
    agglomeration_.restrictField(coarseSources[0], finestResidual, 0);

    if (debug >= 2 && nPreSweeps_)
    {
        Pout<< "Pre-smoothing scaling factors: ";
    }


    // Residual restriction (going to coarser levels)
    for (label leveli = 0; leveli < coarsestLevel; leveli++)
    {
        if (coarseSources.set(leveli + 1))
        {
            // If the optional pre-smoothing sweeps are selected
            // smooth the coarse-grid field for the restriced source
            if (nPreSweeps_)
            {
                coarseCorrFields[leveli] = 0.0;

                smoothers[leveli + 1].smooth
                (
                    coarseCorrFields[leveli],
                    coarseSources[leveli],
                    cmpt,
                    min
                    (
                        nPreSweeps_ +  preSweepsLevelMultiplier_*leveli,
                        maxPreSweeps_
                    )
                );

                scalargpuField ACf
                (
                    const_cast<const scalargpuField&>(scratch1),
                    coarseCorrFields[leveli].size()
                );

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

                // Correct the residual with the new solution
                matrixLevels_[leveli].Amul
                (
                    ACf,
                    coarseCorrFields[leveli],
                    interfaceLevelsBouCoeffs_[leveli],
                    interfaceLevels_[leveli],
                    cmpt
                );

                coarseSources[leveli] -= ACf;
            }

            // Residual is equal to source
            agglomeration_.restrictField
            (
                coarseSources[leveli + 1],
                coarseSources[leveli],
                leveli + 1
            );
        }
    }

    if (debug >= 2 && nPreSweeps_)
    {
//.........这里部分代码省略.........

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
//.........这里部分代码省略.........

void processorFvPatchField<scalar>::updateInterfaceMatrix
(
    scalargpuField& result,
    const scalargpuField&,
    const scalargpuField& coeffs,
    const direction,
    const Pstream::commsTypes commsType
) const
{
    if (this->updatedMatrix())
    {
        return;
    }

    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
    {
        // Fast path.
        if
        (
            outstandingRecvRequest_ >= 0
         && outstandingRecvRequest_ < Pstream::nRequests()
        )
        {
            UPstream::waitRequest(outstandingRecvRequest_);
        }
        // Recv finished so assume sending finished as well.
        outstandingSendRequest_ = -1;
        outstandingRecvRequest_ = -1;

        scalargpuReceiveBuf_ = scalarReceiveBuf_;
        // Consume straight from scalarReceiveBuf_
/*
        forAll(faceCells, elemI)
        {
            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
        }
*/
        matrixPatchOperation(this->patch().index(),
                             result,
                             this->patch().boundaryMesh().mesh().lduAddr(),
                             processorFvPatchScalarFunctor(coeffs.data(),scalargpuReceiveBuf_.data()));
    }
    else
    {
        scalargpuField pnf
        (
            procPatch_.compressedReceive<scalar>(commsType, this->size())()
        );
/*
        forAll(faceCells, elemI)
        {
            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
        }
*/
        matrixPatchOperation(this->patch().index(),
                             result,
                             this->patch().boundaryMesh().mesh().lduAddr(),
                             processorFvPatchScalarFunctor(coeffs.data(),pnf.data()));
    }

    const_cast<processorFvPatchField<scalar>&>(*this).updatedMatrix() = true;
}

void epsilonLowReWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbulence,
    const gpuList<scalar>& cornerWeights,
    const fvPatch& patch,
    scalargpuField& G,
    scalargpuField& epsilon
)
{
    const label patchi = patch.index();

    const scalargpuField& y = turbulence.y()[patchi];

    const scalar Cmu25 = pow025(Cmu_);
    const scalar Cmu75 = pow(Cmu_, 0.75);

    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();

    const tmp<scalargpuField> tnuw = turbulence.nu(patchi);
    const scalargpuField& nuw = tnuw();

    const tmp<scalargpuField> tnutw = turbulence.nut(patchi);
    const scalargpuField& nutw = tnutw();

    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchi];

    const scalargpuField magGradUw(mag(Uw.snGrad()));

    matrixPatchOperation
    (
        patchi,
        epsilon,
        patch.boundaryMesh().mesh().lduAddr(),
        EpsilonLowReCalculateEpsilonFunctor
        (
            yPlusLam_,
            Cmu25,
            Cmu75,
            kappa_,
            cornerWeights.data(),
            y.data(),
            k.getField().data(),
            nuw.data()
        )
    );
	
    matrixPatchOperation
    (
        patchi,
        G,
        patch.boundaryMesh().mesh().lduAddr(),
        EpsilonLowReCalculateGFunctor
        (
            Cmu25,
            kappa_,
            cornerWeights.data(),
            y.data(),
            k.getField().data(),
            nuw.data(),
            nutw.data(),
            magGradUw.data()
        )
    );
}

void Foam::fvMesh::storeOldVol(const scalargpuField& V)
{
    if (curTimeIndex_ < time().timeIndex())
    {
        if (debug)
        {
            Info<< "fvMesh::storeOldVol(const scalarField&) :"
                << " Storing old time volumes since from time " << curTimeIndex_
                << " and time now " << time().timeIndex()
                << " V:" << V.size()
                << endl;
        }


        if (V00Ptr_ && V0Ptr_)
        {
            // Copy V0 into V00 storage
            *V00Ptr_ = *V0Ptr_;
        }

        if (V0Ptr_)
        {
            // Copy V into V0 storage
            V0Ptr_->getField().scalargpuField::operator=(V);
        }
        else
        {
            // Allocate V0 storage, fill with V
            V0Ptr_ = new DimensionedField<scalar, volMesh>
            (
                IOobject
                (
                    "V0",
                    time().timeName(),
                    *this,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                *this,
                dimVolume
            );
            scalargpuField& V0 = (*V0Ptr_).getField();
            // Note: V0 now sized with current mesh, not with (potentially
            //       different size) V.
            V0.setSize(V.size());
            V0 = V;
        }

        curTimeIndex_ = time().timeIndex();

        if (debug)
        {
            Info<< "fvMesh::storeOldVol() :"
                << " Stored old time volumes V0:" << V0Ptr_->size()
                << endl;
            if (V00Ptr_)
            {
                Info<< "fvMesh::storeOldVol() :"
                    << " Stored oldold time volumes V00:" << V00Ptr_->size()
                    << endl;
            }
        }
    }
}

void Foam::GAMGSolver::initVcycle
(
    PtrList<scalargpuField>& coarseCorrFields,
    PtrList<scalargpuField>& coarseSources,
    PtrList<lduMatrix::smoother>& smoothers,
    scalargpuField& scratch1,
    scalargpuField& scratch2
) const
{
    label maxSize = matrix_.diag().size();

    coarseCorrFields.setSize(matrixLevels_.size());
    coarseSources.setSize(matrixLevels_.size());
    smoothers.setSize(matrixLevels_.size() + 1);

    // Create the smoother for the finest level
    smoothers.set
    (
        0,
        lduMatrix::smoother::New
        (
            fieldName_,
            matrix_,
            interfaceBouCoeffs_,
            interfaceIntCoeffs_,
            interfaces_,
            controlDict_
        )
    );

    forAll(matrixLevels_, leveli)
    {
        if (agglomeration_.nCells(leveli) >= 0)
        {
            label nCoarseCells = agglomeration_.nCells(leveli);
            coarseSources.set(leveli,GAMGSolverCache::source(leveli,nCoarseCells));
            //coarseSources.set(leveli, new scalargpuField(nCoarseCells));
        }

        if (matrixLevels_.set(leveli))
        {
            const lduMatrix& mat = matrixLevels_[leveli];

            label nCoarseCells = mat.diag().size();

            maxSize = max(maxSize, nCoarseCells);

            //coarseCorrFields.set(leveli, new scalargpuField(nCoarseCells));
            coarseCorrFields.set(leveli,GAMGSolverCache::corr(leveli,nCoarseCells));

            smoothers.set
            (
                leveli + 1,
                lduMatrix::smoother::New
                (
                    fieldName_,
                    matrixLevels_[leveli],
                    interfaceLevelsBouCoeffs_[leveli],
                    interfaceLevelsIntCoeffs_[leveli],
                    interfaceLevels_[leveli],
                    controlDict_
                )
            );
        }
    }

    if (maxSize > matrix_.diag().size())
    {
        // Allocate some scratch storage
        scratch1.setSize(maxSize);
        scratch2.setSize(maxSize);
    }
}