v0.14.0
OperatorsTester.cpp
Go to the documentation of this file.
1 /**
2  * @file OperatorsTester.cpp
3  * @brief
4  * @version 0.13.2
5  * @date 2022-12-04
6  *
7  * @copyright Copyright (c) 2022
8  *
9  */
10 
11 namespace MoFEM {
12 
13 OperatorsTester::OperatorsTester(const Core &core)
14  : cOre(const_cast<Core &>(core)) {}
15 
16 MoFEMErrorCode
17 OperatorsTester::query_interface(boost::typeindex::type_index type_index,
18  UnknownInterface **iface) const {
19  *iface = const_cast<OperatorsTester *>(this);
20  return 0;
21 }
22 
23 SmartPetscObj<Vec>
24 OperatorsTester::setRandomFields(SmartPetscObj<DM> dm,
25  std::vector<RandomFieldData> random_fields,
26  boost::shared_ptr<Range> ents_ptr) {
27 
28  MoFEM::Interface &m_field = cOre;
29 
30  auto prb_ptr = getProblemPtr(dm);
31 
32  auto get_random_number = [](auto &range) {
33  const auto r = static_cast<double>(std::rand()) / RAND_MAX;
34  return range[0] + r * (range[1] - range[0]);
35  };
36 
37  auto v = createDMVector(dm);
38  double *array;
39  CHKERR VecGetArray(v, &array);
40 
41  int size;
42  VecGetLocalSize(v, &size);
43 
44  MOFEM_LOG_CHANNEL("WORLD");
45 
46  for (auto &rf : random_fields) {
47  const auto field_id = m_field.get_field_bit_number(rf.first);
48  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "OperatorsTester")
49  << "Set random field " << rf.first << " " << static_cast<int>(field_id);
50 
51  if (!ents_ptr) {
52  auto lo = prb_ptr->getNumeredRowDofsPtr()->lower_bound(
53  FieldEntity::getLoBitNumberUId(field_id));
54  const auto hi = prb_ptr->getNumeredRowDofsPtr()->upper_bound(
55  FieldEntity::getHiBitNumberUId(field_id));
56  for (; lo != hi; ++lo) {
57  const auto idx = (*lo)->getPetscLocalDofIdx();
58  if (idx >= 0 && idx < size) {
59  array[idx] = get_random_number(rf.second);
60  }
61  }
62  } else {
63  for (auto pit = ents_ptr->const_pair_begin();
64  pit != ents_ptr->const_pair_end(); ++pit) {
65  auto lo =
66  prb_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>().lower_bound(
67  DofEntity::getLoFieldEntityUId(field_id, pit->first));
68  auto hi =
69  prb_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>().upper_bound(
70  DofEntity::getHiFieldEntityUId(field_id, pit->second));
71  for (; lo != hi; ++lo) {
72  const auto idx = (*lo)->getPetscLocalDofIdx();
73  if (idx >= 0 && idx < size) {
74  array[idx] = get_random_number(rf.second);
75  }
76  }
77  }
78  }
79  }
80 
81  CHKERR VecRestoreArray(v, &array);
82 
83  CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
84  CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
85 
86  return v;
87 };
88 
89 SmartPetscObj<Vec>
90 OperatorsTester::assembleVec(SmartPetscObj<DM> dm, std::string fe_name,
91  boost::shared_ptr<FEMethod> pipeline,
92  SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
93  SmartPetscObj<Vec> delta2_x, double time,
94  double delta_t, CacheTupleWeakPtr cache_ptr) {
95 
96  pipeline->data_ctx = PetscData::CtxSetF | PetscData::CtxSetTime;
97  auto f = createDMVector(dm);
98  pipeline->f = f;
99  pipeline->ts_t = time;
100  pipeline->ts_dt = delta_t;
101 
102  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_SETFUNCTION;
103  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETFUNCTION;
104  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIFUNCTION;
105 
106  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
107  std::ignore = a;
108  std::ignore = v;
109 
110  CHK_THROW_MESSAGE(
111  DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
112  "loc to vec");
113  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
114  "Run assemble pipeline");
115  CHK_THROW_MESSAGE(VecAssemblyBegin(f), "Asmb vector");
116  CHK_THROW_MESSAGE(VecAssemblyEnd(f), "Asmb vector");
117 
118  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE),
119  "Update vec");
120  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE),
121  "Update vec");
122  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD),
123  "Update vec");
124  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD),
125  "Update vec");
126 
127  auto get_norm = [](auto f) {
128  double fnorm;
129  VecNorm(f, NORM_2, &fnorm);
130  return fnorm;
131  };
132  MOFEM_LOG_CHANNEL("WORLD");
133  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Vec norm %3.4e",
134  get_norm(f));
135 
136  return f;
137 }
138 
139 SmartPetscObj<Mat>
140 OperatorsTester::assembleMat(SmartPetscObj<DM> dm, std::string fe_name,
141  boost::shared_ptr<FEMethod> pipeline,
142  SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
143  SmartPetscObj<Vec> delta2_x, double time,
144  double delta_t, CacheTupleWeakPtr cache_ptr) {
145  auto m = createDMMatrix(dm);
146 
147  pipeline->data_ctx =
148  PetscData::CtxSetA | PetscData::CtxSetB | PetscData::CtxSetTime;
149  pipeline->A = m;
150  pipeline->B = m;
151 
152  pipeline->ts_t = time;
153 
154  pipeline->ts_dt = delta_t;
155  pipeline->ts_a = 1. / delta_t;
156  pipeline->ts_aa = 1. / pow(delta_t, 2);
157 
158  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_OPERATORS;
159  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETJACOBIAN;
160  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIJACOBIAN;
161 
162  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
163  std::ignore = a;
164  std::ignore = v;
165 
166  CHK_THROW_MESSAGE(
167  DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
168  "loc to vec");
169  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
170  "Run assemble pipeline");
171  CHK_THROW_MESSAGE(MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
172  CHK_THROW_MESSAGE(MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
173 
174  auto get_norm = [](auto m) {
175  double fnorm;
176  MatNorm(m, NORM_1, &fnorm);
177  return fnorm;
178  };
179  MOFEM_LOG_CHANNEL("WORLD");
180  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Matrix norm %3.4e",
181  get_norm(m));
182 
183  return m;
184 }
185 
186 SmartPetscObj<Vec> OperatorsTester::directionalCentralFiniteDifference(
187  SmartPetscObj<DM> dm, std::string fe_name,
188  boost::shared_ptr<FEMethod> pipeline, SmartPetscObj<Vec> x,
189  SmartPetscObj<Vec> delta_x, SmartPetscObj<Vec> delta2_x,
190  SmartPetscObj<Vec> diff_x, double time, double delta_t, double eps,
191  CacheTupleWeakPtr cache_ptr) {
192 
193  auto axpy = [&](auto v0, auto diff_v, auto eps) {
194  SmartPetscObj<Vec> v;
195  if (v0.use_count()) {
196  v = vectorDuplicate(v0);
197  CHK_THROW_MESSAGE(VecCopy(v0, v), "Cpy");
198  CHK_THROW_MESSAGE(VecAXPY(v, eps, diff_v), "Add");
199  CHK_THROW_MESSAGE(VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD),
200  "ghost");
201  CHK_THROW_MESSAGE(VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD),
202  "ghost");
203  }
204  return v;
205  };
206 
207  auto fm1 = assembleVec(dm, fe_name, pipeline,
208  //
209  axpy(x, diff_x, -eps),
210  //
211  axpy(delta_x, diff_x, -eps),
212  //
213  axpy(delta2_x, diff_x, -eps),
214  //
215  time, delta_t, cache_ptr);
216 
217  auto fp1 = assembleVec(dm, fe_name, pipeline,
218  //
219  axpy(x, diff_x, eps),
220  //
221  axpy(delta_x, diff_x, eps),
222  //
223  axpy(delta2_x, diff_x, eps),
224  //
225  time, delta_t, cache_ptr);
226 
227  auto calculate_derivative = [eps](auto fm1, auto fp1) {
228  int size;
229  VecGetLocalSize(fp1, &size);
230  VecAXPY(fp1, -1, fm1);
231  double *array;
232  VecGetArray(fp1, &array);
233  for (int j = 0; j != size; ++j) {
234  array[j] /= 2 * eps;
235  }
236  VecRestoreArray(fp1, &array);
237  return fp1;
238  };
239 
240  return calculate_derivative(fm1, fp1);
241 }
242 
243 SmartPetscObj<Vec> OperatorsTester::checkCentralFiniteDifference(
244  SmartPetscObj<DM> dm, std::string fe_name,
245  boost::shared_ptr<FEMethod> pipeline_rhs,
246  boost::shared_ptr<FEMethod> pipeline_lhs, SmartPetscObj<Vec> x,
247  SmartPetscObj<Vec> delta_x, SmartPetscObj<Vec> delta2_x,
248  SmartPetscObj<Vec> diff_x, double time, double delta_t, double eps,
249  CacheTupleWeakPtr cache_ptr) {
250 
251  auto fd_diff = directionalCentralFiniteDifference(
252  dm, fe_name, pipeline_rhs, x, delta_x, delta2_x, diff_x, time, delta_t,
253  eps, cache_ptr);
254 
255  auto m = assembleMat(dm, fe_name, pipeline_lhs, x, delta_x, delta2_x, time,
256  delta_t, cache_ptr);
257 
258  auto fm = vectorDuplicate(fd_diff);
259  CHK_THROW_MESSAGE(MatMult(m, diff_x, fm), "Mat mult");
260  CHK_THROW_MESSAGE(VecAXPY(fd_diff, -1., fm), "Add");
261 
262  return fd_diff;
263 }
264 
265 std::pair<SmartPetscObj<Vec>, SmartPetscObj<Vec>>
266 OperatorsTester::setPipelineX(boost::shared_ptr<FEMethod> pipeline,
267  SmartPetscObj<Vec> x, SmartPetscObj<Vec> delta_x,
268  SmartPetscObj<Vec> delta2_x, double delta_t) {
269 
270  // Set dofs vector to finite element instance
271  pipeline->x = x;
272  pipeline->data_ctx |= PetscData::CTX_SET_X;
273 
274  SmartPetscObj<Vec> x_t, x_tt;
275 
276  // Set velocity dofs vector to finite element instance
277 
278  if (delta_x.use_count()) {
279  x_t = vectorDuplicate(x);
280  VecCopy(delta_x, x_t);
281  VecScale(x_t, 1. / delta_t);
282  CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_t, INSERT_VALUES, SCATTER_FORWARD),
283  "ghost");
284  CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_t, INSERT_VALUES, SCATTER_FORWARD),
285  "ghost");
286  pipeline->x_t = x_t;
287  pipeline->data_ctx |= PetscData::CTX_SET_X_T;
288  } else {
289  pipeline->x_t = PETSC_NULL;
290  }
291 
292  // Set acceleration dofs vector to finite element instance
293 
294  if (delta2_x.use_count()) {
295  x_tt = vectorDuplicate(x);
296  VecCopy(delta2_x, x_tt);
297  VecScale(x_tt, 1. / pow(delta_t, 2));
298  CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_tt, INSERT_VALUES, SCATTER_FORWARD),
299  "ghost");
300  CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_tt, INSERT_VALUES, SCATTER_FORWARD),
301  "ghost");
302  pipeline->x_tt = x_tt;
303  pipeline->data_ctx |= PetscData::CTX_SET_X_TT;
304  } else {
305  pipeline->x_tt = PETSC_NULL;
306  }
307 
308  return std::make_pair(x_t, x_tt);
309 }
310 } // namespace MoFEM
MoFEM::OperatorsTester::OperatorsTester
OperatorsTester(const MoFEM::Core &core)
Definition: OperatorsTester.cpp:13
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::PetscData::CtxSetB
static constexpr Switches CtxSetB
Definition: LoopMethods.hpp:38
MoFEM::PetscData::CtxSetA
static constexpr Switches CtxSetA
Definition: LoopMethods.hpp:37
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OperatorsTester::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: OperatorsTester.cpp:17
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::OperatorsTester::setPipelineX
std::pair< SmartPetscObj< Vec >, SmartPetscObj< Vec > > setPipelineX(boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double delta_t)
Set vectors, x, x_t, and x_tt to finite element instance.
Definition: OperatorsTester.cpp:266
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::PetscData::CTX_SET_X_TT
@ CTX_SET_X_TT
Definition: LoopMethods.hpp:29
MoFEM::OperatorsTester::assembleMat
SmartPetscObj< Mat > assembleMat(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
Assemble the left hand side vector.
Definition: OperatorsTester.cpp:140
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::SmartPetscObj::use_count
int use_count() const
Definition: PetscSmartObj.hpp:109
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::OperatorsTester::cOre
MoFEM::Core & cOre
Definition: OperatorsTester.hpp:145
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
MoFEM::PetscData::CtxSetTime
static constexpr Switches CtxSetTime
Definition: LoopMethods.hpp:42
MoFEM::OperatorsTester::directionalCentralFiniteDifference
SmartPetscObj< Vec > directionalCentralFiniteDifference(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Calculate directional directive using finite difference.
Definition: OperatorsTester.cpp:186
MoFEM::DofEntity::getLoFieldEntityUId
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:69
MoFEM::OperatorsTester::setRandomFields
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fields.
Definition: OperatorsTester.cpp:24
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
MoFEM::OperatorsTester::assembleVec
SmartPetscObj< Vec > assembleVec(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
Assemble the right hand side vector.
Definition: OperatorsTester.cpp:90
MoFEM::DofEntity::getHiFieldEntityUId
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:75
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:1044
MoFEM::CacheTupleWeakPtr
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
Definition: FEMultiIndices.hpp:494
MoFEM::PetscData::CTX_SET_X_T
@ CTX_SET_X_T
Definition: LoopMethods.hpp:28
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::PetscData::CtxSetF
static constexpr Switches CtxSetF
Definition: LoopMethods.hpp:36
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::SmartPetscObj< Vec >
MoFEM::PetscData::CTX_SET_X
@ CTX_SET_X
Definition: LoopMethods.hpp:27
MoFEM::OperatorsTester::checkCentralFiniteDifference
SmartPetscObj< Vec > checkCentralFiniteDifference(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline_rhs, boost::shared_ptr< FEMethod > pipeline_lhs, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Check consistency between directional derivative with matrix.
Definition: OperatorsTester.cpp:243
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
Generated by Doxygen 1.8.17 and hosted at University of Glasgow