运行结束后输出密度rho

[wangfei9088]

@李东岳 哈哈，李老师，我不是大佬，我们组里没人用PINN研究化学反应。看论坛里讨论得这么火热，我也想学习学习PINN的，就是还不熟悉。搞这个要向李老师多学习呀！

[尚善若水]

@wangfei9088 请教大佬，参照specieReactionRates的写法，先尝试输出每个反应步的体积平均Qdot(命名sumVQdotRi/V)，计算Qdot时需要物种的比焓h[speciei]，这个如何访问呢？

const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？

bool Foam::functionObjects::reactionSpecieQdots::write()
{
    logFiles::write();

    const label nSpecie = chemistryModel_.nSpecie();
    const label nReaction = chemistryModel_.nReaction();

    // Region volume
    const scalar V = this->V();

    for (label reactioni=0; reactioni<nReaction; reactioni++)
    {
        if (Pstream::master())
        {
            writeTime(file());
            file() << token::TAB << reactioni;
        }

        const PtrList<volScalarField::Internal> RR
        (
            chemistryModel_.reactionRR(reactioni)
        );

        scalar sumVQdotsR = 0;

        for (label speciei=0; speciei<nSpecie; speciei++)
        {
            scalar sumVQdotRi = 0;

            const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？

            if (all())
            {
                sumVQdotRi = fvc::domainIntegrate(RR[speciei]*(-h[speciei])).value();
            }
            else
            {
                sumVQdotRi =
                    gSum
                    (
                        scalarField
                        (
                            fvMeshFunctionObject::mesh_.V()*RR[speciei]*(-h[speciei]),
                            cells()
                        )
                    );
            }

            if (Pstream::master())
            {
                file() << token::TAB << sumVQdotRi/V;
            }

            sumVQdotsR += sumVQdotRi;
        }

        if (Pstream::master())
        {
            file() << token::TAB << sumVQdotsR/V;
        }

        if (Pstream::master())
        {
            file() << nl;
        }
    }

    if (Pstream::master())
    {
        file() << nl << endl;
    }

    return true;
}

[wangfei9088]

@尚善若水 没太看懂const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？这是啥意思。
如果需要得到每个组分的化学焓（生成焓），可以这么写：

const scalar hi = specieThermo_[speciei].Hc(); // 化学焓 J/kg

OpenFOAM版本不同，可能有些许差异，我用的是OpenFOAM-v2012。

[尚善若水]

@wangfei9088 这是我看到基金会版本里面用来算Qdot的，我直接贴上去了。我的问题是specieThermo_这个可以直接访问吗？const label nSpecie = chemistryModel_.nSpecie();可以这么写是前面定义了

我看到specieThermo_是在别的地方定义的。我弄不明白这个。

[wangfei9088]

@尚善若水 看了一下，开发版和of2012差别很大，开发版里不能用specieThermo_了。依据要不要单位，组分的生成焓在开发版里可以这么写：

dimensionedScalar hi = chemistryModel_.thermo().hfi(speciei); //- Enthalpy of formation [J/kg]
或者
scalar hi = chemistryModel_.thermo().hfiValue(speciei); //- Enthalpy of formation [J/kg]

[尚善若水]

@wangfei9088 好的，我仔细看下。谢谢。

[尚善若水]

@wangfei9088 感谢大佬，程序正常输出，正确性还没测试，也许可以和cantera的结果比一下，代码就是上面发的加上根据大佬说的改的。另外，我在想可不可以更进一步，写出每一个反应步每个组分的Qdot，命名为Qdot.Ri.Speciei，类似输出RR.Speciei?

[wangfei9088]

@尚善若水 应该可以吧。RR[speciei][celli]是celli网格里speciei的反应速率，那该网格里该组分的Qdot应该可以这么写：

scalar Qdot-speciei-celli = hi * RR_[speciei][celli];

[欧阳]

@尚善若水

对于你的问题，我最近也在修改，加入以下代码，然后在类定义中添加成员变量就行了。我这边是可以成功输出每个物质对应的Qdot

template<class ReactionThermo, class ThermoType>
Foam::tmp<Foam::volScalarField>
Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
{
    tmp<volScalarField> tQdot
    (
        new volScalarField
        (
            IOobject
            (
                "Qdot",
                this->mesh_.time().timeName(),
                this->mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                IOobject::NO_REGISTER
            ),
            this->mesh_,
            dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
        )
    );
QdotSpecies_;
    if (!QdotSpecies_.size())
    {
        // 首次调用时初始化
        QdotSpecies_.setSize(nSpecie_);
        
        forAll(Y_, i)
        {
            QdotSpecies_.set
            (
                i,
                new volScalarField
                (
                    IOobject
                    (
                        "Qdot." + Y_[i].name(),
                        this->mesh_.time().timeName(),
                        this->mesh_,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE,  
                        IOobject::REGISTER     
                    ),
                    this->mesh_,
                    dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
                )
            );
        }
    }
    
    if (this->chemistry_)
    {
        scalarField& Qdot = tQdot.ref();
        Qdot = 0.0;  
        
        forAll(Y_, i)
        {
            volScalarField& QdotI = QdotSpecies_[i];
            const scalar hi = specieThermo_[i].Hc();
            
            forAll(Qdot, celli)
            {

                QdotI[celli] = -hi*RR_[i][celli];
                

                Qdot[celli] += QdotI[celli];
            }
            
            // 修正边界条件
            QdotSpecies_[i].correctBoundaryConditions();
        }
        
        tQdot.ref().correctBoundaryConditions();
       
    }

    return tQdot;
}

[尚善若水]

@欧阳 多谢分享。这样是不是可以使用writeObjects输出想要的就可以了。我最后是按照wangfei9088大佬建议，仿照Qdot写了一个functionObject，但是后处理的话，会全部输出，占用空间很大。另外我见到的比较多的是按照cantera的处理，输出每步反应里面各个组分的放热和每步反应的总放热，也许你可以参考下。

[欧阳]

@尚善若水 我现在在我自己的代码的基础上，运行的结果是不管用不用writeObjects指定一些物质。所有物质的热量贡献都会输出

[尚善若水]

这里改一下就好了，不然详细机理太占空间了

[欧阳]

@尚善若水 写个functionObject，后处理，这个我没有试过，我之前一直都是在调用模型的原代码上修改的。我去学学，
“按照cantera的处理，输出每步反应里面各个组分的放热和每步反应的总放热”
这个是我现在在实现的，每个化学反应步在每个网格的Qdot，然后一些敏感性ODCantera那边的已经实现了，但是这边一直编译的或多或多出现一些问题

[欧阳]

@尚善若水 大佬有输出过tc吗？我按```
tc[celli] = nReaction*cSum/tc[celli];

StandardChemistryModel.C中的修改出现的结果，tc的值特别的大。反应几乎不发生时，特征时间趋于无穷大。但是这些极大值没有实际意义，反而会造成数值问题。

[尚善若水]

@欧阳 先请教下，组分总放热有什么分析场景吗？我看到较多的是分步的。tc的话你是用的什么燃烧模型？PaSR吗？我以前输出过，我找下

[欧阳]

@尚善若水 喷雾燃烧，当关注一些燃烧的重要时间点，低温反应着火与高温反应着火时间点，还有燃烧衰退时，组分的总放热蛮重要的，尤其当分析出，这个物质的对流，扩散与化学反应哪个占主导后。是的，PaSR

[尚善若水]

@欧阳 我以前是这样输出的，10版本之后好像不行了，我没测试

[欧阳]

@尚善若水

                    IOobject
                    (
                        "Qdot." + Y_[i].name(),
                        this->mesh_.time().timeName(),
                        this->mesh_,
                        IOobject::NO_READ,
                        IOobject::NO_WRITE,  
                        IOobject::REGISTER    
                    ),

这样修改后，我使用

writeObjects
(
    "Qdot.OH" 
);

后，没有输出Qdot.OH

[尚善若水]

@欧阳 先存储一下试试？

//必须得在functions里面定义启用才会计算，然后才能缓存
cacheTemporaryObjects
(
    //RR:methaneReaction:CH4 //speciekReactioni
    //RR:methaneReaction //reactioni

    //Qdot:methaneReaction:CH4
    //Qdot:methaneReaction
);

[欧阳]

@尚善若水 我现在用的是v2306，测试出来的tc值，如果没有反应就特别大