nullSpace.H

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) 2023 PCOpt/NTUA
    Copyright (C) 2023 FOSS GP
-------------------------------------------------------------------------------
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/>.


Class
    Foam::nullSpace

Description
    Update design variables using a null space approach.
    Can handle inequality and bound constraints.

    Reference:
    \verbatim
        Feppon, F., Allaire, G., & Dapogny, C. (2020).
        Null space gradient flows for constrained optimization with
        applications to shape optimization.
        ESAIM: COCV, 26, 90.
        https://doi.org/10.1051/cocv/2020015
    \endverbatim

SourceFiles
    nullSpace.C

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

#ifndef nullSpace_H
#define nullSpace_H

#include "constrainedOptimisationMethod.H"
#include "updateMethod.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                          Class nullSpace Declaration
\*---------------------------------------------------------------------------*/

class nullSpace
:
    public constrainedOptimisationMethod,
    public updateMethod
{
protected:

    // Protected data

        //- Are bound constraints included?
        bool includeBoundConstraints_;

        //- Lagrange multipliers for flow-reated constraints
        scalarField mu_;

        //- Lagrange multipliers for the lower bound constraints
        scalarField l_;

        //- Lagrange multipliers for the upper bound constraints
        scalarField u_;


        // Dual variables

            //- Solve the dual problem?
            bool solveDualProblem_;

            //- Lagrange multipliers of the dual problem for flow-related
            //- constraints
            scalarField dualMu_;

            //- Lagrange multipliers of the dual problem for the lower bound
            //- constraints
            scalarField dualL_;

            //- Lagrange multipliers of the dual problem for the upper bound
            //- constraints
            scalarField dualU_;


        // Constraint subsets

            //- List of saturated or violated constraints
            //  List[0]: flow related constraints
            //  List[1]: lower bound constraints
            //  List[2]: upper bound constraints
            //  Sub-list meaning is the same for all labelListLists that follow
            labelListList iTilda_;

            //- List of saturated or violated constraints (up to epsConstr_)
            labelListList iTildaEps_;

            //- List of constraints that must remain active
            //  Determines the null space update
            labelListList iHat_;

            //- List of constraints the values of which need to be reduced
            //  Determines the range space update
            labelListList iRangeSpace_;


        // Fields holding the updates of the Lagrange multipliers of the primal
        // and dual problems

            scalarField deltaMu_;
            scalarField deltaL_;
            scalarField deltaU_;
            scalarField deltaDualMu_;
            scalarField deltaDualL_;
            scalarField deltaDualU_;


        //- Infinitesimal quantity
        //  Updated during the inner iterations of the dual problem
        scalar eps_;

        //- Maxmimum number of Newton iterations for the dual problem
        label maxNewtonIters_;

        //- Maxmimum number of line search iterations for each Newton iteration
        //- of the dual problem
        label maxLineSearchIters_;

        //- Maxmimum number of CG iterations for obtaining the null space and
        //- range space updates
        label maxCGIters_;

        //- Tolerance of the dual problem
        scalar dualTolerance_;


        // Constant parameters

            //- Value for considering a constraint as marginally active
            //  Used to avoid the frequent change of the active set of
            //  constraints
            scalar epsConstr_;

            //- Value to perturb the design variables by, if all of them
            //- lay on their bounds at the beginning of the optimisation
            scalar epsPerturb_;

            //- Multiplier of the null space update
            scalar aJ_;

            //- Multiplier of the range space update
            scalar aC_;

            //- Change aJ and aC adaptively
            bool adaptiveStep_;

            //- Last cycle on which to reset aJ
            label lastAcceleratedCycle_;

            //- Max change of the design variables, pertaining to the objective
            //- reduction.
            //  If adaptiveStep is set to true, aJ will be set in such a way
            //  to obtain a maxDVChange_ for the design variables, for all
            //  optimisation cycles up to lastAcceleratedCycle_
            autoPtr<scalar> maxDVChange_;

            //- Whether to impose maxDVChange strictly on all optimisation
            //- cycles or just up to the lastAcceleratedCycle
            bool strictMaxDVChange_;

            //- Target reduction of active constraints (range in [0-1])
            //  1: active constraints will become zero (first-order approx)
            //  0: no change in the constraints
            scalar targetConstraintReduction_;

            //- Downplay the importance of the bound constraint reduction
            //- by this ammount
            scalar bcMult_;


    // Protected Member Functions

        //- Update sizes of fields related to the constraints
        void initialise();

        //- Zero the updates computed in the previous optimisation cycle
        void zeroUpdates();


        // Functions related to the solution of the dual subproblem

            //- Solve the dual problem for computing the Lagrange multiplier
            //- using a Newton solver
            void solveDualProblem();

            //- Compute and return residuals based on the current solution
            tmp<scalarField> computeResiduals();

            //- Compute direction for the Newton problem
            void computeNewtonDirection();

            //- Perform line search and return max residual corresponding to
            //- the updated solution
            scalar lineSearch();

            //- Update the current solution using the known direction and the
            //- given step length
            void updateSolution(const scalar step);

            //- Adjust step to satisfy cireteria
            void adjustStep
            (
                scalar& step,
                const scalar value,
                const scalar update
            );


        // Functions related to the computation of the design variables update

            //- Update the constraint subsets
            void updateViolatedConstraintsSubsets();

            //- Update the constraint subsets
            void updateNullAndRangeSpaceSubsets();

            //- Update violated constraints indices (iTilda and iTildaEps)
            void updateViolatedIndices
            (
                const label i,
                const scalarField& constraints
            );

            //- Update constraint indices related to the null and range space
            //- part of the correction
            void updateCorrectionIndices
            (
                const label i,
                const scalarField& LagrangeMults,
                const scalarField& dual
            );

            //- Compute the update of the design variables as a combination of
            //- null space and range space parts
            void correction();

            //- A & v
            //  where A is the Jacobian of all constraints, identified by
            //  the labelList addressings, w.r.t. the active design variables
            tmp<scalarField> Av
            (
                const scalarField& v,
                const labelListList& subsets
            );

            //- A.T & v
            //  where A is the Jacobian of all constraints, identified by
            //  the labelList addressings, w.r.t. the active design variables
            tmp<scalarField> ATv
            (
                const scalarField& v,
                const labelListList& subsets
            );

            //- Collect all constraint values corresponding to the provided
            //- indices
            tmp<scalarField> activeConstraints
            (
                const labelListList& subsets
            );

            //- Computes the parts of ksiJ and ksiC that require a system
            //- solution
            //  Formulated in terms of the flow constraints, which are usually
            //  few
            tmp<scalarField> constraintRelatedUpdate
            (
                const scalarField& rhs,
                const labelListList& subsets
            );

            //- Print statistics on the number of flow- and bound-related
            //- constraints included in the subset
            void statistics
            (
                const labelListList& subset,
                const word& description
            );


private:

    // Private Member Functions

        //- No copy construct
        nullSpace(const nullSpace&) = delete;

        //- No copy assignment
        void operator=(const nullSpace&) = delete;


public:

    //- Runtime type information
    TypeName("nullSpace");


    // Constructors

        //- Construct from components
        nullSpace
        (
            const fvMesh& mesh,
            const dictionary& dict,
            autoPtr<designVariables>& designVars,
            const label nConstraints,
            const word& type
        );


    //- Destructor
    virtual ~nullSpace() = default;


    // Member Functions

        //- Compute design variables correction
        virtual void computeCorrection();

        //- Compute the merit function for line search
        virtual scalar computeMeritFunction();

        //- Directional derivative of the merit function, in the direction of
        //- the correction
        virtual scalar meritFunctionDirectionalDerivative();

        //- Write useful quantities to files
        virtual bool writeData(Ostream& os) const;
};


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

} // End namespace Foam

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

#endif

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