IncrementalOptimization.hpp

Go to the documentation of this file.

/**
 * @file IncrementalOptimization.hpp
 * @brief Physics-independent incremental optimization API
 *
 * The orchestration in this file knows how to minimize an incremental
 * potential whose state is obtained from a nested equilibrium solve.  A
 * concrete dissipation process supplies the control layout, conservative
 * value and gradient, dissipation, acceptance checks, and commit operation
 * through IncrementalOptimizationProblem.
 */
 
#ifndef MOFEM_INCREMENTAL_OPTIMIZATION_HPP
#define MOFEM_INCREMENTAL_OPTIMIZATION_HPP
 
#include <limits>
 
namespace EshelbianPlasticity {
 
enum class EquilibriumSolveStatus { NotEvaluated, Converged, Diverged, Failed };
 
/** Recoverable result of a fixed-control state solve. */
struct EquilibriumSolveReport {
  EquilibriumSolveStatus status = EquilibriumSolveStatus::NotEvaluated;
  MoFEMErrorCode errorCode = 0;
  double equilibriumNorm = std::numeric_limits<double>::max();
};
 
/** Result of one rollback-safe objective evaluation. */
struct IncrementalObjectiveEvaluation {
  double conservativeValue = 0;
  double dissipativeValue = 0;
  /** Null until the corresponding gradient has been requested. */
  SmartPetscObj<Vec> smoothGradient;
  SmartPetscObj<Vec> objectiveGradient;
  double equilibriumNorm = std::numeric_limits<double>::max();
  double smoothGradientNorm = std::numeric_limits<double>::max();
  double objectiveGradientNorm = std::numeric_limits<double>::max();
  EquilibriumSolveStatus stateStatus = EquilibriumSolveStatus::NotEvaluated;
  MoFEMErrorCode stateSolveError = 0;
  PetscInt stateSolveCount = 0;
  /** True when the state and reduced-objective data were reused for the same
   * control. */
  bool cacheHit = false;
};
 
using EquilibriumSolveFunction =
    std::function<MoFEMErrorCode(Vec, EquilibriumSolveReport &)>;
 
struct IncrementalOptimizationContext;
 
/**
 * Physics contract consumed by the generic incremental optimizer.
 *
 * Plasticity, cohesive fracture, and contact/friction implement this
 * interface without adding their variables or constitutive assumptions to
 * IncrementalOptimizationContext.
 */
struct IncrementalOptimizationProblem {
  virtual ~IncrementalOptimizationProblem() = default;
 
  virtual std::string getName() const = 0;
  virtual MoFEMErrorCode validate() = 0;
 
  /** Parent DM containing every optimization-control field. */
  virtual DM getControlDM() const = 0;
 
  /** Create the problem's committed and candidate control vectors. */
  virtual MoFEMErrorCode createControlData(
      SmartPetscObj<Vec> &committed_control,
      SmartPetscObj<Vec> &candidate_control) = 0;
 
  /** Add physics-specific rows, columns, and elements to constraintDM. */
  virtual MoFEMErrorCode configureConstraintDM(DM constraint_dm) = 0;
 
  /** Set initial multipliers using constraint-DM row numbering. */
  virtual MoFEMErrorCode
  initialiseInequalityMultipliers(DM constraint_dm, Vec multipliers) = 0;
 
  /** Define the zero-increment/reference point for a new load step. */
  virtual MoFEMErrorCode initialiseReferenceControl(
      Vec committed_control, Vec reference_control) = 0;
 
  virtual MoFEMErrorCode solveState(Vec state,
                                    EquilibriumSolveReport &report) = 0;
  virtual MoFEMErrorCode restoreState(Vec state) = 0;
  virtual MoFEMErrorCode evaluateConservativeValue(double &value) = 0;
 
  /** Assemble the reduced conservative gradient. */
  virtual MoFEMErrorCode
  assembleConservativeGradient(Vec smooth_gradient) = 0;
  virtual MoFEMErrorCode evaluateDissipation(Vec control, double &value) = 0;
  virtual MoFEMErrorCode assembleObjectiveGradient(
      Vec control, Vec smooth_gradient, Vec objective_gradient) = 0;
  virtual MoFEMErrorCode evaluateInequalityConstraints(Vec control,
                                                        Vec constraints) = 0;
  virtual MoFEMErrorCode assembleInequalityJacobian(
      Vec control, Vec smooth_gradient, Mat jacobian) = 0;
 
  /** Assemble a finite generalized row without a usable smooth gradient. */
  virtual MoFEMErrorCode
  assembleInequalityJacobianWithoutSmoothGradient(Vec control,
                                                   Mat jacobian) = 0;
 
  /** Physics-specific primal/dual checks using both DM layouts. */
  virtual MoFEMErrorCode validateSolution(
      DM constraint_dm, Vec control, Vec smooth_gradient,
      Vec inequality_multipliers, PetscReal gradient_tolerance,
      PetscReal constraint_tolerance, std::string &diagnostics) = 0;
 
  /**
   * Commit physics history and all next-step control vectors atomically.
   * `solution` and `equilibrated_state` are read-only accepted inputs; a
   * successful implementation may modify only physics history and the three
   * control outputs.
   */
  virtual MoFEMErrorCode commit(Vec solution, Vec equilibrated_state,
                                Vec committed_control, Vec candidate_control,
                                Vec reference_control) = 0;
};
 
/**
 * @brief Generic transaction/cache data; contains no constitutive state.
 *
 * Terminology used below:
 *
 * - **Control**: variables selected by TAO.  For plasticity its tensor part is
 *   `Delta H^p`, the increment of `H^p = log(F^p)`, which is the logarithmic
 *   plastic strain in the symmetric, zero-plastic-spin model.  Its scalar part
 *   is the epigraph coordinate `Delta kappa`; committed `kappa_n` remains in
 *   physics-owned mesh storage sharing the same field layout.
 * - **Mechanical state**: equilibrium unknowns, such as displacement, solved
 *   at fixed control by the nested TS/SNES solve.
 * - **Trial**: one control proposed by TAO and evaluated without changing
 *   committed material history.
 * - **Rejected trial**: a trial for which the nested state solve fails or its
 *   residual exceeds `equilibriumTolerance`; its changes are discarded.
 * - **Committed data**: the accepted rollback point and physics-owned material
 *   history; trial evaluations must not modify them.
 * - **Rollback**: restore the accepted control and last converged mechanical
 *   state after rejecting a trial.
 * - **Reference**: the zero-increment control at the beginning of a load step.
 * - **Baseline**: the equilibrated state and conservative value at the
 *   reference control; the latter is subtracted to form the energy increment.
 * - **Last valid state**: the most recent equilibrated trial state, used as the
 *   starting point for the next nested solve and as its recovery state.
 * - **Cached evaluation**: results retained for an identical repeated control,
 *   avoiding another equilibrium solve.
 * - **Commit**: after final validation, atomically update material history and
 *   make the accepted state the starting point of the next load step.
 */
struct IncrementalOptimizationContext {
  /** Physics-specific callbacks used by the generic optimiser. */
  boost::shared_ptr<IncrementalOptimizationProblem> problem;
 
  /** Last accepted control; source for rollback after a rejected trial. */
  SmartPetscObj<Vec> committedControl;
  /** Current TAO trial control, synchronized before physics evaluation. */
  SmartPetscObj<Vec> candidateControl;
  /** Reduced conservative gradient evaluated at candidateControl. */
  SmartPetscObj<Vec> gradient;
  /** Complete objective gradient, including local incremental resistance. */
  SmartPetscObj<Vec> objectiveGradient;
  /** Rectangular problem mapping controls to inequality constraints. */
  SmartPetscObj<DM> constraintDM;
  /** Values of the TAO inequalities, using the convention c_i >= 0. */
  SmartPetscObj<Vec> inequalityConstraints;
  /** Jacobian of the inequality map with respect to the control. */
  SmartPetscObj<Mat> inequalityJacobian;
 
  /** Mechanical state equilibrated for candidateControl. */
  SmartPetscObj<Vec> equilibratedState;
  /** Control defining zero increment at the start of the load step. */
  SmartPetscObj<Vec> referenceControl;
  /** Equilibrated-state snapshot corresponding to referenceControl. */
  SmartPetscObj<Vec> baselineState;
  /** Latest converged state, used to warm-start or restore trial solves. */
  SmartPetscObj<Vec> lastValidState;
  /** Control associated with cachedState and cachedEvaluation. */
  SmartPetscObj<Vec> cachedControl;
  /** Equilibrated state associated with cachedControl. */
  SmartPetscObj<Vec> cachedState;
  /** Conservative value at referenceControl, subtracted from the objective. */
  double baselineConservativeValue = 0;
  /** Values and optional conservative gradient for cachedControl/cachedState.
   * The complete objective gradient is a mutable TAO workspace, not cached. */
  IncrementalObjectiveEvaluation cachedEvaluation;
 
  /** Fixed-control state solve callback. */
  EquilibriumSolveFunction stateSolve;
  /** Number of nested equilibrium solves performed in this context. */
  PetscInt stateSolveCount = 0;
  /** Number of TAO trials rejected because equilibrium was not obtained. */
  PetscInt rejectedTrialCount = 0;
  /** Accumulated wall time spent in nested equilibrium solves. */
  double equilibriumSeconds = 0;
  /** Accumulated wall time spent assembling conservative gradients. */
  double gradientAssemblySeconds = 0;
  /** Maximum accepted residual norm from a nested equilibrium solve. */
  double equilibriumTolerance = 1e-8;
  /** Whether the reference objective value and state are initialized. */
  bool baselineValid = false;
  /** Whether cachedControl exactly identifies a reusable evaluation. */
  bool cacheValid = false;
 
  MoFEMErrorCode rollbackTrial();
};
 
boost::shared_ptr<IncrementalOptimizationContext>
createIncrementalOptimizationContext(
    boost::shared_ptr<IncrementalOptimizationProblem> problem,
    SmartPetscObj<Vec> state);
 
/** Equilibrate and evaluate values without assembling an adjoint. */
MoFEMErrorCode evaluateIncrementalObjective(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, IncrementalObjectiveEvaluation &evaluation);
 
/** Reuse the value transaction and assemble its conservative gradient. */
MoFEMErrorCode evaluateIncrementalObjectiveAndGradient(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, IncrementalObjectiveEvaluation &evaluation);
 
/** Objective/gradient transaction used by TAO and directional tests.
 * A null objective_gradient requests only the objective value. */
MoFEMErrorCode evaluateIncrementalOptimizationTAOObjectiveAndGradient(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, PetscReal &objective, Vec objective_gradient,
    IncrementalObjectiveEvaluation &evaluation);
 
MoFEMErrorCode evaluateIncrementalOptimizationInequalityConstraints(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, Vec constraints);
 
/** Assemble a constraint Jacobian from the matching objective transaction. */
MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobianFromEvaluation(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, const IncrementalObjectiveEvaluation &evaluation,
    Mat jacobian);
 
MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobian(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec control, Mat jacobian);
 
/** Run the generic outer TAO solve and ask the problem to commit its state. */
MoFEMErrorCode solveIncrementalOptimizationTAO(
    const boost::shared_ptr<IncrementalOptimizationContext> &context,
    Vec state);
 
} // namespace EshelbianPlasticity
 
#endif // MOFEM_INCREMENTAL_OPTIMIZATION_HPP