Loading...
Searching...
No Matches
ITHACAerror.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12
13License
14 This file is part of ITHACA-FV
15
16 ITHACA-FV is free software: you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 ITHACA-FV is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
28\*---------------------------------------------------------------------------*/
29
30#include "ITHACAerror.H"
31
34
35// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
36
37namespace ITHACAutilities
38{
39
40template<class Type, template<class> class PatchField, class GeoMesh>
41double errorFrobRel(GeometricField<Type, PatchField, GeoMesh>& field1,
42 GeometricField<Type, PatchField, GeoMesh>& field2, List<label>* labels)
43{
44 double err;
45 autoPtr<GeometricField<Type, PatchField, GeoMesh >> errField;
46 autoPtr<GeometricField<Type, PatchField, GeoMesh >> field1_S;
47 autoPtr<GeometricField<Type, PatchField, GeoMesh >> field2_S;
48 autoPtr<fvMeshSubset> submesh;
49
50 if (labels != NULL)
51 {
52 submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(field1.mesh()));
53#if OPENFOAM >= 1812
54 submesh->setCellSubset(* labels);
55#else
56 submesh->setLargeCellSubset(* labels);
57#endif
58 GeometricField<Type, PatchField, GeoMesh> field1tmp(submesh->interpolate(
59 field1));
60 GeometricField<Type, PatchField, GeoMesh> field2tmp(submesh->interpolate(
61 field2));
62 field1_S = autoPtr<GeometricField<Type, PatchField, GeoMesh >>
63 (new GeometricField<Type, PatchField, GeoMesh>(field1tmp.clone()));
64 field2_S = autoPtr<GeometricField<Type, PatchField, GeoMesh >>
65 (new GeometricField<Type, PatchField, GeoMesh>(field2tmp.clone()));
66 }
67 else
68 {
69 field1_S = autoPtr<GeometricField<Type, PatchField, GeoMesh >>
70 (new GeometricField<Type, PatchField, GeoMesh>(field1));
71 field2_S = autoPtr<GeometricField<Type, PatchField, GeoMesh >>
72 (new GeometricField<Type, PatchField, GeoMesh>(field2));
73 }
74
75 errField = autoPtr<GeometricField<Type, PatchField, GeoMesh >>
76 (new GeometricField<Type, PatchField, GeoMesh>(field1_S() - field2_S()));
77
78 if (frobNorm(field1) <= 1e-6)
79 {
80 err = 0;
81 }
82 else
83 {
84 err = frobNorm(errField()) / frobNorm(field1_S());
85 }
86
87 return err;
88}
89
90
91template double errorFrobRel(GeometricField<scalar, fvPatchField, volMesh>&
92 field1,
93 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
94template double errorFrobRel(GeometricField<vector, fvPatchField, volMesh>&
95 field1,
96 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
97
98template double errorFrobRel(GeometricField<scalar, fvsPatchField, surfaceMesh>
99 &
100 field1,
101 GeometricField<scalar, fvsPatchField, surfaceMesh>& field2,
102 List<label>* labels);
103
104
105template<typename T>
106double errorLinfRel(GeometricField<T, fvPatchField, volMesh>& field1,
107 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
108{
109 double err;
110 autoPtr<GeometricField<T, fvPatchField, volMesh >> errField;
111 autoPtr<GeometricField<T, fvPatchField, volMesh >> field1_S;
112 autoPtr<GeometricField<T, fvPatchField, volMesh >> field2_S;
113 autoPtr<fvMeshSubset> submesh;
114
115 if (labels != NULL)
116 {
117 submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(field1.mesh()));
118#if OPENFOAM >= 1812
119 submesh->setCellSubset(* labels);
120#else
121 submesh->setLargeCellSubset(* labels);
122#endif
123 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
124 field1));
125 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
126 field2));
127 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
128 (new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
129 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
130 (new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
131 }
132 else
133 {
134 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
135 (new GeometricField<T, fvPatchField, volMesh>(field1));
136 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
137 (new GeometricField<T, fvPatchField, volMesh>(field2));
138 }
139
140 errField = autoPtr<GeometricField<T, fvPatchField, volMesh >>
141 (new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
142
143 if (LinfNorm(field1_S()) <= 1e-6)
144 {
145 err = 0;
146 }
147 else
148 {
149 err = LinfNorm(errField()) / LinfNorm(field1_S());
150 }
151
152 return err;
153}
154
155
156template double errorLinfRel(GeometricField<scalar, fvPatchField, volMesh>&
157 field1,
158 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
159template double errorLinfRel(GeometricField<vector, fvPatchField, volMesh>&
160 field1,
161 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
162
163template<>
164double errorL2Abs(GeometricField<vector, fvPatchField, volMesh>& field1,
165 GeometricField<vector, fvPatchField, volMesh>& field2, volScalarField& Volumes)
166
167{
168 volScalarField diffFields2 = (((field1 - field2) & (field1 - field2)) *
169 Volumes).ref();
170 double err = Foam::sqrt(gSum(diffFields2));
171 return err;
172}
173
174template<>
175double errorL2Abs(GeometricField<scalar, fvPatchField, volMesh>& field1,
176 GeometricField<scalar, fvPatchField, volMesh>& field2, volScalarField& Volumes)
177
178{
179 volScalarField diffFields2 = (((field1 - field2) * (field1 - field2)) *
180 Volumes).ref();
181 double err = Foam::sqrt(gSum(diffFields2));
182 return err;
183}
184
185template<typename T>
186double errorL2Abs(GeometricField<T, fvPatchField, volMesh>& field1,
187 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
188{
189 autoPtr<GeometricField<T, fvPatchField, volMesh >> errField;
190 autoPtr<GeometricField<T, fvPatchField, volMesh >> field1_S;
191 autoPtr<GeometricField<T, fvPatchField, volMesh >> field2_S;
192 autoPtr<fvMeshSubset> submesh;
193
194 if (labels != NULL)
195 {
196 submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(field1.mesh()));
197#if OPENFOAM >= 1812
198 submesh->setCellSubset(* labels);
199#else
200 submesh->setLargeCellSubset(* labels);
201#endif
202 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
203 field1));
204 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
205 field2));
206 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
207 (new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
208 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
209 (new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
210 }
211 else
212 {
213 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
214 (new GeometricField<T, fvPatchField, volMesh>(field1));
215 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
216 (new GeometricField<T, fvPatchField, volMesh>(field2));
217 }
218
219 errField = autoPtr<GeometricField<T, fvPatchField, volMesh >>
220 (new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
221 double err = L2Norm(errField());
222 return err;
223}
224
225template double errorL2Abs(GeometricField<scalar, fvPatchField, volMesh>&
226 field1,
227 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
228template double errorL2Abs(GeometricField<vector, fvPatchField, volMesh>&
229 field1,
230 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
231
232
233template<class T, template<class> class PatchField, class GeoMesh>
234Eigen::MatrixXd errorFrobRel(PtrList<GeometricField<T, PatchField, GeoMesh >>&
235 fields1,
236 PtrList<GeometricField<T, PatchField, GeoMesh >>& fields2, List<label>* labels)
237{
238 Eigen::VectorXd err;
239
240 if (fields1.size() != fields2.size())
241 {
242 Info << "The two fields do not have the same size, code will abort" << endl;
243 exit(0);
244 }
245
246 err.resize(fields1.size(), 1);
247
248 for (label k = 0; k < fields1.size(); k++)
249 {
250 err(k, 0) = errorFrobRel(fields1[k], fields2[k], labels);
251 Info << " Error is " << err[k] << endl;
252 }
253
254 return err;
255}
256
257template Eigen::MatrixXd errorFrobRel(
258 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields1,
259 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields2,
260 List<label>* labels);
261template Eigen::MatrixXd errorFrobRel(
262 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields1,
263 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields2,
264 List<label>* labels);
265template Eigen::MatrixXd errorFrobRel(
266 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& fields1,
267 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& fields2,
268 List<label>* labels);
269
270
271template<typename T>
272Eigen::MatrixXd errorL2Abs(
273 PtrList<GeometricField<T, fvPatchField, volMesh >>& fields1,
274 PtrList<GeometricField<T, fvPatchField, volMesh >>& fields2,
275 PtrList<volScalarField>& Volumes)
276{
277 M_Assert(fields1.size() == fields2.size(),
278 "The two fields do not have the same size, code will abort");
279 M_Assert(fields1.size() == Volumes.size(),
280 "The volumes field and the two solution fields do not have the same size, code will abort");
281 Eigen::VectorXd err;
282 err.resize(fields1.size(), 1);
283
284 for (label k = 0; k < fields1.size(); k++)
285 {
286 err(k, 0) = errorL2Abs(fields1[k], fields2[k], Volumes[k]);
287 Info << " Error is " << err[k] << endl;
288 }
289
290 return err;
291}
292
293template Eigen::MatrixXd errorL2Abs(
294 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields1,
295 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields2,
296 PtrList<volScalarField>& Volumes);
297template Eigen::MatrixXd errorL2Abs(
298 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields1,
299 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields2,
300 PtrList<volScalarField>& Volumes);
301
302template<typename T>
303Eigen::MatrixXd errorL2Abs(PtrList<GeometricField<T, fvPatchField, volMesh >>&
304 fields1,
305 PtrList<GeometricField<T, fvPatchField, volMesh >>& fields2,
306 List<label>* labels)
307{
308 Eigen::VectorXd err;
309
310 if (fields1.size() != fields2.size())
311 {
312 Info << "The two fields do not have the same size, code will abort" << endl;
313 exit(0);
314 }
315
316 err.resize(fields1.size(), 1);
317
318 for (label k = 0; k < fields1.size(); k++)
319 {
320 err(k, 0) = errorL2Abs(fields1[k], fields2[k], labels);
321 Info << " Error is " << err[k] << endl;
322 }
323
324 return err;
325}
326
327template Eigen::MatrixXd errorL2Abs(
328 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields1,
329 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields2,
330 List<label>* labels);
331template Eigen::MatrixXd errorL2Abs(
332 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields1,
333 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields2,
334 List<label>* labels);
335
336template<typename T>
337double errorL2Rel(GeometricField<T, fvPatchField, volMesh>& field1,
338 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
339{
340 double err;
341 autoPtr<GeometricField<T, fvPatchField, volMesh >> errField;
342 autoPtr<GeometricField<T, fvPatchField, volMesh >> field1_S;
343 autoPtr<GeometricField<T, fvPatchField, volMesh >> field2_S;
344 autoPtr<fvMeshSubset> submesh;
345
346 if (labels != NULL)
347 {
348 submesh = autoPtr<fvMeshSubset>(new fvMeshSubset(field1.mesh()));
349#if OPENFOAM >= 1812
350 submesh->setCellSubset(* labels);
351#else
352 submesh->setLargeCellSubset(* labels);
353#endif
354 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
355 field1));
356 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
357 field2));
358 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
359 (new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
360 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
361 (new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
362 }
363 else
364 {
365 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
366 (new GeometricField<T, fvPatchField, volMesh>(field1));
367 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh >>
368 (new GeometricField<T, fvPatchField, volMesh>(field2));
369 }
370
371 errField = autoPtr<GeometricField<T, fvPatchField, volMesh >>
372 (new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
373
374 if (L2Norm(field1) <= 1e-6)
375 {
376 err = 0;
377 }
378 else
379 {
380 err = L2Norm(errField()) / L2Norm(
381 field1_S());
382 }
383
384 return err;
385}
386
387template double errorL2Rel(GeometricField<scalar, fvPatchField, volMesh>&
388 field1,
389 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
390template double errorL2Rel(GeometricField<vector, fvPatchField, volMesh>&
391 field1,
392 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
393
394template<typename T>
395Eigen::MatrixXd errorL2Rel(PtrList<GeometricField<T, fvPatchField, volMesh >>&
396 fields1, PtrList<GeometricField<T, fvPatchField, volMesh >>& fields2,
397 List<label>* labels)
398{
399 Eigen::VectorXd err;
400
401 if (fields1.size() != fields2.size())
402 {
403 Info << "The two fields do not have the same size, code will abort" << endl;
404 exit(0);
405 }
406
407 err.resize(fields1.size(), 1);
408
409 for (label k = 0; k < fields1.size(); k++)
410 {
411 err(k, 0) = errorL2Rel(fields1[k], fields2[k], labels);
412 Info << " Error is " << err[k] << endl;
413 }
414
415 return err;
416}
417
418template Eigen::MatrixXd errorL2Rel(
419 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields1,
420 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& fields2,
421 List<label>* labels);
422template Eigen::MatrixXd errorL2Rel(
423 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields1,
424 PtrList<GeometricField<vector, fvPatchField, volMesh >>& fields2,
425 List<label>* labels);
426
427
428}
ITHACAerror.H
Header file of the ITHACAerror file.
ITHACAutilities
Namespace to implement some useful assign operation of OF fields.
Definition ITHACAassign.C:28
ITHACAutilities::errorLinfRel
double errorLinfRel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in Linf norm.
Definition ITHACAerror.C:106
ITHACAutilities::errorL2Rel
double errorL2Rel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in L2 norm.
Definition ITHACAerror.C:337
ITHACAutilities::frobNorm
double frobNorm(GeometricField< Type, PatchField, GeoMesh > &field)
Evaluate the Frobenius norm of a field.
Definition ITHACAnorm.C:355
ITHACAutilities::errorFrobRel
double errorFrobRel(GeometricField< Type, PatchField, GeoMesh > &field1, GeometricField< Type, PatchField, GeoMesh > &field2, List< label > *labels)
Computes the relative error between two Fields in the Frobenius norm.
Definition ITHACAerror.C:41
ITHACAutilities::L2Norm
double L2Norm(const T v)
Compute the L2 norm of v.
Definition ITHACAnorm.C:300
Generated by   doxygen 1.16.1