PlasticFlowPreconditionerTests.hpp

Go to the documentation of this file.

#ifndef MOFEM_PLASTIC_FLOW_PRECONDITIONER_TESTS_HPP
#define MOFEM_PLASTIC_FLOW_PRECONDITIONER_TESTS_HPP
 
#include <LogarithmicDirectionalDerivatives.hpp>
 
namespace EshelbianPlasticity {
namespace {
 
/** Flow-only matrices on the control DM leave the latent rows empty. Put a
 * unit diagonal there so the flow block can be factorised and applied. */
MoFEMErrorCode addLatentIdentityRows(EshelbianCore &ep, DM control_dm,
                                     Mat matrix) {
  MoFEMFunctionBegin;
  IS raw_is = nullptr;
  CHKERR DMMoFEMGetFieldIS(control_dm, RowColData::ROW,
                           ep.plasticLatentField.c_str(), &raw_is);
  SmartPetscObj<IS> latent_is(raw_is);
  PetscInt size = 0;
  const PetscInt *indices = nullptr;
  CHKERR ISGetLocalSize(latent_is, &size);
  CHKERR ISGetIndices(latent_is, &indices);
  for (PetscInt i = 0; i != size; ++i)
    CHKERR MatSetValue(matrix, indices[i], indices[i], 1., ADD_VALUES);
  CHKERR ISRestoreIndices(latent_is, &indices);
  CHKERR MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
  MoFEMFunctionReturn(0);
}
 
/** Check the Gauss-point local curvature against shell/gradient differences
 * and verify numerical factor refresh at changing flow and hardening. */
MoFEMErrorCode testLocalFlowPreconditioner(
    EshelbianCore &ep, DM flow_dm, Vec flow_template, DM latent_dm, Vec latent,
    Vec latent_centre,
    PlasticIncrementalOptimizationInternal::PlasticLocalLawParameters
        parameters,
    Vec direction) {
  using namespace PlasticIncrementalOptimizationInternal;
  MoFEMFunctionBegin;
  Mat raw_matrix = nullptr;
  CHKERR DMCreateMatrix(flow_dm, &raw_matrix);
  SmartPetscObj<Mat> matrix(raw_matrix);
  auto flow = vectorDuplicate(flow_template);
  auto action = vectorDuplicate(flow);
  auto reference = vectorDuplicate(flow);
  auto trial = vectorDuplicate(flow);
  auto plus = vectorDuplicate(flow);
  auto minus = vectorDuplicate(flow);
  auto zero_latent = vectorDuplicate(latent);
  auto solution = vectorDuplicate(flow);
  auto solve_residual = vectorDuplicate(flow);
  KSP raw_ksp = nullptr;
  CHKERR KSPCreate(ep.mField.get_comm(), &raw_ksp);
  SmartPetscObj<KSP> ksp(raw_ksp);
  CHKERR KSPSetOptionsPrefix(ksp, "incremental_optimization_tensor_");
  PC pc = nullptr;
  CHKERR KSPGetPC(ksp, &pc);
  CHKERR PCSetType(pc, PCBJACOBI);
  CHKERR KSPSetFromOptions(ksp);
  CHKERR VecZeroEntries(zero_latent);
  parameters.regularisationEpsilon = 1e-4;
  PetscReal direction_norm = 0;
  CHKERR VecNorm(direction, NORM_2, &direction_norm);
  const double step = 1e-7 / std::max(1., direction_norm);
  for (const double hardening : {0., .3}) {
    parameters.hardeningModulus = hardening;
    for (const double flow_scale : {0., .01, .1}) {
      CHKERR VecCopy(direction, flow);
      CHKERR VecShift(flow, .13);
      CHKERR VecScale(flow, flow_scale);
      CHKERR assemblePlasticReconstructedLocalFlowPreconditioner(
          ep, flow_dm, flow, latent_dm, latent, latent_dm, latent_centre,
          parameters, matrix);
      CHKERR MatMult(matrix, direction, action);
      CHKERR VecZeroEntries(reference);
      CHKERR assemblePlasticReconstructedLocalFlowHessianAction(
          ep, flow_dm, flow, latent_dm, latent, latent_dm, latent_centre,
          direction, zero_latent, parameters, reference);
      PetscReal norm = 0, action_error = 0;
      CHKERR VecNorm(reference, NORM_2, &norm);
      CHKERR VecAXPY(reference, -1., action);
      CHKERR VecNorm(reference, NORM_2, &action_error);
      if (!std::isfinite(action_error) || action_error > 1e-11 * (1. + norm))
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Local matrix/action mismatch %g, norm %g, "
                "flow scale %g, hardening %g",
                action_error, norm, flow_scale, hardening);
 
      const auto gradient = [&](double shift, Vec out) {
        MoFEMFunctionBegin;
        CHKERR VecWAXPY(trial, shift, direction, flow);
        CHKERR VecZeroEntries(out);
        CHKERR assemblePlasticReconstructedLocalFlowGradient(
            ep, flow_dm, trial, latent_dm, latent, latent_dm, latent_centre,
            parameters, out);
        MoFEMFunctionReturn(0);
      };
      CHKERR gradient(step, plus);
      CHKERR gradient(-step, minus);
      CHKERR VecAXPY(plus, -1., minus);
      CHKERR VecScale(plus, .5 / step);
      CHKERR VecAXPY(plus, -1., action);
      PetscReal fd_error = 0;
      CHKERR VecNorm(plus, NORM_2, &fd_error);
      if (!std::isfinite(fd_error) || fd_error > 1e-8 + 1e-5 * norm)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Local matrix gradient difference error %g, "
                "norm %g, flow scale %g, hardening %g",
                fd_error, norm, flow_scale, hardening);
      MOFEM_LOG("EP", Sev::inform)
          << "Local flow preconditioner verified: flow scale " << flow_scale
          << ", hardening " << hardening << ", action error " << action_error
          << ", FD error " << fd_error;
 
      // Reuse one Pmat and one PC across changing flow and hardening. Check
      // the actual inverse action, so stale LU factors fail this regression.
      CHKERR addLatentIdentityRows(ep, flow_dm, matrix);
      CHKERR KSPSetOperators(ksp, matrix, matrix);
      CHKERR KSPSetUp(ksp);
      CHKERR PCApply(pc, direction, solution);
      CHKERR MatMult(matrix, solution, solve_residual);
      CHKERR VecAXPY(solve_residual, -1., direction);
      PetscReal solve_error = 0;
      CHKERR VecNorm(solve_residual, NORM_2, &solve_error);
      if (!std::isfinite(solve_error) ||
          solve_error > 1e-8 * std::max(direction_norm, PETSC_SMALL))
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Updated flow preconditioner inverse failed: "
                "relative residual %g, flow scale %g, hardening %g",
                solve_error / std::max(direction_norm, PETSC_SMALL), flow_scale,
                hardening);
      MOFEM_LOG("EP", Sev::inform)
          << "Flow preconditioner factors verified: flow scale " << flow_scale
          << ", hardening " << hardening << ", relative inverse residual "
          << solve_error / std::max(direction_norm, PETSC_SMALL);
    }
  }
  MoFEMFunctionReturn(0);
}
 
using TestMaterialActionByElement = std::map<EntityHandle, MatrixDouble>;
 
/** Independent fixed-F differences are a test oracle only. Production reads
 * Ue directly and never reconstructs the perturbed elastic deformation. */
struct OpTestLocalMaterialFlowAction : VolUserDataOperator {
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  OpTestLocalMaterialFlowAction(
      boost::shared_ptr<DataAtIntegrationPts> data,
      boost::shared_ptr<MatrixDouble> direction,
      boost::shared_ptr<TestMaterialActionByElement> actions)
      : VolUserDataOperator(NOSPACE, OPSPACE), data(data), direction(direction),
        actions(actions) {}
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEXES(3, i, j, k);
    FTensor::Index<'a', 5> a;
    FTensor::Index<'b', 5> b;
    const int count = getGaussPts().size2();
    if (count <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Local material assembly test requires multiple Gauss points");
    auto hp = data->getFTensorPlasticH(count);
    auto ue = data->getFTensorStretch(count);
    auto he = data->getFTensorLogStretch(count);
    auto v = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        *direction, count)();
    auto out = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
        (*actions)[getFEEntityHandle()], count)();
    const auto function = [](const FTensor::Tensor2_symmetric<double, 3> &h,
                             EigenMatrix::Fun<double> f) {
      FTENSOR_INDEXES(3, m, n);
      FTensor::Tensor1<double, 3> values;
      FTensor::Tensor2<double, 3, 3> vectors;
      vectors(m, n) = h(m, n);
      CHK_THROW_MESSAGE(computeEigenValuesSymmetric(vectors, values),
                        "Diagonalise independent material action test tensor");
      return EigenMatrix::getMat(values, vectors, f);
    };
    const auto &basis = Tensor2SymmetricDeviatorBasis::getBasis();
    for (int gg = 0; gg != count; ++gg) {
      FTensor::Tensor2_symmetric<double, 3> h, u, elastic_h;
      h(i, j) = hp(i, j);
      u(i, j) = ue(i, j);
      elastic_h(i, j) = he(i, j);
      const auto fp = function(h, [](double x) { return std::exp(x); });
      FTensor::Tensor2<double, 3, 3> anchor;
      anchor(i, j) = u(i, k) * fp(k, j);
      FTensor::Tensor2<double, 5, 5> map;
      constexpr double step = 2e-6;
      for (int col = 0; col != 5; ++col) {
        FTensor::Tensor1<double, 5> unit, coordinates;
        unit(a) = 0.;
        unit(col) = 1.;
        FTensor::Tensor2_symmetric<double, 3> tensor_v, difference;
        tensor_v(i, j) = basis(i, j, a) * unit(a);
        const auto evaluate = [&](double shift) {
          FTensor::Tensor2_symmetric<double, 3> trial, metric;
          trial(i, j) = h(i, j) + shift * tensor_v(i, j);
          const auto inverse =
              function(trial, [](double x) { return std::exp(-x); });
          FTensor::Tensor2<double, 3, 3> fe, cauchy;
          fe(i, j) = anchor(i, k) * inverse(k, j);
          cauchy(i, j) = fe(k, i) * fe(k, j);
          metric(i, j) = (cauchy(i, j) || cauchy(j, i)) / 2.;
          return function(metric, [](double x) { return .5 * std::log(x); });
        };
        const auto plus = evaluate(step), minus = evaluate(-step);
        difference(i, j) = (plus(i, j) - minus(i, j)) / (2. * step);
        coordinates(a) = basis(i, j, a) * difference(i, j);
        for (int row = 0; row != 5; ++row)
          map(row, col) = coordinates(row);
      }
      FTensor::Tensor2<double, 5, 5> tangent;
      CHKERR data->physicsPtr->evaluateElasticLogarithmicMaterialHessian(
          getFEEntityHandle(), gg, *data, elastic_h, tangent);
      FTensor::Tensor1<double, 5> direction_h, stress, source;
      direction_h(a) = map(a, b) * v(b);
      stress(a) = tangent(a, b) * direction_h(b);
      source(a) = map(b, a) * stress(b);
      out(a) = source(a);
      // The zero-adjoint material response independently validates the pure
      // tangent API and packed shear conventions, for both material adapters.
      LogarithmicMaterialCurvature oracle;
      oracle.stretch(i, j) = elastic_h(i, j);
      oracle.direction(i, j) = basis(i, j, a) * direction_h(a);
      oracle.adjoint(i, j) = 0.;
      oracle.stressAdjoint(a) = 0.;
      oracle.stressDirection(a) = 0.;
      CHKERR data->physicsPtr->evaluateHelmholtzCurvature(getFEEntityHandle(),
                                                          gg, *data, oracle);
      const FTensor::Tensor2_symmetric<double, 3> covector(
          oracle.stretchResult(0), .5 * oracle.stretchResult(1),
          .5 * oracle.stretchResult(2), oracle.stretchResult(3),
          .5 * oracle.stretchResult(4), oracle.stretchResult(5));
      FTensor::Tensor1<double, 5> error;
      error(a) = stress(a) - basis(i, j, a) * covector(i, j);
      if (std::sqrt(error(a) * error(a)) >
          1e-9 * (1. + std::sqrt(stress(a) * stress(a))))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Pure material Hessian differs from zero-adjoint response");
      ++hp;
      ++ue;
      ++he;
      ++v;
      ++out;
    }
    MoFEMFunctionReturn(0);
  }
  boost::shared_ptr<DataAtIntegrationPts> data;
  boost::shared_ptr<MatrixDouble> direction;
  boost::shared_ptr<TestMaterialActionByElement> actions;
};
 
struct OpTestLoadLocalMaterialAction : VolUserDataOperator {
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  OpTestLoadLocalMaterialAction(
      boost::shared_ptr<TestMaterialActionByElement> actions,
      boost::shared_ptr<MatrixDouble> values)
      : VolUserDataOperator(NOSPACE, OPSPACE), actions(actions),
        values(values) {}
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTensor::Index<'a', 5> a;
    const int count = getGaussPts().size2();
    auto input = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        actions->at(getFEEntityHandle()), count)();
    auto output = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
        *values, count)();
    for (int gg = 0; gg != count; ++gg) {
      output(a) = input(a);
      ++input;
      ++output;
    }
    MoFEMFunctionReturn(0);
  }
  boost::shared_ptr<TestMaterialActionByElement> actions;
  boost::shared_ptr<MatrixDouble> values;
};
 
MoFEMErrorCode assembleTestLocalMaterialFlowAction(EshelbianCore &ep,
                                                   DM flow_dm, Vec direction,
                                                   Vec action) {
  MoFEMFunctionBegin;
  const int order =
      2 * (std::max({ep.spaceOrder, ep.materialH1Order, ep.plasticFlowOrder,
                     ep.plasticLatentOrder}) +
           1);
  const auto plastic_only = [&ep](FEMethod *method) {
    return ep.plasticVolumes->find(method->getFEEntityHandle()) !=
           ep.plasticVolumes->end();
  };
  auto actions = boost::make_shared<TestMaterialActionByElement>();
  auto values = boost::make_shared<MatrixDouble>();
  auto mechanical =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  mechanical->ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, mechanical);
  mechanical->getRuleHook = [order](int, int, int) { return order; };
  mechanical->exeTestHook = plastic_only;
  mechanical->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<5>(
      SmartPetscObj<DM>(flow_dm, true), VolUserDataOperator::OPCOL,
      ep.plasticFlowField, values, SmartPetscObj<Vec>(direction, true), MBTET));
  mechanical->getOpPtrVector().push_back(
      new OpTestLocalMaterialFlowAction(ep.dataAtPts, values, actions));
  TetPolynomialBase::switchCacheBaseOn<HDIV>({mechanical.get()});
  const auto error =
      DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, mechanical);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({mechanical.get()});
  CHKERR error;
  auto flow = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  flow->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  flow->ksp_f = action;
  flow->getRuleHook = mechanical->getRuleHook;
  flow->exeTestHook = plastic_only;
  CHKERR MoFEM::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      flow->getOpPtrVector(), {NOSPACE}, ep.materialH1Positions);
  flow->getOpPtrVector().push_back(
      new OpTestLoadLocalMaterialAction(actions, values));
  using Form = FormsIntegrators<VolUserDataOperator>::Assembly<A>::LinearForm<
      GAUSS>::OpBaseTimesVector<1, 5, 5>;
  flow->getOpPtrVector().push_back(new Form(ep.plasticFlowField, values));
  CHKERR VecZeroEntries(action);
  CHKERR DMoFEMLoopFiniteElements(flow_dm, ep.elementVolumeName, flow);
  CHKERR VecAssemblyBegin(action);
  CHKERR VecAssemblyEnd(action);
  CHKERR VecGhostUpdateBegin(action, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(action, ADD_VALUES, SCATTER_REVERSE);
  MoFEMFunctionReturn(0);
}
 
/** Independent fixed-state/fixed-adjoint gradient difference, including
 * higher-order basis integration. No perturbed equilibrium or adjoint solve. */
MoFEMErrorCode testFrozenFlowPreconditioner(
    EshelbianCore &ep,
    PlasticIncrementalOptimizationInternal::
        ReducedBulkMaterialEnergyHessianWorkspace &workspace,
    TS ts, DM flow_dm, Vec base_flow, Vec base_state, Vec direction) {
  using namespace PlasticIncrementalOptimizationInternal;
  MoFEMFunctionBegin;
  Mat raw_matrix = nullptr;
  CHKERR DMCreateMatrix(flow_dm, &raw_matrix);
  SmartPetscObj<Mat> matrix(raw_matrix);
  KSP raw_ksp = nullptr;
  CHKERR KSPCreate(ep.mField.get_comm(), &raw_ksp);
  SmartPetscObj<KSP> ksp(raw_ksp);
  CHKERR KSPSetOptionsPrefix(ksp, "incremental_optimization_tensor_");
  PC pc = nullptr;
  CHKERR KSPGetPC(ksp, &pc);
  CHKERR PCSetType(pc, PCBJACOBI);
  CHKERR KSPSetFromOptions(ksp);
  auto current_flow = vectorDuplicate(base_flow);
  for (const double flow_shift : {0., .03, .12}) {
    CHKERR VecWAXPY(current_flow, flow_shift, direction, base_flow);
    CHKERR prepareReducedBulkMaterialEnergyHessian(
        ep, ts, flow_dm, current_flow, base_state, workspace);
    CHKERR MatZeroEntries(matrix);
    CHKERR assembleReducedBulkMaterialEnergyFlowPreconditioner(ep, workspace,
                                                               matrix);
    auto action = vectorDuplicate(base_flow);
    CHKERR MatMult(matrix, direction, action);
    auto material_action = vectorDuplicate(base_flow);
    CHKERR assembleTestLocalMaterialFlowAction(ep, flow_dm, direction,
                                               material_action);
    auto trial = vectorDuplicate(base_flow);
    auto plus = vectorDuplicate(base_flow);
    auto minus = vectorDuplicate(base_flow);
    const auto set_flow = [&](Vec flow) {
      MoFEMFunctionBegin;
      CHKERR setReducedHessianControl(ep, flow_dm, flow);
      CHKERR setStateOnMesh(ep, base_state);
      MoFEMFunctionReturn(0);
    };
    constexpr double step = 1e-4;
    const auto gradient = [&](double shift, Vec out) {
      MoFEMFunctionBegin;
      CHKERR VecWAXPY(trial, shift, direction, current_flow);
      CHKERR set_flow(trial);
      CHKERR assemblePlasticAdjointContribution(
          ep, flow_dm, getReducedBulkMaterialEnergyBaseAdjoint(workspace), out);
      MoFEMFunctionReturn(0);
    };
    const auto plus_error = gradient(step, plus);
    const auto minus_error = plus_error ? 0 : gradient(-step, minus);
    const auto restore_error = set_flow(current_flow);
    CHKERR plus_error;
    CHKERR minus_error;
    CHKERR restore_error;
    CHKERR VecAXPY(plus, -1., minus);
    CHKERR VecScale(plus, .5 / step);
    CHKERR VecAXPY(plus, 1., material_action);
    PetscReal reference_norm = 0;
    CHKERR VecNorm(plus, NORM_2, &reference_norm);
    CHKERR VecAXPY(plus, -1., action);
    PetscReal error = 0;
    CHKERR VecNorm(plus, NORM_2, &error);
    if (!std::isfinite(error) || error > 1e-8 + 1e-5 * reference_norm)
      SETERRQ(
          ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
          "Frozen flow plus local material curvature failed: error %g, norm %g",
          error, reference_norm);
 
    auto second = vectorDuplicate(direction);
    CHKERR VecCopy(direction, second);
    CHKERR VecShift(second, .17);
    auto second_action = vectorDuplicate(action);
    CHKERR MatMult(matrix, second, second_action);
    PetscScalar uv = 0, vu = 0;
    CHKERR VecDot(direction, second_action, &uv);
    CHKERR VecDot(second, action, &vu);
    const double symmetry_error = std::abs(PetscRealPart(uv - vu));
    if (symmetry_error > 1e-10 * (1. + std::max(std::abs(PetscRealPart(uv)),
                                                std::abs(PetscRealPart(vu)))))
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Frozen flow curvature symmetry error %g", symmetry_error);
    MOFEM_LOG("EP", Sev::inform)
        << "Frozen flow plus local material preconditioner verified: FD error "
        << error << ", norm " << reference_norm << ", symmetry error "
        << symmetry_error;
    auto solution = vectorDuplicate(base_flow);
    auto residual = vectorDuplicate(base_flow);
    CHKERR addLatentIdentityRows(ep, flow_dm, matrix);
    CHKERR KSPSetOperators(ksp, matrix, matrix);
    CHKERR KSPSetUp(ksp);
    CHKERR PCApply(pc, direction, solution);
    CHKERR MatMult(matrix, solution, residual);
    CHKERR VecAXPY(residual, -1., direction);
    PetscReal inverse_error = 0., rhs_norm = 0.;
    CHKERR VecNorm(residual, NORM_2, &inverse_error);
    CHKERR VecNorm(direction, NORM_2, &rhs_norm);
    if (!std::isfinite(inverse_error) || inverse_error > 1e-8 * rhs_norm)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Local material Pmat factor refresh failed: relative residual %g",
              inverse_error / rhs_norm);
    MOFEM_LOG("EP", Sev::inform)
        << "Local material Pmat factors verified: flow shift " << flow_shift
        << ", relative inverse residual " << inverse_error / rhs_norm;
  }
  CHKERR prepareReducedBulkMaterialEnergyHessian(ep, ts, flow_dm, base_flow,
                                                 base_state, workspace);
  MoFEMFunctionReturn(0);
}
 
} // namespace
} // namespace EshelbianPlasticity
#endif