v0.14.0
Schur.hpp
Go to the documentation of this file.
1 /**
2  * @file Schur.hpp
3  * @brief Assemble Schur complement
4  * @date 2023-02-02
5  *
6  * @copyright Copyright (c) 2023
7  *
8  * To create nested system of Schur complements, you push sequence of operator,
9  * to set up data on entities, and then assemble complements.
10  *
11  * \note Try separate floating points operations from book keeping. Also, align
12  * memory that blocks follow floating point operations.
13  *
14  */
15 
16 #ifndef __SCHUR_HPP__
17 #define __SCHUR_HPP__
18 
19 namespace MoFEM {
20 
21 constexpr const char MoFEM_BLOCK_MAT[] = "mofem_block_mat";
22 
23 /**
24  * @brief Structure to register events for Schur block assembly and solver
25  */
26 struct SchurEvents {
27  static PetscLogEvent MOFEM_EVENT_schurMatSetValues;
28  static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd;
30  static PetscLogEvent MOFEM_EVENT_BlockStructureMult;
31  static PetscLogEvent MOFEM_EVENT_BlockStructureSolve;
32  static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols;
33  SchurEvents();
34 };
35 
36 struct SchurElemMats;
37 struct SchurElemMatsBlock;
39 
41 
42  OpSchurAssembleBase() = delete;
43 
44  using MatSetValuesRaw = boost::function<MoFEMErrorCode(
45  Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n,
46  const PetscInt idxn[], const PetscScalar v[], InsertMode addv)>;
48 
49 private:
51 };
52 
54 
55 /**
56  * @brief Construct a new Op Schur Assemble End object
57  *
58  * @param fields_name list of fields
59  * @param field_ents list of entities on which schur complement is applied
60  * (can be empty)
61  * @param sequence_of_aos list of maps from base problem to Schur complement
62  * matrix
63  * @param sequence_of_mats list of Schur complement matrices
64  * @param sym_schur true if Schur complement is symmetric
65  * @param symm_op true if block diagonal is symmetric
66  */
68  std::vector<std::string> fields_name,
69  std::vector<boost::shared_ptr<Range>> field_ents,
70  std::vector<SmartPetscObj<AO>> sequence_of_aos,
71  std::vector<SmartPetscObj<Mat>> sequence_of_mats,
72  std::vector<bool> sym_schur, bool symm_op,
73  boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
74 
75 /**
76  * @brief Construct a new Op Schur Assemble End object
77  *
78  * @param fields_name list of fields
79  * @param field_ents list of entities on which schur complement is applied
80  * (can be empty)
81  * @param sequence_of_aos list of maps from base problem to Schur complement
82  * matrix
83  * @param sequence_of_mats list of Schur complement matrices
84  * @param sym_schur true if Schur complement is symmetric
85  * @param diag_eps add epsilon on diagonal of inverted matrix
86  * @param symm_op true if block diagonal is symmetric
87  */
89  std::vector<std::string> fields_name,
90  std::vector<boost::shared_ptr<Range>> field_ents,
91  std::vector<SmartPetscObj<AO>> sequence_of_aos,
92  std::vector<SmartPetscObj<Mat>> sequence_of_mats,
93  std::vector<bool> sym_schur, std::vector<double> diag_eps, bool symm_op,
94  boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
95 
96 using SchurFieldPair = std::pair<std::string, std::string>;
97 
98 using SchurFEOpsFEandFields = std::vector<
99 
100  std::pair<std::string, std::vector<SchurFieldPair>>
101 
102  >;
103 
104 /**
105  * @brief Create a Mat Diag Blocks object
106  *
107  * @return Mat
108  */
109 boost::shared_ptr<BlockStructure> createBlockMatStructure(
110 
111  DM dm, //< dm
112  SchurFEOpsFEandFields schur_fe_op_vec //< block elements
113 
114 );
115 
116 using SchurShellMatData =
117  std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>;
118 
119 /**
120  * @brief Create a Schur Mat object
121  *
122  * @param dm
123  * @param data
124  * @return std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>
125  */
126 SchurShellMatData createBlockMat(DM dm, boost::shared_ptr<BlockStructure> data);
127 
128 using NestSchurData = std::tuple<
129 
130  std::array<SmartPetscObj<Mat>, 4>,
131  std::array<boost::shared_ptr<BlockStructure>, 4>,
132  boost::shared_ptr<BlockStructure>,
133  std::pair<SmartPetscObj<IS>, SmartPetscObj<IS>>
134 
135  >;
136 
137 /**
138  * @brief Get the Schur Nest Mat Array object
139  *
140  * @param dms schur dm, and block dm
141  * @param block mat A data
142  * @param fields_name list of fields
143  * @param field_ents list of entities on which schur complement is applied
144  * @param add_preconditioner_block add block for preconditioner
145  * @return boost::shared_ptr<NestSchurData>
146  */
147 boost::shared_ptr<NestSchurData>
149  boost::shared_ptr<BlockStructure> block_mat_data,
150 
151  std::vector<std::string> fields_name, //< a00 fields
152  std::vector<boost::shared_ptr<Range>> field_ents, //< a00 ents
153  bool add_preconditioner_block = false
154 
155 );
156 
157 /**
158  * @brief Switch preconditioner
159  *
160  * @param block_mat_data
161  * @return MoFEMErrorCode
162  */
164 schurSwitchPreconditioner(boost::shared_ptr<BlockStructure> block_mat_data);
165 
166 /**
167  * @brief Save block matrix as a mesh
168  *
169  * @param block_mat_data
170  * @param filename
171  * @return MoFEMErrorCode
172  */
174 schurSaveBlockMesh(boost::shared_ptr<BlockStructure> block_mat_data,
175  std::string filename);
176 
177 /**
178  * @brief Create a Mat Diag Blocks object
179  *
180  * \code {.cpp}
181  *
182  * auto [nested_mat, nested_data_ptr] = createSchurNestedMatrix(
183  *
184  * getNestSchurData(
185  *
186  * {schur_dm, block_dm}, shell_data,
187  *
188  * {"TENSOR"}, {nullptr}
189  *
190  * )
191  *
192  * );
193  *
194  * \endcode
195  *
196  * @return Mat
197  */
198 std::pair<SmartPetscObj<Mat>, boost::shared_ptr<NestSchurData>>
199 createSchurNestedMatrix(boost::shared_ptr<NestSchurData> schur_net_data_ptr);
200 
201 template <>
202 MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr<NestSchurData>);
203 
204 /**
205  * @brief Set PC for Schur block
206  *
207  * @param pc
208  * @return MoFEMErrorCode
209  */
211 
213  SchurBackendMatSetValuesPtr() = delete;
214  using MatSetValuesPtr = boost::function<MoFEMErrorCode(
215  Mat M, const EntitiesFieldData::EntData &row_data,
216  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
217  InsertMode iora)>;
218  static MatSetValuesPtr matSetValuesPtr; ///< backend assembly function
219  static MatSetValuesPtr
220  matSetValuesBlockPtr; ///< backend assembly block mat function
221  static MatSetValuesPtr
222  matSetValuesPreconditionedBlockPtr; ///< backend assembly block
223  ///< preconditioner mat function
224 };
225 
226 template <>
229  const EntitiesFieldData::EntData &col_data,
230  const MatrixDouble &mat, InsertMode iora);
231 
232 template <>
233 inline MoFEMErrorCode
235  const VectorDouble &nf, InsertMode iora) {
236  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
237 }
238 
239 template <>
240 MoFEMErrorCode MatSetValues<AssemblyTypeSelector<SCHUR>>(
241  Mat M, const EntitiesFieldData::EntData &row_data,
242  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
243  InsertMode iora);
244 
245 template <>
246 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<SCHUR>>(
247  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
248  InsertMode iora) {
249  return VecSetValues<SchurElemMats>(V, data, nf, iora);
250 }
251 
252 /***
253  * @brief Specialization of MatSetValues for BlockStructure
254  */
255 template <>
258  const EntitiesFieldData::EntData &col_data,
259  const MatrixDouble &mat, InsertMode iora);
260 
261 /***
262  * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_MAT>
263  */
264 template <>
265 inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
266  Mat M, const EntitiesFieldData::EntData &row_data,
267  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
268  InsertMode iora) {
269  return MatSetValues<BlockStructure>(M, row_data, col_data, mat, iora);
270 }
271 
272 /***
273  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_MAT>
274  */
275 template <>
276 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
277  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
278  InsertMode iora) {
279  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
280 }
281 
282 /***
283  * @brief Specialization of MatSetValues for SchurElemMatsBlock
284  */
285 template <>
288  const EntitiesFieldData::EntData &row_data,
289  const EntitiesFieldData::EntData &col_data,
290  const MatrixDouble &mat, InsertMode iora);
291 
292 /***
293  * @brief Specialization of VecSetValues for SchurElemMatsBlock
294  */
295 template <>
296 inline MoFEMErrorCode
298  const VectorDouble &nf, InsertMode iora) {
299  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
300 }
301 
302 /***
303  * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
304  */
305 template <>
306 inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
307  Mat M, const EntitiesFieldData::EntData &row_data,
308  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
309  InsertMode iora) {
310  return MatSetValues<SchurElemMatsBlock>(M, row_data, col_data, mat, iora);
311 }
312 
313 /***
314  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
315  */
316 template <>
317 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
318  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
319  InsertMode iora) {
320  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
321 }
322 
323 /***
324  * @brief Specialization of MatSetValues for SchurElemMatsPreconditionedBlock
325  */
326 template <>
328  Mat M, const EntitiesFieldData::EntData &row_data,
329  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
330  InsertMode iora);
331 
332 /***
333  * @brief Specialisation of MatSetValues for
334  * AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>
335  */
336 template <>
337 inline MoFEMErrorCode
338 MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
339  Mat M, const EntitiesFieldData::EntData &row_data,
340  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
341  InsertMode iora) {
342  return MatSetValues<SchurElemMatsPreconditionedBlock>(M, row_data, col_data,
343  mat, iora);
344 }
345 
346 /***
347  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
348  */
349 template <>
350 inline MoFEMErrorCode
351 VecSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
352  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
353  InsertMode iora) {
354  return ::VecSetValues(V, data.getIndices().size(),
355  &*data.getIndices().begin(), &*nf.begin(), iora);
356 }
357 
358 } // namespace MoFEM
359 
360 #endif //__SCHUR_HPP__
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::MatSetValues< BlockStructure >
MoFEMErrorCode MatSetValues< BlockStructure >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2548
MoFEM::setSchurMatSolvePC
MoFEMErrorCode setSchurMatSolvePC(SmartPetscObj< PC > pc)
Set PC for Schur block.
Definition: Schur.cpp:2476
MoFEM::MoFEM_BLOCK_MAT
constexpr const char MoFEM_BLOCK_MAT[]
Definition: Schur.hpp:21
MoFEM::SchurFieldPair
std::pair< std::string, std::string > SchurFieldPair
Definition: Schur.hpp:96
MoFEM::OpSchurAssembleBase::OpSchurAssembleBase
OpSchurAssembleBase()=delete
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::OpSchurAssembleBase::MatSetValuesRaw
boost::function< MoFEMErrorCode(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar v[], InsertMode addv)> MatSetValuesRaw
Definition: Schur.hpp:46
MoFEM::SchurBackendMatSetValuesPtr::MatSetValuesPtr
boost::function< MoFEMErrorCode(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)> MatSetValuesPtr
Definition: Schur.hpp:217
MoFEM::SchurEvents::SchurEvents
SchurEvents()
Definition: Schur.cpp:297
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureMult
static PetscLogEvent MOFEM_EVENT_BlockStructureMult
Definition: Schur.hpp:30
MoFEM::SchurEvents::MOFEM_EVENT_schurMatSetValues
static PetscLogEvent MOFEM_EVENT_schurMatSetValues
Definition: Schur.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::VecSetValues< SchurElemMatsBlock >
MoFEMErrorCode VecSetValues< SchurElemMatsBlock >(Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf, InsertMode iora)
Definition: Schur.hpp:297
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, std::vector< SmartPetscObj< AO >> sequence_of_aos, std::vector< SmartPetscObj< Mat >> sequence_of_mats, std::vector< bool > sym_schur, bool symm_op, boost::shared_ptr< BlockStructure > diag_blocks)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2442
MoFEM::createBlockMat
SchurShellMatData createBlockMat(DM dm, boost::shared_ptr< BlockStructure > data)
Create a Schur Mat object.
Definition: Schur.cpp:1387
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureSolve
static PetscLogEvent MOFEM_EVENT_BlockStructureSolve
Definition: Schur.hpp:31
MoFEM::schurSwitchPreconditioner
MoFEMErrorCode schurSwitchPreconditioner(boost::shared_ptr< BlockStructure > block_mat_data)
Switch preconditioner.
Definition: Schur.cpp:2620
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::SchurBackendMatSetValuesPtr::SchurBackendMatSetValuesPtr
SchurBackendMatSetValuesPtr()=delete
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2437
MoFEM::SchurEvents::MOFEM_EVENT_opSchurAssembleEnd
static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd
Definition: Schur.hpp:28
MoFEM::SchurElemMatsBlock
Definition: Schur.cpp:2554
MoFEM::createSchurNestedMatrix
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< NestSchurData > > createSchurNestedMatrix(boost::shared_ptr< NestSchurData > schur_net_data_ptr)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:2409
MoFEM::MatSetValues< SchurElemMatsPreconditionedBlock >
MoFEMErrorCode MatSetValues< SchurElemMatsPreconditionedBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2611
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesPtr
static MatSetValuesPtr matSetValuesPtr
backend assembly function
Definition: Schur.hpp:218
MoFEM::SchurElemMatsPreconditionedBlock
Definition: Schur.cpp:2587
MoFEM::SchurBackendMatSetValuesPtr
Definition: Schur.hpp:212
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
MoFEM::OpSchurAssembleBase
Definition: Schur.hpp:40
convert.n
n
Definition: convert.py:82
MoFEM::SchurElemMats
Schur complement data storage.
Definition: Schur.cpp:122
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesPreconditionedBlockPtr
static MatSetValuesPtr matSetValuesPreconditionedBlockPtr
Definition: Schur.hpp:222
MoFEM::MatSetValues< SchurElemMats >
MoFEMErrorCode MatSetValues< SchurElemMats >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2504
MoFEM::OpSchurAssembleBase::matSetValuesSchurRaw
static MatSetValuesRaw matSetValuesSchurRaw
Definition: Schur.hpp:47
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::SchurEvents
Structure to register events for Schur block assembly and solver.
Definition: Schur.hpp:26
MoFEM::SchurShellMatData
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< BlockStructure > > SchurShellMatData
Definition: Schur.hpp:117
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1551
MoFEM::VecSetValues< SchurElemMats >
MoFEMErrorCode VecSetValues< SchurElemMats >(Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf, InsertMode iora)
Definition: Schur.hpp:234
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1084
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesBlockPtr
static MatSetValuesPtr matSetValuesBlockPtr
backend assembly block mat function
Definition: Schur.hpp:220
MoFEM::SmartPetscObj< AO >
MoFEM::MatSetValues< SchurElemMatsBlock >
MoFEMErrorCode MatSetValues< SchurElemMatsBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2580
MoFEM::NestSchurData
std::tuple< std::array< SmartPetscObj< Mat >, 4 >, std::array< boost::shared_ptr< BlockStructure >, 4 >, boost::shared_ptr< BlockStructure >, std::pair< SmartPetscObj< IS >, SmartPetscObj< IS > > > NestSchurData
Definition: Schur.hpp:135
MoFEM::SchurEvents::MOFEM_EVENT_zeroRowsAndCols
static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols
Definition: Schur.hpp:32
MoFEM::schurSaveBlockMesh
MoFEMErrorCode schurSaveBlockMesh(boost::shared_ptr< BlockStructure > block_mat_data, std::string filename)
Save block matrix as a mesh.
Definition: Schur.cpp:2635
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureSetValues
static PetscLogEvent MOFEM_EVENT_BlockStructureSetValues
Definition: Schur.hpp:29
MoFEM::SchurFEOpsFEandFields
std::vector< std::pair< std::string, std::vector< SchurFieldPair > > > SchurFEOpsFEandFields
Definition: Schur.hpp:102
MoFEM::VecSetValues< EssentialBcStorage >
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
Definition: FormsIntegrators.cpp:76
MoFEM::getNestSchurData
boost::shared_ptr< NestSchurData > getNestSchurData(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:1991
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1942