auxiliary_logarithmic_stress_output_atom.cpp

Go to the documentation of this file.

/**
 * @file auxiliary_logarithmic_stress_output_atom.cpp
 * @brief Check production output meanings and diagnostic state preservation.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <MatrixFunction.hpp>
#include <NeoHookeanLogarithmicMaterial.hpp>
using namespace EshelbianPlasticity;
 
#include <array>
#include <cmath>
 
namespace {
 
using Material = NeoHookeanLogarithmicMaterial;
using Coordinates = Material::Coordinates;
 
MoFEMErrorCode checkNear(const double actual, const double expected,
                         const char *name, const double tolerance = 2.e-6) {
  MoFEMFunctionBegin;
  if (!std::isfinite(actual) ||
      std::abs(actual - expected) >
          tolerance * std::max(1., std::abs(expected)))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "%s: expected %.16e, obtained %.16e", name, expected, actual);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setState(EshelbianCore &ep, Vec state, const Coordinates &t_d,
                        const double theta, const Coordinates &t_stress) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(state);
  const Problem *problem = nullptr;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  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 || dof->getDofOrder() != 0)
      continue;
    const auto &name = dof->getName();
    const auto component = dof->getDofCoeffIdx();
    double value = 0.;
    if (name == ep.logDeviator)
      value = t_d(component);
    else if (name == ep.logJacobian)
      value = theta;
    else if (name == ep.auxiliaryLogStress)
      value = t_stress(component);
    else if (name == ep.rotAxis)
      value = rotation.at(component);
    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);
}
 
MoFEMErrorCode readTag(moab::Interface &mesh, const EntityHandle vertex,
                       const std::string &name, const int size,
                       double *values) {
  MoFEMFunctionBegin;
  Tag tag;
  CHKERR mesh.tag_get_handle(name.c_str(), tag);
  int length;
  CHKERR mesh.tag_get_length(tag, length);
  if (length != size)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Output tag %s has length %d instead of %d", name.c_str(), length,
            size);
  CHKERR mesh.tag_get_data(tag, &vertex, 1, values);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkTensor(moab::Interface &mesh, const EntityHandle vertex,
                           const std::string &name,
                           const Material::SymmetricTensor &t_expected) {
  MoFEMFunctionBegin;
  std::array<double, 9> values;
  CHKERR readTag(mesh, vertex, name, values.size(), values.data());
  auto t_actual = getFTensor2FromPtr<3, 3>(values.data());
  FTENSOR_INDEXES(3, i, j);
  FTensor::Tensor2<double, 3, 3> t_error;
  t_error(i, j) = t_actual(i, j) - t_expected(i, j);
  CHKERR checkNear(std::sqrt(t_error(i, j) * t_error(i, j)), 0., name.c_str());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkVectorUnchanged(Vec before, Vec after,
                                    const char *description) {
  MoFEMFunctionBegin;
  auto difference = vectorDuplicate(after);
  CHKERR VecWAXPY(difference, -1., before, after);
  PetscReal norm;
  CHKERR VecNorm(difference, NORM_INFINITY, &norm);
  CHKERR checkNear(norm, 0., description, 0.);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkOutput(const std::string &file, const Coordinates &t_d,
                           const double theta, const Coordinates &t_stress,
                           const bool variation) {
  MoFEMFunctionBegin;
  const Material::Parameters parameters{1.7, 8.5};
  Material::InverseState inverse;
  Material::VolumeState volume;
  double deviatoric_energy, gap;
  CHKERR Material::evaluateDeviator(parameters, t_d, deviatoric_energy);
  CHKERR Material::evaluateVolume(parameters, theta, volume);
  CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
  CHKERR Material::evaluateFenchelGap(parameters, t_d, t_stress, inverse, gap);
  FTensor::Index<'a', 5> a;
  FTENSOR_INDEXES(3, i, j, k);
  Coordinates t_copy_error;
  t_copy_error(a) = t_d(a) - inverse.tMaterialDeviator(a);
  const double copy_error = std::sqrt(t_copy_error(a) * t_copy_error(a));
  const double physical_energy = deviatoric_energy + volume.energy;
  const double mixed_energy =
      t_stress(a) * t_d(a) - inverse.conjugateEnergy + volume.energy;
  if (gap <= .01 || copy_error <= .01)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Output test requires distinct material and kinematic copies");
  const auto t_d_tensor = Material::getTensor(t_d);
  const auto t_stress_tensor = Material::getTensor(t_stress);
  FTensor::Tensor2<double, 3, 3> t_commutator;
  t_commutator(i, j) = t_d_tensor(i, k) * t_stress_tensor(k, j) -
                       t_stress_tensor(i, k) * t_d_tensor(k, j);
  if (std::sqrt(t_commutator(i, j) * t_commutator(i, j)) <= .01)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Output test requires noncoaxial D and Td");
  const auto t_log_stretch = Material::getTensor(t_d, theta);
  FTensor::Tensor1<double, 3> t_eigenvalues;
  FTensor::Tensor2<double, 3, 3> t_eigenvectors;
  t_eigenvectors(i, j) = t_log_stretch(i, j);
  CHKERR computeEigenValuesSymmetric(t_eigenvectors, t_eigenvalues);
  const Material::SymmetricTensor t_stretch =
      EigenMatrix::getMat(t_eigenvalues, t_eigenvectors,
                          [](double value) { return std::exp(value); });
 
  moab::Core output;
  CHKERR output.load_file(file.c_str());
  for (const auto *name : {"FenchelGap", "PointwiseCopyMismatch"}) {
    Tag tag;
    if (output.tag_get_handle(name, tag) != MB_TAG_NOT_FOUND)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Unexpected material diagnostic output tag %s", name);
  }
  for (const auto *prefix :
       {"D_", "Td_", "D_m_", "VarD_", "VarTd_", "ResD_", "ResTd_"})
    for (int c = 0; c != 5; ++c) {
      const std::string name = prefix + std::to_string(c);
      Tag tag;
      if (output.tag_get_handle(name.c_str(), tag) != MB_TAG_NOT_FOUND)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Unexpected scalar coordinate output tag %s", name.c_str());
    }
  Range vertices;
  CHKERR output.get_entities_by_type(0, MBVERTEX, vertices);
  if (vertices.empty())
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Production output %s contains no vertices", file.c_str());
  for (const auto vertex : vertices) {
    CHKERR checkTensor(output, vertex, "D", t_d_tensor);
    CHKERR checkTensor(output, vertex, "Td", t_stress_tensor);
    CHKERR checkTensor(output, vertex, "D_m",
                       Material::getTensor(inverse.tMaterialDeviator));
    CHKERR checkTensor(output, vertex, "LogSpatialStretch", t_log_stretch);
    CHKERR checkTensor(output, vertex, "SpatialStretch", t_stretch);
    const std::array<std::pair<const char *, double>, 6> scalar_values{
        {{"theta", theta},
         {"J", volume.jacobian},
         {"PhysicalEnergy", physical_energy},
         {"MixedStoredEnergy", mixed_energy},
         {"Vartheta", .25 * theta},
         {"Restheta", -.75 * theta}}};
    double saved_physical = 0., saved_mixed = 0.;
    for (const auto &[name, expected] : scalar_values) {
      if (std::string(name) == "Vartheta" && !variation)
        continue;
      double actual;
      CHKERR readTag(output, vertex, name, 1, &actual);
      CHKERR checkNear(actual, expected, name);
      if (std::string(name) == "PhysicalEnergy")
        saved_physical = actual;
      else if (std::string(name) == "MixedStoredEnergy")
        saved_mixed = actual;
    }
    CHKERR checkNear(saved_physical - saved_mixed, gap,
                     "PhysicalEnergy - MixedStoredEnergy = FenchelGap", 6.e-6);
    Coordinates t_residual;
    t_residual(a) = -.75 * t_d(a);
    CHKERR checkTensor(output, vertex, "ResD", Material::getTensor(t_residual));
    t_residual(a) = -.75 * t_stress(a);
    CHKERR checkTensor(output, vertex, "ResTd",
                       Material::getTensor(t_residual));
    if (variation) {
      Coordinates t_variation;
      t_variation(a) = .25 * t_d(a);
      CHKERR checkTensor(output, vertex, "VarD",
                         Material::getTensor(t_variation));
      CHKERR checkTensor(output, vertex, "VarLogSpatialStretch",
                         Material::getTensor(t_variation, .25 * theta));
      t_variation(a) = .25 * t_stress(a);
      CHKERR checkTensor(output, vertex, "VarTd",
                         Material::getTensor(t_variation));
    }
  }
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "Checked %d production output vertices in %s; "
                     "gap %.6e, pointwise copy mismatch %.6e\n",
                     static_cast<int>(vertices.size()), file.c_str(), gap,
                     copy_error);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runOutputAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  if (ep.mField.get_comm_size() != 1)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "The output inspection atom runs on one rank");
  CHKERR ep.setElasticElementOps(1);
  const Coordinates t_d(.12, -.08, .06, -.04, .09);
  const Coordinates t_stress(-.8, .4, .5, .3, -.7);
  const double theta = .13;
  auto state = createDMVector(ep.dmElastic);
  CHKERR setState(ep, state, t_d, theta, t_stress);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, state, INSERT_VALUES,
                                 SCATTER_REVERSE);
  auto variation = vectorDuplicate(state);
  CHKERR VecCopy(state, variation);
  CHKERR VecScale(variation, .25);
  CHKERR VecGhostUpdateBegin(variation, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(variation, INSERT_VALUES, SCATTER_FORWARD);
  auto residual = vectorDuplicate(state);
  CHKERR VecCopy(state, residual);
  CHKERR VecScale(residual, -.75);
  CHKERR VecGhostUpdateBegin(residual, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(residual, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecCopy(state, ep.solTSStep);
  CHKERR VecScale(ep.solTSStep, -.5);
  CHKERR VecGhostUpdateBegin(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
  auto previous = vectorDuplicate(ep.solTSStep);
  CHKERR VecCopy(ep.solTSStep, previous);
  auto mesh_before = createDMVector(ep.dM);
  auto mesh_after = vectorDuplicate(mesh_before);
  CHKERR DMoFEMMeshToLocalVector(ep.dM, mesh_before, INSERT_VALUES,
                                 SCATTER_FORWARD);
 
  // Leave nonzero rate data in the shared material scratch. Mesh-only output
  // must evaluate the zero-rate inverse, independent of the preceding solve.
  boost::shared_ptr<VolumeElementForcesAndSourcesCore> rhs, lhs;
  CHKERR ep.setVolumeElementOps(1, true, false, rhs, lhs);
  ep.mField.getInterface<OperatorsTester>()->assembleVec(
      ep.dmElastic, ep.elementVolumeName, rhs, state, variation, nullptr, 1.,
      .4, {});
 
  const std::string skin_file = "auxiliary_output_skin.h5m";
  CHKERR ep.postProcessResults(1, skin_file, residual, variation);
  CHKERR DMoFEMMeshToLocalVector(ep.dM, mesh_after, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR checkVectorUnchanged(mesh_before, mesh_after,
                              "Skin output changed mesh coefficients");
  CHKERR checkVectorUnchanged(previous, ep.solTSStep,
                              "Skin output changed previous coefficients");
  CHKERR checkOutput(skin_file, t_d, theta, t_stress, true);
 
  const std::string skeleton_file = "auxiliary_output_skeleton.h5m";
  CHKERR ep.postProcessSkeletonResults(1, skeleton_file, residual);
  CHKERR DMoFEMMeshToLocalVector(ep.dM, mesh_after, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR checkVectorUnchanged(mesh_before, mesh_after,
                              "Skeleton output changed mesh coefficients");
  CHKERR checkVectorUnchanged(previous, ep.solTSStep,
                              "Skeleton output changed previous coefficients");
  CHKERR checkOutput(skeleton_file, t_d, theta, t_stress, false);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Auxiliary logarithmic stress output atom passed\n");
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] = "Check production auxiliary material diagnostic output.\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(runOutputAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}