IncrementalMechanics.cpp

Go to the documentation of this file.

/**
 * @file IncrementalMechanics.cpp
 * @brief Mechanical equilibrium shared by the incremental-optimisation physics
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <IncrementalMechanics.hpp>
#include <PlasticIncrementalOptimization.hpp>
#include <impl/PlasticIncrementalOptimizationInternal.hpp>
 
namespace EshelbianPlasticity {
 
using namespace PlasticIncrementalOptimizationInternal;
 
IncrementalMechanics::IncrementalMechanics(EshelbianCore &ep,
                                           SmartPetscObj<TS> ts)
    : ep(ep), timeSolver(std::move(ts)) {}
 
MoFEMErrorCode IncrementalMechanics::publishTrialControl() {
  MoFEMFunctionBegin;
  if (!ep.dmIncrementalOptimization || !ep.incrementalTrialControl)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Incremental-optimisation control storage is unavailable");
  CHKERR VecGhostUpdateBegin(ep.incrementalTrialControl, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ep.incrementalTrialControl, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmIncrementalOptimization,
                                 ep.incrementalTrialControl, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::ROW);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::stateControlBlocksEqual(Vec first,
                                                             Vec second,
                                                             PetscBool &equal) {
  MoFEMFunctionBegin;
  equal = PETSC_TRUE;
  for (const auto &field : stateControlFields) {
    IS raw_is = nullptr;
    CHKERR DMMoFEMGetFieldIS(ep.dmIncrementalOptimization, RowColData::ROW,
                             field.c_str(), &raw_is);
    SmartPetscObj<IS> field_is(raw_is);
    Vec first_block = nullptr, second_block = nullptr;
    CHKERR VecGetSubVector(first, field_is, &first_block);
    CHKERR VecGetSubVector(second, field_is, &second_block);
    PetscBool block_equal = PETSC_FALSE;
    CHKERR VecEqual(first_block, second_block, &block_equal);
    CHKERR VecRestoreSubVector(second, field_is, &second_block);
    CHKERR VecRestoreSubVector(first, field_is, &first_block);
    if (!block_equal) {
      equal = PETSC_FALSE;
      break;
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::zeroStateControlBlocks(Vec control) {
  MoFEMFunctionBegin;
  for (const auto &field : stateControlFields) {
    IS raw_is = nullptr;
    CHKERR DMMoFEMGetFieldIS(ep.dmIncrementalOptimization, RowColData::ROW,
                             field.c_str(), &raw_is);
    SmartPetscObj<IS> field_is(raw_is);
    Vec block = nullptr;
    CHKERR VecGetSubVector(control, field_is, &block);
    CHKERR VecZeroEntries(block);
    CHKERR VecRestoreSubVector(control, field_is, &block);
  }
  CHKERR VecGhostUpdateBegin(control, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(control, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
IncrementalMechanics::snapshotHistory(std::vector<FieldSnapshot> &snapshot) {
  MoFEMFunctionBegin;
  snapshot.clear();
  snapshot.resize(historyFields.size());
  for (std::size_t f = 0; f != historyFields.size(); ++f) {
    auto &field_snapshot = snapshot[f];
    auto save = [&](boost::shared_ptr<FieldEntity> field_entity) {
      MoFEMFunctionBegin;
      auto data = field_entity->getEntFieldData();
      field_snapshot.emplace(field_entity->getEnt(),
                             std::vector<double>(data.begin(), data.end()));
      MoFEMFunctionReturn(0);
    };
    CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
        save, historyFields[f]);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
IncrementalMechanics::solveState(Vec state, EquilibriumSolveReport &report) {
  MoFEMFunctionBegin;
  CHKERR publishTrialControl();
  std::vector<FieldSnapshot> history;
  CHKERR snapshotHistory(history);
  PetscBool same_control = PETSC_FALSE;
  if (cacheValid)
    CHKERR stateControlBlocksEqual(ep.incrementalTrialControl, cachedControl,
                                   same_control);
  int local_changed = !cacheValid || !same_control ||
                      cachedPhysicalTime != ep.currentPhysicalTime ||
                      history != cachedHistory;
  int changed = 0;
  CHKERR MPI_Allreduce(&local_changed, &changed, 1, MPI_INT, MPI_MAX,
                       ep.mField.get_comm());
  if (!changed) {
    CHKERR VecCopy(cachedState, state);
    CHKERR setStateOnMesh(ep, state);
    report = cachedReport;
    report.mechanicalSolveCount = 0;
    MoFEMFunctionReturnHot(0);
  }
  // Only converged states are cached. A rejected trial keeps the previous
  // cache usable while the optimiser rolls the control back.
  gradientValid = false;
  ++mechanicalSolves;
  CHKERR solveEquilibriumStateTS(ep, timeSolver, state, report);
  report.mechanicalSolveCount = 1;
  if (report.status == EquilibriumSolveStatus::Converged) {
    if (!cachedControl)
      cachedControl = vectorDuplicate(ep.incrementalTrialControl);
    if (!cachedState)
      cachedState = vectorDuplicate(state);
    CHKERR VecCopy(ep.incrementalTrialControl, cachedControl);
    CHKERR VecCopy(state, cachedState);
    cachedHistory = std::move(history);
    cachedPhysicalTime = ep.currentPhysicalTime;
    cachedReport = report;
    cacheValid = true;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::fixedDataAreCurrent(bool &current) {
  MoFEMFunctionBegin;
  current = true;
  if (!cacheValid)
    MoFEMFunctionReturnHot(0);
  std::vector<FieldSnapshot> history;
  CHKERR snapshotHistory(history);
  int local_changed =
      cachedPhysicalTime != ep.currentPhysicalTime || history != cachedHistory;
  int changed = 0;
  CHKERR MPI_Allreduce(&local_changed, &changed, 1, MPI_INT, MPI_MAX,
                       ep.mField.get_comm());
  current = !changed;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::restoreState(Vec state) {
  MoFEMFunctionBegin;
  CHKERR publishTrialControl();
  CHKERR setStateOnMesh(ep, state);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::evaluateBaseline(Vec state,
                                                      double &value) {
  MoFEMFunctionBegin;
  // The baseline is the equilibrium without any increment. Zero the control
  // blocks that enter equilibrium, solve, and put the trial control back so
  // the current trial is not lost.
  auto trial_before = vectorDuplicate(ep.incrementalTrialControl);
  CHKERR VecCopy(ep.incrementalTrialControl, trial_before);
  CHKERR zeroStateControlBlocks(ep.incrementalTrialControl);
  EquilibriumSolveReport report;
  const MoFEMErrorCode solve_error =
      solveEquilibriumStateTS(ep, timeSolver, state, report);
  MoFEMErrorCode value_error = 0;
  if (!solve_error && report.status == EquilibriumSolveStatus::Converged)
    value_error = evaluateEquilibratedMechanicalValue(ep, timeSolver, value);
  CHKERR VecCopy(trial_before, ep.incrementalTrialControl);
  CHKERR VecGhostUpdateBegin(ep.incrementalTrialControl, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ep.incrementalTrialControl, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR solve_error;
  if (report.status != EquilibriumSolveStatus::Converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
            "Zero-increment baseline did not equilibrate (residual %g)",
            report.equilibriumNorm);
  CHKERR value_error;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::evaluateConservativeValue(double &value) {
  return evaluateEquilibratedMechanicalValue(ep, timeSolver, value);
}
 
MoFEMErrorCode
IncrementalMechanics::refreshConservativeGradient(DM control_dm) {
  MoFEMFunctionBegin;
  if (!cacheValid || !cachedState)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "The conservative gradient needs an equilibrated state");
  if (!cachedGradient)
    cachedGradient = createDMVector(control_dm);
  CHKERR VecZeroEntries(cachedGradient);
  if (supportsHessianAction()) {
    // Keep the base adjoint; prepareHessian() freezes it for the shell.
    if (!gradientWorkspace)
      gradientWorkspace = createReducedBulkMaterialEnergyHessianWorkspace();
    CHKERR prepareReducedBulkMaterialEnergyHessian(
        ep, timeSolver, control_dm, ep.incrementalTrialControl, cachedState,
        *gradientWorkspace);
    CHKERR assemblePlasticAdjointContribution(
        ep, control_dm,
        getReducedBulkMaterialEnergyBaseAdjoint(*gradientWorkspace),
        cachedGradient);
  } else {
    CHKERR assembleReducedBulkMaterialEnergyGradient(ep, timeSolver, control_dm,
                                                     cachedGradient);
  }
  gradientValid = true;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
IncrementalMechanics::assembleConservativeGradient(DM control_dm,
                                                   Vec gradient) {
  MoFEMFunctionBegin;
  if (!gradientValid)
    CHKERR refreshConservativeGradient(control_dm);
  CHKERR VecCopy(cachedGradient, gradient);
  MoFEMFunctionReturn(0);
}
 
bool IncrementalMechanics::supportsHessianAction() const {
  return ep.physicalEquations &&
         ep.physicalEquations->providesHelmholtzCurvature();
}
 
MoFEMErrorCode IncrementalMechanics::prepareHessian(DM control_dm,
                                                    Vec control) {
  MoFEMFunctionBegin;
  if (!supportsHessianAction())
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "The selected material does not provide the curvature of the "
            "reduced mechanical energy");
  if (!gradientValid)
    CHKERR refreshConservativeGradient(control_dm);
  if (!hessianWorkspace)
    hessianWorkspace = createReducedBulkMaterialEnergyHessianWorkspace();
  CHKERR freezeReducedBulkMaterialEnergyHessian(*gradientWorkspace, control,
                                                *hessianWorkspace);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::applyHessian(Vec direction, Vec action) {
  MoFEMFunctionBegin;
  if (!hessianWorkspace)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Mechanical Hessian action requested before prepareHessian()");
  CHKERR assembleReducedBulkMaterialEnergyHessianAction(
      ep, timeSolver, *hessianWorkspace, direction, action);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode IncrementalMechanics::assembleHessianPreconditioner(Mat metric) {
  MoFEMFunctionBegin;
  if (!hessianWorkspace)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Mechanical metric block requested before prepareHessian()");
  CHKERR assembleReducedBulkMaterialEnergyFlowPreconditioner(
      ep, *hessianWorkspace, metric, false);
  MoFEMFunctionReturn(0);
}
 
PetscInt IncrementalMechanics::getDirectionalSolveCount() const {
  return hessianWorkspace ? getReducedBulkMaterialEnergyDirectionalSolveCount(
                                *hessianWorkspace)
                          : 0;
}
 
} // namespace EshelbianPlasticity