void Foam::chtPimpleApp::Manager::write() const
{
    Region_FLUID::Storage& storage = regions().region_FLUID().storage();


    const volScalarField& T     = storage.thermo().T();
    const volScalarField& kappa = storage.Kappa();

    forAll(T.boundaryField(), patchI)
    {
        const fvPatchScalarField& Tpatch = T.boundaryField()[patchI];

        if (isType<solidWallMixedTemperatureCoupledFvPatchScalarField>(Tpatch))
        {
            const fvPatchScalarField& Tpatch = T.boundaryField()[patchI];

            // Thermal conductivity at this patch
            const scalarField& kappaPatch = kappa.boundaryField()[patchI];

            // Compute per-face heat flux: q = kappa * (gradT · n) * faceArea
            scalarField heatFluxPatch = kappaPatch * Tpatch.snGrad() * Tpatch.patch().magSf();

            surfaceScalarField heatFlux
            (
                IOobject
                (
                    "heatFlux",       
                    this->time().timeName(),  
                    T.mesh(),       
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                T.mesh(),            
                dimensionedScalar("zero", dimensionSet(1,0,-3,0,0,0,0), 0.0) // initialize all faces to 0
            );

            heatFlux.boundaryField()[patchI] = heatFluxPatch;

            heatFlux.write();
        }
    }
}