OperatorsTester.cpp

Go to the documentation of this file.

/**
 * @file OperatorsTester.cpp
 * @brief
 * @version 0.13.2
 * @date 2022-12-04
 *
 * @copyright Copyright (c) 2022
 *
 */
 
namespace MoFEM {
 
OperatorsTester::OperatorsTester(const Core &core)
    : cOre(const_cast<Core &>(core)) {}
 
MoFEMErrorCode
OperatorsTester::query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
  *iface = const_cast<OperatorsTester *>(this);
  return 0;
}
 
SmartPetscObj<Vec> OperatorsTester::setRandomFields(
    SmartPetscObj<DM> dm, std::vector<RandomFieldData> random_fields,
    boost::shared_ptr<Range> ents_ptr, RowColData r) {
 
  MoFEM::Interface &m_field = cOre;
 
  auto prb_ptr = getProblemPtr(dm);
  std::mt19937 rng(0x4d584431u);
 
  auto get_random_number = [&rng](auto &range) {
    std::uniform_real_distribution<double> dist(range[0], range[1]);
    return dist(rng);
  };
 
  auto v = createDMVector(dm, r);
  double *array;
  CHKERR VecGetArray(v, &array);
 
  int size;
  CHKERR VecGetLocalSize(v, &size);
 
  MOFEM_LOG_CHANNEL("WORLD");
 
  for (auto &rf : random_fields) {
    const auto field_id = m_field.get_field_bit_number(rf.first);
    MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "OperatorsTester")
        << "Set random field " << rf.first << " " << static_cast<int>(field_id);
 
    auto get_numbered_dofs = [&]() {
      switch (r) {
      case RowColData::ROW:
        return prb_ptr->getNumeredRowDofsPtr();
      case RowColData::COL:
        return prb_ptr->getNumeredColDofsPtr();
      default:
        throw std::runtime_error("Invalid RowColData");
      }
    };
 
    if (!ents_ptr) {
      auto lo = get_numbered_dofs()->lower_bound(
          FieldEntity::getLoBitNumberUId(field_id));
      const auto hi = get_numbered_dofs()->upper_bound(
          FieldEntity::getHiBitNumberUId(field_id));
      for (; lo != hi; ++lo) {
        const auto idx = (*lo)->getPetscLocalDofIdx();
        if (idx >= 0 && idx < size) {
          array[idx] = get_random_number(rf.second);
        }
      }
    } else {
      for (auto pit = ents_ptr->const_pair_begin();
           pit != ents_ptr->const_pair_end(); ++pit) {
        auto lo = get_numbered_dofs()->get<Unique_mi_tag>().lower_bound(
            DofEntity::getLoFieldEntityUId(field_id, pit->first));
        auto hi = get_numbered_dofs()->get<Unique_mi_tag>().upper_bound(
            DofEntity::getHiFieldEntityUId(field_id, pit->second));
        for (; lo != hi; ++lo) {
          const auto idx = (*lo)->getPetscLocalDofIdx();
          if (idx >= 0 && idx < size) {
            array[idx] = get_random_number(rf.second);
          }
        }
      }
    }
  }
 
  CHKERR VecRestoreArray(v, &array);
 
  CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
 
  return v;
};
 
SmartPetscObj<Vec>
OperatorsTester::assembleVec(SmartPetscObj<DM> dm, std::string fe_name,
                             boost::shared_ptr<FEMethod> pipeline,
                             SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
                             SmartPetscObj<Vec> delta2_x, double time,
                             double delta_t, CacheTupleWeakPtr cache_ptr) {
 
  pipeline->data_ctx = PetscData::CtxSetF | PetscData::CtxSetTime;
  auto f = createDMVector(dm);
  pipeline->f = f;
  pipeline->ts_t = time;
  pipeline->ts_dt = delta_t;
 
  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_SETFUNCTION;
  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETFUNCTION;
  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIFUNCTION;
 
  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
  std::ignore = a;
  std::ignore = v;
 
  CHK_THROW_MESSAGE(
      DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
      "loc to vec");
  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
                    "Run assemble pipeline");
  CHK_THROW_MESSAGE(VecAssemblyBegin(f), "Asmb vector");
  CHK_THROW_MESSAGE(VecAssemblyEnd(f), "Asmb vector");
 
  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD),
                    "Update vec");
 
  auto get_norm = [](auto f) {
    double fnorm;
    VecNorm(f, NORM_2, &fnorm);
    return fnorm;
  };
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Vec norm %3.4e",
                      get_norm(f));
 
  return f;
}
 
SmartPetscObj<Mat>
OperatorsTester::assembleMat(SmartPetscObj<DM> dm, std::string fe_name,
                             boost::shared_ptr<FEMethod> pipeline,
                             SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
                             SmartPetscObj<Vec> delta2_x, double time,
                             double delta_t, CacheTupleWeakPtr cache_ptr) {
  auto m = createDMMatrix(dm);
 
  pipeline->data_ctx =
      PetscData::CtxSetA | PetscData::CtxSetB | PetscData::CtxSetTime;
  pipeline->A = m;
  pipeline->B = m;
 
  pipeline->ts_t = time;
 
  pipeline->ts_dt = delta_t;
  pipeline->ts_a = 1. / delta_t;
  pipeline->ts_aa = 1. / pow(delta_t, 2);
 
  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_OPERATORS;
  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETJACOBIAN;
  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIJACOBIAN;
 
  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
  std::ignore = a;
  std::ignore = v;
 
  CHK_THROW_MESSAGE(
      DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
      "loc to vec");
  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
                    "Run assemble pipeline");
  CHK_THROW_MESSAGE(MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
  CHK_THROW_MESSAGE(MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
 
  auto get_norm = [](auto m) {
    double fnorm;
    MatNorm(m, NORM_1, &fnorm);
    return fnorm;
  };
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Matrix norm %3.4e",
                      get_norm(m));
 
  return m;
}
 
SmartPetscObj<Vec> OperatorsTester::directionalCentralFiniteDifference(
    SmartPetscObj<DM> dm, std::string fe_name,
    boost::shared_ptr<FEMethod> pipeline, SmartPetscObj<Vec> x,
    SmartPetscObj<Vec> delta_x, SmartPetscObj<Vec> delta2_x,
    SmartPetscObj<Vec> diff_x, double time, double delta_t, double eps,
    CacheTupleWeakPtr cache_ptr) {
 
  auto axpy = [&](auto v0, auto diff_v, auto eps) {
    SmartPetscObj<Vec> v;
    if (v0.use_count()) {
      v = vectorDuplicate(v0);
      CHK_THROW_MESSAGE(VecCopy(v0, v), "Cpy");
      CHK_THROW_MESSAGE(VecAXPY(v, eps, diff_v), "Add");
      CHK_THROW_MESSAGE(VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD),
                        "ghost");
      CHK_THROW_MESSAGE(VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD),
                        "ghost");
    }
   return v;
  };
 
  auto fm1 = assembleVec(dm, fe_name, pipeline,
                         //
                         axpy(x, diff_x, -eps),
                         //
                         axpy(delta_x, diff_x, -eps),
                         //
                         axpy(delta2_x, diff_x, -eps),
                         //
                         time, delta_t, cache_ptr);
 
  auto fp1 = assembleVec(dm, fe_name, pipeline,
                         //
                         axpy(x, diff_x, eps),
                         //
                         axpy(delta_x, diff_x, eps),
                         //
                         axpy(delta2_x, diff_x, eps),
                         //
                         time, delta_t, cache_ptr);
 
  auto calculate_derivative = [eps](auto fm1, auto fp1) {
    int size;
    VecGetLocalSize(fp1, &size);
    VecAXPY(fp1, -1, fm1);
    double *array;
    VecGetArray(fp1, &array);
    for (int j = 0; j != size; ++j) {
      array[j] /= 2 * eps;
    }
    VecRestoreArray(fp1, &array);
    return fp1;
  };
 
  return calculate_derivative(fm1, fp1);
}
 
SmartPetscObj<Vec> OperatorsTester::checkCentralFiniteDifference(
    SmartPetscObj<DM> dm, std::string fe_name,
    boost::shared_ptr<FEMethod> pipeline_rhs,
    boost::shared_ptr<FEMethod> pipeline_lhs, SmartPetscObj<Vec> x,
    SmartPetscObj<Vec> delta_x, SmartPetscObj<Vec> delta2_x,
    SmartPetscObj<Vec> diff_x, double time, double delta_t, double eps,
    CacheTupleWeakPtr cache_ptr) {
 
  auto fd_diff = directionalCentralFiniteDifference(
      dm, fe_name, pipeline_rhs, x, delta_x, delta2_x, diff_x, time, delta_t,
      eps, cache_ptr);
 
  auto m = assembleMat(dm, fe_name, pipeline_lhs, x, delta_x, delta2_x, time,
                       delta_t, cache_ptr);
 
  auto fm = vectorDuplicate(fd_diff);
  CHK_THROW_MESSAGE(MatMult(m, diff_x, fm), "Mat mult");
  CHK_THROW_MESSAGE(VecAXPY(fd_diff, -1., fm), "Add");
 
  return fd_diff;
}
 
MoFEMErrorCode OperatorsTester::checkScalarCentralFiniteDifference(
    Vec x, Vec direction, Vec gradient, const double eps,
    const ScalarFunction &evaluate, PetscReal &finite_difference,
    PetscReal &assembled_derivative) {
  MoFEMFunctionBegin;
  if (!(eps > 0) || !evaluate)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Scalar directional-derivative test requires a positive step "
            "and evaluator");
  auto x_plus = vectorDuplicate(x);
  auto x_minus = vectorDuplicate(x);
  CHKERR VecWAXPY(x_plus, eps, direction, x);
  CHKERR VecWAXPY(x_minus, -eps, direction, x);
  CHKERR VecGhostUpdateBegin(x_plus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x_plus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(x_minus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x_minus, INSERT_VALUES, SCATTER_FORWARD);
  PetscReal value_plus = 0;
  PetscReal value_minus = 0;
  CHKERR evaluate(x_plus, value_plus);
  CHKERR evaluate(x_minus, value_minus);
  finite_difference = (value_plus - value_minus) / (2 * eps);
  PetscScalar derivative = 0;
  CHKERR VecDot(gradient, direction, &derivative);
  assembled_derivative = PetscRealPart(derivative);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OperatorsTester::checkVectorCentralFiniteDifference(
    Vec x, Vec direction, Vec function_layout, Mat jacobian, const double eps,
    const VectorFunction &evaluate, Vec error) {
  MoFEMFunctionBegin;
  if (!(eps > 0) || !evaluate)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Vector directional-derivative test requires a positive step "
            "and evaluator");
  auto x_plus = vectorDuplicate(x);
  auto x_minus = vectorDuplicate(x);
  auto value_plus = vectorDuplicate(function_layout);
  auto value_minus = vectorDuplicate(function_layout);
  auto assembled = vectorDuplicate(function_layout);
  CHKERR VecWAXPY(x_plus, eps, direction, x);
  CHKERR VecWAXPY(x_minus, -eps, direction, x);
  CHKERR VecGhostUpdateBegin(x_plus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x_plus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(x_minus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x_minus, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR evaluate(x_plus, value_plus);
  CHKERR evaluate(x_minus, value_minus);
  CHKERR VecCopy(value_plus, error);
  CHKERR VecAXPY(error, -1., value_minus);
  CHKERR VecScale(error, 1. / (2 * eps));
  CHKERR MatMult(jacobian, direction, assembled);
  CHKERR VecAXPY(error, -1., assembled);
  MoFEMFunctionReturn(0);
}
 
std::pair<SmartPetscObj<Vec>, SmartPetscObj<Vec>>
OperatorsTester::setPipelineX(boost::shared_ptr<FEMethod> pipeline,
                              SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
                              SmartPetscObj<Vec> delta2_x, double delta_t) {
 
  // Set dofs vector to finite element instance
  pipeline->x = x;
  pipeline->data_ctx |= PetscData::CTX_SET_X;
 
  SmartPetscObj<Vec> x_t, x_tt;
 
  // Set velocity dofs vector to finite element instance
 
  if (delta_x.use_count()) {
    x_t = vectorDuplicate(x);
    VecCopy(delta_x, x_t);
    VecScale(x_t, 1. / delta_t);
    CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_t, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_t, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    pipeline->x_t = x_t;
    pipeline->data_ctx |= PetscData::CTX_SET_X_T;
  } else {
    pipeline->x_t = PETSC_NULLPTR;
  }
 
  // Set acceleration dofs vector to finite element instance
 
  if (delta2_x.use_count()) {
    x_tt = vectorDuplicate(x);
    VecCopy(delta2_x, x_tt);
    VecScale(x_tt, 1. / pow(delta_t, 2));
    CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_tt, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_tt, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    pipeline->x_tt = x_tt;
    pipeline->data_ctx |= PetscData::CTX_SET_X_TT;
  } else {
    pipeline->x_tt = PETSC_NULLPTR;
  }
 
  return std::make_pair(x_t, x_tt);
}
} // namespace MoFEM