plastic_convergence_atom.cpp

Go to the documentation of this file.

/** @file plastic_convergence_atom.cpp
 * @brief Check scaled joint gradient and hardening-sensitive slack acceptance.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <PlasticIncrementalOptimization.hpp>
using namespace EshelbianPlasticity;
 
namespace {
 
/** Give each scaled gradient block the same prescribed Euclidean norm. */
MoFEMErrorCode setGradientBlocks(PlasticIncrementalOptimizationProblem &problem,
                                 Vec gradient, const double norm) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(gradient);
  for (const auto &[field, scale] : problem.getControlFieldScales()) {
    IS raw_is = nullptr;
    CHKERR DMMoFEMGetFieldIS(problem.getControlDM(), RowColData::ROW,
                             field.c_str(), &raw_is);
    SmartPetscObj<IS> field_is(raw_is);
    PetscInt size = 0;
    CHKERR ISGetSize(field_is, &size);
    if (!size)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Gradient fixture requires both control fields");
    Vec block = nullptr;
    CHKERR VecGetSubVector(gradient, field_is, &block);
    const auto error = VecSet(block, norm / (scale * std::sqrt(double(size))));
    CHKERR VecRestoreSubVector(gradient, field_is, &block);
    CHKERR error;
  }
  MoFEMFunctionReturn(0);
}
 
/** Only the leading Legendre coefficient carries a constant latent value. */
MoFEMErrorCode setSlack(EshelbianCore &ep, Vec control, const double slack) {
  MoFEMFunctionBegin;
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmIncrementalOptimization, &problem);
  PetscInt size = 0;
  CHKERR VecGetLocalSize(control, &size);
  PetscScalar *values = nullptr;
  CHKERR VecGetArray(control, &values);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscLocalDofIdx();
    if (dof->getName() == ep.plasticLatentField && index >= 0 && index < size)
      values[index] = dof->getEntDofIdx() == 0 ? std::log(slack) : 0.;
  }
  CHKERR VecRestoreArray(control, &values);
  CHKERR VecGhostUpdateBegin(control, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(control, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runConvergenceAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  auto ts = createTS(ep.mField.get_comm());
  CHKERR TSSetDM(ts, ep.dmElastic);
  auto state = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(state);
  auto context = createPlasticIncrementalOptimizationContext(ep, ts, state);
  auto &problem = getPlasticIncrementalOptimizationProblem(context);
  auto gradient = vectorDuplicate(context->candidateControl);
  constexpr double tolerance = 1.e-3;
  double norm = 0;
  std::string diagnostics;
  CHKERR setGradientBlocks(problem, gradient, 0.6 * tolerance);
  CHKERR problem.updateGradientDiagnostics(gradient, norm);
  if (std::abs(norm - std::sqrt(2.) * 0.6 * tolerance) > 1.e-12 * tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Joint gradient norm has incorrect scaling or MPI ownership");
  CHKERR problem.validateSolution(context->candidateControl, gradient,
                                  tolerance, diagnostics);
  // Each scaled block is below tolerance, but their joint norm exceeds it.
  CHKERR setGradientBlocks(problem, gradient, 0.8 * tolerance);
  CHKERR PetscPushErrorHandler(PetscReturnErrorHandler, nullptr);
  const auto rejection = problem.validateSolution(
      context->candidateControl, gradient, tolerance, diagnostics);
  CHKERR PetscPopErrorHandler();
  if (!rejection)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Unresolved joint gradient was accepted");
 
  double slack_tolerance = 0;
  CHKERR PetscOptionsGetReal(
      nullptr, nullptr, "-incremental_optimization_plastic_slack_tolerance",
      &slack_tolerance, nullptr);
  const double factor = std::max(1., problem.getIsotropicHardeningModulus() /
                                         problem.getInitialYieldStress());
  bool converged = false;
  CHKERR setSlack(ep, context->candidateControl,
                  0.5 * slack_tolerance / factor);
  CHKERR problem.assessProximalConvergence(context->candidateControl, converged,
                                           diagnostics);
  if (!converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Resolved physical slack was rejected");
  CHKERR setSlack(ep, context->candidateControl, 2. * slack_tolerance / factor);
  CHKERR problem.assessProximalConvergence(context->candidateControl, converged,
                                           diagnostics);
  if (converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Unresolved strain/hardening slack was accepted");
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Joint gradient and combined slack acceptance verified\n");
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] = "Verify the two plastic acceptance criteria.\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
  auto core_log = logging::core::get();
  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
  LogManager::setLog("EP");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
  LogManager::setLog("EPSELF");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
  LogManager::setLog("EPSYNC");
  try {
    CHKERR runAuxiliaryLogarithmicStressTest(runConvergenceAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}