neohookean_direct_logarithmic_atom.cpp

Go to the documentation of this file.

/**
 * @file neohookean_direct_logarithmic_atom.cpp
 * @brief Direct six-coefficient constitutive law and coupled Jacobian checks.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
using namespace EshelbianPlasticity;
 
namespace {
 
using VolumeElement = VolumeElementForcesAndSourcesCore;
constexpr double c10 = 1.7;
constexpr double bulk = 8.5;
constexpr double deltaTime = .5;
 
SmartPetscObj<IS> getFieldIS(EshelbianCore &ep, const std::string &name) {
  IS indices;
  CHK_THROW_MESSAGE(
      DMMoFEMGetFieldIS(ep.dmElastic, RowColData::ROW, name.c_str(), &indices),
      "Get direct material test field indices");
  return SmartPetscObj<IS>(indices);
}
 
MoFEMErrorCode getFieldNorm(Vec vector, IS field, PetscReal &norm) {
  MoFEMFunctionBegin;
  Vec subvector;
  CHKERR VecGetSubVector(vector, field, &subvector);
  CHKERR VecNorm(subvector, NORM_INFINITY, &norm);
  CHKERR VecRestoreSubVector(vector, field, &subvector);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setState(EshelbianCore &ep, Vec state) {
  MoFEMFunctionBegin;
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  const std::array<double, 6> h{{.12, .035, -.021, -.045, .027, .018}};
  const std::array<double, 3> rotation{{.14, -.11, .07}};
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (index < first || index >= last)
      continue;
    double value = .01 * std::sin(.37 * (index + 1));
    if (dof->getName() == ep.stretchTensor) {
      if (dof->getDofOrder() == 0)
        value = h.at(dof->getDofCoeffIdx());
      else
        value *= .2 / (1 + dof->getDofOrder());
    } else if (dof->getName() == ep.rotAxis && dof->getDofOrder() == 0) {
      value = rotation.at(dof->getDofCoeffIdx());
    }
    CHKERR VecSetValue(state, index, value, INSERT_VALUES);
  }
  CHKERR VecAssemblyBegin(state);
  CHKERR VecAssemblyEnd(state);
  CHKERR VecGhostUpdateBegin(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
struct Audit {
  double noncoaxial = 0.;
  double rateGradient = 0.;
  double energyError = 0.;
  PetscInt points = 0;
};
 
struct OpAuditState : VolumeElement::UserDataOperator {
  OpAuditState(boost::shared_ptr<DataAtIntegrationPts> data,
               boost::shared_ptr<Audit> audit)
      : VolumeElement::UserDataOperator(NOSPACE, OPSPACE), dataPtr(data),
        auditPtr(audit) {}
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEXES(3, i, j, k);
    FTENSOR_INDEX(6, L);
    const int points = getGaussPts().size2();
    if (points <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Direct material atom requires multiple Gauss points");
    auto t_h = dataPtr->getFTensorLogStretch(points);
    auto t_p = dataPtr->getFTensorAdjointPdstretch(points);
    auto t_gradient = dataPtr->getFTensorGradLogStretchDot(points);
    FTensor::Tensor2_symmetric<double, 3> t_symmetric_p;
    FTensor::Tensor2<double, 3, 3> t_commutator;
    for (int gg = 0; gg != points; ++gg) {
      t_symmetric_p(i, j) = (t_p(i, j) || t_p(j, i)) / 2.;
      t_commutator(i, j) =
          t_h(i, k) * t_symmetric_p(k, j) - t_symmetric_p(i, k) * t_h(k, j);
      auditPtr->noncoaxial =
          std::max(auditPtr->noncoaxial,
                   std::sqrt(t_commutator(i, j) * t_commutator(i, j)));
      auditPtr->rateGradient =
          std::max(auditPtr->rateGradient,
                   std::sqrt(t_gradient(L, i) * t_gradient(L, i)));
      ++t_h;
      ++t_p;
      ++t_gradient;
    }
    auditPtr->points += points;
    MoFEMFunctionReturn(0);
  }
  boost::shared_ptr<DataAtIntegrationPts> dataPtr;
  boost::shared_ptr<Audit> auditPtr;
};
 
struct OpReferenceStress : VolumeElement::UserDataOperator {
  OpReferenceStress(boost::shared_ptr<DataAtIntegrationPts> data,
                    MatrixPtr stress, boost::shared_ptr<Audit> audit)
      : VolumeElement::UserDataOperator(NOSPACE, OPSPACE), dataPtr(data),
        stressPtr(stress), auditPtr(audit) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEXES(3, i, j, k);
    FTENSOR_INDEX(6, L);
    constexpr auto t_identity = FTensor::Kronecker_Delta<int>();
    const auto t_packed_basis = FTensor::SymmLTensor<double, 3>();
    const int points = getGaussPts().size2();
    auto t_h = dataPtr->getFTensorLogStretch(points);
    auto t_energy = getFTensor0FromVec(dataPtr->energyAtPts);
    auto t_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<6, -1, DL>, DL>::size(
            *stressPtr, points)();
    FTensor::Tensor2<double, 3, 3> t_argument, t_exponential, t_term, t_product;
    FTensor::Tensor2<double, 3, 3> t_material_stress;
    for (int gg = 0; gg != points; ++gg) {
      const double theta = t_h(i, i);
      t_argument(i, j) = 2. * t_h(i, j) - (2. * theta / 3.) * t_identity(i, j);
      const double argument_norm =
          std::sqrt(t_argument(i, j) * t_argument(i, j));
      if (!(argument_norm < .75))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Test state exceeds matrix exponential reference bound");
      // At ||2D|| < .75, 24 Taylor terms bound the omitted series by 1e-27.
      // This reference does not call the material kernel or its spectral map.
      t_exponential(i, j) = t_identity(i, j);
      t_term(i, j) = t_identity(i, j);
      for (int order = 1; order <= 24; ++order) {
        t_product(i, j) = t_term(i, k) * t_argument(k, j);
        t_term(i, j) = t_product(i, j) / order;
        t_exponential(i, j) += t_term(i, j);
      }
      const double trace_exponential = t_exponential(i, i);
      const double jacobian = std::exp(theta);
      t_material_stress(i, j) =
          2. * c10 *
              (t_exponential(i, j) -
               (trace_exponential / 3.) * t_identity(i, j)) +
          (bulk * jacobian * (jacobian - 1.)) * t_identity(i, j);
      t_stress(L) = t_packed_basis(i, j, L) * t_material_stress(i, j);
      const double energy = c10 * (trace_exponential - 3.) +
                            .5 * bulk * (jacobian - 1.) * (jacobian - 1.);
      const double error = std::abs(t_energy - energy);
      if (!std::isfinite(error))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Nonfinite direct material energy");
      auditPtr->energyError = std::max(auditPtr->energyError, error);
      ++t_h;
      ++t_energy;
      ++t_stress;
    }
    MoFEMFunctionReturn(0);
  }
  boost::shared_ptr<DataAtIntegrationPts> dataPtr;
  MatrixPtr stressPtr;
  boost::shared_ptr<Audit> auditPtr;
};
 
MoFEMErrorCode checkExternalStrain(EshelbianCore &ep, SmartPetscObj<Vec> state,
                                   SmartPetscObj<Vec> zero) {
  MoFEMFunctionBegin;
  Range entities;
  CHKERR ep.mField.get_moab().get_entities_by_type(0, MBTET, entities);
  auto blocks = boost::make_shared<ExternalStrainVec>();
  blocks->emplace_back("EXTERNALSTRAIN_FIRST", std::vector<double>{.07, 2.},
                       entities);
  blocks->emplace_back("EXTERNALSTRAIN_SECOND", std::vector<double>{-.01, 5.},
                       entities);
  blocks->emplace_back("EXTERNALSTRAIN_OTHER", std::vector<double>{1., 20.},
                       Range());
  const auto first_scale = [](double time) { return 1. + 2. * time; };
  const auto second_scale = [](double time) { return time - .3; };
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scaling{
      {"EXTERNALSTRAIN_FIRST",
       boost::make_shared<TimeScale>("", false, first_scale)},
      {"EXTERNALSTRAIN_SECOND",
       boost::make_shared<TimeScale>("", false, second_scale)}};
 
  auto production = boost::make_shared<VolumeElement>(ep.mField);
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, production);
  production->getOpPtrVector().push_back(
      ep.physicalEquations->returnOpSpatialPhysicalExternalStrain(
          ep.stretchTensor, ep.dataAtPts, blocks, scaling));
  auto reference = boost::make_shared<VolumeElement>(ep.mField);
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, reference);
  auto pressure = boost::make_shared<double>(0.);
  using OpReferenceLoad = FormsIntegrators<VolUserDataOperator>::Assembly<
      A>::LinearForm<GAUSS>::OpSource<1, 6>;
  reference->getOpPtrVector().push_back(
      new OpReferenceLoad(ep.stretchTensor, [pressure](double, double, double) {
        return FTensor::Tensor1<double, 6>{-*pressure, 0., 0.,
                                           -*pressure, 0., -*pressure};
      }));
  auto *tester = ep.mField.getInterface<OperatorsTester>();
  for (const double time : {.2, .7}) {
    const double load_time = ep.physicalTimeFlg ? ep.currentPhysicalTime : time;
    *pressure = 3. * (2. * .07 * first_scale(load_time) +
                      5. * -.01 * second_scale(load_time));
    auto assembled =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, production,
                            state, zero, zero, time, deltaTime, {});
    auto expected =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, reference,
                            state, zero, zero, time, deltaTime, {});
    PetscReal norm, error;
    CHKERR VecNorm(expected, NORM_INFINITY, &norm);
    CHKERR VecAXPY(assembled, -1., expected);
    CHKERR VecNorm(assembled, NORM_INFINITY, &error);
    if (!(norm > 0.) || !std::isfinite(error) ||
        error > 1.e-12 * std::max(1., norm))
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "External-strain residual error %g at time %g (norm %g)", error,
              time, norm);
    CHKERR PetscPrintf(ep.mField.get_comm(),
                       "Direct external strain: time %.2f, pressure %.6f, "
                       "assembled error %.3e\n",
                       time, *pressure, error);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runDirectAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  if (ep.physicalEquations->getFeatures().test(
          PhysicalEquations::AUXILIARY_LOGARITHMIC_STRESS))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Direct atom requires the ordinary six-coefficient field");
  auto state = createDMVector(ep.dmElastic);
  auto rates = createDMVector(ep.dmElastic);
  auto zero = createDMVector(ep.dmElastic);
  CHKERR setState(ep, state);
  // OperatorsTester takes increments and divides by deltaTime internally.
  CHKERR VecCopy(state, rates);
  CHKERR VecZeroEntries(zero);
  CHKERR VecCopy(state, ep.solTSStep);
  CHKERR VecGhostUpdateBegin(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
  boost::shared_ptr<VolumeElement> rhs, lhs;
  CHKERR ep.setContactElementRhsOps(ep.contactTreeRhs);
  CHKERR ep.setVolumeElementOps(1, true, false, rhs, lhs);
  auto audit = boost::make_shared<Audit>();
  rhs->getOpPtrVector().push_back(new OpAuditState(ep.dataAtPts, audit));
  auto *tester = ep.mField.getInterface<OperatorsTester>();
  auto jacobian = tester->assembleMat(ep.dmElastic, ep.elementVolumeName, lhs,
                                      state, rates, zero, 1., deltaTime, {});
  auto residual = tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs,
                                      state, rates, zero, 1., deltaTime, {});
  const std::vector<std::string> fields{ep.stretchTensor, ep.piolaStress,
                                        ep.bubbleField,   ep.rotAxis,
                                        ep.spatialL2Disp, ep.hybridSpatialDisp};
  std::map<std::string, SmartPetscObj<IS>> indices;
  for (const auto &field : fields) {
    auto field_is = getFieldIS(ep, field);
    PetscInt size;
    CHKERR ISGetSize(field_is, &size);
    if (!size)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Direct test field %s is empty", field.c_str());
    indices.emplace(field, field_is);
  }
  auto trial = vectorDuplicate(state);
  auto trial_rates = vectorDuplicate(state);
  auto direction = vectorDuplicate(state);
  auto mixed = vectorDuplicate(state);
  auto error = vectorDuplicate(residual);
  auto action = vectorDuplicate(residual);
  CHKERR VecZeroEntries(mixed);
  const OperatorsTester::VectorFunction evaluate = [&](Vec input, Vec output) {
    MoFEMFunctionBegin;
    CHKERR VecCopy(input, trial);
    CHKERR VecGhostUpdateBegin(trial, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(trial, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecCopy(trial, trial_rates);
    auto value =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs, trial,
                            trial_rates, zero, 1., deltaTime, {});
    CHKERR VecCopy(value, output);
    MoFEMFunctionReturn(0);
  };
  auto check_direction = [&](const std::string &name) {
    MoFEMFunctionBegin;
    CHKERR tester->checkVectorCentralFiniteDifference(
        state, direction, residual, jacobian, 2.e-6, evaluate, error);
    CHKERR MatMult(jacobian, direction, action);
    PetscReal maximum = 0.;
    for (const auto &field : fields) {
      PetscReal row_error, row_action;
      CHKERR getFieldNorm(error, indices.at(field), row_error);
      CHKERR getFieldNorm(action, indices.at(field), row_action);
      const double tolerance = 2.e-8 + 2.e-6 * row_action;
      if (!std::isfinite(row_error) || row_error > tolerance)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Direct Jacobian row %s, direction %s: error %g exceeds %g",
                field.c_str(), name.c_str(), row_error, tolerance);
      maximum = std::max(maximum, row_error);
    }
    CHKERR PetscPrintf(ep.mField.get_comm(),
                       "Direct Jacobian direction %s: error %.3e\n",
                       name.c_str(), maximum);
    MoFEMFunctionReturn(0);
  };
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  for (const auto &field : fields) {
    const int components = field == ep.stretchTensor ? 6 : 1;
    for (int component = 0; component != components; ++component) {
      CHKERR VecZeroEntries(direction);
      for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
        const auto index = dof->getPetscGlobalDofIdx();
        if (index < first || index >= last || dof->getName() != field ||
            (components == 6 && dof->getDofCoeffIdx() != component))
          continue;
        CHKERR VecSetValue(direction, index, std::sin(.71 * (index + 1)) + .3,
                           INSERT_VALUES);
      }
      CHKERR VecAssemblyBegin(direction);
      CHKERR VecAssemblyEnd(direction);
      PetscReal norm;
      CHKERR VecNormalize(direction, &norm);
      if (!(norm > 0.))
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Empty direct material direction");
      CHKERR VecGhostUpdateBegin(direction, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(direction, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecAXPY(mixed, 1., direction);
      CHKERR check_direction(field + ":" + std::to_string(component));
    }
  }
  CHKERR VecCopy(mixed, direction);
  CHKERR VecNormalize(direction, nullptr);
  CHKERR VecGhostUpdateBegin(direction, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(direction, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR check_direction("mixed");
 
  // Zero P and all rates isolate the direct material row of the production
  // residual. Retain the varying, non-diagonal logarithmic field.
  CHKERR VecCopy(state, trial);
  for (const auto &field : {ep.piolaStress, ep.bubbleField}) {
    Vec subvector;
    CHKERR VecGetSubVector(trial, indices.at(field), &subvector);
    CHKERR VecZeroEntries(subvector);
    CHKERR VecRestoreSubVector(trial, indices.at(field), &subvector);
  }
  CHKERR VecGhostUpdateBegin(trial, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(trial, INSERT_VALUES, SCATTER_FORWARD);
  auto direct = tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs,
                                    trial, zero, zero, 1., deltaTime, {});
  auto reference = boost::make_shared<VolumeElement>(ep.mField);
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, reference);
  reference->getOpPtrVector().push_back(
      ep.physicalEquations->returnOpCalculateHelmholtzFreeEnergy(ep.dataAtPts,
                                                                 nullptr));
  auto reference_stress = boost::make_shared<MatrixDouble>();
  reference->getOpPtrVector().push_back(
      new OpReferenceStress(ep.dataAtPts, reference_stress, audit));
  using OpReferenceResidual = FormsIntegrators<VolUserDataOperator>::Assembly<
      A>::LinearForm<GAUSS>::OpBaseTimesVector<1, 6, -1>;
  reference->getOpPtrVector().push_back(
      new OpReferenceResidual(ep.stretchTensor, reference_stress));
  auto independent =
      tester->assembleVec(ep.dmElastic, ep.elementVolumeName, reference, trial,
                          zero, zero, 1., deltaTime, {});
  CHKERR VecAXPY(independent, -1., direct);
  PetscReal stress_error, stress_norm;
  CHKERR getFieldNorm(independent, indices.at(ep.stretchTensor), stress_error);
  CHKERR getFieldNorm(direct, indices.at(ep.stretchTensor), stress_norm);
  if (!std::isfinite(stress_error) ||
      stress_error > 1.e-11 * std::max(1., stress_norm))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Independent direct material residual error %g, norm %g",
            stress_error, stress_norm);
  std::array<double, 3> local{audit->noncoaxial, audit->rateGradient,
                              audit->energyError},
      global;
  PetscInt points;
  CHKERR MPI_Allreduce(local.data(), global.data(), 3, MPI_DOUBLE, MPI_MAX,
                       ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->points, &points, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  if (!points || !(global[0] > 1.e-7) || !(global[1] > 1.e-7) ||
      !std::isfinite(global[2]) || global[2] > 1.e-12)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Direct reference coverage failed: noncoaxial %g, rate gradient "
            "%g, energy error %g",
            global[0], global[1], global[2]);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Direct logarithmic material passed: stress error %.3e, "
                     "energy error %.3e, noncoaxial %.3e, rate gradient %.3e\n",
                     stress_error, global[2], global[0], global[1]);
  CHKERR checkExternalStrain(ep, state, zero);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] =
    "Check direct logarithmic Neo-Hookean material equations.\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(runDirectAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}