[ChasingNeutrons]

Hi guys, I'm trying to create a volScalarField for a square matrix of values. Specifically, in a nuclear reactor physics context, I'm trying to make some convenient data structure to hold the scattering cross-section matrix which is a property of a given material and for g energy groups will have g^2 entries.

Having had great success and finding much convenience in storing other properties in the form of PtrList<volScalarField> my first inclination is to attempt to extend this method to a second dimension. For example, for the property sigmaT I've written in createFields.H:

PtrList<volScalarField> sigmaT(energyGroups);
forAll(sigmaT, i)
{
    sigmaT.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "sigmaT"+i,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar("", dimensionSet(0,-1,0,0,0,0,0), 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
    );
}

For my matrix of scattering values therefore I've tried:

PtrList<volScalarField> sigmaTo(energyGroups);
forAll(sigmaTo, energyJ)
{
    sigmaTo.set
    (
        energyJ,
        new volScalarField
        (
            IOobject
            (
                "sigmaTo"+Foam::name(energyJ),
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar("", dimensionSet(0,-1,0,0,0,0,0), 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
    );
}


PtrList<PtrList<volScalarField> > sigmaFrom(energyGroups);
forAll(sigmaFrom,energyI)
{    
    sigmaFrom.set(energyI,sigmaTo);
}

Unfortunately, this results in the following error:

In file included from /opt/openfoam30/src/OpenFOAM/lnInclude/tmp.H:139:0,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/PtrListI.H:27,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/PtrList.H:436,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/List.C:30,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/List.H:259,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/labelList.H:48,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/UPstream.H:42,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/Pstream.H:42,
                 from /opt/openfoam30/src/OpenFOAM/lnInclude/parRun.H:35,
                 from /opt/openfoam30/src/finiteVolume/lnInclude/fvCFD.H:4,
                 from MOCMEMRsolver.C:35:
/opt/openfoam30/src/OpenFOAM/lnInclude/tmpI.H: In instantiation of ‘Foam::tmp<T>::~tmp() [with T = Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >]’:
createFields.H:241:31:   required from here
/opt/openfoam30/src/OpenFOAM/lnInclude/tmpI.H:111:9: error: ‘class Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ has no member named ‘okToDelete’
         if (ptr_->okToDelete())
         ^
/opt/openfoam30/src/OpenFOAM/lnInclude/tmpI.H:118:13: error: ‘class Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ has no member named ‘operator--’
             ptr_->operator--();
             ^
/opt/openfoam30/src/OpenFOAM/lnInclude/tmpI.H: In instantiation of ‘T* Foam::tmp<T>::ptr() const [with T = Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >]’:
/opt/openfoam30/src/OpenFOAM/lnInclude/PtrListI.H:140:47:   required from ‘Foam::autoPtr<T> Foam::PtrList<T>::set(Foam::label, const Foam::tmp<T>&) [with T = Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >; Foam::label = int]’
createFields.H:241:31:   required from here
/opt/openfoam30/src/OpenFOAM/lnInclude/tmpI.H:162:10: error: ‘class Foam::PtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ has no member named ‘resetRefCount’
          ptr->resetRefCount();
          ^

My suspicion is that I've done something very silly by not generating energyGroup different instances of sigmaTo. I ask your patience as I'm still far from competent in either C++ or OpenFOAM! If anyone can suggest a relatively painless fix it would be very much appreciated.

Cheers.

[marupio]

Look at the VectorN library in foam-extend. It has been used specifically for nuclear reactions before. See Ivor Clifford's work.

[ChasingNeutrons]

Thanks David - I'm actually modifying and gutting one of the descendants of Ivor Clifford's work, a coupled thermal-hydraulics, neutronics, mechanics code called GeN-Foam (although I'm only interested in the neutronics at the moment). I'm trying to stick quite close to its format in the hope of potentially being able to replace its neutronics solver with my own in the future, hence the use of PtrLists.

Unfortunately, it seems to be the case that its implementation of scattering isn't entirely compatible with the method I'm using (or, at least, it would be quite inefficient to use it). If I don't have any luck soon, however, then I'll definitely take a look at foam-extend. Cheers.

[ChasingNeutrons]

Managed to get it working exactly as desired. More than anything the problem was getting the right syntax:

Info<<"sigmaS\n" << endl;
PtrList<PtrList<volScalarField> > sigmaS(energyGroups);
forAll(sigmaS, energyI)
{
    sigmaS.set(energyI, new PtrList<volScalarField>(energyGroups));
    forAll(sigmaS[energyI], energyJ)
    {
        sigmaS[energyI].set
        (
            energyJ,
            new volScalarField
            (
                IOobject
                (
                    "sigmaS"+Foam::name(energyI)+Foam::name(energyJ),
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                dimensionedScalar("", dimensionSet(0,-1,0,0,0,0,0), 0.0),
                zeroGradientFvPatchScalarField::typeName
            )
        );
    }
}

I was helped along by a previous post here: http://www.cfd-online.com/Forums/ope...ist-usage.html

[Jerryfan]

Hi,


This is a very good example of using PtrList. I have a simpler question in this case, why do you choose PtrList<PtrList<>> instead of ListList<>? Is there any reason that we cann't use the normal list?