SnesCtx.cpp

Go to the documentation of this file.

/**
 * @brief Context for PETSc SNES, i.e. nonlinear solver
 * @author Anonymous contributor under MIT license
 *
 */
 
namespace MoFEM {
 
struct SnesCtx::SnesCtxImpl {
 
  using PairNameFEMethodPtr = MoFEM::PairNameFEMethodPtr;
  using FEMethodsSequence = MoFEM::FEMethodsSequence;
  using BasicMethodsSequence = MoFEM::BasicMethodsSequence;
  using HookFunction = SnesCtx::HookFunction;
 
  SnesCtxImpl(Interface &m_field, std::string problem_name);
 
  virtual ~SnesCtxImpl() = default;
 
  FEMethodsSequence &getSetOperators() { return loopsOperator; }
 
  FEMethodsSequence &getComputeRhs() { return loopsRhs; }
 
  BasicMethodsSequence &getPreProcComputeRhs() { return preProcessRhs; }
 
  BasicMethodsSequence &getPostProcComputeRhs() { return postProcessRhs; }
 
  BasicMethodsSequence &getPreProcSetOperators() { return preProcessOperator; }
 
  BasicMethodsSequence &getPostProcSetOperators() {
    return postProcessOperator;
  }
 
  FEMethodsSequence &getLoadTangent() { return loopsLoadTangent; }
 
  BasicMethodsSequence &getPreProcLoadTangent() {
    return preProcessLoadTangent;
  }
 
  BasicMethodsSequence &getPostProcLoadTangent() {
    return postProcessLoadTangent;
  }
 
  MoFEMErrorCode copyLoops(const SnesCtxImpl &snes_ctx);
 
  MoFEMErrorCode clearLoops();
 
  HookFunction &getLoadTangentHook() { return loadTangentHook; }
 
  HookFunction &getRhsHook() { return rhsHook; }
 
  friend PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx);
 
  friend PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx);
 
  friend MoFEMErrorCode SNESMoFEMSetAssmblyType(SNES snes,
                                                MatAssemblyType type);
  friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh);
 
  friend MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes,
                                                 MatAssemblyType type);
 
  friend MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh);
 
  friend MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its,
                                               PetscReal fgnorm,
                                               SnesCtx *snes_ctx);
 
  friend MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its,
                                               PetscReal fgnorm,
                                               SnesCtx *snes_ctx);
 
  friend PetscErrorCode SnesLoadTangent(SNES snes, Vec u, Vec F, void *ctx);
 
protected:
  PetscLogEvent MOFEM_EVENT_SnesRhs; ///< Log events to assemble residual
  PetscLogEvent MOFEM_EVENT_SnesMat; ///< Log events to assemble tangent matrix
  PetscLogEvent MOFEM_EVENT_SnesLoadTangent; ///< Log events to assemble tangent
                                             ///< of load vector
 
  MoFEM::Interface &mField; ///< database Interface
  moab::Interface &moab;    ///< moab Interface
 
  std::string problemName; ///< problem name
 
  FEMethodsSequence loopsOperator; ///< Sequence of finite elements instances
                                   ///< assembling tangent matrix
  FEMethodsSequence loopsRhs;      ///< Sequence of finite elements instances
                                   ///< assembling residual vector
  FEMethodsSequence
      loopsLoadTangent; ///< Sequence of finite elements instances assembling
                        ///< tangent of load vector in arc length control
 
  BasicMethodsSequence preProcessOperator;  ///< Sequence of methods run before
                                            ///< tangent matrix is assembled
  BasicMethodsSequence postProcessOperator; ///< Sequence of methods run after
                                            ///< tangent matrix is assembled
  BasicMethodsSequence
      preProcessLoadTangent; ///< Sequence of methods run before
                             ///< tangent of load vector in arc
                             ///< length control is assembled
  BasicMethodsSequence
      preProcessRhs; ///< Sequence of methods run before residual is assembled
  BasicMethodsSequence
      postProcessRhs; ///< Sequence of methods run after residual is assembled
  BasicMethodsSequence
      postProcessLoadTangent; ///< Sequence of methods run after tangent of load
                              ///< vector in arc length control is assembled
 
  MoFEMTypes bH =
      MF_EXIST; ///< If set to MF_EXIST check if element exist, default MF_EXIST
  bool zeroPreCondMatrixB = true; ///< If true zero matrix, otherwise user need
                                  ///< to do it, default true
  bool vErify = false;            ///< If true verify vector
  MatAssemblyType typeOfAssembly =
      MAT_FINAL_ASSEMBLY; ///< type of assembly at the end
 
  HookFunction loadTangentHook; ///< Hook function called at the end of load
                                ///< tangent evaluation
 
  HookFunction rhsHook; ///< Hook function called at the end of residual
                        ///< evaluation
};
 
SnesCtx::SnesCtx(Interface &m_field, std::string problem_name)
    : snesCtxImpl(boost::movelib::make_unique<SnesCtx::SnesCtxImpl>(
          m_field, problem_name)) {}
 
SnesCtx::~SnesCtx() = default;
 
FEMethodsSequence &SnesCtx::getSetOperators() {
  return snesCtxImpl->getSetOperators();
}
 
FEMethodsSequence &SnesCtx::getComputeRhs() {
  return snesCtxImpl->getComputeRhs();
}
 
BasicMethodsSequence &SnesCtx::getPreProcComputeRhs() {
  return snesCtxImpl->getPreProcComputeRhs();
}
 
BasicMethodsSequence &SnesCtx::getPostProcComputeRhs() {
  return snesCtxImpl->getPostProcComputeRhs();
}
 
BasicMethodsSequence &SnesCtx::getPreProcSetOperators() {
  return snesCtxImpl->getPreProcSetOperators();
}
 
BasicMethodsSequence &SnesCtx::getPostProcSetOperators() {
  return snesCtxImpl->getPostProcSetOperators();
}
 
FEMethodsSequence &SnesCtx::getLoadTangent() {
  return snesCtxImpl->getLoadTangent();
}
 
BasicMethodsSequence &SnesCtx::getPreProcLoadTangent() {
  return snesCtxImpl->getPreProcLoadTangent();
}
 
BasicMethodsSequence &SnesCtx::getPostProcLoadTangent() {
  return snesCtxImpl->getPostProcLoadTangent();
}
 
MoFEMErrorCode SnesCtx::copyLoops(const SnesCtx &snes_ctx) {
  return snesCtxImpl->copyLoops(*snes_ctx.snesCtxImpl);
}
 
MoFEMErrorCode SnesCtx::clearLoops() { return snesCtxImpl->clearLoops(); }
 
SnesCtx::HookFunction &SnesCtx::getLoadTangentHook() {
  return snesCtxImpl->getLoadTangentHook();
}
 
SnesCtx::HookFunction &SnesCtx::getRhsHook() {
  return snesCtxImpl->getRhsHook();
}
 
SnesCtx::SnesCtxImpl::SnesCtxImpl(Interface &m_field, std::string problem_name)
    : mField(m_field), moab(m_field.get_moab()), problemName(problem_name) {
  PetscLogEventRegister("LoopSNESRhs", 0, &MOFEM_EVENT_SnesRhs);
  PetscLogEventRegister("LoopSNESMat", 0, &MOFEM_EVENT_SnesMat);
  PetscLogEventRegister("LoopSNESLoadTangent", 0, &MOFEM_EVENT_SnesLoadTangent);
  if (!LogManager::checkIfChannelExist("SNES_WORLD")) {
    auto core_log = logging::core::get();
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "SNES_WORLD"));
    LogManager::setLog("SNES_WORLD");
    MOFEM_LOG_TAG("SNES_WORLD", "SNES");
  }
}
 
MoFEMErrorCode SnesCtx::SnesCtxImpl::copyLoops(const SnesCtxImpl &snes_ctx) {
  MoFEMFunctionBeginHot;
  loopsOperator = snes_ctx.loopsOperator;
  loopsRhs = snes_ctx.loopsRhs;
  preProcessOperator = snes_ctx.preProcessOperator;
  postProcessOperator = snes_ctx.postProcessOperator;
  preProcessRhs = snes_ctx.preProcessRhs;
  postProcessRhs = snes_ctx.postProcessRhs;
  loopsLoadTangent = snes_ctx.loopsLoadTangent;
  preProcessLoadTangent = snes_ctx.preProcessLoadTangent;
  postProcessLoadTangent = snes_ctx.postProcessLoadTangent;
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode SnesCtx::SnesCtxImpl::clearLoops() {
  MoFEMFunctionBeginHot;
  loopsOperator.clear();
  loopsRhs.clear();
  preProcessOperator.clear();
  postProcessOperator.clear();
  preProcessRhs.clear();
  postProcessRhs.clear();
  loopsLoadTangent.clear();
  preProcessLoadTangent.clear();
  postProcessLoadTangent.clear();
  MoFEMFunctionReturnHot(0);
}
 
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx) {
  auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
  MoFEMFunctionBegin;
  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
  Vec dx;
  CHKERR SNESGetSolutionUpdate(snes, &dx);
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
  if (snes_ctx->vErify) {
    // Verify finite elements, check for not a number
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
    MPI_Comm comm = PetscObjectComm((PetscObject)f);
    PetscSynchronizedPrintf(comm, "SNES Verify x\n");
    const Problem *prb_ptr;
    CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
    CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
        prb_ptr, COL, x);
  }
  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
      snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
 
  auto zero_ghost_vec = [](Vec g) {
    MoFEMFunctionBegin;
    Vec l;
    CHKERR VecGhostGetLocalForm(g, &l);
    double *a;
    CHKERR VecGetArray(l, &a);
    int s;
    CHKERR VecGetLocalSize(l, &s);
    for (int i = 0; i != s; ++i)
      a[i] = 0;
    CHKERR VecRestoreArray(l, &a);
    CHKERR VecGhostRestoreLocalForm(g, &l);
    MoFEMFunctionReturn(0);
  };
  CHKERR zero_ghost_vec(f);
 
  auto vec_assemble_switch = boost::movelib::make_unique<bool>(true);
  auto cache_ptr = boost::make_shared<CacheTuple>();
  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
                                                 cache_ptr);
 
  auto set = [&](auto &fe) {
    MoFEMFunctionBeginHot;
    fe.snes = snes;
    fe.snes_x = x;
    fe.snes_dx = dx;
    fe.snes_f = f;
    fe.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
    fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
    fe.data_ctx = PetscData::CtxSetF | PetscData::CtxSetX | PetscData::CtxSetDX;
 
    CHKERR SNESGetKSP(snes, &fe.ksp);
 
    // If SNES is of type SNESNEWTONAL then get load parameter, will act as a
    // psuo time
#if PETSC_VERSION_GE(3, 22, 0)
    PetscBool same;
    CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
    if (same) {
      CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
      fe.data_ctx |= PetscData::CtxSetTime;
    }
#endif
 
    fe.cacheWeakPtr = cache_ptr;
    MoFEMFunctionReturnHot(0);
  };
 
  auto unset = [&](auto &fe) {
    fe.snes_ctx = SnesMethod::CTX_SNESNONE;
    fe.ksp_ctx = KspMethod::CTX_KSPNONE;
    fe.data_ctx = PetscData::CtxSetNone;
  };
 
  for (auto &bit : snes_ctx->preProcessRhs) {
    bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_preProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
  }
 
  for (auto &lit : snes_ctx->loopsRhs) {
    lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*lit.second);
    CHKERR snes_ctx->mField.loop_finite_elements(
        snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
        cache_ptr);
    unset(*lit.second);
    if (snes_ctx->vErify) {
      // Verify finite elements, check for not a number
      CHKERR VecAssemblyBegin(f);
      CHKERR VecAssemblyEnd(f);
      MPI_Comm comm = PetscObjectComm((PetscObject)f);
      PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
                              lit.first.c_str());
      const Problem *prb_ptr;
      CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
      CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
          prb_ptr, ROW, f);
    }
 
    vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
  }
 
  for (auto &bit : snes_ctx->postProcessRhs) {
    bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_postProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
  }
 
  if (*(vec_assemble_switch)) {
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
  }
 
  if (snes_ctx->rhsHook) {
    CHKERR snes_ctx->rhsHook(snes, x, f, cache_ptr, ctx);
  }
 
  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}
 
PetscErrorCode SnesLoadTangent(SNES snes, Vec x, Vec f, void *ctx) {
  auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
  MoFEMFunctionBegin;
  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
  Vec dx;
  CHKERR SNESGetSolutionUpdate(snes, &dx);
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
  if (snes_ctx->vErify) {
    // Verify finite elements, check for not a number
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
    MPI_Comm comm = PetscObjectComm((PetscObject)f);
    PetscSynchronizedPrintf(comm, "SNES Verify x\n");
    const Problem *prb_ptr;
    CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
    CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
        prb_ptr, COL, x);
  }
  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
      snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
 
  auto zero_ghost_vec = [](Vec g) {
    MoFEMFunctionBegin;
    Vec l;
    CHKERR VecGhostGetLocalForm(g, &l);
    double *a;
    CHKERR VecGetArray(l, &a);
    int s;
    CHKERR VecGetLocalSize(l, &s);
    for (int i = 0; i != s; ++i)
      a[i] = 0;
    CHKERR VecRestoreArray(l, &a);
    CHKERR VecGhostRestoreLocalForm(g, &l);
    MoFEMFunctionReturn(0);
  };
  CHKERR zero_ghost_vec(f);
 
  auto vec_assemble_switch = boost::movelib::make_unique<bool>(true);
  auto cache_ptr = boost::make_shared<CacheTuple>();
  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
                                                 cache_ptr);
 
  auto set = [&](auto &fe) {
    MoFEMFunctionBeginHot;
    fe.snes = snes;
    fe.snes_x = x;
    fe.snes_dx = dx;
    fe.snes_f = f;
    fe.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
    fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
    fe.data_ctx = PetscData::CtxSetF | PetscData::CtxSetX | PetscData::CtxSetDX;
 
    CHKERR SNESGetKSP(snes, &fe.ksp);
 
    // If SNES is of type SNESNEWTONAL then get load parameter, will act as a
    // pseudo time
#if PETSC_VERSION_GE(3, 22, 0)
    PetscBool same;
    CHKERR PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &same);
    if (same) {
      CHKERR SNESNewtonALGetLoadParameter(snes, &fe.ts_t);
      fe.data_ctx |= PetscData::CtxSetTime;
    }
#endif
 
    fe.cacheWeakPtr = cache_ptr;
    MoFEMFunctionReturnHot(0);
  };
 
  auto unset = [&](auto &fe) {
    fe.snes_ctx = SnesMethod::CTX_SNESNONE;
    fe.ksp_ctx = KspMethod::CTX_KSPNONE;
    fe.data_ctx = PetscData::CtxSetNone;
  };
 
  for (auto &bit : snes_ctx->preProcessLoadTangent) {
    bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_preProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
  }
 
  for (auto &lit : snes_ctx->loopsLoadTangent) {
    lit.second->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*lit.second);
    CHKERR snes_ctx->mField.loop_finite_elements(
        snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
        cache_ptr);
    unset(*lit.second);
    if (snes_ctx->vErify) {
      // Verify finite elements, check for not a number
      CHKERR VecAssemblyBegin(f);
      CHKERR VecAssemblyEnd(f);
      MPI_Comm comm = PetscObjectComm((PetscObject)f);
      PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
                              lit.first.c_str());
      const Problem *prb_ptr;
      CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
      CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
          prb_ptr, ROW, f);
    }
 
    vec_assemble_switch = boost::move(lit.second->vecAssembleSwitch);
  }
 
  for (auto &bit : snes_ctx->postProcessLoadTangent) {
    bit->vecAssembleSwitch = boost::move(vec_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_postProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    vec_assemble_switch = boost::move(bit->vecAssembleSwitch);
  }
 
  if (*(vec_assemble_switch)) {
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
  }
 
  if (snes_ctx->loadTangentHook) {
    CHKERR snes_ctx->loadTangentHook(snes, x, f, cache_ptr, ctx);
  }
 
  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesLoadTangent, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}
 
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx) {
  auto *snes_ctx = static_cast<SnesCtx *>(ctx)->snesCtxImpl.get();
  MoFEMFunctionBegin;
  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
 
  auto mat_assemble_switch = boost::movelib::make_unique<bool>(true);
  auto cache_ptr = boost::make_shared<CacheTuple>();
  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
                                                 cache_ptr);
  Vec dx;
  CHKERR SNESGetSolutionUpdate(snes, &dx);
 
  auto set = [&](auto &fe) {
    MoFEMFunctionBeginHot;
    fe.snes = snes;
    fe.snes_x = x;
    fe.snes_dx = dx;
    fe.snes_A = A;
    fe.snes_B = B;
    fe.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
    fe.ksp_ctx = KspMethod::CTX_OPERATORS;
    fe.data_ctx = PetscData::CtxSetA | PetscData::CtxSetB | PetscData::CtxSetX |
                  PetscData::CtxSetDX;
 
    CHKERR SNESGetKSP(snes, &fe.ksp);
 
    fe.cacheWeakPtr = cache_ptr;
    MoFEMFunctionReturnHot(0);
  };
 
  auto unset = [&](auto &fe) {
    fe.snes_ctx = SnesMethod::CTX_SNESNONE;
    fe.ksp_ctx = KspMethod::CTX_KSPNONE;
    fe.data_ctx = PetscData::CtxSetNone;
  };
 
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
      snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
 
  if (*mat_assemble_switch) {
    if (snes_ctx->zeroPreCondMatrixB)
      CHKERR MatZeroEntries(B);
  }
  for (auto &bit : snes_ctx->preProcessOperator) {
    bit->matAssembleSwitch = boost::move(mat_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_preProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    mat_assemble_switch = boost::move(bit->matAssembleSwitch);
  }
 
  for (auto &lit : snes_ctx->loopsOperator) {
    lit.second->matAssembleSwitch = boost::move(mat_assemble_switch);
    CHKERR set(*lit.second);
    CHKERR snes_ctx->mField.loop_finite_elements(
        snes_ctx->problemName, lit.first, *(lit.second), nullptr, snes_ctx->bH,
        cache_ptr);
    unset(*lit.second);
    mat_assemble_switch = boost::move(lit.second->matAssembleSwitch);
  }
 
  for (auto &bit : snes_ctx->postProcessOperator) {
    bit->matAssembleSwitch = boost::move(mat_assemble_switch);
    CHKERR set(*bit);
    CHKERR snes_ctx->mField.problem_basic_method_postProcess(
        snes_ctx->problemName, *bit);
    unset(*bit);
    mat_assemble_switch = boost::move(bit->matAssembleSwitch);
  }
 
  if (*mat_assemble_switch) {
    CHKERR MatAssemblyBegin(B, snes_ctx->typeOfAssembly);
    CHKERR MatAssemblyEnd(B, snes_ctx->typeOfAssembly);
  }
  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type) {
 
  auto get_snes_ctx = [](SNES snes) {
    SnesCtx *snes_ctx;
    CHK_THROW_MESSAGE(SNESGetApplicationContext(snes, &snes_ctx),
                      "Cannot get SNES context");
    return snes_ctx->snesCtxImpl.get();
  };
 
  MoFEMFunctionBegin;
  get_snes_ctx(snes)->typeOfAssembly = type;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh) {
 
  auto get_snes_ctx = [](SNES snes) {
    SnesCtx *snes_ctx;
    CHK_THROW_MESSAGE(SNESGetApplicationContext(snes, &snes_ctx),
                      "Cannot get SNES context");
    return snes_ctx->snesCtxImpl.get();
  };
 
  MoFEMFunctionBegin;
  get_snes_ctx(snes)->bH = bh;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm,
                                      SnesCtx *ctx) {
  MoFEMFunctionBegin;
  auto snes_ctx = ctx->snesCtxImpl.get();
  auto &m_field = snes_ctx->mField;
  auto problem_ptr = m_field.get_problem(snes_ctx->problemName);
  auto fields_ptr = m_field.get_fields();
  auto dofs = problem_ptr->numeredRowDofsPtr;
 
  std::vector<double> lnorms(fields_ptr->size(), 0),
      norms(fields_ptr->size(), 0);
 
  Vec res;
  CHKERR SNESGetFunction(snes, &res, NULL, NULL);
 
  const double *r;
  CHKERR VecGetArrayRead(res, &r);
  {
    int f = 0;
    for (auto fi : *fields_ptr) {
      const auto lo_uid = FieldEntity::getLoBitNumberUId(fi->bitNumber);
      const auto hi_uid = FieldEntity::getHiBitNumberUId(fi->bitNumber);
      const auto hi = dofs->get<Unique_mi_tag>().upper_bound(hi_uid);
      for (auto lo = dofs->get<Unique_mi_tag>().lower_bound(lo_uid); lo != hi;
           ++lo) {
        const DofIdx loc_idx = (*lo)->getPetscLocalDofIdx();
        if (loc_idx >= 0 && loc_idx < problem_ptr->nbLocDofsRow) {
          lnorms[f] += PetscRealPart(PetscSqr(r[loc_idx]));
        }
      }
      ++f;
    }
  }
  CHKERR VecRestoreArrayRead(res, &r);
 
  MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
                 MPIU_SUM, PetscObjectComm((PetscObject)snes));
 
  {
    MOFEM_LOG("SNES_WORLD", Sev::inform)
        << its << " Function norm " << boost::format("%14.12e") % (double)fgnorm
        << " [";
    int f = 0;
    for (auto fi : *fields_ptr) {
      MOFEM_LOG("SNES_WORLD", Sev::inform)
          << its << "\t Field norm "
          << boost::format("%14.12e") % (double)PetscSqrtReal(norms[f]) << " "
          << fi->getName();
      ++f;
    }
    MOFEM_LOG("SNES_WORLD", Sev::inform) << its << "]";
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm,
                                      SnesCtx *ctx) {
  MoFEMFunctionBegin;
  auto snes_ctx = ctx->snesCtxImpl.get();
  auto &m_field = snes_ctx->mField;
  auto problem_ptr = m_field.get_problem(snes_ctx->problemName);
  auto fields_ptr = m_field.get_fields();
  auto dofs = problem_ptr->numeredRowDofsPtr;
 
  std::vector<double> lnorms(3 * fields_ptr->size(), 0),
      norms(3 * fields_ptr->size(), 0);
 
  Vec x_update, res;
  if (its == 0)
    CHKERR SNESGetSolution(snes, &x_update);
  else
    CHKERR SNESGetSolutionUpdate(snes, &x_update);
  CHKERR SNESGetFunction(snes, &res, PETSC_NULLPTR, PETSC_NULLPTR);
 
  const double *x, *r;
  ;
  CHKERR VecGetArrayRead(res, &r);
  CHKERR VecGetArrayRead(x_update, &x);
 
  {
    int f = 0;
    for (auto fi : *fields_ptr) {
      const auto lo_uid = FieldEntity::getLoBitNumberUId(fi->bitNumber);
      const auto hi_uid = FieldEntity::getHiBitNumberUId(fi->bitNumber);
      const auto hi = dofs->get<Unique_mi_tag>().upper_bound(hi_uid);
      for (auto lo = dofs->get<Unique_mi_tag>().lower_bound(lo_uid); lo != hi;
           ++lo) {
        const DofIdx loc_idx = (*lo)->getPetscLocalDofIdx();
        if (loc_idx >= 0 && loc_idx < problem_ptr->nbLocDofsRow) {
          lnorms[3 * f] += std::abs(x[loc_idx] * r[loc_idx]);
          lnorms[3 * f + 1] += r[loc_idx] * r[loc_idx];
          lnorms[3 * f + 2] += x[loc_idx] * x[loc_idx];
        }
      }
      ++f;
    }
  }
  CHKERR VecRestoreArrayRead(res, &r);
  CHKERR VecRestoreArrayRead(x_update, &x);
 
  MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
                 MPIU_SUM, PetscObjectComm((PetscObject)snes));
 
  double energy_norm = 0, residual_norm = 0, increment_norm = 0;
  for (int f = 0; f != fields_ptr->size(); ++f) {
    energy_norm += norms[3 * f];
    residual_norm += norms[3 * f + 1];
    increment_norm += norms[3 * f + 2];
  }
  energy_norm = std::sqrt(energy_norm);
  residual_norm = std::sqrt(residual_norm);
  increment_norm = std::sqrt(increment_norm);
 
  {
    MOFEM_LOG("SNES_WORLD", Sev::inform)
        << its << " Function norm Enrgy/Residual/Increment"
        << boost::format(" %14.12e ") % energy_norm
        << boost::format(" %14.12e ") % residual_norm
        << boost::format(" %14.12e ") % increment_norm << "[";
    int f = 0;
    for (auto fi : *fields_ptr) {
      MOFEM_LOG("SNES_WORLD", Sev::inform)
          << its << "\t Energy/Residual/Increment norm "
          << boost::format("%14.12e") % std::sqrt(norms[3 * f]) << "\t"
          << boost::format("%14.12e") % std::sqrt(norms[3 * f + 1]) << " "
          << boost::format("%14.12e") % std::sqrt(norms[3 * f + 2]) << "\t "
          << fi->getName();
      ++f;
    }
    MOFEM_LOG("SNES_WORLD", Sev::inform) << its << " ]";
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace MoFEM