TSElasticPreAndPostStep.cpp

Go to the documentation of this file.

/** @file TSElasticPreAndPostStep.cpp
 * @brief Elastic TS pre/post-step callbacks and incremental kinematics reset.
 * @date 2023-05-13
 *
 * @license{This project is released under the MIT License.}
 *
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/lambert_w.hpp>
 
#include <TSElasticPreAndPostStep.hpp>
#include <array>
#include <numbers>
 
namespace EshelbianPlasticity {
 
namespace {
 
/** Shorten only the constant L2 coefficient block; higher modes are untouched.
 * Used in the elastic TS pre-step, outside nonlinear/adjoint evaluations.
 */
struct OpShortenConstantRotation : VolUserDataOperator {
  OpShortenConstantRotation(const std::string &field, Vec solution,
                            PetscInt &count)
      : VolUserDataOperator(field, OPROW), solutionVec(solution),
        changedCount(count) {
    doEntities.fill(false);
    doEntities[MBTET] = true;
  }
  MoFEMErrorCode doWork(int, EntityType, EntData &data) override;
  Vec solutionVec;
  PetscInt &changedCount;
};
 
MoFEMErrorCode OpShortenConstantRotation::doWork(int, EntityType,
                                                 EntData &data) {
  MoFEMFunctionBegin;
  if (!data.getIndices().size())
    MoFEMFunctionReturnHot(0);
  if (data.getIndices().size() < 3 || data.getFieldData().size() < 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Incomplete constant rotation block");
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(solutionVec, &first, &last);
  const auto *indices = data.getIndices().data().data();
  if (indices[0] < first || indices[0] >= last)
    MoFEMFunctionReturnHot(0);
  if (indices[1] < first || indices[2] >= last)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Constant rotation block crosses PETSc ownership ranges");
  const double phi0 = data.getN()(0, 0);
  auto coefficient = data.getFTensor1FieldData<3>();
  const double angle = std::abs(phi0) * coefficient.l2();
  if (!std::isfinite(angle))
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Nonfinite constant rotation");
  constexpr double pi = std::numbers::pi;
  // Stay away from both the pi branch switch and the 2*pi rank loss.
  if (angle < 1.5 * pi)
    MoFEMFunctionReturnHot(0);
  const double turns = std::round(angle / (2 * pi));
  const double factor = 1. - 2 * pi * turns / angle;
  std::array<double, 3> values;
  auto shortened = getFTensor1FromPtr<3>(values.data());
  FTENSOR_INDEX(3, i);
  shortened(i) = factor * coefficient(i);
  CHKERR VecSetValues(solutionVec, 3, indices, values.data(), INSERT_VALUES);
  ++changedCount;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode resetIncrementalKinematics(EshelbianCore *ep_ptr, TS ts) {
  MoFEMFunctionBegin;
 
  Vec solution;
  CHKERR TSGetSolution(ts, &solution);
  const auto shorten_rotation = [&]() {
    MoFEMFunctionBegin;
    // The product rotation is invariant under a constant-mode branch change.
    // Raw-coordinate mass penalties/viscosity are not invariant.
    if (ep_ptr->rotSelector != LARGE_ROT || ep_ptr->alphaR != 0. ||
        ep_ptr->alphaViscousR != 0.)
      MoFEMFunctionReturnHot(0);
    CHKERR VecGhostUpdateBegin(solution, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(solution, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(ep_ptr->dmElastic, solution, INSERT_VALUES,
                                   SCATTER_REVERSE);
    auto fe =
        boost::make_shared<VolumeElementForcesAndSourcesCore>(ep_ptr->mField);
    fe->getUserPolynomialBase() =
        boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
    fe->getRuleHook = [](int, int, int) { return 1; };
    CHKERR AddHOOps<3, 3, 3>::add(fe->getOpPtrVector(), {L2},
                                  ep_ptr->materialH1Positions,
                                  ep_ptr->frontAdjEdges);
    PetscInt local_count = 0, global_count;
    fe->getOpPtrVector().push_back(
        new OpShortenConstantRotation(ep_ptr->rotAxis, solution, local_count));
    CHKERR DMoFEMLoopFiniteElements(ep_ptr->dmElastic,
                                    ep_ptr->elementVolumeName, fe);
    CHKERR VecAssemblyBegin(solution);
    CHKERR VecAssemblyEnd(solution);
    CHKERR MPI_Allreduce(&local_count, &global_count, 1, MPIU_INT, MPI_SUM,
                         ep_ptr->mField.get_comm());
    if (global_count) {
      CHKERR VecGhostUpdateBegin(solution, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(solution, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToLocalVector(ep_ptr->dmElastic, solution, INSERT_VALUES,
                                     SCATTER_REVERSE);
      MOFEM_LOG("EP", Sev::inform) << "Shortened constant L2 rotation on "
                                   << global_count << " elements";
    }
    MoFEMFunctionReturn(0);
  };
 
  const auto zero_fields = [&]() {
    MoFEMFunctionBegin;
    auto zero_field = [&](const std::string &field_name) {
      MoFEMFunctionBegin;
      auto is_mng = ep_ptr->mField.getInterface<ISManager>();
      SmartPetscObj<IS> is;
      CHKERR is_mng->isCreateProblemFieldAndRankLocal(
          "ELASTIC_PROBLEM", ROW, field_name, 0, MAX_DOFS_ON_ENTITY, is);
      const int *index_ptr;
      CHKERR ISGetIndices(is, &index_ptr);
      int size;
      CHKERR ISGetLocalSize(is, &size);
      double *array;
      CHKERR VecGetArray(solution, &array);
      for (int ii = 0; ii != size; ++ii) {
        array[index_ptr[ii]] = 0;
      }
      CHKERR VecRestoreArray(solution, &array);
      CHKERR ISRestoreIndices(is, &index_ptr);
      MoFEMFunctionReturn(0);
    };
 
    MOFEM_LOG("EP", Sev::verbose)
        << "Resetting incremental stretch and rotation after H1 transfer";
    CHKERR zero_field(ep_ptr->stretchTensor);
    CHKERR zero_field(ep_ptr->rotAxis);
    CHKERR VecGhostUpdateBegin(solution, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(solution, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  switch (ep_ptr->gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR shorten_rotation();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
  case SMALL_ROT:
    CHKERR zero_fields();
    break;
  default:
    SETERRQ(
        ep_ptr->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
        "Unsupported gradient approximator for incremental kinematics reset");
  }
 
  MoFEMFunctionReturn(0);
}
 
constexpr char elasticPreStepKey[] = "ep_elastic_pre_step";
struct ElasticPreStepContext {
  EshelbianCore *ep;
  bool h1Update;
  bool reportedCoordinatePenalty = false;
};
 
/** Preserve optional H1 history bookkeeping before resetting coordinates. */
MoFEMErrorCode elasticPreStep(TS ts) {
  MoFEMFunctionBegin;
  PetscContainer container = nullptr;
  CHKERR PetscObjectQuery(getPetscObject(ts), elasticPreStepKey,
                          reinterpret_cast<PetscObject *>(&container));
  if (!container)
    SETERRQ(PetscObjectComm(getPetscObject(ts)), MOFEM_DATA_INCONSISTENCY,
            "Missing elastic TS pre-step context");
  void *pointer = nullptr;
  CHKERR PetscContainerGetPointer(container, &pointer);
  auto &context = *static_cast<ElasticPreStepContext *>(pointer);
  if (context.h1Update)
    CHKERR TSElasticPreAndPostStep::preStepFun(ts);
  if (context.ep->gradApproximator == NO_H1_CONFIGURATION) {
    if (context.ep->rotSelector == LARGE_ROT &&
        (context.ep->alphaR != 0. || context.ep->alphaViscousR != 0.) &&
        !context.reportedCoordinatePenalty) {
      MOFEM_LOG("EP", Sev::inform)
          << "Constant rotation shortening disabled: alpha_r or "
             "viscosity_alpha_r acts on raw rotation coordinates";
      context.reportedCoordinatePenalty = true;
    }
    CHKERR resetIncrementalKinematics(context.ep, ts);
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
MoFEMErrorCode TSElasticPreAndPostStep::setPreStep(EshelbianCore *ep_ptr, TS ts,
                                             bool h1_update) {
  MoFEMFunctionBegin;
  PetscContainer container = nullptr;
  CHKERR PetscObjectQuery(getPetscObject(ts), elasticPreStepKey,
                          reinterpret_cast<PetscObject *>(&container));
  if (!container) {
    CHKERR PetscContainerCreate(ep_ptr->mField.get_comm(), &container);
    auto *context = new ElasticPreStepContext{ep_ptr, h1_update};
    CHKERR PetscContainerSetPointer(container, context);
    CHKERR PetscContainerSetUserDestroy(container, [](void *pointer) {
      MoFEMFunctionBegin;
      delete static_cast<ElasticPreStepContext *>(pointer);
      MoFEMFunctionReturn(0);
    });
    CHKERR PetscObjectCompose(getPetscObject(ts), elasticPreStepKey,
                              getPetscObject(container));
    CHKERR PetscContainerDestroy(&container);
  } else {
    void *pointer = nullptr;
    CHKERR PetscContainerGetPointer(container, &pointer);
    auto &context = *static_cast<ElasticPreStepContext *>(pointer);
    context.ep = ep_ptr;
    context.h1Update = h1_update;
  }
  CHKERR TSSetPreStep(ts, elasticPreStep);
  MoFEMFunctionReturn(0);
}
 
TSElasticPreAndPostStep::ProjectionScope::ProjectionScope(EshelbianCore *ep_ptr)
    : comm(ep_ptr->mField.get_comm()) {
  // Callbacks share one active projection context, as in the existing solvers.
  CHK_THROW_MESSAGE(epPtr != nullptr ? MOFEM_DATA_INCONSISTENCY : 0,
                    "Elastic projection is already active");
  try {
    CHK_THROW_MESSAGE(postStepInitialise(ep_ptr),
                      "Initialise elastic spatial projection");
    active = true;
  } catch (...) {
    CHKERRABORT(comm, postStepDestroy());
    throw;
  }
}
 
TSElasticPreAndPostStep::ProjectionScope::~ProjectionScope() {
  CHKERRABORT(comm, close());
}
 
MoFEMErrorCode TSElasticPreAndPostStep::ProjectionScope::close() {
  MoFEMFunctionBegin;
  if (active) {
    CHKERR postStepDestroy();
    active = false;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
TSElasticPreAndPostStep::postStepInitialise(EshelbianCore *ep_ptr) {
  MoFEMFunctionBegin;
 
  epPtr = ep_ptr;
 
  // H1 displacement may be excluded from dmElastic; use the full parent DM.
  dmPrjSpatial = createDM(ep_ptr->mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, ep_ptr->dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, ep_ptr->spatialH1Disp);
  CHKERR DMMoFEMAddElement(dmPrjSpatial, ep_ptr->elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  auto create_post_step_ksp = [&]() {
    auto ksp = createKSP(epPtr->mField.get_comm());
 
    auto set_up = [&]() {
      MoFEMFunctionBegin;
      using DomainEle = VolumeElementForcesAndSourcesCore;
      using DomainEleOp = DomainEle::UserDataOperator;
      using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
          PETSC>::BiLinearForm<GAUSS>::OpMass<1, 3>;
      using OpRhs = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
          GAUSS>::OpBaseTimesVector<1, 3, 3>;
      auto fe_lhs = boost::make_shared<DomainEle>(ep_ptr->mField);
      auto fe_rhs = boost::make_shared<DomainEle>(ep_ptr->mField);
 
      fe_lhs->getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
          new CGGUserPolynomialBase(nullptr, true));
      fe_rhs->getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
          new CGGUserPolynomialBase(nullptr, true));
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          fe_lhs->getOpPtrVector(), {L2}, ep_ptr->materialH1Positions,
          ep_ptr->frontAdjEdges);
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          fe_rhs->getOpPtrVector(), {L2}, ep_ptr->materialH1Positions,
          ep_ptr->frontAdjEdges);
 
      fe_lhs->getOpPtrVector().push_back(
          new OpDomainMass(ep_ptr->spatialH1Disp, ep_ptr->spatialH1Disp));
      auto w_ptr = boost::make_shared<MatrixDouble>();
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>(ep_ptr->spatialL2Disp, w_ptr));
      fe_rhs->getOpPtrVector().push_back(
          new OpRhs(ep_ptr->spatialH1Disp, w_ptr));
 
      CHKERR DMMoFEMKSPSetComputeOperators(
          dmPrjSpatial, ep_ptr->elementVolumeName, fe_lhs, nullptr, nullptr);
      CHKERR DMMoFEMKSPSetComputeRHS(dmPrjSpatial, ep_ptr->elementVolumeName,
                                     fe_rhs, nullptr, nullptr);
 
      CHKERR KSPAppendOptionsPrefix(ksp, "prjspatial_");
      CHKERR KSPSetFromOptions(ksp);
      CHKERR KSPSetDM(ksp, dmPrjSpatial);
      CHKERR KSPSetUp(ksp);
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(set_up(), "set up");
 
    return ksp;
  };
 
  prjKsp = create_post_step_ksp();
  prjD = createDMVector(dmPrjSpatial);
  prjF = vectorDuplicate(prjD);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode TSElasticPreAndPostStep::postStepDestroy() {
  MoFEMFunctionBegin;
  prjKsp.reset();
  prjD.reset();
  prjF.reset();
  if (dmPrjSpatial) {
    // Break DM/pipeline reference cycles while the parent problem is alive.
    CHKERR DMMoFEMClearDMCtx(dmPrjSpatial);
    dmPrjSpatial.reset();
  }
  epPtr = nullptr;
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode TSElasticPreAndPostStep::preStepFun(TS ts) {
  MoFEMFunctionBegin;
  MOFEM_LOG("EP", Sev::inform) << "Pre step";
 
  double time;
  CHKERR TSGetTime(ts, &time);
 
  if (EshelbianCore::crackingOn) {
 
    auto debug_crack = [&]() {
      MoFEMFunctionBegin;
 
      PetscBool debug_crack_mesh = PETSC_FALSE;
      CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_crack_mesh",
                                 &debug_crack_mesh, PETSC_NULLPTR);
      if (debug_crack_mesh) {
 
        auto append_meshsets_from_block = [](MoFEM::Interface &m_field,
                                             const std::string &block_name,
                                             const int dim,
                                             std::vector<EntityHandle> &r) {
          MoFEMFunctionBegin;
          auto mesh_mng = m_field.getInterface<MeshsetsManager>();
          auto bcs = mesh_mng->getCubitMeshsetPtr(
 
              std::regex((boost::format("%s(.*)") % block_name).str())
 
          );
 
          for (auto bc : bcs) {
            Range ents;
            CHKERR m_field.get_moab().get_entities_by_dimension(
                bc->getMeshset(), dim, ents, true);
            if (!ents.empty()) {
              r.push_back(bc->getMeshset());
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        if (!epPtr->mField.get_comm_rank() &&
            !EshelbianCore::noCrackExtension) {
          if (epPtr->maxMovedFaces) {
            CHKERR epPtr->mField.getInterface<MeshsetsManager>()
                ->addEntitiesToMeshset(BLOCKSET, epPtr->addCrackMeshsetId,
                                       *epPtr->maxMovedFaces);
          }
 
          auto meshset_ptr = get_temp_meshset_ptr(epPtr->mField.get_moab());
          Range tets;
          CHKERR epPtr->mField.get_moab().get_entities_by_dimension(
              *meshset_ptr, 3, tets);
          CHKERR epPtr->mField.get_moab().add_entities(*meshset_ptr, tets);
 
          std::vector<EntityHandle> meshsets;
          meshsets.push_back(*meshset_ptr);
 
          CHKERR append_meshsets_from_block(epPtr->mField, "CRACK", 2,
                                            meshsets);
          CHKERR append_meshsets_from_block(epPtr->mField, "FRONT", 1,
                                            meshsets);
          CHKERR append_meshsets_from_block(epPtr->mField, "EDGE", 1,
                                            meshsets);
 
          int time_step = 0;
          CHKERR TSGetStepNumber(ts, &time_step);
          std::string file_name =
              "crack_meshsets_" + std::to_string(time_step) + ".h5m";
          CHKERR epPtr->mField.get_moab().write_file(file_name.c_str(), "MOAB",
                                                     nullptr, meshsets.data(),
                                                     static_cast<int>(meshsets.size()));
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
      if (EshelbianCore::crackingStartTime < time) {
        CHKERR epPtr->setNewFrontCoordinates();
        CHKERR debug_crack();
        CHKERR epPtr->addCrackSurfaces();
      }
      CHKERR epPtr->projectGeometry(0, time);
      CHKERR epPtr->projectMaterialTags();
    }
 
  Vec T;
  CHKERR TSGetSolution(ts, &T);
 
  CHKERR VecCopy(T, epPtr->solTSStep);
  CHKERR VecGhostUpdateBegin(epPtr->solTSStep, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(epPtr->solTSStep, INSERT_VALUES, SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode TSElasticPreAndPostStep::postStepFun(TS ts) {
  MoFEMFunctionBegin;
 
  double time;
  CHKERR TSGetTime(ts, &time);
 
  MOFEM_LOG("EP", Sev::inform) << "Solve H1 post-step";
  CHKERR VecZeroEntries(prjF);
  CHKERR KSPSolve(prjKsp, prjF, prjD);
  CHKERR VecGhostUpdateBegin(prjD, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(prjD, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dmPrjSpatial, prjD, INSERT_VALUES,
                                 SCATTER_REVERSE);
 
  // The projected H1 field now carries the accepted deformation. Reset the
  // incremental stretch and rotation before crack calculations and TS
  // monitors so that the accepted increment is not applied a second time.
  CHKERR resetIncrementalKinematics(epPtr, ts);
 
  if (EshelbianCore::crackingOn) {
    switch (EshelbianCore::energyReleaseSelector) {
    case GRIFFITH_FORCE:
    case GRIFFITH_SKELETON:
      MOFEM_LOG("EP", Sev::inform) << "Calculate Griffith force";
      CHKERR epPtr->calculateFaceMaterialForce(1, ts);
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Energy release selector not handled");
    };
    if (time - EshelbianCore::crackingStartTime >
        std::numeric_limits<double>::epsilon()) {
      MOFEM_LOG("EP", Sev::inform) << "Calculate faces orientation";
 
      EshelbianCore::potentialCrackArrest = false;
      CHKERR epPtr->calculateOrientation(1, true);
    }
    CHKERR epPtr->calculateCrackArea(epPtr->currentCrackAreaPtr);
  }
 
  // Clear history after rebasing incremental kinematics onto the H1 field.
  // With no H1 configuration, retain it until the next pre-step copy.
  if (EshelbianCore::gradApproximator != NO_H1_CONFIGURATION)
    CHKERR VecZeroEntries(epPtr->solTSStep);
 
  PetscBool local_stop = PETSC_FALSE;
  PetscBool global_stop = PETSC_FALSE;
  switch (EshelbianCore::solverType) {
  case EshelbianCore::LoadFactor:
    if (!epPtr->mField.get_comm_rank() && EshelbianCore::crackingOn) {
      const double area_increase =
          *epPtr->currentCrackAreaPtr - epPtr->oldCrackArea;
      const bool crack_arrest_stops =
          EshelbianCore::potentialCrackArrest &&
          !EshelbianCore::propagateUnderCompression;
      if (area_increase > EshelbianCore::maxCrackExtension ||
          epPtr->avgGriffithsEnergy < epPtr->griffithEnergy ||
          crack_arrest_stops) {
        MOFEM_LOG_C("EP", Sev::warning,
                    "Stopping crack growth. Area increase: %3.12e, "
                    "average Griffiths energy: %3.12e, critical Griffiths "
                    "energy: %3.12e",
                    area_increase, epPtr->avgGriffithsEnergy,
                    epPtr->griffithEnergy);
        local_stop = PETSC_TRUE;
      }
    }
    MPI_Bcast(&local_stop, 1, MPIU_BOOL, 0, PETSC_COMM_WORLD);
    global_stop = local_stop;
 
    if (global_stop) {
      TSSetConvergedReason(ts, TS_CONVERGED_USER);
      EshelbianCore::noCrackExtension = true;
      EshelbianCore::nbStepsNoCrackExtension += 1;
    } else {
      EshelbianCore::noCrackExtension = false;
      EshelbianCore::nbStepsNoCrackExtension = 0;
    }
    break;
  default:
    break;
  }
 
  MoFEMFunctionReturn(0);
};
} // namespace EshelbianPlasticity