tmp类使用问题，承接上次的帖子问题

[悬浮的猪]

根据各位老师的建议这次我将代码直接贴上来，补充说明问题，有老师说我的错误提示应该在tmp的使用不当上面，但是我现在对具体出错的地方还是不甚明白。我之前对于tmp的理解停留在一个智能的指针分配管理者上，知道它能减少内存的调配，提高效率，但实际所知并不多，对它的使用也停留在仿照源代码湍流模型上的类似公式。这次的问题我直接贴上来，承接上次帖子贴出来的问题：

tmp<fvScalarMatrix> epstEqn
(
    fvm::ddt(epsilont_)
  + fvm::div(phi_, epsilont_)
  - fvm::laplacian(alphaEff(), epsilont_) 
 ==
     Cp1_*Gt*epsilont_/kt_
  + Cp2_*G*epsilont_/k_
  - Cd1_*(sqr(epsilont_)/kt_)
  - Cd2_*((epsilont_*epsilon_)/k_)
);

epstEqn.ref().relax();
epstEqn.ref().boundaryManipulate(epsilont_.boundaryFieldRef());
solve(epstEqn);   
bound(epsilont_, epsilontMin_);
tmp<fvScalarMatrix> ktEqn
(
    fvm::ddt(kt_)
  + fvm::div(phi_, kt_)
  - fvm::laplacian(alphaEff(), kt_)
 ==
    - Gt//速度产生项
  - fvm::Sp(2*epsilont_/kt_, kt_)
  
);

ktEqn.ref().relax();
solve(ktEqn);
bound(kt_, ktMin_);

correctalphat();

[悬浮的猪]

仅仅运行一步直接报错，这是具体的错误提示：

Time = 1

 smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0667125, No Iterations 3
 smoothSolver:  Solving for Uy, Initial residual = 1, Final residual =     0.0804578, No Iterations 3
 GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0500194, No Iterations 4
 time step continuity errors : sum local = 0.879188, global = -0.272208, cumulative = -0.272208
 smoothSolver:  Solving for T, Initial residual = 1, Final residual = 0.0597181, No Iterations 2
smoothSolver:  Solving for epsilon, Initial residual = 0.418304, Final residual = 0.0304297, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.00011942, No Iterations 1
smoothSolver:  Solving for epsilont, Initial residual = 1, Final residual = 0.0345744, No Iterations 2
bounding epsilont, min: -1800.64 max: 30630.5 average: 1357.93


  --> FOAM FATAL ERROR: 
 tmp<N4Foam14GeometricFieldIdNS_12fvPatchFieldENS_7volMeshEEE> deallocated

From function const T& Foam::tmp<T>::operator()() const [with T = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]
in file /opt/openfoam6/src/OpenFOAM/lnInclude/tmpI.H at line 278.

[风大仙]

@悬浮的猪 请问楼主解决了么？我修改湍流模型也遇到了一样的问题

[李东岳]

这个问题是某个tmp类型，已经销毁掉了，但是又被继续调用了。tmp类型使用的使用可以后面加()试一下

[风大仙]

@李东岳 李老师我加了（）还是不行啊。。。纠结两天了。。。。
顺手贴下问题：主要是在of org版本上添加kOmegaSSTIDDES模型，但是org版本只有kOmegaSSTDES,所以我把openfoam +版本的kOmegaSSTIDDES移植过来，涉及到长度尺度dTilda()函数，目前跟露珠遇到的问题一样，编译通过但计算一步就崩，报错信息：

--> FOAM FATAL ERROR: 
tmp<N4Foam14GeometricFieldINS_6TensorIdEENS_12fvPatchFieldENS_7volMeshEEE> deallocated

    From function const T& Foam::tmp<T>::operator()() const [with T = Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh>]
    in file /home/fxr/OpenFOAM/OpenFOAM-7/src/OpenFOAM/lnInclude/tmpI.H at line 278.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >::operator()() const at ??:?
#3  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::mag<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#4  Foam::LESModels::kOmegaSSTIDDESTr<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::correct() at ??:?
#5  ? in "/home/fxr/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt64Opt/bin/pimpleFoam"
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? in "/home/fxr/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt64Opt/bin/pimpleFoam"
Aborted (core dumped)

dTilda代码：

template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSSTIDDESTr<BasicTurbulenceModel>::dTilda
(
    const volScalarField& magGradU,
    const volScalarField& CDES
) const
{
    const volScalarField& k = this->k_;
    const volScalarField& omega = this->omega_;

    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
    const volScalarField lLES(CDES*this->delta());

    const volScalarField alpha(this->alpha());
    const volScalarField expTerm(exp(sqr(alpha)));

    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
    tmp<volScalarField> fe1 =
        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;

    const volScalarField fdTilda(max(1 - fdt(magGradU), fB));

   
    return max
    (
        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
        dimensionedScalar("SMALL", dimLength, SMALL)
    );
}

后面k方程中包含了dTilda的调用：

 // Turbulent kinetic energy equation
    volScalarField CDES(this->CDES(F1));
    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
     ==
        alpha()*rho()*Pk(G)
      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
      - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(tgradU), CDES)(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

各位老师有没有什么建议啊？？？谢谢了

[李东岳]

你最好debug到具体某一行（比如这一行去掉就没问题，这一行放进去就有问题）

[风大仙]

@李东岳 谢谢，按照您的方法试了，就是这行的问题，一去掉就可以算

- fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(tgradU), CDES)(), k_)

[李东岳]

tgradU是不是tmp，前面怎么用的

[风大仙]

@李东岳 李老师，前面有个epsilonByk函数，

template<class BasicTurbulenceModel>
tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
(
    const volScalarField& F1,
    const volTensorField& gradU
) const
{
    volScalarField CDES(this->CDES(F1));
    return sqrt(this->k_())/dTilda(mag(gradU), CDES)()();
}

但是我把tgradU换成gradU,编译报错：

kOmegaSSTIDDESTr.C: In member function ‘virtual void Foam::LESModels::kOmegaSSTIDDESTr<BasicTurbulenceModel>::correct()’:
kOmegaSSTIDDESTr.C:1248:59: error: ‘gradU’ was not declared in this scope
       - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(gradU), CDES)(), k_)
                                                           ^~~~~
kOmegaSSTIDDESTr.C:1248:59: note: suggested alternative: ‘tgradU’
       - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(gradU), CDES)(), k_)
                                                           ^~~~~
                                                           tgradU
/home/fxr/OpenFOAM/OpenFOAM-7/wmake/rules/General/transform:25: recipe for target 'Make/linux64GccDPInt64Opt/IDDESturbulentTransportModels.o' failed
make: *** [Make/linux64GccDPInt64Opt/IDDESturbulentTransportModels.o] Error 1

[李东岳]

前面的tgradU怎么用的，把涉及到tgradU的代码贴一下

[风大仙]

@李东岳 感谢李老师，前面整个correct函数是这样的

template<class BasicTurbulenceModel>
void kOmegaSSTIDDESTr<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    volScalarField& nut = this->nut_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    BasicTurbulenceModel::correct();

    volScalarField::Internal divU
    (
        fvc::div(fvc::absolute(this->phi(), U))()()
    );

    tmp<volTensorField> tgradU = fvc::grad(U);
    volScalarField S2(2*magSqr(symm(tgradU())));
    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());
    volScalarField::Internal G(this->GName(), nut()*GbyNu);
    tgradU.clear();

    // Update omega and G at the wall
    omega_.boundaryFieldRef().updateCoeffs();

    volScalarField CDkOmega
    (
        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
    );

    volScalarField F1(this->F1(CDkOmega));
    volScalarField F23(this->F23());

    {
        volScalarField::Internal gamma(this->gamma(F1));
        volScalarField::Internal beta(this->beta(F1));

        // Turbulent frequency equation
        tmp<fvScalarMatrix> omegaEqn
        (
            fvm::ddt(alpha, rho, omega_)
          + fvm::div(alphaRhoPhi, omega_)
          - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
         ==
            alpha()*rho()*gamma
           *min
            (
                GbyNu,
                (c1_/a1_)*betaStar_*omega_()
               *max(a1_*omega_(), b1_*F23()*sqrt(S2()))
            )
          - fvm::SuSp((2.0/3.0)*alpha()*rho()*gamma*divU, omega_)
          - fvm::Sp(alpha()*rho()*beta*omega_(), omega_)
          - fvm::SuSp
            (
                alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
                omega_
            )
          + Qsas(S2(), gamma, beta)
          + omegaSource()
          + fvOptions(alpha, rho, omega_)
        );

        omegaEqn.ref().relax();
        fvOptions.constrain(omegaEqn.ref());
        omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
        solve(omegaEqn);
        fvOptions.correct(omega_);
        bound(omega_, this->omegaMin_);
    }
k方程。。。

[李东岳]

• 你把tgradU.clear();去掉，
• 然后使用的时候用tgradU()()

[风大仙]

@李东岳 非常感谢李老师不厌其烦的指导！！tmp的问题已经解决了，不过又出现了量纲的问题，我自己再看看。