sensitivitySurfaceIncompressible.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2007-2021 PCOpt/NTUA
    Copyright (C) 2013-2021 FOSS GP
    Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "sensitivitySurfaceIncompressible.H"
#include "incompressibleAdjointSolver.H"
#include "primitivePatchInterpolation.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
#include "faMatrices.H"
#include "famSup.H"
#include "famLaplacian.H"
#include "volSurfaceMapping.H"
#include "fixedValueFaPatchFields.H"
#include "zeroGradientFaPatchFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace incompressible
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

defineTypeNameAndDebug(sensitivitySurface, 1);
addToRunTimeSelectionTable
(
    adjointSensitivity,
    sensitivitySurface,
    dictionary
);

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

void sensitivitySurface::addGeometricSens()
{
    if (includeObjective_)
    {
        // Grab objective refs
        PtrList<objective>& functions
            (objectiveManager_.getObjectiveFunctions());
        // Compute sens for all points in parameterized patches.
        // Interfacing points will be accumulated later
        autoPtr<pointBoundaryVectorField> pointSensdSdb
        (
            createZeroBoundaryPointFieldPtr<vector>(mesh_)
        );
        autoPtr<pointBoundaryVectorField> pointSensdndb
        (
            createZeroBoundaryPointFieldPtr<vector>(mesh_)
        );
        // Geometric (or "direct") sensitivities are better computed directly
        // on the points
        for (const label patchI : sensitivityPatchIDs_)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            const vectorField nf(patch.nf());

            // point sens result for patch
            vectorField& patchdSdb = pointSensdSdb()[patchI];
            vectorField& patchdndb = pointSensdndb()[patchI];

            vectorField dSdbMultiplierTot(patch.size(), Zero);
            vectorField dndbMultiplierTot(patch.size(), Zero);
            for (auto& fun : functions)
            {
                dSdbMultiplierTot += fun.weight()*fun.dSdbMultiplier(patchI);
                dndbMultiplierTot += fun.weight()*fun.dndbMultiplier(patchI);
            }
            // Correspondence of local point addressing to global point
            // addressing
            const labelList& meshPoints = patch.patch().meshPoints();
            //  List with mesh faces. Global addressing
            const faceList& faces = mesh_.faces();
            //  Each local patch point belongs to these local patch faces
            //  (local numbering)
            const labelListList& patchPointFaces = patch.patch().pointFaces();
            //  index of first face in patch
            const label patchStartIndex = patch.start();
            //  geometry differentiation engine
            deltaBoundary dBoundary(mesh_);
            //  Loop over patch points.
            //  Collect contributions from each boundary face this point
            //  belongs to
            forAll(meshPoints, ppI)
            {
                const labelList& pointFaces = patchPointFaces[ppI];
                for (label localFaceIndex : pointFaces)
                {
                    label globalFaceIndex = patchStartIndex + localFaceIndex;
                    const face& faceI = faces[globalFaceIndex];
                    // Point coordinates. All indices in global numbering
                    const pointField p(faceI.points(mesh_.points()));
                    tensorField p_d(faceI.size(), Zero);
                    forAll(faceI, facePointI)
                    {
                        if (faceI[facePointI] == meshPoints[ppI])
                        {
                            p_d[facePointI] = tensor::I;
                        }
                    }
                    const tensorField deltaNormals
                    (
                        dBoundary.makeFaceCentresAndAreas_d(p, p_d)
                    );

                    // Element [1] is the variation in the (dimensional) normal
                    const tensor& deltaSf = deltaNormals[1];
                    patchdSdb[ppI] +=
                        dSdbMultiplierTot[localFaceIndex] & deltaSf;

                    // Element [2] is the variation in the unit normal
                    const tensor& deltaNf = deltaNormals[2];
                    patchdndb[ppI] +=
                        dndbMultiplierTot[localFaceIndex] & deltaNf;
                }
            }
        }
        // Do parallel communications to avoid wrong values at processor
        // boundaries
        vectorField dSdbGlobal(mesh_.nPoints(), Zero);
        vectorField dndbGlobal(mesh_.nPoints(), Zero);
        for (const label patchI : sensitivityPatchIDs_)
        {
            const labelList& meshPoints =
                mesh_.boundaryMesh()[patchI].meshPoints();
            forAll(meshPoints, ppI)
            {
                const label globaPointI = meshPoints[ppI];
                dSdbGlobal[globaPointI] += pointSensdSdb()[patchI][ppI];
                dndbGlobal[globaPointI] += pointSensdndb()[patchI][ppI];
            }
        }
        // Accumulate over processors
        syncTools::syncPointList
        (
            mesh_, dSdbGlobal, plusEqOp<vector>(), vector::zero
        );
        syncTools::syncPointList
        (
            mesh_, dndbGlobal, plusEqOp<vector>(), vector::zero
        );
        // Transfer back to local fields and map to faces
        for (const label patchI : sensitivityPatchIDs_)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            const labelList& meshPoints = patch.patch().meshPoints();
            const scalarField& magSf = patch.magSf();
            pointSensdSdb()[patchI].map(dSdbGlobal, meshPoints);
            pointSensdndb()[patchI].map(dndbGlobal, meshPoints);
            // Map dSf/dx and dnf/dx term from points to faces.  Ideally, all
            // sensitivities should be computed at points rather than faces.
            PrimitivePatchInterpolation<polyPatch> patchInter(patch.patch());
            vectorField dSdbFace
            (
                patchInter.pointToFaceInterpolate(pointSensdSdb()[patchI])
            );
            // dSdb already contains the face area. Divide with it to make it
            // compatible with the rest of the terms
            dSdbFace /= magSf;

            tmp<vectorField> tdndbFace =
                patchInter.pointToFaceInterpolate(pointSensdndb()[patchI]);
            const vectorField& dndbFace = tdndbFace();

            // Add to sensitivity fields
            wallFaceSensVecPtr_()[patchI] += dSdbFace + dndbFace;
        }
    }
}


void sensitivitySurface::setSuffixName()
{
    word suffix(dict().getOrDefault<word>("suffix", word::null));
    // Determine suffix for fields holding the sens
    if (includeMeshMovement_)
    {
        shapeSensitivitiesBase::setSuffix
        (
            adjointVars_.solverName() + "ESI" + suffix
        );
    }
    else
    {
        shapeSensitivitiesBase::setSuffix
        (
            adjointVars_.solverName() + "SI" + suffix
        );
    }
}


void sensitivitySurface::smoothSensitivities()
{
    // Read in parameters
    const label iters(dict().getOrDefault<label>("iters", 500));
    const scalar tolerance(dict().getOrDefault<scalar>("tolerance", 1.e-06));
    autoPtr<faMesh> aMeshPtr(nullptr);

    IOobject faceLabels
    (
        "faceLabels",
        mesh_.time().findInstance
        (
            mesh_.dbDir()/faMesh::meshSubDir,
            "faceLabels",
            IOobject::READ_IF_PRESENT
        ),
        faMesh::meshSubDir,
        mesh_,
        IOobject::READ_IF_PRESENT,
        IOobject::NO_WRITE
    );

    // If the faMesh already exists, read it
    if (faceLabels.typeHeaderOk<labelIOList>(false))
    {
        Info<< "Reading the already constructed faMesh" << endl;
        aMeshPtr.reset(new faMesh(mesh_));
    }
    else
    {
        // Dictionary used to construct the faMesh
        dictionary faMeshDefinition;

        IOobject faMeshDefinitionDict
        (
            "faMeshDefinition",
            mesh_.time().caseSystem(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        );

        // If the faMeshDefinitionDict exists, use it to construct the mesh
        if (faMeshDefinitionDict.typeHeaderOk<IOdictionary>(false))
        {
            Info<< "Reading faMeshDefinition from system " << endl;
            faMeshDefinition = IOdictionary(faMeshDefinitionDict);
        }
        // Otherwise, faMesh is generated from all patches on which we compute
        // sensitivities
        else
        {
            Info<< "Constructing faMeshDefinition from sensitivity patches"
                << endl;
            wordList polyMeshPatches(sensitivityPatchIDs_.size());
            label i(0);
            for (const label patchID : sensitivityPatchIDs_)
            {
                polyMeshPatches[i++] = mesh_.boundary()[patchID].name();
            }
            faMeshDefinition.add<wordList>("polyMeshPatches", polyMeshPatches);
            (void)faMeshDefinition.subDictOrAdd("boundary");
        }

        // Construct faMesh
        aMeshPtr.reset(new faMesh(mesh_, faMeshDefinition));
    }
    faMesh& aMesh = aMeshPtr.ref();

    // Physical radius of the smoothing, provided either directly or computed
    // based on the average 'length' of boundary faces
    const scalar Rphysical
        (dict().getOrDefault<scalar>("radius", computeRadius(aMesh)));
    DebugInfo
        << "Physical radius of the sensitivity smoothing "
        << Rphysical << nl << endl;

    // Radius used as the diffusivity in the Helmholtz filter, computed as a
    // function of the physical radius
    const dimensionedScalar RpdeSqr
    (
        "RpdeSqr", dimArea, sqr(Rphysical/(2.*::sqrt(3.)))
    );

    dimensionedScalar one("1", dimless, 1.);

    // Mapping engine
    const volSurfaceMapping vsm(aMesh);

    // Source term in faMatrix needs to be an areaField
    areaScalarField sens
    (
        IOobject
        (
            "sens",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        aMesh,
        dimensionedScalar(dimless, Zero),
        faPatchFieldBase::zeroGradientType()
    );

    // Copy sensitivities to area field
    sens.primitiveFieldRef() =
        vsm.mapToSurface<scalar>(wallFaceSensNormalPtr_());

    // Initialisation of the smoothed sensitivities field based on the original
    // sensitivities
    areaScalarField smoothedSens("smoothedSens", sens);
    for (label iter = 0; iter < iters; ++iter)
    {
        Info<< "Sensitivity smoothing iteration " << iter << endl;

        faScalarMatrix smoothEqn
        (
            fam::Sp(one, smoothedSens)
          - fam::laplacian(RpdeSqr, smoothedSens)
         ==
            sens
        );

        smoothEqn.relax();

        const scalar residual(mag(smoothEqn.solve().initialResidual()));

        DebugInfo
            << "Max smoothSens " << gMax(mag(smoothedSens)()) << endl;

        // Print execution time
        mesh_.time().printExecutionTime(Info);

        // Check convergence
        if (residual < tolerance)
        {
            Info<< "\n***Reached smoothing equation convergence limit, "
                   "iteration " << iter << "***\n\n";
            break;
        }
    }

    // Field used to write the smoothed sensitivity field to file
    volScalarField volSmoothedSens
    (
        IOobject
        (
            "smoothedSurfaceSens" + surfaceFieldSuffix_,
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar(dimless, Zero)
    );

    // Transfer result back to volField and write
    vsm.mapToVolume(smoothedSens, volSmoothedSens.boundaryFieldRef());
    volSmoothedSens.write();
}


scalar sensitivitySurface::computeRadius(const faMesh& aMesh)
{
    scalar averageArea(gAverage(aMesh.S().field()));
    const Vector<label>& geometricD = mesh_.geometricD();
    const boundBox& bounds = mesh_.bounds();
    forAll(geometricD, iDir)
    {
        if (geometricD[iDir] == -1)
        {
            averageArea /= bounds.span()[iDir];
        }
    }
    scalar mult = dict().getOrDefault<scalar>("meanRadiusMultiplier", 10);

    return mult*pow(averageArea, scalar(1)/scalar(mesh_.nGeometricD() - 1));
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

sensitivitySurface::sensitivitySurface
(
    const fvMesh& mesh,
    const dictionary& dict,
    incompressibleAdjointSolver& adjointSolver
)
:
    adjointSensitivity(mesh, dict, adjointSolver),
    shapeSensitivitiesBase(mesh, dict),
    includeSurfaceArea_(false),
    includePressureTerm_(false),
    includeGradStressTerm_(false),
    includeTransposeStresses_(false),
    useSnGradInTranposeStresses_(false),
    includeDivTerm_(false),
    includeDistance_(false),
    includeMeshMovement_(false),
    includeObjective_(false),
    writeGeometricInfo_(false),
    smoothSensitivities_(false),
    eikonalSolver_(nullptr),
    meshMovementSolver_(nullptr),

    nfOnPatchPtr_(nullptr),
    SfOnPatchPtr_(nullptr),
    CfOnPatchPtr_(nullptr)
{
    read();
    setSuffixName();

    // Allocate boundary field pointer
    wallFaceSensVecPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
    wallFaceSensNormalPtr_.reset(createZeroBoundaryPtr<scalar>(mesh_));
    wallFaceSensNormalVecPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));

    // Allocate fields to contain geometric info
    if (writeGeometricInfo_)
    {
        nfOnPatchPtr_.reset
        (
            new volVectorField
            (
                IOobject
                (
                    "nfOnPatch",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                vector::zero
            )
        );

        SfOnPatchPtr_.reset
        (
            new volVectorField
            (
                IOobject
                (
                    "SfOnPatch",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                vector::zero
            )
        );

        CfOnPatchPtr_.reset
        (
            new volVectorField
            (
                IOobject
                (
                    "CfOnPatch",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                vector::zero
            )
        );
    }

    // Allocate appropriate space for the sensitivity field
    computeDerivativesSize();
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void sensitivitySurface::read()
{
    includeSurfaceArea_ =
        dict().getOrDefault<bool>("includeSurfaceArea", true);
    includePressureTerm_ =
        dict().getOrDefault<bool>("includePressure", true);
    includeGradStressTerm_ =
        dict().getOrDefault<bool>("includeGradStressTerm", true);
    includeTransposeStresses_ =
        dict().getOrDefault<bool>("includeTransposeStresses", true);
    useSnGradInTranposeStresses_ =
        dict().getOrDefault<bool>("useSnGradInTranposeStresses", false);
    includeDivTerm_ = dict().getOrDefault<bool>("includeDivTerm", false);
    includeDistance_ =
        dict().getOrDefault<bool>
        (
            "includeDistance",
            adjointVars_.adjointTurbulence().ref().includeDistance()
        );
    includeMeshMovement_ =
        dict().getOrDefault<bool>("includeMeshMovement", true);
    includeObjective_ =
        dict().getOrDefault<bool>("includeObjectiveContribution", true);
    writeGeometricInfo_ =
        dict().getOrDefault<bool>("writeGeometricInfo", false);
    smoothSensitivities_ =
        dict().getOrDefault<bool>("smoothSensitivities", false);

    // Allocate new solvers if necessary
    if (includeDistance_ && !eikonalSolver_)
    {
        eikonalSolver_.reset
        (
            new adjointEikonalSolver
            (
                mesh_,
                dict_,
                primalVars_.RASModelVariables(),
                adjointVars_,
                sensitivityPatchIDs_
            )
        );
    }
    if (includeMeshMovement_ && !meshMovementSolver_)
    {
        meshMovementSolver_.reset
        (
            new adjointMeshMovementSolver
            (
                mesh_,
                dict_,
                *this,
                sensitivityPatchIDs_,
                eikonalSolver_
            )
        );
    }
}


bool sensitivitySurface::readDict(const dictionary& dict)
{
    if (sensitivity::readDict(dict))
    {
        if (eikonalSolver_)
        {
            eikonalSolver_().readDict(dict);
        }

        if (meshMovementSolver_)
        {
            meshMovementSolver_().readDict(dict);
        }

        return true;
    }

    return false;
}


void sensitivitySurface::computeDerivativesSize()
{
    label nFaces(0);
    for (const label patchI : sensitivityPatchIDs_)
    {
        nFaces += mesh_.boundary()[patchI].size();
    }
    derivatives_.setSize(nFaces);
}


void sensitivitySurface::accumulateIntegrand(const scalar dt)
{
    // Grab references
    const volScalarField& p = primalVars_.p();
    const volVectorField& U = primalVars_.U();

    const volScalarField& pa = adjointVars_.pa();
    const volVectorField& Ua = adjointVars_.Ua();
    autoPtr<incompressibleAdjoint::adjointRASModel>& adjointTurbulence =
        adjointVars_.adjointTurbulence();

    // Accumulate source for additional post-processing PDEs, if necessary
    if (includeDistance_)
    {
        eikonalSolver_->accumulateIntegrand(dt);
    }

    if (includeMeshMovement_)
    {
        meshMovementSolver_->accumulateIntegrand(dt);
    }

    // Terms from the adjoint turbulence model
    const boundaryVectorField& adjointTMsensitivities =
        adjointTurbulence->wallShapeSensitivities();

    DebugInfo
        << "    Calculating adjoint sensitivity. " << endl;

    tmp<volScalarField> tnuEff = adjointTurbulence->nuEff();
    const volScalarField& nuEff = tnuEff.ref();

    // Sensitivities do not include locale surface area by default.
    // The part of the sensitivities that multiplies dxFace/db follows

    // Deal with the stress part first since it's the most awkward in terms
    // of memory managment
    if (includeGradStressTerm_)
    {
        // Terms corresponding to contributions from converting delta
        // to thetas are added through the corresponding adjoint
        // boundary conditions instead of grabbing contributions from
        // the objective function.  Useful to have a unified
        // formulation for low- and high-re meshes

        tmp<volVectorField> tgradp = fvc::grad(p);
        const volVectorField& gradp = tgradp.ref();
        for (const label patchI : sensitivityPatchIDs_)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            tmp<vectorField> tnf = patch.nf();
            const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
            wallFaceSensVecPtr_()[patchI] -=
                (Uab & tnf)*gradp.boundaryField()[patchI]*dt;
        }
        tgradp.clear();

        // We only need to modify the boundaryField of gradU locally.
        // If grad(U) is cached then
        // a. The .ref() call fails since the tmp is initialised from a
        //    const ref
        // b. we would be changing grad(U) for all other places in the code
        //    that need it
        // So, always allocate new memory and avoid registering the new field
        tmp<volTensorField> tgradU =
            volTensorField::New("gradULocal", fvc::grad(U));
        volTensorField::Boundary& gradUbf = tgradU.ref().boundaryFieldRef();

        // Explicitly correct the boundary gradient to get rid of the
        // tangential component
        forAll(mesh_.boundary(), patchI)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            if (isA<wallFvPatch>(patch))
            {
                tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
                gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
            }
        }

        tmp<volSymmTensorField> tstress = nuEff*twoSymm(tgradU);
        const volSymmTensorField& stress = tstress.cref();
        autoPtr<volVectorField> ptemp
            (Foam::createZeroFieldPtr<vector>(mesh_, "temp", sqr(dimVelocity)));
        volVectorField& temp = ptemp.ref();
        for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
        {
            unzipRow(stress, idir, temp);
            volTensorField gradStressDir(fvc::grad(temp));
            for (const label patchI : sensitivityPatchIDs_)
            {
                const fvPatch& patch = mesh_.boundary()[patchI];
                tmp<vectorField> tnf = patch.nf();
                const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
                wallFaceSensVecPtr_()[patchI] +=
                    (
                        Uab.component(idir)
                       *(gradStressDir.boundaryField()[patchI] & tnf)
                    )*dt;
            }
        }
    }

    // Transpose part of the adjoint stresses
    // Dealt with separately to deallocate gradUa as soon as possible
    if (includeTransposeStresses_)
    {
        tmp<volTensorField> tgradUa = fvc::grad(Ua);
        const volTensorField::Boundary& gradUabf =
            tgradUa.cref().boundaryField();

        for (const label patchI : sensitivityPatchIDs_)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            tmp<vectorField> tnf = patch.nf();
            const vectorField& nf = tnf();
            vectorField gradUaNf
                (
                    useSnGradInTranposeStresses_
                  ? (Ua.boundaryField()[patchI].snGrad() & nf)*nf
                  : (gradUabf[patchI] & nf)
                );
            wallFaceSensVecPtr_()[patchI] -=
                nuEff.boundaryField()[patchI]
               *(gradUaNf & U.boundaryField()[patchI].snGrad())*nf;
        }
    }

    for (const label patchI : sensitivityPatchIDs_)
    {
        const fvPatch& patch = mesh_.boundary()[patchI];
        tmp<vectorField> tnf = patch.nf();
        const vectorField& nf = tnf();

        // Includes spurious tangential gradU part. Deprecated
        /*
        vectorField stressAndPressureTerm =
              (
                - (
                       Ua.boundaryField()[patchI].snGrad()
                    + (gradUa.boundaryField()[patchI] & nf)
                  ) * nuEff.boundaryField()[patchI]
                + pa.boundaryField()[patchI] *nf
              ) & gradU.boundaryField()[patchI].T();
        */

        // Adjoint stress term
        vectorField stressTerm
        (
          - (
                Ua.boundaryField()[patchI].snGrad()
              & U.boundaryField()[patchI].snGrad()
            )
          * nuEff.boundaryField()[patchI]
          * nf
        );

        if (includeDivTerm_)
        {
            stressTerm +=
                scalar(1./3.)*nuEff.boundaryField()[patchI]
              * (
                    ((Ua.boundaryField()[patchI].snGrad() &nf)*nf)
                  & U.boundaryField()[patchI].snGrad()
                )
              * nf;
        }

        // Adjoint pressure terms
        vectorField pressureTerm(patch.size(), Zero);
        if (includePressureTerm_)
        {
            pressureTerm =
            (
                (nf*pa.boundaryField()[patchI])
                & U.boundaryField()[patchI].snGrad()
            )* nf;
        }

        PtrList<objective>& functions
            (objectiveManager_.getObjectiveFunctions());

        // Term from objectives including x directly (e.g. moments)
        vectorField dxdbMultiplierTot(pressureTerm.size(), Zero);
        if (includeObjective_)
        {
            forAll(functions, funcI)
            {
                dxdbMultiplierTot +=
                    functions[funcI].weight()
                  * (
                        functions[funcI].dxdbDirectMultiplier(patchI)
                    );
            }
        }

        // Fill in sensitivity fields
        wallFaceSensVecPtr_()[patchI] +=
        (
            stressTerm
          + pressureTerm
          + adjointTMsensitivities[patchI]
          + dxdbMultiplierTot
        )*dt;
    }

    // Add terms from physics other than the typical incompressible flow eqns
    adjointSolver_.additionalSensitivityMapTerms
        (wallFaceSensVecPtr_(), sensitivityPatchIDs_, dt);

    // Add the sensitivity part corresponding to changes of the normal vector
    // Computed at points and mapped to faces
    addGeometricSens();
}


void sensitivitySurface::assembleSensitivities()
{
    // Update geometric fields for use by external users
    if (writeGeometricInfo_)
    {
        for (const label patchI : sensitivityPatchIDs_)
        {
            const fvPatch& patch = mesh_.boundary()[patchI];
            tmp<vectorField> tnf = patch.nf();
            const vectorField& nf = tnf();
            const vectorField& Sf = patch.Sf();
            const vectorField& Cf = patch.Cf();

            nfOnPatchPtr_().boundaryFieldRef()[patchI] = nf;
            SfOnPatchPtr_().boundaryFieldRef()[patchI] = Sf;
            CfOnPatchPtr_().boundaryFieldRef()[patchI] = Cf;
        }
    }

    // Solve extra equations if necessary
    // Solved using accumulated sources over time
    autoPtr<boundaryVectorField> distanceSensPtr(nullptr);
    if (includeDistance_)
    {
        eikonalSolver_->solve();
        distanceSensPtr.reset(createZeroBoundaryPtr<vector>(mesh_));
        const boundaryVectorField& sens =
            eikonalSolver_->distanceSensitivities();
        for (const label patchI : sensitivityPatchIDs_)
        {
            distanceSensPtr()[patchI] = sens[patchI];
        }
    }

    autoPtr<boundaryVectorField> meshMovementSensPtr(nullptr);
    if (includeMeshMovement_)
    {
        meshMovementSolver_->solve();
        meshMovementSensPtr.reset(createZeroBoundaryPtr<vector>(mesh_));
        const boundaryVectorField& sens =
            meshMovementSolver_->meshMovementSensitivities();
        for (const label patchI : sensitivityPatchIDs_)
        {
            meshMovementSensPtr()[patchI] = sens[patchI];
        }
    }


    // Project to normal face vector
    label nPassedFaces(0);
    for (const label patchI : sensitivityPatchIDs_)
    {
        const fvPatch& patch = mesh_.boundary()[patchI];
        tmp<vectorField> tnf(patch.nf());
        const vectorField& nf = tnf();

        // Distance related terms
        if (includeDistance_)
        {
            wallFaceSensVecPtr_()[patchI] += distanceSensPtr()[patchI];
        }

        // Mesh movement related terms
        if (includeMeshMovement_)
        {
            wallFaceSensVecPtr_()[patchI] += meshMovementSensPtr()[patchI];
        }

        if (includeSurfaceArea_)
        {
            wallFaceSensVecPtr_()[patchI] *= patch.magSf();
        }

        wallFaceSensNormalPtr_()[patchI] = wallFaceSensVecPtr_()[patchI] & nf;
        wallFaceSensNormalVecPtr_()[patchI] =
            wallFaceSensNormalPtr_()[patchI] * nf;

        forAll(patch, fI)
        {
            derivatives_[nPassedFaces + fI]
                = wallFaceSensNormalPtr_()[patchI][fI];
        }
        nPassedFaces += patch.size();
    }

    // Smooth sensitivities if needed
    if (smoothSensitivities_)
    {
        smoothSensitivities();
    }
}


void sensitivitySurface::clearSensitivities()
{
    // Reset terms in post-processing PDEs
    if (includeDistance_)
    {
        eikonalSolver_->reset();
    }
    if (includeMeshMovement_)
    {
        meshMovementSolver_->reset();
    }
    // Reset sensitivity fields
    adjointSensitivity::clearSensitivities();
    shapeSensitivitiesBase::clearSensitivities();
}


autoPtr<adjointEikonalSolver>& sensitivitySurface::getAdjointEikonalSolver()
{
    return eikonalSolver_;
}


void sensitivitySurface::write(const word& baseName)
{
    setSuffixName();
    adjointSensitivity::write();
    shapeSensitivitiesBase::write();

    if (writeGeometricInfo_)
    {
        nfOnPatchPtr_().write();
        SfOnPatchPtr_().write();
        CfOnPatchPtr_().write();
    }
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam
} // End namespace incompressible

// ************************************************************************* //