v0.15.0
Loading...
Searching...
No Matches
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
11namespace MoFEM {
12
14 : cOre(const_cast<Core &>(core)) {}
15
17OperatorsTester::query_interface(boost::typeindex::type_index type_index,
18 UnknownInterface **iface) const {
19 *iface = const_cast<OperatorsTester *>(this);
20 return 0;
21}
22
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(
54 const auto hi = prb_ptr->getNumeredRowDofsPtr()->upper_bound(
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
91 boost::shared_ptr<FEMethod> pipeline,
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;
105
106 auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
107 std::ignore = a;
108 std::ignore = v;
109
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
141 boost::shared_ptr<FEMethod> pipeline,
143 SmartPetscObj<Vec> delta2_x, double time,
144 double delta_t, CacheTupleWeakPtr cache_ptr) {
145 auto m = createDMMatrix(dm);
146
147 pipeline->data_ctx =
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;
161
162 auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
163 std::ignore = a;
164 std::ignore = v;
165
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
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) {
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
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
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
265std::pair<SmartPetscObj<Vec>, SmartPetscObj<Vec>>
266OperatorsTester::setPipelineX(boost::shared_ptr<FEMethod> pipeline,
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_NULLPTR;
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_NULLPTR;
306 }
307
308 return std::make_pair(x_t, x_tt);
309}
310} // namespace MoFEM
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
constexpr double a
static const double eps
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define CHKERR
Inline error check.
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:514
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition DMMoFEM.hpp:1047
FTensor::Index< 'm', 3 > m
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Core (interface) class.
Definition Core.hpp:82
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
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.
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fields.
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.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
OperatorsTester(const MoFEM::Core &core)
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.
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.
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.
static constexpr Switches CtxSetA
static constexpr Switches CtxSetF
static constexpr Switches CtxSetB
static constexpr Switches CtxSetTime
intrusive_ptr for managing petsc objects
base class for all interface classes