Inheritance diagram for MoFEM::ConstrainMatrixCtxImpl:

Collaboration diagram for MoFEM::ConstrainMatrixCtxImpl:

Public Member Functions
	ConstrainMatrixCtxImpl (MPI_Comm comm, DM x_dm, DM y_dm)

	ConstrainMatrixCtxImpl (MPI_Comm comm, VecScatter scatter)

	~ConstrainMatrixCtxImpl () override=default

MoFEMErrorCode	setConstrainMatrix (Mat constrain_matrix, Mat inverse_cct) override

Mat	getConstrainMatrix () const override

Mat	getInverseCCT () const override

MoFEMErrorCode	setSystemMatrix (Mat system_matrix) override

Mat	getSystemMatrix () const override

MoFEMErrorCode	initializeQorP (Vec x) override

MoFEMErrorCode	initializeQTKQ () override

MoFEMErrorCode	recalculateCT () override

MoFEMErrorCode	recalculateCTC () override

MoFEMErrorCode	applyQ (Vec x, Vec f) override

MoFEMErrorCode	applyP (Vec x, Vec f) override

MoFEMErrorCode	applyR (Vec x, Vec f) override

MoFEMErrorCode	applyRT (Vec x, Vec f) override

MoFEMErrorCode	applyCTC_QTKQ (Vec x, Vec f) override

Public Member Functions inherited from MoFEM::ConstrainMatrixCtx
virtual	~ConstrainMatrixCtx ()=default

Private Member Functions
MoFEMErrorCode	invalidateProjectionData ()

MoFEMErrorCode	invalidateSystemData ()

MoFEMErrorCode	validatePetscObject (PetscObject object, const char *object_name) const

MoFEMErrorCode	validateVectorLayout (Vec vec, PetscInt local_size, PetscInt global_size, const char *vector_name) const

MoFEMErrorCode	validateInverseCCTLayout (Mat constrain_matrix, Mat inverse_cct) const

MoFEMErrorCode	validateSystemMatrixLayout () const

MoFEMErrorCode	applyInverseCCT (Vec x, Vec f, bool transpose)

MoFEMErrorCode	applyQImpl (Vec x, Vec f, bool transpose_inverse)

Private Attributes
MPI_Comm	comm

SmartPetscObj< DM >	xDm

SmartPetscObj< DM >	yDm

SmartPetscObj< Mat >	C

SmartPetscObj< Mat >	CT

SmartPetscObj< Mat >	inverseCCT

SmartPetscObj< Mat >	CTC

SmartPetscObj< Mat >	K

SmartPetscObj< Vec >	Cx

SmartPetscObj< Vec >	CCTm1Cx

SmartPetscObj< Vec >	CTCCTm1Cx

SmartPetscObj< Vec >	CTCx

SmartPetscObj< Vec >	X

SmartPetscObj< Vec >	Qx

SmartPetscObj< Vec >	KQx

SmartPetscObj< VecScatter >	scatter

bool	projectionInitialized = false

bool	qtkqInitialized = false

bool	buildScatterFromDms = true

PetscInt	fullLocalSize = 0

PetscInt	fullGlobalSize = 0

Detailed Description

Definition at line 115 of file ConstrainMatrixCtx.cpp.

Constructor & Destructor Documentation

◆ ConstrainMatrixCtxImpl() [1/2]

MoFEM::ConstrainMatrixCtxImpl::ConstrainMatrixCtxImpl	(	MPI_Comm	comm,
		DM	x_dm,
		DM	y_dm
	)

Definition at line 180 of file ConstrainMatrixCtx.cpp.

    : comm(comm), xDm(x_dm, true), yDm(y_dm, true) {
  getProjectionEvents();
}

◆ ConstrainMatrixCtxImpl() [2/2]

MoFEM::ConstrainMatrixCtxImpl::ConstrainMatrixCtxImpl	(	MPI_Comm	comm,
		VecScatter	scatter
	)

Definition at line 185 of file ConstrainMatrixCtx.cpp.

    : comm(comm), scatter(scatter, true), buildScatterFromDms(false) {
  getProjectionEvents();
}

◆ ~ConstrainMatrixCtxImpl()

MoFEM::ConstrainMatrixCtxImpl::~ConstrainMatrixCtxImpl ( )

overridedefault

Member Function Documentation

◆ applyCTC_QTKQ()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyCTC_QTKQ	(	Vec	x,
		Vec	f
	)

overridevirtual

Apply C^T C + Q^T K Q.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 658 of file ConstrainMatrixCtx.cpp.

                                                                 {
  MoFEMFunctionBegin;
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionCTCPlusQTKQ, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionCTCPlusQTKQ, comm);
 
  CHKERR initializeQorP(x);
  CHKERR initializeQTKQ();
  CHKERR validateSystemMatrixLayout();
  CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize,
                              "C^T C + Q^T K Q input");
  CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize,
                              "C^T C + Q^T K Q output");
 
  // Consume x into CTCx before Qx/f below can overwrite it, so that x and f
  // may safely alias the same Vec.
  CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(CTC, X, CTCx);
 
  CHKERR applyQImpl(x, Qx, false);
  CHKERR MatMult(K, Qx, KQx);
  CHKERR applyQImpl(KQx, f, true);
  CHKERR VecScatterBegin(scatter, CTCx, f, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecScatterEnd(scatter, CTCx, f, ADD_VALUES, SCATTER_REVERSE);
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ applyInverseCCT()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyInverseCCT	(	Vec	x,
		Vec	f,
		bool	transpose
	)

private

Definition at line 525 of file ConstrainMatrixCtx.cpp.

                                                                       {
  MoFEMFunctionBeginHot;
  if (transpose) {
    PetscBool has_mult_transpose;
    CHKERR MatHasOperation(inverseCCT, MATOP_MULT_TRANSPOSE,
                           &has_mult_transpose);
    if (!has_mult_transpose) {
      MOFEM_LOG("SELF", Sev::error)
          << "Inverse C C^T action does not implement MatMultTranspose";
      SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
              "Inverse C C^T action does not implement MatMultTranspose");
    }
    CHKERR MatMultTranspose(inverseCCT, x, f);
  } else {
    CHKERR MatMult(inverseCCT, x, f);
  }
  MoFEMFunctionReturnHot(0);
}

◆ applyP()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyP	(	Vec	x,
		Vec	f
	)

overridevirtual

Apply P = C^T S C for the supplied inverse action S.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 575 of file ConstrainMatrixCtx.cpp.

                                                          {
  MoFEMFunctionBegin;
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionP, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionP, comm);
 
  CHKERR initializeQorP(x);
  CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "P input");
  CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "P output");
  CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(C, X, Cx);
  CHKERR applyInverseCCT(Cx, CCTm1Cx, false);
  CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
  CHKERR VecZeroEntries(f);
  CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
  CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ applyQ()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyQ	(	Vec	x,
		Vec	f
	)

overridevirtual

Apply Q = I - C^T S C for the supplied inverse action S.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 563 of file ConstrainMatrixCtx.cpp.

                                                          {
  MoFEMFunctionBegin;
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionQ, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionQ, comm);
 
  CHKERR applyQImpl(x, f, false);
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ applyQImpl()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyQImpl	(	Vec	x,
		Vec	f,
		bool	transpose_inverse
	)

private

Definition at line 545 of file ConstrainMatrixCtx.cpp.

                                                                          {
  MoFEMFunctionBegin;
  CHKERR initializeQorP(x);
  CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "Q input");
  CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "Q output");
  CHKERR VecCopy(x, f);
  CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(C, X, Cx);
  CHKERR applyInverseCCT(Cx, CCTm1Cx, transpose_inverse);
  CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
  CHKERR VecScale(CTCCTm1Cx, -1);
  CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, ADD_VALUES, SCATTER_REVERSE);
  MoFEMFunctionReturn(0);
}

◆ applyR()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyR	(	Vec	x,
		Vec	f
	)

overridevirtual

Apply R = C^T S for the supplied inverse action S.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 599 of file ConstrainMatrixCtx.cpp.

                                                          {
  MoFEMFunctionBegin;
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionR, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionR, comm);
 
  if (!projectionInitialized) {
    MOFEM_LOG("SELF", Sev::error)
        << "R requires initializeQorP() or a previously created Q/P matrix";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Projection context is not initialized");
  }
  PetscInt constraint_local_size;
  PetscInt constraint_global_size;
  CHKERR VecGetLocalSize(Cx, &constraint_local_size);
  CHKERR VecGetSize(Cx, &constraint_global_size);
  CHKERR validateVectorLayout(x, constraint_local_size, constraint_global_size,
                              "R input");
  CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "R output");
  CHKERR applyInverseCCT(x, CCTm1Cx, false);
  CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
  CHKERR VecZeroEntries(f);
  CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
  CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ applyRT()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::applyRT	(	Vec	x,
		Vec	f
	)

overridevirtual

Apply R^T = S^T C for the supplied inverse action S.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 630 of file ConstrainMatrixCtx.cpp.

                                                           {
  MoFEMFunctionBegin;
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionRT, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionRT, comm);
 
  if (!projectionInitialized) {
    MOFEM_LOG("SELF", Sev::error)
        << "R^T requires initializeQorP() or a previously created Q/P matrix";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Projection context is not initialized");
  }
  PetscInt constraint_local_size;
  PetscInt constraint_global_size;
  CHKERR VecGetLocalSize(Cx, &constraint_local_size);
  CHKERR VecGetSize(Cx, &constraint_global_size);
  CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "R^T input");
  CHKERR validateVectorLayout(f, constraint_local_size, constraint_global_size,
                              "R^T output");
  CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(C, X, Cx);
  CHKERR applyInverseCCT(Cx, f, true);
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ getConstrainMatrix()

Mat MoFEM::ConstrainMatrixCtxImpl::getConstrainMatrix ( ) const

overridevirtual

Get the constraint matrix C.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 355 of file ConstrainMatrixCtx.cpp.

355{ return C.get(); }

◆ getInverseCCT()

Mat MoFEM::ConstrainMatrixCtxImpl::getInverseCCT ( ) const

overridevirtual

Get the supplied inverse action for C C^T.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 357 of file ConstrainMatrixCtx.cpp.

357{ return inverseCCT.get(); }

◆ getSystemMatrix()

Mat MoFEM::ConstrainMatrixCtxImpl::getSystemMatrix ( ) const

overridevirtual

Get the matrix K used by C^T C + Q^T K Q.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 367 of file ConstrainMatrixCtx.cpp.

367{ return K.get(); }

◆ initializeQorP()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::initializeQorP ( Vec x )

overridevirtual

Initialize matrices, vectors, and the scatter used by Q, P, R, and R^T. The operation is idempotent until the C/inverse-action pair is replaced.

Parameters

x	Vector defining the full problem layout.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 369 of file ConstrainMatrixCtx.cpp.

                                                           {
  MoFEMFunctionBegin;
  if (projectionInitialized)
    MoFEMFunctionReturnHot(0);
 
  CHKERR validateInverseCCTLayout(C, inverseCCT);
  CHKERR validatePetscObject(getPetscObject(x), "full vector x");
  if (buildScatterFromDms) {
    CHKERR validatePetscObject(getPetscObject(xDm.get()), "full problem DM");
    CHKERR validatePetscObject(getPetscObject(yDm.get()), "constraint DM");
  } else {
    CHKERR validatePetscObject(getPetscObject(scatter.get()),
                               "constraint scatter");
  }
 
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionInit, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionInit, comm);
 
  Mat raw_ct;
  CHKERR MatTranspose(C, MAT_INITIAL_MATRIX, &raw_ct);
  SmartPetscObj<Mat> new_ct(raw_ct);
 
  Vec raw_x;
  Vec raw_cx;
  CHKERR MatCreateVecs(C, &raw_x, &raw_cx);
  SmartPetscObj<Vec> new_x(raw_x);
  SmartPetscObj<Vec> new_cx(raw_cx);
 
  Vec raw_cctm1_cx;
  CHKERR VecDuplicate(new_cx, &raw_cctm1_cx);
  SmartPetscObj<Vec> new_cctm1_cx(raw_cctm1_cx);
 
  Vec raw_ct_cctm1_cx;
  CHKERR VecDuplicate(new_x, &raw_ct_cctm1_cx);
  SmartPetscObj<Vec> new_ct_cctm1_cx(raw_ct_cctm1_cx);
 
  SmartPetscObj<VecScatter> new_scatter = scatter;
  if (buildScatterFromDms) {
    MoFEM::Interface *x_m_field_ptr;
    MoFEM::Interface *y_m_field_ptr;
    CHKERR DMoFEMGetInterfacePtr(xDm, &x_m_field_ptr);
    CHKERR DMoFEMGetInterfacePtr(yDm, &y_m_field_ptr);
    if (x_m_field_ptr != y_m_field_ptr) {
      MOFEM_LOG("SELF", Sev::error)
          << "The full problem DM and constraint DM belong to different MoFEM "
             "interfaces";
      SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
              "DMs must belong to the same MoFEM interface");
    }
 
    const MoFEM::Problem *x_problem_ptr;
    const MoFEM::Problem *y_problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(xDm, &x_problem_ptr);
    CHKERR DMMoFEMGetProblemPtr(yDm, &y_problem_ptr);
    CHKERR x_m_field_ptr->getInterface<VecManager>()->vecScatterCreate(
        x, x_problem_ptr->getName(), ROW, new_x, y_problem_ptr->getName(), COL,
        new_scatter);
  }
 
  PetscInt new_full_local_size;
  PetscInt new_full_global_size;
  CHKERR VecGetLocalSize(x, &new_full_local_size);
  CHKERR VecGetSize(x, &new_full_global_size);
 
  PetscInt c_rows;
  PetscInt c_columns;
  PetscInt inverse_cct_rows;
  PetscInt inverse_cct_columns;
  CHKERR MatGetSize(C, &c_rows, &c_columns);
  CHKERR MatGetSize(inverseCCT, &inverse_cct_rows, &inverse_cct_columns);
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Constraint projection initialized: C size (" << c_rows << ", "
      << c_columns << "), inverse C C^T action size (" << inverse_cct_rows
      << ", " << inverse_cct_columns << "), full vector size "
      << new_full_global_size;
 
  CT = new_ct;
  X = new_x;
  Cx = new_cx;
  CCTm1Cx = new_cctm1_cx;
  CTCCTm1Cx = new_ct_cctm1_cx;
  scatter = new_scatter;
  fullLocalSize = new_full_local_size;
  fullGlobalSize = new_full_global_size;
  projectionInitialized = true;
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ initializeQTKQ()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::initializeQTKQ ( )

overridevirtual

Initialize matrices and work vectors used by C^T C + Q^T K Q. The operation is idempotent until C or K is replaced.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 472 of file ConstrainMatrixCtx.cpp.

                                                      {
  MoFEMFunctionBegin;
  if (qtkqInitialized)
    MoFEMFunctionReturnHot(0);
 
  CHKERR validatePetscObject(getPetscObject(C.get()), "constraint matrix C");
  CHKERR validateSystemMatrixLayout();
 
  const auto &events = getProjectionEvents();
  CHKERR PetscLogEventBegin(events.projectionInit, 0, 0, 0, 0);
  LogEventScope event_scope(events.projectionInit, comm);
 
  Mat raw_ctc;
  CHKERR MatTransposeMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &raw_ctc);
  SmartPetscObj<Mat> new_ctc(raw_ctc);
 
  Vec raw_qx;
  Vec raw_kqx;
  CHKERR MatCreateVecs(K, &raw_qx, &raw_kqx);
  SmartPetscObj<Vec> new_qx(raw_qx);
  SmartPetscObj<Vec> new_kqx(raw_kqx);
 
  Vec raw_ctcx;
  CHKERR MatCreateVecs(new_ctc, nullptr, &raw_ctcx);
  SmartPetscObj<Vec> new_ctcx(raw_ctcx);
 
  PetscInt ctc_rows;
  PetscInt ctc_columns;
  CHKERR MatGetSize(new_ctc, &ctc_rows, &ctc_columns);
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Initialized C^T C + Q^T K Q work data: C^T C size (" << ctc_rows
      << ", " << ctc_columns << ")";
 
  CTC = new_ctc;
  Qx = new_qx;
  KQx = new_kqx;
  CTCx = new_ctcx;
  qtkqInitialized = true;
 
  CHKERR event_scope.end();
  MoFEMFunctionReturn(0);
}

◆ invalidateProjectionData()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::invalidateProjectionData ( )

private

Definition at line 329 of file ConstrainMatrixCtx.cpp.

                                                                {
  MoFEMFunctionBeginHot;
  projectionInitialized = false;
  CT.reset();
  Cx.reset();
  CCTm1Cx.reset();
  CTCCTm1Cx.reset();
  X.reset();
  fullLocalSize = 0;
  fullGlobalSize = 0;
  if (buildScatterFromDms)
    scatter.reset();
  CHKERR invalidateSystemData();
  MoFEMFunctionReturnHot(0);
}

◆ invalidateSystemData()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::invalidateSystemData ( )

private

Definition at line 319 of file ConstrainMatrixCtx.cpp.

                                                            {
  MoFEMFunctionBeginHot;
  qtkqInitialized = false;
  CTC.reset();
  CTCx.reset();
  Qx.reset();
  KQx.reset();
  MoFEMFunctionReturnHot(0);
}

◆ recalculateCT()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::recalculateCT ( )

overridevirtual

Recalculate C^T after an in-place value change in C. The layout and nonzero structure of C must be unchanged. The caller must update the supplied inverse action in the same transaction. Projection MATSHELL objects already created from this context do not receive a PETSc state increment; recreate them before passing the changed operator to a consumer which caches by matrix state.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 460 of file ConstrainMatrixCtx.cpp.

                                                     {
  MoFEMFunctionBegin;
  if (!projectionInitialized)
    MoFEMFunctionReturnHot(0);
 
  CHKERR validateInverseCCTLayout(C, inverseCCT);
  Mat raw_mat = CT.get();
  CHKERR MatTranspose(C, MAT_REUSE_MATRIX, &raw_mat);
  MOFEM_LOG("WORLD", Sev::verbose) << "Recalculated C^T";
  MoFEMFunctionReturn(0);
}

◆ recalculateCTC()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::recalculateCTC ( )

overridevirtual

Recalculate C^T C after an in-place value change in C. The layout and nonzero structure of C must be unchanged.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 515 of file ConstrainMatrixCtx.cpp.

                                                      {
  MoFEMFunctionBegin;
  if (!qtkqInitialized)
    MoFEMFunctionReturnHot(0);
  Mat raw_mat = CTC.get();
  CHKERR MatTransposeMatMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &raw_mat);
  MOFEM_LOG("WORLD", Sev::verbose) << "Recalculated C^T C";
  MoFEMFunctionReturn(0);
}

◆ setConstrainMatrix()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::setConstrainMatrix	(	Mat	constrain_matrix,
		Mat	inverse_cct
	)

overridevirtual

Set the constraint matrix C and its row-space inverse action atomically.

inverse_cct represents (C C^T)^-1 on the selected constraint space, or the corresponding embedded generalized inverse for a masked full space. It must have the row layout of C in both dimensions and support MatMult(). MatMultTranspose() is additionally required by applyRT() and applyCTC_QTKQ(). Both PETSc objects are retained.

All data derived from the previous pair are invalidated. For in-place value changes which preserve layouts and sparsity, update inverse_cct in the same transaction and call recalculateCT() and recalculateCTC().

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 345 of file ConstrainMatrixCtx.cpp.

                                                                           {
  MoFEMFunctionBegin;
  CHKERR validateInverseCCTLayout(constrain_matrix, inverse_cct);
  C = SmartPetscObj<Mat>(constrain_matrix, true);
  inverseCCT = SmartPetscObj<Mat>(inverse_cct, true);
  CHKERR invalidateProjectionData();
  MoFEMFunctionReturn(0);
}

◆ setSystemMatrix()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::setSystemMatrix ( Mat system_matrix )

overridevirtual

Set the matrix K used by C^T C + Q^T K Q and retain a PETSc reference. Work vectors derived from the previous K are invalidated.

Implements MoFEM::ConstrainMatrixCtx.

Definition at line 359 of file ConstrainMatrixCtx.cpp.

                                                                        {
  MoFEMFunctionBegin;
  CHKERR validatePetscObject(getPetscObject(system_matrix), "system matrix K");
  K = SmartPetscObj<Mat>(system_matrix, true);
  CHKERR invalidateSystemData();
  MoFEMFunctionReturn(0);
}

◆ validateInverseCCTLayout()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::validateInverseCCTLayout	(	Mat	constrain_matrix,
		Mat	inverse_cct
	)		const

private

Definition at line 236 of file ConstrainMatrixCtx.cpp.

                                                                        {
  MoFEMFunctionBegin;
  CHKERR validatePetscObject(getPetscObject(constrain_matrix),
                             "constraint matrix C");
  CHKERR validatePetscObject(getPetscObject(inverse_cct),
                             "inverse C C^T action");
 
  PetscInt c_local_rows;
  PetscInt c_local_columns;
  PetscInt c_global_rows;
  PetscInt c_global_columns;
  PetscInt inverse_local_rows;
  PetscInt inverse_local_columns;
  PetscInt inverse_global_rows;
  PetscInt inverse_global_columns;
  CHKERR MatGetLocalSize(constrain_matrix, &c_local_rows, &c_local_columns);
  CHKERR MatGetSize(constrain_matrix, &c_global_rows, &c_global_columns);
  CHKERR MatGetLocalSize(inverse_cct, &inverse_local_rows,
                         &inverse_local_columns);
  CHKERR MatGetSize(inverse_cct, &inverse_global_rows, &inverse_global_columns);
  if (inverse_local_rows != c_local_rows ||
      inverse_local_columns != c_local_rows ||
      inverse_global_rows != c_global_rows ||
      inverse_global_columns != c_global_rows) {
    MOFEM_LOG("SELF", Sev::error)
        << "Inverse C C^T layout, local (" << inverse_local_rows << ", "
        << inverse_local_columns << "), global (" << inverse_global_rows << ", "
        << inverse_global_columns
        << "), does not match the row layout of C, local " << c_local_rows
        << ", global " << c_global_rows;
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Inverse C C^T action has an incompatible layout");
  }
 
  PetscBool has_mult;
  CHKERR MatHasOperation(inverse_cct, MATOP_MULT, &has_mult);
  if (!has_mult)
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Inverse C C^T action does not implement MatMult");
 
  PetscObject associated_c = nullptr;
  CHKERR PetscObjectQuery(getPetscObject(inverse_cct),
                          inverseCCTConstraintMatrixKey, &associated_c);
  if (associated_c && associated_c != getPetscObject(constrain_matrix)) {
    MOFEM_LOG("SELF", Sev::error)
        << "KSP-backed inverse C C^T action was created for a different "
           "constraint matrix";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Inverse C C^T action is associated with a different C");
  }
  MoFEMFunctionReturn(0);
}

◆ validatePetscObject()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::validatePetscObject	(	PetscObject	object,
		const char *	object_name
	)		const

private

Definition at line 192 of file ConstrainMatrixCtx.cpp.

                                                                           {
  MoFEMFunctionBeginHot;
  if (!object) {
    MOFEM_LOG("SELF", Sev::error)
        << "Null PETSc object supplied for " << object_name;
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY, "Null PETSc object");
  }
 
  MPI_Comm object_comm;
  CHKERR PetscObjectGetComm(object, &object_comm);
  int relation;
  CHKERR MPI_Comm_compare(comm, object_comm, &relation);
  if (relation > MPI_CONGRUENT) {
    MOFEM_LOG("SELF", Sev::error)
        << object_name << " uses a communicator incompatible with the context";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "PETSc object uses an incompatible communicator");
  }
  MoFEMFunctionReturnHot(0);
}

◆ validateSystemMatrixLayout()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::validateSystemMatrixLayout ( ) const

private

Definition at line 290 of file ConstrainMatrixCtx.cpp.

                                                                        {
  MoFEMFunctionBegin;
  CHKERR validatePetscObject(getPetscObject(K.get()), "system matrix K");
  PetscInt local_rows;
  PetscInt local_columns;
  PetscInt global_rows;
  PetscInt global_columns;
  CHKERR MatGetLocalSize(K.get(), &local_rows, &local_columns);
  CHKERR MatGetSize(K.get(), &global_rows, &global_columns);
  if (local_rows != local_columns || global_rows != global_columns) {
    MOFEM_LOG("SELF", Sev::error)
        << "System matrix K must be square, local size (" << local_rows << ", "
        << local_columns << "), global size (" << global_rows << ", "
        << global_columns << ")";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "System matrix K must have matching row and column layouts");
  }
  if (projectionInitialized &&
      (local_rows != fullLocalSize || global_rows != fullGlobalSize)) {
    MOFEM_LOG("SELF", Sev::error)
        << "System matrix K layout (" << local_rows << ", " << global_rows
        << ") does not match the full vector layout (" << fullLocalSize << ", "
        << fullGlobalSize << ")";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "System matrix K is incompatible with the full vector layout");
  }
  MoFEMFunctionReturn(0);
}

◆ validateVectorLayout()

MoFEMErrorCode MoFEM::ConstrainMatrixCtxImpl::validateVectorLayout	(	Vec	vec,
		PetscInt	local_size,
		PetscInt	global_size,
		const char *	vector_name
	)		const

private

Definition at line 215 of file ConstrainMatrixCtx.cpp.

                                                                            {
  MoFEMFunctionBeginHot;
  CHKERR validatePetscObject(getPetscObject(vec), vector_name);
  PetscInt actual_local_size;
  PetscInt actual_global_size;
  CHKERR VecGetLocalSize(vec, &actual_local_size);
  CHKERR VecGetSize(vec, &actual_global_size);
  if (actual_local_size != local_size || actual_global_size != global_size) {
    MOFEM_LOG("SELF", Sev::error)
        << vector_name << " layout is (" << actual_local_size << ", "
        << actual_global_size << "), expected (" << local_size << ", "
        << global_size << ")";
    SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
            "Vector layout is incompatible with the projection context");
  }
  MoFEMFunctionReturnHot(0);
}

Member Data Documentation

◆ buildScatterFromDms

bool MoFEM::ConstrainMatrixCtxImpl::buildScatterFromDms = true

private

Definition at line 174 of file ConstrainMatrixCtx.cpp.

◆ C

SmartPetscObj<Mat> MoFEM::ConstrainMatrixCtxImpl::C

private

Definition at line 158 of file ConstrainMatrixCtx.cpp.

◆ CCTm1Cx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::CCTm1Cx

private

Definition at line 164 of file ConstrainMatrixCtx.cpp.

◆ comm

MPI_Comm MoFEM::ConstrainMatrixCtxImpl::comm

private

Definition at line 154 of file ConstrainMatrixCtx.cpp.

◆ CT

SmartPetscObj<Mat> MoFEM::ConstrainMatrixCtxImpl::CT

private

Definition at line 159 of file ConstrainMatrixCtx.cpp.

◆ CTC

SmartPetscObj<Mat> MoFEM::ConstrainMatrixCtxImpl::CTC

private

Definition at line 161 of file ConstrainMatrixCtx.cpp.

◆ CTCCTm1Cx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::CTCCTm1Cx

private

Definition at line 165 of file ConstrainMatrixCtx.cpp.

◆ CTCx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::CTCx

private

Definition at line 166 of file ConstrainMatrixCtx.cpp.

◆ Cx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::Cx

private

Definition at line 163 of file ConstrainMatrixCtx.cpp.

◆ fullGlobalSize

PetscInt MoFEM::ConstrainMatrixCtxImpl::fullGlobalSize = 0

private

Definition at line 177 of file ConstrainMatrixCtx.cpp.

◆ fullLocalSize

PetscInt MoFEM::ConstrainMatrixCtxImpl::fullLocalSize = 0

private

Definition at line 176 of file ConstrainMatrixCtx.cpp.

◆ inverseCCT

SmartPetscObj<Mat> MoFEM::ConstrainMatrixCtxImpl::inverseCCT

private

Definition at line 160 of file ConstrainMatrixCtx.cpp.

◆ K

SmartPetscObj<Mat> MoFEM::ConstrainMatrixCtxImpl::K

private

Definition at line 162 of file ConstrainMatrixCtx.cpp.

◆ KQx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::KQx

private

Definition at line 169 of file ConstrainMatrixCtx.cpp.

◆ projectionInitialized

bool MoFEM::ConstrainMatrixCtxImpl::projectionInitialized = false

private

Definition at line 172 of file ConstrainMatrixCtx.cpp.

◆ qtkqInitialized

bool MoFEM::ConstrainMatrixCtxImpl::qtkqInitialized = false

private

Definition at line 173 of file ConstrainMatrixCtx.cpp.

◆ Qx

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::Qx

private

Definition at line 168 of file ConstrainMatrixCtx.cpp.

◆ scatter

SmartPetscObj<VecScatter> MoFEM::ConstrainMatrixCtxImpl::scatter

private

Definition at line 170 of file ConstrainMatrixCtx.cpp.

◆ X

SmartPetscObj<Vec> MoFEM::ConstrainMatrixCtxImpl::X

private

Definition at line 167 of file ConstrainMatrixCtx.cpp.

◆ xDm

SmartPetscObj<DM> MoFEM::ConstrainMatrixCtxImpl::xDm

private

Definition at line 155 of file ConstrainMatrixCtx.cpp.

◆ yDm

SmartPetscObj<DM> MoFEM::ConstrainMatrixCtxImpl::yDm

private

Definition at line 156 of file ConstrainMatrixCtx.cpp.

The documentation for this struct was generated from the following file:

src/petsc/impl/ConstrainMatrixCtx.cpp

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ConstrainMatrixCtxImpl() [1/2]

◆ ConstrainMatrixCtxImpl() [2/2]

◆ ~ConstrainMatrixCtxImpl()

Member Function Documentation

◆ applyCTC_QTKQ()

◆ applyInverseCCT()

◆ applyP()

◆ applyQ()

◆ applyQImpl()

◆ applyR()

◆ applyRT()

◆ getConstrainMatrix()

◆ getInverseCCT()

◆ getSystemMatrix()

◆ initializeQorP()

◆ initializeQTKQ()

◆ invalidateProjectionData()

◆ invalidateSystemData()

◆ recalculateCT()

◆ recalculateCTC()

◆ setConstrainMatrix()

◆ setSystemMatrix()

◆ validateInverseCCTLayout()

◆ validatePetscObject()

◆ validateSystemMatrixLayout()

◆ validateVectorLayout()

Member Data Documentation

◆ buildScatterFromDms

◆ C

◆ CCTm1Cx

◆ comm

◆ CT

◆ CTC

◆ CTCCTm1Cx

◆ CTCx

◆ Cx

◆ fullGlobalSize

◆ fullLocalSize

◆ inverseCCT

◆ K

◆ KQx

◆ projectionInitialized

◆ qtkqInitialized

◆ Qx

◆ scatter

◆ X

◆ xDm

◆ yDm