[Shiqi Meng]

The EEqn.H is

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he)
      + fvm::div(phi, he,"div(phi,H)")
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p(),
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(D_H, he)//Alia
     ==
        rho*(U&g)
      + Qdot
      + coalParcels.Sh(he)
      + radiation->Sh(thermo, he)
      + fvOptions(rho, he)
    );

where the return of source term coalParcels.Sh(he) is

           return
                hsTrans()/Vdt
              - fvm::SuSp(hsCoeff()/(Cp*Vdt), hs)
              + hsCoeff()/(Cp*Vdt)*hs;

the first and third terms are volScalarField, the second one is a fvmScalarMatrix.

How they calculate?