C++ scalarField::size方法代码示例

Foam::scalar checkUniformTimeStep(const scalarField& t)
{
    // check that a uniform time step has been applied
    scalar deltaT = -1.0;
    if (t.size() > 1)
    {
        for (label i = 1; i < t.size(); i++)
        {
            scalar dT = t[i] - t[i-1];
            if (deltaT < 0)
            {
                deltaT = dT;
            }

            if (mag(deltaT - dT) > SMALL)
            {
                FatalErrorIn("checkUniformTimeStep(const scalarField&)")
                    << "Unable to process data with a variable time step"
                    << exit(FatalError);
            }
        }
    }
    else
    {
        FatalErrorIn("checkUniformTimeStep(const scalarField&)")
            << "Unable to create FFT with a single value"
            << exit(FatalError);
    }

    return deltaT;
}

Foam::tmp<Foam::scalarField> Foam::GAMGInterface::agglomerateCoeffs
(
    const scalarField& fineCoeffs
) const
{
    tmp<scalarField> tcoarseCoeffs(new scalarField(size(), Zero));
    scalarField& coarseCoeffs = tcoarseCoeffs.ref();

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

    forAll(faceRestrictAddressing_, ffi)
    {
        coarseCoeffs[faceRestrictAddressing_[ffi]] += fineCoeffs[ffi];
    }

    return tcoarseCoeffs;
}

void Foam::dx<Type>::writeDXData
(
    const pointField& points,
    const scalarField& values,
    Ostream& os
) const
{
    // Write data
    os  << "object 3 class array type float rank 0 items "
        << values.size()
        << " data follows" << nl;

    forAll(values, elemI)
    {
        os << float(values[elemI]) << nl;
    }

    if (values.size() == points.size())
    {
        os  << nl << "attribute \"dep\" string \"positions\""
            << nl << nl;
    }
    else
    {
        os  << nl << "attribute \"dep\" string \"connections\""
            << nl << nl;
    }
}

int main(int argc, char *argv[])
{
    argList::noParallel();
    #include "addDictOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createFields.H"

    Info<< "Reading data file" << endl;
    CSV<scalar> pData("pressure", dict, "Data");

    // time history data
    const scalarField t(pData.x());

    // pressure data
    const scalarField p(pData.y());

    if (t.size() < N)
    {
        FatalErrorIn(args.executable())
            << "Block size N = " << N
            << " is larger than number of data = " << t.size()
            << exit(FatalError);
    }

    Info<< "    read " << t.size() << " values" << nl << endl;


    Info<< "Creating noise FFT" << endl;
    noiseFFT nfft(checkUniformTimeStep(t), p);

    nfft -= pRef;

    fileName baseFileName(pData.fName().lessExt());

    graph Pf(nfft.RMSmeanPf(N, min(nfft.size()/N, nw)));
    Info<< "    Creating graph for " << Pf.title() << endl;
    Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat);

    graph Lf(nfft.Lf(Pf));
    Info<< "    Creating graph for " << Lf.title() << endl;
    Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat);

    graph Ldelta(nfft.Ldelta(Lf, f1, fU));
    Info<< "    Creating graph for " << Ldelta.title() << endl;
    Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat);

    graph Pdelta(nfft.Pdelta(Pf, f1, fU));
    Info<< "    Creating graph for " << Pdelta.title() << endl;
    Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat);

    Info<< nl << "End\n" << endl;

    return 0;
}

Foam::scalar Foam::lduMatrix::solver::normFactor
(
    const scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
    scalarField wA(x.size());
    scalarField tmpField(x.size());

    matrix_.Amul(wA, x, interfaceBouCoeffs_, interfaces_, cmpt);

    return normFactor(x, b, wA, tmpField);
}

void Foam::GAMGPreconditioner::precondition
(
    scalarField& wA,
    const scalarField& rA,
    const direction cmpt
) const
{
    wA = 0.0;
    scalarField AwA(wA.size());
    scalarField finestCorrection(wA.size());
    scalarField finestResidual(rA);

    // Create coarse grid correction fields
    PtrList<scalarField> coarseCorrX;

    // Create coarse grid sources
    PtrList<scalarField> coarseB;

    // Create the smoothers for all levels
    PtrList<lduSmoother> smoothers;

    // Initialise the above data structures
    GAMG_.initVcycle(coarseCorrX, coarseB, smoothers);

    for (label cycle=0; cycle<nVcycles_; cycle++)
    {
        GAMG_.Vcycle
        (
            smoothers,
            wA,
            rA,
            AwA,
            finestCorrection,
            finestResidual,
            coarseCorrX,
            coarseB,
            cmpt
        );

        if (cycle < nVcycles_ - 1)
        {
            // Calculate finest level residual field
            matrix_.Amul(AwA, wA, coupleBouCoeffs_, interfaces_, cmpt);
            finestResidual = rA;
            finestResidual -= AwA;
        }
    }
}

Foam::scalar Foam::lduMatrix::solver::normFactor
(
    const scalarField& x,
    const scalarField& b,
    const scalarField& Ax,
    scalarField& tmpField,
    const direction cmpt
) const
{
    // Calculate A dot reference value of x
//     matrix_.sumA(tmpField, coupleBouCoeffs_, interfaces_);
//     tmpField *= gAverage(x);

    // Calculate normalisation factor using full multiplication
    // with mean value.  HJ, 5/Nov/2007
    scalar xRef = gAverage(x);
    matrix_.Amul
    (
        tmpField,
        scalarField(x.size(), xRef),
        coupleBouCoeffs_,
        interfaces_,
        cmpt
    );

    return gSum(mag(Ax - tmpField) + mag(b - tmpField)) + matrix_.small_;

    // At convergence this simpler method is equivalent to the above
    // return 2*gSumMag(b) + matrix_.small_;
}

Foam::lduSolverPerformance Foam::gmresSolver::solve
(
    scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
    // Prepare solver performance
    lduSolverPerformance solverPerf(typeName, fieldName());

    scalarField wA(x.size());
    scalarField rA(x.size());

    // Calculate initial residual
    matrix_.Amul(wA, x, coupleBouCoeffs_, interfaces_, cmpt);

    // Use rA as scratch space when calculating the normalisation factor
    scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);

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

    // Calculate residual
    forAll (rA, i)
    {
        rA[i] = b[i] - wA[i];
    }

void Foam::DILUPreconditioner::calcReciprocalD
(
    scalarField& rD,
    const lduMatrix& matrix
)
{
    scalar* __restrict__ rDPtr = rD.begin();

    const label* const __restrict__ uPtr =
        matrix.lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr =
        matrix.lduAddr().lowerAddr().begin();

    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
    const scalar* const __restrict__ lowerPtr = matrix.lower().begin();

    register label nFaces = matrix.upper().size();
    for (register label face=0; face<nFaces; face++)
    {
        rDPtr[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rDPtr[lPtr[face]];
    }


    // Calculate the reciprocal of the preconditioned diagonal
    register label nCells = rD.size();

    for (register label cell=0; cell<nCells; cell++)
    {
        rDPtr[cell] = 1.0/rDPtr[cell];
    }
}

void Foam::DILUPreconditioner::calcReciprocalD
(
    scalarField& rD,
    const lduMatrix& matrix
)
{
    scalar* __restrict__ rDPtr = rD.begin();

    const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();

    const scalar* const __restrict__ upperPtr = matrix.upper().begin();
    const scalar* const __restrict__ lowerPtr = matrix.lower().begin();

    label nFaces = matrix.upper().size();
    for (label face=0; face<nFaces; face++)
    {
        Pout<< "Adapting diagonal for cell:" << uPtr[face]
            << " contributions from cell " << lPtr[face]
            << " from " << rDPtr[uPtr[face]];

        rDPtr[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rDPtr[lPtr[face]];

        Pout<< " to " << rDPtr[uPtr[face]] << endl;
    }


    // Calculate the reciprocal of the preconditioned diagonal
    label nCells = rD.size();

    for (label cell=0; cell<nCells; cell++)
    {
        rDPtr[cell] = 1.0/rDPtr[cell];
    }
}

void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
(
    const scalarField&,
    scalarField& result,
    const lduMatrix&,
    const scalarField& coeffs,
    const direction cmpt,
    const Pstream::commsTypes commsType,
    const bool switchToLhs
) const
{
    scalarField pnf
    (
        procInterface_.compressedReceive<scalar>(commsType, coeffs.size())
    );
    transformCoupleField(pnf, cmpt);

    const unallocLabelList& faceCells = procInterface_.faceCells();

    if (switchToLhs)
    {
        forAll(faceCells, elemI)
        {
            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
        }
    }

Foam::lduSolverPerformance Foam::cufflink_DiagPCG_Parallel::solve
(
    scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
   
	if (!Pstream::parRun()) {
		FatalErrorIn("cufflink:\nMulti-GPU Solver cannot be run serially.  Choose the serial version of the solver")<< "Fatal Error"<< abort(FatalError);
	}

	OFSolverPerformance OFSP;//container for solver performance

	#include "../CFL_Headers/getGPUStorage.H"

	// --- Setup class containing solver performance data
    lduSolverPerformance solverPerf("cufflink_DiagPCG_Parallel",fieldName());

    register label nCells = x.size();
    register label NFaces = matrix().lower().size();
    OFSP.nCells = nCells;
    OFSP.nFaces = NFaces;
    OFSP.debugCusp = false;
    
    if (lduMatrix::debug >= 2) {
    	OFSP.debugCusp = true;
    }

//interface gathering

	//used when launching kernals in machines with multiple gpus
	int gpusPerMachine = readInt(dict().lookup("gpusPerMachine"));

	//get all the interface information
	#include "../CFL_Headers/getInterfaces.H"

	OFInterfaces.gpusPerMachine = gpusPerMachine;

	if(!OFInterfaces.checkGPUCount(OFInterfaces.gpusPerMachine)) {
		FatalErrorIn("checkGPUCount:\nThe number of GPUs per machine is not correct.\nChange the gpusPerMachine variable in fvSolution")<< "Fatal Error"<< abort(FatalError);
	}
//end interface gathering

	#include "../CFL_Headers/initializeCusp.H"

	//set the device
	cudaError_t cudareturn; 
	cudareturn = cudaSetDevice(Pstream::myProcNo() % (gpusPerMachine));
	if (cudareturn == cudaErrorInvalidDevice){ perror("cudaSetDevice returned  cudaErrorInvalidDevice");  } 

	CFL_DiagPCG_Parallel(&CES, &OFSP, &OFInterfaces);//call the CUDA code to solve the system
	
	thrust::copy(CES.X.begin(),CES.X.end(),x.begin());//copy the x vector back to Openfoam
	solverPerf.initialResidual() = OFSP.iRes;//return initial residual
	solverPerf.finalResidual() = OFSP.fRes;//return final residual
	solverPerf.nIterations() = OFSP.nIterations;//return the number of iterations
	//solverPerf.checkConvergence(tolerance(), relTolerance()) ; what should be passed here and does this waste time?
    return solverPerf;
}

Foam::scalarField Foam::CompositionModel<CloudType>::X
(
    const label phaseI,
    const scalarField& Y
) const
{
    const phaseProperties& props = phaseProps_[phaseI];
    scalarField X(Y.size(), 0.0);
    scalar WInv = 0.0;
    switch (props.phase())
    {
        case phaseProperties::GAS:
        {
            forAll(Y, i)
            {
                label gid = props.globalIds()[i];
                WInv += Y[i]/thermo_.carrier().W(gid);
                X[i] = Y[i]/thermo_.carrier().W(gid);
            }
            break;
        }
        case phaseProperties::LIQUID:
        {
            forAll(Y, i)
            {
                label gid = props.globalIds()[i];
                WInv += Y[i]/thermo_.liquids().properties()[gid].W();
                X[i] += Y[i]/thermo_.liquids().properties()[gid].W();
            }
            break;
        }

Foam::lduSolverPerformance Foam::bicgStabSolver::solve
(
    scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
    // Prepare solver performance
    lduSolverPerformance solverPerf(typeName, fieldName());

    scalarField p(x.size());
    scalarField r(x.size());

    // Calculate initial residual
    matrix_.Amul(p, x, coupleBouCoeffs_, interfaces_, cmpt);

    scalar normFactor = this->normFactor(x, b, p, r, cmpt);

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

    // Calculate residual
    forAll (r, i)
    {
        r[i] = b[i] - p[i];
    }

Foam::lduSolverPerformance Foam::cufflink_CG::solve
(
    scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
	OFSolverPerformance OFSP;//container for solver performance

	#include "../CFL_Headers/getGPUStorage.H"

	// --- Setup class containing solver performance data
    lduSolverPerformance solverPerf("cufflink_CG",fieldName());

    register label nCells = x.size();
    register label NFaces = matrix().lower().size();
    OFSP.nCells = nCells;
    OFSP.nFaces = NFaces;
    OFSP.debugCusp = false;
    
    if (lduMatrix::debug >= 2) {
         OFSP.debugCusp = true;
    }

	#include "../CFL_Headers/initializeCusp.H"

	CFL_CG(&CES, &OFSP);//call the CUDA code to solve the system
	
	thrust::copy(CES.X.begin(),CES.X.end(),x.begin());//copy the x vector back to Openfoam
	solverPerf.initialResidual() = OFSP.iRes;//return initial residual
	solverPerf.finalResidual() = OFSP.fRes;//return final residual
	solverPerf.nIterations() = OFSP.nIterations;//return the number of iterations
	//solverPerf.checkConvergence(tolerance(), relTolerance()) ; what should be passed here and does this waste time?
    return solverPerf;
}