6#undef PETSC_VERSION_RELEASE 
    7#define PETSC_VERSION_RELEASE 1 
    9#if PETSC_VERSION_GE(3, 6, 0) 
   10#include <petsc/private/petscimpl.h> 
   12#include <petsc-private/petscimpl.h> 
   15#if PETSC_VERSION_GE(3, 6, 0) 
   16#include <petsc/private/dmimpl.h>  
   19#include <petsc-private/dmimpl.h>   
   20#include <petsc-private/vecimpl.h>  
   34  template <InsertMode MODE>
 
   37                                    MatSORType flag, PetscReal shift,
 
   38                                    PetscInt its, PetscInt lits, Vec x);
 
 
   86  boost::ptr_vector<PCMGSubMatrixCtx>
 
 
  143template <InsertMode MODE>
 
  147  CHKERR MatShellGetContext(
a, &void_ctx);
 
  151  CHKERR VecScatterBegin(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
 
  152  CHKERR VecScatterEnd(ctx->
sCat, x, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
 
  154  CHKERR VecScatterBegin(ctx->
sCat, ctx->
F, f, MODE, SCATTER_FORWARD);
 
  155  CHKERR VecScatterEnd(ctx->
sCat, ctx->
F, f, MODE, SCATTER_FORWARD);
 
 
  161  return sub_mat_mult_generic<INSERT_VALUES>(
a, x, f);
 
 
  165  return sub_mat_mult_generic<ADD_VALUES>(
a, x, f);
 
 
  169                           PetscReal shift, PetscInt its, PetscInt lits,
 
  176  CHKERR MatShellGetContext(mat, &void_ctx);
 
  180  CHKERR VecScatterBegin(ctx->
sCat, b, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
 
  181  CHKERR VecScatterEnd(ctx->
sCat, b, ctx->
X, INSERT_VALUES, SCATTER_REVERSE);
 
  182  CHKERR MatSOR(ctx->
A, ctx->
X, 
omega, flag, shift, its, lits, ctx->
F);
 
  183  CHKERR VecScatterBegin(ctx->
sCat, ctx->
F, x, INSERT_VALUES, SCATTER_FORWARD);
 
  184  CHKERR VecScatterEnd(ctx->
sCat, ctx->
F, x, INSERT_VALUES, SCATTER_FORWARD);
 
 
  193  CHKERRABORT(PETSC_COMM_WORLD, 
ierr);
 
 
  214#define GET_DM_FIELD(DM)                                                       \ 
  216      static_cast<DMCtx *>(DM->data)->getInterface<DMMGViaApproxOrdersCtx>();  \ 
 
  220  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
 
  228                                                       bool create_sub_matrix,
 
  230  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  240      CHKERR AOApplicationToPetscIS(dm_field->aO, is2);
 
  242    if (create_sub_matrix) {
 
  248        CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
 
  251                              &(dm_field->shellMatrixCtxPtr.back()),
 
  253        CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT,
 
  255        CHKERR MatShellSetOperation(sub_a_raw, MATOP_MULT_ADD,
 
  257        CHKERR MatShellSetOperation(sub_a_raw, MATOP_SOR,
 
  259        dm_field->kspOperators.emplace_back(
 
  263#if PETSC_VERSION_GE(3, 8, 0) 
  264        CHKERR MatCreateSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
 
  266        CHKERR MatGetSubMatrix(
A, is2, is2, MAT_INITIAL_MATRIX, &sub_a_raw);
 
  268        dm_field->kspOperators.emplace_back(
 
  277  PetscInfo(dm, 
"Push back IS to DMMGViaApproxOrders\n");
 
 
  282  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  285  if (dm_field->coarseningIS.back()) {
 
  286    dm_field->coarseningIS.pop_back();
 
  288  dm_field->kspOperators.pop_back();
 
  289  PetscInfo(dm, 
"Pop back IS to DMMGViaApproxOrders\n");
 
 
  294  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  297  CHKERR dm_field->destroyCoarseningIS();
 
  298  PetscInfo(dm, 
"Clear DMs data structures\n");
 
 
  314  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  319    ((
DMCtxImpl *)(dm->data))->incrementReference();
 
  330#if PETSC_VERSION_GE(3, 17, 0) 
  331  PetscInfo(dm, 
"Create DMMGViaApproxOrders reference = %d\n",
 
  334  PetscInfo1(dm, 
"Create DMMGViaApproxOrders reference = %d\n",
 
 
  341  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  344  CHKERR dm_field->destroyCoarseningIS();
 
 
  351  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  355  int leveldown = dm->leveldown;
 
  357  if (dm_field->kspOperators.empty()) {
 
  361    CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
 
  362    if (dm_field->kspOperators.empty()) {
 
  364              "data inconsistency, operator can not be set");
 
  366    if (
static_cast<int>(dm_field->kspOperators.size()) < leveldown) {
 
  368              "data inconsistency, no IS for that level");
 
  370    *M = dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown];
 
  371    CHKERR PetscObjectReference((PetscObject)*M);
 
  375#if PETSC_VERSION_GE(3, 17, 0) 
  376  PetscInfo(dm, 
"Create Matrix DMMGViaApproxOrders leveldown = %d\n",
 
  379  PetscInfo1(dm, 
"Create Matrix DMMGViaApproxOrders leveldown = %d\n",
 
 
  387  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  390  CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
 
  391  CHKERR DMCreate(comm, dmc);
 
  392  (*dmc)->data = dm->data;
 
  394  CHKERR DMGetType(dm, &type);
 
  395  CHKERR DMSetType(*dmc, type);
 
  397#if PETSC_VERSION_GE(3, 17, 0) 
  398  PetscInfo(dm, 
"Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
 
  400  PetscInfo1(dm, 
"Coarsen DMMGViaApproxOrders leveldown = %d\n", dm->leveldown);
 
 
  407  PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
 
  408  PetscValidHeaderSpecific(dm2, DM_CLASSID, 1);
 
  412  CHKERR PetscObjectGetComm((PetscObject)dm1, &comm);
 
  419  int dm_down_leveldown = dm_down->leveldown;
 
  420  int dm_up_leveldown = dm_up->leveldown;
 
  421#if PETSC_VERSION_GE(3, 17, 0) 
  423            "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d " 
  424            "dm2_leveldown = %d\n",
 
  425            dm_down_leveldown, dm_up_leveldown);
 
  428             "Create interpolation DMMGViaApproxOrders dm1_leveldown = %d " 
  429             "dm2_leveldown = %d\n",
 
  430             dm_down_leveldown, dm_up_leveldown);
 
  437    if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_down_leveldown) {
 
  440    is_down = dm_field->coarseningIS[dm_field->coarseningIS.size() - 1 -
 
  442    CHKERR ISGetSize(is_down, &M);
 
  443    CHKERR ISGetLocalSize(is_down, &
m);
 
  448    if (
static_cast<int>(dm_field->coarseningIS.size()) < dm_up_leveldown) {
 
  453            ->coarseningIS[dm_field->coarseningIS.size() - 1 - dm_up_leveldown];
 
  455    CHKERR ISGetLocalSize(is_up, &
n);
 
  461  CHKERR MatCreate(comm, mat);
 
  463  CHKERR MatSetType(*mat, MATMPIAIJ);
 
  464  CHKERR MatMPIAIJSetPreallocation(*mat, 1, PETSC_NULLPTR, 0, PETSC_NULLPTR);
 
  467  PetscLayout rmap, cmap;
 
  468  CHKERR MatGetLayouts(*mat, &rmap, &cmap);
 
  469  int rstart, rend, cstart, cend;
 
  470  CHKERR PetscLayoutGetRange(rmap, &rstart, &rend);
 
  471  CHKERR PetscLayoutGetRange(cmap, &cstart, &cend);
 
  473  const int *row_indices_ptr, *col_indices_ptr;
 
  474  CHKERR ISGetIndices(is_down, &row_indices_ptr);
 
  475  CHKERR ISGetIndices(is_up, &col_indices_ptr);
 
  477  map<int, int> idx_map;
 
  478  for (
int ii = 0; ii < 
m; ii++) {
 
  479    idx_map[row_indices_ptr[ii]] = rstart + ii;
 
  482  CHKERR MatZeroEntries(*mat);
 
  484  for (
int jj = 0; jj < 
n; jj++) {
 
  485    map<int, int>::iterator mit = idx_map.find(col_indices_ptr[jj]);
 
  486    if (mit != idx_map.end()) {
 
  487      CHKERR MatSetValue(*mat, mit->second, cstart + jj, 1, INSERT_VALUES);
 
  491  CHKERR ISRestoreIndices(is_down, &row_indices_ptr);
 
  492  CHKERR ISRestoreIndices(is_up, &col_indices_ptr);
 
  494  CHKERR MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
 
  495  CHKERR MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
 
  498    *vec = PETSC_NULLPTR;
 
 
  505  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
 
  507  int leveldown = dm->leveldown;
 
  509  if (dm_field->kspOperators.empty()) {
 
  512#if PETSC_VERSION_GE(3, 5, 3) 
  514        dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
 
  518        dm_field->kspOperators[dm_field->kspOperators.size() - 1 - leveldown],
 
  523#if PETSC_VERSION_GE(3, 17, 0) 
  524  PetscInfo(dm, 
"Create global vector DMMGViaApproxOrders leveldown = %d\n",
 
  527  PetscInfo1(dm, 
"Create global vector DMMGViaApproxOrders leveldown = %d\n",
 
 
  535  PetscOptionsBegin(PETSC_COMM_WORLD, 
"",
 
  536                           "MOFEM Multi-Grid (Orders) pre-conditioner", 
"none");
 
  538  CHKERR PetscOptionsInt(
"-mofem_mg_levels", 
"nb levels of multi-grid solver",
 
  540  CHKERR PetscOptionsInt(
"-mofem_mg_coarse_order",
 
  541                         "approximation order of coarse level", 
"", 2,
 
  543  CHKERR PetscOptionsInt(
"-mofem_mg_order_at_last_level", 
"order at last level",
 
  545  CHKERR PetscOptionsInt(
"-mofem_mg_verbose", 
"nb levels of multi-grid solver",
 
  547  PetscBool shell_sub_a = 
shellSubA ? PETSC_TRUE : PETSC_FALSE;
 
  548  CHKERR PetscOptionsBool(
"-mofem_mg_shell_a", 
"use shell matrix as sub matrix",
 
  549                          "", shell_sub_a, &shell_sub_a, NULL);
 
 
  571                    ->getPetscGlobalDofIdx();
 
  577  string problem_name = problem_ptr->
getName();
 
  579                                              order_at_next_level, is);
 
 
  596  CHKERR PetscObjectGetComm((PetscObject)
dM, &comm);
 
  598  MOFEM_LOG_C(
"PCMGViaApproximationOrders", Sev::inform, 
"set MG levels %u\n",
 
  604  std::vector<IS> is_vec(
nbLevels + 1);
 
  607  for (
int kk = 0; kk < 
nbLevels; kk++) {
 
  611    CHKERR ISGetSize(is_vec[kk], &is_glob_size[kk]);
 
  612    CHKERR ISGetLocalSize(is_vec[kk], &is_loc_size[kk]);
 
  615                "Nb. dofs at level [ %d ] global %u local %d\n", kk,
 
  616                is_glob_size[kk], is_loc_size[kk]);
 
  619    if (is_glob_size[kk] == 0) {
 
  624  for (
int kk = 0; kk != 
nbLevels; kk++) {
 
  625    if (kk == 
nbLevels - 1 && use_mat_a) {
 
  642  for (
unsigned int kk = 0; kk < is_vec.size(); kk++) {
 
 
  652boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx>
 
  654  return boost::make_shared<PCMGSetUpViaApproxOrdersCtx>(dm, 
A, use_shell_mat);
 
 
  658    PC pc, boost::shared_ptr<PCMGSetUpViaApproxOrdersCtx> ctx, 
int verb) {
 
  660  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
 
  663  auto core_log = logging::core::get();
 
  665                                            "PCMGViaApproximationOrders"));
 
  667  MOFEM_LOG_TAG(
"PCMGViaApproximationOrders", 
"PCMGViaApproximationOrders");
 
  669  MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::verbose)
 
  670      << 
"Setup PCMGSetUpViaApproxOrders";
 
  673  CHKERR ctx->buildProlongationOperator(
true, verb);
 
  675#if PETSC_VERSION_GE(3, 8, 0) 
  676  CHKERR PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
 
  678  CHKERR PCMGSetGalerkin(pc, PETSC_FALSE);
 
  681  CHKERR PCMGSetLevels(pc, ctx->nbLevels, NULL);
 
  683  MOFEM_LOG(
"PCMGViaApproximationOrders", Sev::noisy)
 
  684      << 
"Setup PCMGSetUpViaApproxOrders <- end";
 
 
Implementation of DM context. You should not use it directly.
Discrete manager interface for MoFEM.
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
header of multi-grid solver for p- adaptivity
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr double omega
Save field DOFS on vertices/tags.
MoFEMErrorCode DMMGViaApproxOrdersPopBackCoarseningIS(DM dm)
Pop IS form MG via approximation orders.
PetscErrorCode DMDestroy_MoFEM(DM dm)
Destroys dm with MoFEM data structure.
MoFEMErrorCode DMMGViaApproxOrdersPushBackCoarseningIS(DM dm, IS is, Mat A, bool create_sub_matrix, bool shell_sub_a)
Push back coarsening level to MG via approximation orders.
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
MoFEMErrorCode DMCoarsen_MGViaApproxOrders(DM dm, MPI_Comm comm, DM *dmc)
Coarsen DM.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
MoFEMErrorCode DMMGViaApproxOrdersClearCoarseningIS(DM dm)
Clear approximation orders.
MoFEMErrorCode DMCreateMatrix_MGViaApproxOrders(DM dm, Mat *M)
Create matrix for Multi-Grid via approximation orders.
PetscErrorCode DMSetOperators_MoFEM(DM dm)
Set operators for MoFEM dm.
MoFEMErrorCode isCreateProblemOrder(const std::string problem_name, RowColData rc, int min_order, int max_order, IS *is) const
create IS for given order range (collective)
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
const double n
refractive index of diffusive medium
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
MoFEMErrorCode DMCreateInterpolation_MGViaApproxOrders(DM dm1, DM dm2, Mat *mat, Vec *vec)
Create interpolation matrix between data managers dm1 and dm2.
static MoFEMErrorCode ksp_set_operators(KSP ksp, Mat A, Mat B, void *ctx)
MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
MoFEMErrorCode DMCreateGlobalVector_MGViaApproxOrders(DM dm, Vec *g)
Create global vector for DMGViaApproxOrders.
MoFEMErrorCode DMMGViaApproxOrdersSetAO(DM dm, AO ao)
Set DM ordering.
auto matCreateVecs(Mat mat)
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
MoFEMErrorCode DMCreate_MGViaApproxOrders(DM dm)
Create DM data structure for Multi-Grid via approximation orders.
PetscErrorCode DMDestroy_MGViaApproxOrders(DM dm)
Destroy DM.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
MoFEMErrorCode sub_mat_mult_add(Mat a, Vec x, Vec f)
MoFEMErrorCode sub_mat_mult(Mat a, Vec x, Vec f)
FTensor::Index< 'm', 3 > m
Structure for DM for multi-grid via approximation orders.
std::vector< SmartPetscObj< Mat > > kspOperators
Get KSP operators.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
SmartPetscObj< DM > coarsenDM
std::vector< SmartPetscObj< IS > > coarseningIS
Coarsening IS.
virtual ~DMMGViaApproxOrdersCtx()
MoFEMErrorCode destroyCoarseningIS()
boost::ptr_vector< PCMGSubMatrixCtx > shellMatrixCtxPtr
Shell sub-matrix context.
Deprecated interface functions.
Section manager is used to create indexes and sections.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
int orderAtLastLevel
set maximal evaluated order
DM dM
Distributed mesh manager.
virtual ~PCMGSetUpViaApproxOrdersCtx()=default
virtual MoFEMErrorCode destroyIsAtLevel(int kk, IS *is)
Destroy IS if internally created.
int coarseOrder
approximation order of coarse level
virtual MoFEMErrorCode createIsAtLevel(int kk, IS *is)
Set IS for levels.
virtual MoFEMErrorCode getOptions()
get options from line command
PCMGSetUpViaApproxOrdersCtx(DM dm, Mat a, bool shell_sub_a)
virtual MoFEMErrorCode buildProlongationOperator(bool use_mat_a, int verb=0)
Set up data structures for MG.
Mat A
Matrix at fine level.
int nbLevels
number of multi-grid levels
friend MoFEMErrorCode sub_mat_mult_generic(Mat a, Vec x, Vec f)
virtual ~PCMGSubMatrixCtx()=default
friend MoFEMErrorCode sub_mat_sor(Mat mat, Vec b, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
PCMGSubMatrixCtx(Mat a, IS is)
PetscLogEvent MOFEM_EVENT_sor
MoFEMErrorCode initData(Vec x)
PetscLogEvent MOFEM_EVENT_mult
SmartPetscObj< VecScatter > sCat
keeps basic data about problem
DofIdx getNbLocalDofsRow() const
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
intrusive_ptr for managing petsc objects
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.