auxiliary_logarithmic_stress_jacobian_atom.cpp

Go to the documentation of this file.

/**
 * @file auxiliary_logarithmic_stress_jacobian_atom.cpp
 * @brief Check the production coupled auxiliary residual and Jacobian.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <AuxiliaryLogarithmicStressOperators.hpp>
using namespace EshelbianPlasticity;
 
namespace {
 
using VolumeElement = VolumeElementForcesAndSourcesCore;
using Material = NeoHookeanLogarithmicMaterial;
 
SmartPetscObj<IS> getFieldIS(EshelbianCore &ep, const std::string &name) {
  IS raw_is;
  CHK_THROW_MESSAGE(
      DMMoFEMGetFieldIS(ep.dmElastic, RowColData::ROW, name.c_str(), &raw_is),
      "Get test field indices");
  return SmartPetscObj<IS>(raw_is);
}
 
MoFEMErrorCode setState(EshelbianCore &ep, Vec state) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(state);
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  const std::array<double, 5> deviator{{.12, -.08, .06, -.04, .09}};
  const std::array<double, 5> stress{{-.21, .12, .31, .19, -.11}};
  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;
    const auto &name = dof->getName();
    const auto component = dof->getDofCoeffIdx();
    double value = .01 * std::sin(.37 * (index + 1));
    if (name == ep.logDeviator || name == ep.logJacobian ||
        name == ep.auxiliaryLogStress || name == ep.rotAxis) {
      if (dof->getDofOrder() == 0) {
        if (name == ep.logDeviator)
          value = deviator.at(component);
        else if (name == ep.logJacobian)
          value = .09;
        else if (name == ep.auxiliaryLogStress)
          value = stress.at(component);
        else
          value = rotation.at(component);
      } else {
        value *= .2 / (1 + dof->getDofOrder());
      }
    }
    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 setPlasticHistory(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  if (!ep.plasticVolume)
    MoFEMFunctionReturnHot(0);
  double history_norm = .2;
  CHKERR PetscOptionsGetReal(nullptr, nullptr,
                             "-auxiliary_plastic_history_norm", &history_norm,
                             nullptr);
  if (!std::isfinite(history_norm) || history_norm <= 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic Jacobian fixture requires positive history norm");
  const FTensor::Tensor1<double, 5> t_history(.13, -.29, .17, .31, -.11);
  FTensor::Index<'A', 5> A;
  const double scale = history_norm / std::sqrt(t_history(A) * t_history(A));
  auto set_history = [&](boost::shared_ptr<FieldEntity> field_entity) {
    MoFEMFunctionBegin;
    auto values = field_entity->getEntFieldData();
    if (values.size() != 5)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Plastic Jacobian fixture requires five P0 history coordinates");
    auto t_values = getFTensor1FromPtr<5>(&values[0]);
    t_values(A) = scale * t_history(A);
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      set_history, ep.plasticHField, ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
struct GaussAudit {
  PetscInt elements = 0;
  double noncoaxial = 0;
  double variation = 0;
  double plasticHistory = 0;
  double plasticStretch = 0;
  double plasticNoncoaxial = 0;
};
 
struct OpAuditGaussState : public VolumeElement::UserDataOperator {
  OpAuditGaussState(boost::shared_ptr<DataAtIntegrationPts> data,
                    boost::shared_ptr<GaussAudit> audit,
                    const Material::Parameters parameters, const double alpha_u)
      : VolumeElement::UserDataOperator(NOSPACE, OPSPACE), dataPtr(data),
        auditPtr(audit), parameters(parameters), alphaU(alpha_u) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    if (getFEMethod()->snes) {
      PetscBool domain_error;
      CHKERR SNESGetFunctionDomainError(getFEMethod()->snes, &domain_error);
      if (domain_error)
        MoFEMFunctionReturnHot(0);
    }
    const int points = getGaussPts().size2();
    if (points <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Coupled Jacobian test requires multiple Gauss points");
    if (alphaU != 0)
      CHKERR checkMaterialRate(points);
    auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        *dataPtr->auxiliaryData->logDeviator, points)();
    auto t_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
            *dataPtr->auxiliaryData->stress, points)();
    auto t_h_p = dataPtr->getFTensorPlasticH(points);
    auto t_f_p = dataPtr->getFTensorPlasticF(points);
    FTensor::Index<'A', 5> A;
    FTENSOR_INDEXES(3, i, j, k);
    FTensor::Tensor1<double, 5> t_first, t_difference;
    t_first(A) = t_d(A);
    for (int gg = 0; gg != points; ++gg) {
      const auto t_d_tensor = Tensor2SymmetricDeviatorBasis::getTensor(t_d);
      const auto t_stress_tensor =
          Tensor2SymmetricDeviatorBasis::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);
      auditPtr->noncoaxial =
          std::max(auditPtr->noncoaxial,
                   std::sqrt(t_commutator(i, j) * t_commutator(i, j)));
      t_difference(A) = t_d(A) - t_first(A);
      auditPtr->variation = std::max(
          auditPtr->variation, std::sqrt(t_difference(A) * t_difference(A)));
      if (std::abs(t_h_p(i, i)) > 1.e-12 ||
          std::abs(determinantTensor3by3(t_f_p) - 1.) > 1.e-12)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Plastic history must remain trace-free and isochoric");
      auditPtr->plasticHistory = std::max(auditPtr->plasticHistory,
                                          std::sqrt(t_h_p(i, j) * t_h_p(i, j)));
      const auto t_identity = FTensor::Kronecker_Delta<double>();
      FTensor::Tensor2<double, 3, 3> t_plastic_difference;
      t_plastic_difference(i, j) = t_f_p(i, j) - t_identity(i, j);
      auditPtr->plasticStretch = std::max(
          auditPtr->plasticStretch,
          std::sqrt(t_plastic_difference(i, j) * t_plastic_difference(i, j)));
      t_commutator(i, j) =
          t_h_p(i, k) * t_d_tensor(k, j) - t_d_tensor(i, k) * t_h_p(k, j);
      const double history_deviator_commutator =
          std::sqrt(t_commutator(i, j) * t_commutator(i, j));
      t_commutator(i, j) = t_h_p(i, k) * t_stress_tensor(k, j) -
                           t_stress_tensor(i, k) * t_h_p(k, j);
      auditPtr->plasticNoncoaxial = std::max(
          auditPtr->plasticNoncoaxial,
          std::min(history_deviator_commutator,
                   std::sqrt(t_commutator(i, j) * t_commutator(i, j))));
      ++t_d;
      ++t_stress;
      ++t_h_p;
      ++t_f_p;
    }
    ++auditPtr->elements;
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode checkMaterialRate(const int points) {
    MoFEMFunctionBegin;
    const auto &fields = *dataPtr->auxiliaryData;
    const auto &material = *dataPtr->auxiliaryMaterialData;
    using CoordinateView = GetFTensor1FromMatType<5, -1, DL>;
    auto get_coordinates = [&](const auto &values) {
      return MatrixSizeHelper<CoordinateView, DL>::get(*values, points)();
    };
    auto t_d = get_coordinates(fields.logDeviator);
    auto t_rate = get_coordinates(fields.logDeviatorDot);
    auto t_stress = get_coordinates(fields.stress);
    auto t_elastic_stress = get_coordinates(material.elasticStress);
    auto t_dm = get_coordinates(material.materialDeviator);
    auto t_energy =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
            *material.energy, points)();
    VectorDouble mixed_energy;
    CHKERR evaluateAuxiliaryMixedStoredEnergy(fields, material, mixed_energy);
    auto t_mixed = getFTensor0FromVec(mixed_energy);
    FTensor::Index<'A', 5> A;
    for (int gg = 0; gg != points; ++gg) {
      Material::Coordinates t_material_deviator, t_forward_stress, t_error;
      t_material_deviator(A) = t_dm(A);
      double material_energy;
      CHKERR Material::evaluateDeviator(parameters, t_material_deviator,
                                        material_energy, &t_forward_stress);
      t_error(A) = t_stress(A) - alphaU * t_rate(A) - t_forward_stress(A);
      const double constitutive_error = std::sqrt(t_error(A) * t_error(A));
      t_error(A) = t_elastic_stress(A) - t_forward_stress(A);
      const double elastic_error = std::sqrt(t_error(A) * t_error(A));
      const double expected_energy =
          material_energy + t_forward_stress(A) * (t_d(A) - t_dm(A)) +
          t_energy(AuxiliaryLogarithmicStressMaterialData::VOLUMETRIC);
      const double energy_error = std::abs(t_mixed - expected_energy);
      if (!std::isfinite(constitutive_error) || constitutive_error > 1.e-10 ||
          !std::isfinite(elastic_error) || elastic_error > 1.e-10 ||
          !std::isfinite(energy_error) || energy_error > 1.e-10)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Auxiliary rate material check: constitutive error %g, "
                "elastic stress error %g, stored energy error %g",
                constitutive_error, elastic_error, energy_error);
      ++t_d;
      ++t_rate;
      ++t_stress;
      ++t_elastic_stress;
      ++t_dm;
      ++t_energy;
      ++t_mixed;
    }
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<DataAtIntegrationPts> dataPtr;
  boost::shared_ptr<GaussAudit> auditPtr;
  Material::Parameters parameters;
  double alphaU;
};
 
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 zeroField(Vec vector, IS field) {
  MoFEMFunctionBegin;
  Vec subvector;
  CHKERR VecGetSubVector(vector, field, &subvector);
  CHKERR VecZeroEntries(subvector);
  CHKERR VecRestoreSubVector(vector, field, &subvector);
  CHKERR VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkCrossBlock(EshelbianCore &ep, Mat jacobian,
                               const std::string &row,
                               const std::string &column,
                               const bool symmetric = true) {
  MoFEMFunctionBegin;
  auto row_is = getFieldIS(ep, row);
  auto column_is = getFieldIS(ep, column);
  Mat raw_forward, raw_reverse, raw_transpose;
  CHKERR MatCreateSubMatrix(jacobian, row_is, column_is, MAT_INITIAL_MATRIX,
                            &raw_forward);
  SmartPetscObj<Mat> forward(raw_forward);
  CHKERR MatCreateSubMatrix(jacobian, column_is, row_is, MAT_INITIAL_MATRIX,
                            &raw_reverse);
  SmartPetscObj<Mat> reverse(raw_reverse);
  CHKERR MatTranspose(reverse, MAT_INITIAL_MATRIX, &raw_transpose);
  SmartPetscObj<Mat> transpose(raw_transpose);
  PetscReal norm, error;
  CHKERR MatNorm(forward, NORM_INFINITY, &norm);
  CHKERR MatAXPY(forward, -1., transpose, DIFFERENT_NONZERO_PATTERN);
  CHKERR MatNorm(forward, NORM_INFINITY, &error);
  if (!std::isfinite(norm) || norm <= 1.e-13 || !std::isfinite(error) ||
      (symmetric ? error > 1.e-12 * std::max(1., norm)
                 : error <= 1.e-12 * std::max(1., norm)))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Cross block %s/%s: norm %g, transpose difference %g, symmetric %d",
            row.c_str(), column.c_str(), norm, error,
            static_cast<int>(symmetric));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkDomainRecovery(EshelbianCore &ep,
                                   boost::shared_ptr<VolumeElement> rhs,
                                   Vec state, Vec rates, Vec accelerations,
                                   const double dt) {
  MoFEMFunctionBegin;
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  TsCtx ts_context(ep.mField, problem->getName());
  ts_context.getLoopsIFunction().emplace_back(ep.elementVolumeName, rhs);
  auto ts = createTS(ep.mField.get_comm());
  CHKERR TSSetTimeStep(ts, dt);
  CHKERR TSSetI2Function(ts, nullptr, TsSetI2Function, &ts_context);
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  struct ResidualContext {
    TS ts;
    Vec rates, accelerations;
  } residual_context{ts, rates, accelerations};
  auto input = vectorDuplicate(state);
  auto residual = vectorDuplicate(state);
  auto reference = vectorDuplicate(state);
  // Include PETSc's SNES input lock and MoFEM's TS assembly lifecycle, which
  // OperatorsTester alone does not exercise. I2Function also covers rates.
  CHKERR SNESSetFunction(
      snes, residual,
      [](SNES, Vec x, Vec f, void *ctx) -> PetscErrorCode {
        auto &data = *static_cast<ResidualContext *>(ctx);
        return TSComputeI2Function(data.ts, 1., x, data.rates,
                                   data.accelerations, f);
      },
      &residual_context);
  CHKERR VecCopy(state, input);
  CHKERR SNESComputeFunction(snes, input, reference);
  PetscReal reference_norm;
  CHKERR VecNorm(reference, NORM_2, &reference_norm);
  if (!std::isfinite(reference_norm))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Domain-recovery reference residual is not finite");
 
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(input, &first, &last);
  PetscInt injected = 0;
  if (ep.mField.get_comm_rank() == 0) {
    for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
      const auto index = dof->getPetscGlobalDofIdx();
      if (index >= first && index < last && dof->getName() == ep.logJacobian &&
          dof->getDofOrder() == 0) {
        // exp(theta) is finite, but the volumetric energy exp(2 theta)
        // overflows. Only one owned cell is outside the material domain.
        CHKERR VecSetValue(input, index, 400., INSERT_VALUES);
        injected = 1;
        break;
      }
    }
  }
  PetscInt total_injected;
  CHKERR MPI_Allreduce(&injected, &total_injected, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  if (total_injected != 1)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Domain-recovery fixture requires one owned logarithmic volume");
  CHKERR VecAssemblyBegin(input);
  CHKERR VecAssemblyEnd(input);
  CHKERR SNESComputeFunction(snes, input, residual);
  PetscBool domain_error;
  CHKERR SNESGetFunctionDomainError(snes, &domain_error);
  PetscInt local_domain = domain_error ? 1 : 0, domains;
  CHKERR MPI_Allreduce(&local_domain, &domains, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  PetscReal rejected_norm;
  CHKERR VecNorm(residual, NORM_2, &rejected_norm);
  if (!domains || std::isfinite(rejected_norm))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Nonrepresentable volume must flag the SNES residual domain");
  for (Vec vector : {input.get(), residual.get()}) {
    PetscInt lock;
    CHKERR VecLockGet(vector, &lock);
    if (lock)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Rejected SNES residual left a vector locked");
  }
 
  CHKERR VecCopy(state, input);
  CHKERR SNESComputeFunction(snes, input, residual);
  CHKERR SNESGetFunctionDomainError(snes, &domain_error);
  local_domain = domain_error ? 1 : 0;
  CHKERR MPI_Allreduce(&local_domain, &domains, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  CHKERR VecAXPY(residual, -1., reference);
  PetscReal error;
  CHKERR VecNorm(residual, NORM_2, &error);
  if (domains || !std::isfinite(error) ||
      error > 1.e-12 * std::max(1., reference_norm))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "SNES residual did not recover after domain rejection: error %g",
            error);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Auxiliary residual domain recovery passed: error %.3e\n",
                     error);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runJacobianAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  CHKERR setPlasticHistory(ep);
  auto state = createDMVector(ep.dmElastic);
  auto rates = createDMVector(ep.dmElastic);
  auto accelerations = createDMVector(ep.dmElastic);
  CHKERR setState(ep, state);
  // OperatorsTester takes increments and scales them by dt and dt squared.
  // Nonzero histories exercise the shared mechanical rate-data producers.
  constexpr double dt = .4;
  CHKERR VecCopy(state, rates);
  CHKERR VecScale(rates, .3);
  CHKERR VecCopy(state, accelerations);
  CHKERR VecScale(accelerations, -.2);
  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<GaussAudit>();
  Material::Parameters parameters{};
  CHKERR PetscOptionsGetReal(nullptr, nullptr, "-neo_hookean_c10",
                             &parameters.c10, nullptr);
  CHKERR PetscOptionsGetReal(nullptr, nullptr, "-neo_hookean_K",
                             &parameters.bulkModulus, nullptr);
  CHKERR Material::validateParameters(parameters);
  rhs->getOpPtrVector().push_back(
      new OpAuditGaussState(ep.dataAtPts, audit, parameters, ep.alphaU));
  auto *tester = ep.mField.getInterface<OperatorsTester>();
  auto jacobian = tester->assembleMat(ep.dmElastic, ep.elementVolumeName, lhs,
                                      state, rates, accelerations, 1., dt, {});
  auto residual = tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs,
                                      state, rates, accelerations, 1., dt, {});
 
  std::vector<std::string> fields{ep.logDeviator,        ep.logJacobian,
                                  ep.auxiliaryLogStress, ep.piolaStress,
                                  ep.bubbleField,        ep.rotAxis,
                                  ep.spatialL2Disp,      ep.hybridSpatialDisp};
  std::map<std::string, SmartPetscObj<IS>> field_indices;
  for (const auto &name : fields) {
    auto indices = getFieldIS(ep, name);
    PetscInt size;
    CHKERR ISGetSize(indices, &size);
    if (!size)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Missing coupled test field %s", name.c_str());
    field_indices.emplace(name, indices);
  }
 
  auto zero_rates = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(zero_rates);
  auto rate_residual =
      tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs, state,
                          zero_rates, zero_rates, 1., dt, {});
  CHKERR VecAYPX(rate_residual, -1., residual);
  for (const auto &[field, active] :
       {std::pair{ep.spatialL2Disp, ep.alphaW != 0 || ep.alphaRho != 0},
        std::pair{ep.rotAxis,
                  ep.alphaViscousR != 0 || ep.alphaViscousOmega != 0},
        std::pair{ep.auxiliaryLogStress, ep.alphaU != 0},
        std::pair{ep.logDeviator, false}, std::pair{ep.logJacobian, false}}) {
    PetscReal norm;
    CHKERR getFieldNorm(rate_residual, field_indices.at(field), norm);
    if (!std::isfinite(norm) || (active ? norm <= 1.e-12 : norm > 1.e-12))
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Rate residual %s: norm %g, active %d", field.c_str(), norm,
              static_cast<int>(active));
  }
 
  if (ep.alphaU != 0) {
    auto zero_deviator_rates = vectorDuplicate(rates);
    CHKERR VecCopy(rates, zero_deviator_rates);
    CHKERR zeroField(zero_deviator_rates, field_indices.at(ep.logDeviator));
    auto zero_deviator_residual =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs, state,
                            zero_deviator_rates, accelerations, 1., dt, {});
    auto deviator_rate_effect = vectorDuplicate(residual);
    CHKERR VecWAXPY(deviator_rate_effect, -1., zero_deviator_residual,
                    residual);
    for (const auto &field : fields) {
      PetscReal norm;
      CHKERR getFieldNorm(deviator_rate_effect, field_indices.at(field), norm);
      const bool active = field == ep.auxiliaryLogStress;
      if (!std::isfinite(norm) || (active ? norm <= 1.e-12 : norm > 1.e-12))
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Deviator rate must affect only Td: row %s, norm %g",
                field.c_str(), norm);
    }
 
    // All other rates, including log J, stay nonzero. Removing alpha_u must
    // give the identical residual once the D rate vanishes.
    const auto alpha_u = ep.alphaU;
    ep.alphaU = 0;
    boost::shared_ptr<VolumeElement> elastic_rhs, elastic_lhs;
    CHKERR ep.setVolumeElementOps(1, true, false, elastic_rhs, elastic_lhs);
    ep.alphaU = alpha_u;
    auto elastic_residual =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, elastic_rhs,
                            state, rates, accelerations, 1., dt, {});
    CHKERR VecAXPY(elastic_residual, -1., zero_deviator_residual);
    PetscReal error;
    CHKERR VecNorm(elastic_residual, NORM_INFINITY, &error);
    if (!std::isfinite(error) || error > 1.e-12)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Zero D rate must recover alpha_u=0 residual: error %g", error);
    CHKERR PetscPrintf(ep.mField.get_comm(),
                       "Auxiliary deviator rate isolation and zero-rate "
                       "recovery passed: error %.3e\n",
                       error);
  }
 
  auto trial = vectorDuplicate(state);
  auto trial_rates = vectorDuplicate(state);
  auto trial_accelerations = vectorDuplicate(state);
  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);
    // Perturb both increments consistently with the assembled TS shifts.
    CHKERR VecWAXPY(trial_rates, -1., state, trial);
    CHKERR VecCopy(trial_rates, trial_accelerations);
    CHKERR VecAXPY(trial_rates, 1., rates);
    CHKERR VecAXPY(trial_accelerations, 1., accelerations);
    auto value =
        tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs, trial,
                            trial_rates, trial_accelerations, 1., dt, {});
    CHKERR VecCopy(value, output);
    MoFEMFunctionReturn(0);
  };
  auto direction = vectorDuplicate(state);
  auto mixed = vectorDuplicate(state);
  auto error = vectorDuplicate(residual);
  auto action = vectorDuplicate(residual);
  CHKERR VecZeroEntries(mixed);
  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 max_error = 0;
    for (const auto &row : fields) {
      PetscReal row_error, row_action;
      CHKERR getFieldNorm(error, field_indices.at(row), row_error);
      CHKERR getFieldNorm(action, field_indices.at(row), 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,
                "Jacobian row %s, direction %s: FD error %g exceeds %g "
                "(Jacobian action %g)",
                row.c_str(), name.c_str(), row_error, tolerance, row_action);
      max_error = std::max(max_error, row_error);
    }
    CHKERR PetscPrintf(ep.mField.get_comm(),
                       "Auxiliary coupled Jacobian direction %s: error %.3e\n",
                       name.c_str(), max_error);
    MoFEMFunctionReturn(0);
  };
  for (const auto &name : fields) {
    CHKERR VecZeroEntries(direction);
    Vec subvector;
    auto indices = field_indices.at(name);
    CHKERR VecGetSubVector(direction, indices, &subvector);
    PetscInt first, last;
    CHKERR VecGetOwnershipRange(subvector, &first, &last);
    PetscScalar *values;
    CHKERR VecGetArray(subvector, &values);
    for (PetscInt index = first; index != last; ++index)
      values[index - first] = std::sin(.71 * (index + 1)) + .3;
    CHKERR VecRestoreArray(subvector, &values);
    CHKERR VecRestoreSubVector(direction, indices, &subvector);
    PetscReal norm;
    CHKERR VecNormalize(direction, &norm);
    if (!(norm > 0.))
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Zero test direction for %s", name.c_str());
    CHKERR VecGhostUpdateBegin(direction, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(direction, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAXPY(mixed, 1., direction);
    CHKERR check_direction(name);
  }
  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");
 
  for (const auto &material : {ep.logDeviator, ep.logJacobian})
    for (const auto &mechanical : {ep.piolaStress, ep.bubbleField, ep.rotAxis})
      CHKERR checkCrossBlock(ep, jacobian, material, mechanical);
  CHKERR checkCrossBlock(ep, jacobian, ep.logDeviator, ep.logJacobian);
  CHKERR checkCrossBlock(ep, jacobian, ep.logDeviator, ep.auxiliaryLogStress,
                         ep.alphaU == 0);
 
  // At zero total Piola stress there is no DD material Hessian: f(D) has
  // already been replaced by Td:D-f*(Td), and only stress work contributes DD.
  CHKERR VecCopy(state, trial);
  CHKERR zeroField(trial, field_indices.at(ep.piolaStress));
  CHKERR zeroField(trial, field_indices.at(ep.bubbleField));
  auto zero_stress_jacobian =
      tester->assembleMat(ep.dmElastic, ep.elementVolumeName, lhs, trial, rates,
                          accelerations, 1., dt, {});
  Mat raw_dd;
  CHKERR MatCreateSubMatrix(
      zero_stress_jacobian, field_indices.at(ep.logDeviator),
      field_indices.at(ep.logDeviator), MAT_INITIAL_MATRIX, &raw_dd);
  SmartPetscObj<Mat> dd(raw_dd);
  PetscReal dd_norm;
  CHKERR MatNorm(dd, NORM_INFINITY, &dd_norm);
  if (!std::isfinite(dd_norm) || dd_norm > 1.e-13)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Zero-stress DD block must vanish; norm %g", dd_norm);
 
  PetscInt elements;
  double noncoaxial, variation;
  CHKERR MPI_Allreduce(&audit->elements, &elements, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->noncoaxial, &noncoaxial, 1, MPI_DOUBLE, MPI_MAX,
                       ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->variation, &variation, 1, MPI_DOUBLE, MPI_MAX,
                       ep.mField.get_comm());
  if (!elements || !(noncoaxial > 1.e-4) || !(variation > 1.e-7))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Insufficient Gauss-state coverage: elements %d, commutator %g, "
            "variation %g",
            static_cast<int>(elements), noncoaxial, variation);
  if (ep.plasticVolume) {
    const std::array<double, 3> local{
        audit->plasticHistory, audit->plasticStretch, audit->plasticNoncoaxial};
    std::array<double, 3> global;
    CHKERR MPI_Allreduce(local.data(), global.data(), 3, MPI_DOUBLE, MPI_MAX,
                         ep.mField.get_comm());
    if (global[0] < 1.e-3 || global[1] < 1.e-3 || global[2] < 1.e-4)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Insufficient plastic state coverage: Hp norm %g, Fp-I norm %g, "
              "minimum Hp/D and Hp/Td commutator %g",
              global[0], global[1], global[2]);
    CHKERR PetscPrintf(ep.mField.get_comm(),
                       "Auxiliary plastic Jacobian: Hp norm %.3e, "
                       "Fp-I norm %.3e, commutator %.3e\n",
                       global[0], global[1], global[2]);
  }
  CHKERR checkDomainRecovery(ep, rhs, state, rates, accelerations, dt);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Auxiliary coupled Jacobian atom passed: "
                     "commutator %.3e, Gauss variation %.3e\n",
                     noncoaxial, variation);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] = "Verify the full auxiliary logarithmic stress Jacobian.\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(runJacobianAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}