3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
resultevaluation3d.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
19
26#ifndef DUMUX_BENCHMARKRESULT_HH
27#define DUMUX_BENCHMARKRESULT_HH
28
29#include <dumux/porousmediumflow/sequential/onemodelproblem.hh>
30
31namespace Dumux
32{
33
37struct BenchmarkResult
38{
39private:
40 template<int dim>
41 struct FaceLayout
42 {
43 bool contains (Dune::GeometryType gt)
44 {
45 return gt.dim() == dim-1;
46 }
47 };
48
49public:
50 double relativeL2Error;
51 double ergrad;
52 double ervell2;
53 double uMinExact;
54 double uMaxExact;
55 double uMin;
56 double uMax;
57 double flux0;
58 double flux1;
59 double fluy0;
60 double fluy1;
61 double fluz0;
62 double fluz1;
63 double sumf;
64 double sumflux;
65 double exactflux0;
66 double exactflux1;
67 double exactfluy0;
68 double exactfluy1;
69 double exactfluz0;
70 double exactfluz1;
71 double errflx0;
72 double errflx1;
73 double errfly0;
74 double errfly1;
75 double erflm;
76 double ener1;
77 double ener2;
78 double eren;
79 double uMean;
80
81 template<class Grid, class ProblemType, class SolutionType>
82 void evaluate(const Grid& grid, ProblemType& problem,
83 SolutionType& solution, bool pureNeumann = false)
84 {
85 using Entity = typename Grid::Traits::template Codim<0>::Entity;
86 using Geometry = typename Entity::Geometry;
87 using GV = typename Grid::LevelGridView;
88 using IS = typename GV::IndexSet;
89 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
90 using ct = typename Grid::ctype;
91
92 enum{dim = Grid::dimension};
93
94 using JacobianInverseTransposed = typename Geometry::JacobianInverseTransposed;
95 using ReferenceElements = Dune::ReferenceElements<ct, dim>;
96 const GV& gridview(grid.levelGridView(grid.maxLevel()));
97 const IS& indexset(gridview.indexSet());
98 Mapper elementMapper(gridview, Dune::mcmgElementLayout());
99 Mapper faceMapper(gridview, Dune::mcmgLayout(Dune::Codim<1>()));
100 SolutionType& exactSol(grid.levelGridView(grid.maxLevel()));
101
102 uMean = 0;
103 double domainVolume = 0;
104 for (const auto& element : elements(gridview))
105 {
106 // get volume
107 double volume = element.geometry().volume();
108
109 // cell index
110 int indexi = elementMapper.index(element);
111
112 // get approximate solution value
113 uMean += volume*(problem.variables().cellData(indexi).globalPressure());
114
115 // add to domainVolume
116 domainVolume += volume;
117 }
118 uMean /= domainVolume;
119
120 if (pureNeumann) {
121 for (int i = 0; i < (int) problem.variables().pressure().size(); i++)
122 {
123 double press = problem.variables().cellData(i).globalPressure() - uMean;
124 problem.variables().cellData(i).setGlobalPressure(press);
125 }
126 }
127
128
129 uMinExact = 1e100;
130 uMaxExact = -1e100;
131 uMin = 1e100;
132 uMax = -1e100;
133 flux0 = 0;
134 flux1 = 0;
135 fluy0 = 0;
136 fluy1 = 0;
137 fluz0 = 0;
138 fluz1 = 0;
139 sumf = 0;
140 sumflux = 0;
141 exactflux0 = 0;
142 exactflux1 = 0;
143 exactfluy0 = 0;
144 exactfluy1 = 0;
145 exactfluz0 = 0;
146 exactfluz1 = 0;
147 erflm = 0;
148 ener1 = 0;
149 ener2 = 0;
150 double numerator = 0;
151 double denominator = 0;
152 double numeratorGrad = 0;
153 double denominatorGrad = 0;
154 double numeratorFlux = 0;
155 double denominatorFlux = 0;
156 bool exactOutput = true;
157 for (const auto& element : elements(gridview))
158 {
159 // element geometry
160 const Geometry& geometry = element.geometry();
161
162 // cell geometry type
163 Dune::GeometryType gt = geometry.type();
164
165 // cell center in reference element
166 const Dune::FieldVector<ct,dim>&
167 local = ReferenceElements::general(gt).position(0,0);
168
169 // get global coordinate of cell center
170 Dune::FieldVector<ct,dim> globalPos = geometry.global(local);
171
172 // get exact solution value
173 double exactValue = problem.exact(globalPos);
174
175 // cell index
176 int indexi = elementMapper.index(element);
177
178 // evaluate exact solution vector
179 exactSol[indexi] = exactValue;
180
181 // get approximate solution value
182 double approximateValue = problem.variables().cellData(indexi).globalPressure();
183
184 // update uMinExact and uMaxExact
185 using std::min;
186 using std::max;
187 uMinExact = min(uMinExact, exactValue);
188 uMaxExact = max(uMaxExact, exactValue);
189
190 // update uMin and uMax
191 uMin = min(uMin, approximateValue);
192 uMax = max(uMax, approximateValue);
193
194 // cell volume, assume linear map here
195 double volume = geometry.volume();
196
197 // update sumf
198 sumf += volume*(problem.source(globalPos, element, local)[0]);
199
200 // get the absolute permeability
201 Dune::FieldMatrix<double,dim,dim> K = problem.spatialParams().K(globalPos, element, local);
202
203 numerator += volume*(exactValue - approximateValue)*(exactValue - approximateValue);
204 denominator += volume*exactValue*exactValue;
205
206 int i = -1;
207 Dune::FieldVector<ct,2*dim> fluxVector;
208 Dune::FieldVector<ct,dim> exactGradient;
209 for (const auto& intersection : intersections(gridview, element))
210 {
211 // local number of facet
212 i = intersection.indexInInside();
213
214 // global number of face
215 int faceIndex = faceMapper.template map<1>(element, i);
216
217 Dune::FieldVector<double,dim> faceGlobal = intersection.geometry().center();
218 double faceVol = intersection.geometry().volume();
219
220 // get normal vector
221 Dune::FieldVector<double,dim> unitOuterNormal = intersection.centerUnitOuterNormal();
222
223 // get the approximate solution on the face
224 double approximateFace = (*solution.pressTrace)[faceIndex];
225
226 // get the exact gradient
227 exactGradient = problem.exactGrad(faceGlobal);
228
229 // get the negative exact velocity
230 Dune::FieldVector<double,dim> KGrad(0);
231 K.umv(exactGradient, KGrad);
232
233 // calculate the exact normal velocity
234 double exactFlux = KGrad*unitOuterNormal;
235
236 // get the approximate normalvelocity
237 double approximateFlux = solution.normalVelocity[indexi][i];
238
239 // calculate the difference in the normal velocity
240 double fluxDiff = exactFlux + approximateFlux;
241
242 // update mean value error
243 using std::abs;
244 erflm = max(erflm, abs(fluxDiff));
245
246 numeratorFlux += volume*fluxDiff*fluxDiff;
247 denominatorFlux += volume*exactFlux*exactFlux;
248
249 // calculate the fluxes through the element faces
250 exactFlux *= faceVol;
251 approximateFlux *= faceVol;
252 fluxVector[i] = approximateFlux;
253
254 //if (is.boundary()) {
255 if (!intersection.neighbor()) {
256 if (abs(faceGlobal[2]) < 1e-6) {
257 fluz0 += approximateFlux;
258 exactfluz0 += exactFlux;
259 ener2 += -approximateFlux*approximateFace;
260 }
261 else if (abs(faceGlobal[2] - 1.0) < 1e-6) {
262 fluz1 += approximateFlux;
263 exactfluz1 += exactFlux;
264 ener2 += -approximateFlux*approximateFace;
265 }
266 if (abs(faceGlobal[1]) < 1e-6) {
267 fluy0 += approximateFlux;
268 exactfluy0 += exactFlux;
269 ener2 += -approximateFlux*approximateFace;
270 }
271 else if (abs(faceGlobal[1] - 1.0) < 1e-6) {
272 fluy1 += approximateFlux;
273 exactfluy1 += exactFlux;
274 ener2 += -approximateFlux*approximateFace;
275 }
276 else if (faceGlobal[0] < 1e-6) {
277 flux0 += approximateFlux;
278 exactflux0 += exactFlux;
279 ener2 += -approximateFlux*approximateFace;
280 }
281 else if (abs(faceGlobal[0] - 1.0) < 1e-6) {
282 flux1 += approximateFlux;
283 exactflux1 += exactFlux;
284 ener2 += -approximateFlux*approximateFace;
285 }
286 }
287 }
288
289 // calculate velocity on reference element as the Raviart-Thomas-0
290 // interpolant of the fluxes
291 Dune::FieldVector<double, dim> refVelocity;
292 // simplices
293 if (geometry.type().isSimplex()) {
294 for (int dimIdx = 0; dimIdx < dim; dimIdx++)
295 {
296 refVelocity[dimIdx] = -fluxVector[dim - 1 - dimIdx];
297 for (int fIdx = 0; fIdx < dim + 1; fIdx++)
298 {
299 refVelocity[dimIdx] += fluxVector[fIdx]/(dim + 1);
300 }
301 }
302 }
303 // cubes
304 else if (geometry.type().isCube()){
305 for (int j = 0; j < dim; j++)
306 refVelocity[j] = 0.5 * (fluxVector[2*j + 1] - fluxVector[2*j]);
307 }
308 // 3D prism and pyramids
309 else {
310 DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
311 }
312
313 // get the transposed Jacobian of the element mapping
314 const JacobianInverseTransposed& jacobianInv = geometry.jacobianInverseTransposed(local);
315 JacobianInverseTransposed jacobianT(jacobianInv);
316 jacobianT.invert();
317
318 // calculate the element velocity by the Piola transformation
319 Dune::FieldVector<ct,dim> elementVelocity(0);
320 jacobianT.umtv(refVelocity, elementVelocity);
321 elementVelocity /= geometry.integrationElement(local);
322
323 // get the approximate gradient
324 Dune::FieldVector<ct,dim> approximateGradient;
325 K.solve(approximateGradient, elementVelocity);
326
327 // get the exact gradient
328 exactGradient = problem.exactGrad(globalPos);
329
330 // the difference between exact and approximate gradient
331 Dune::FieldVector<ct,dim> gradDiff(exactGradient);
332 gradDiff += approximateGradient;
333
334 // add to energy
335 ener1 += volume*(approximateGradient*elementVelocity);
336
337 numeratorGrad += volume*(gradDiff*gradDiff);
338 denominatorGrad += volume*(exactGradient*exactGradient);
339 }
340
341 using std::sqrt;
342 using std::abs;
343 using std::max;
344 relativeL2Error = sqrt(numerator/denominator);
345 ergrad = sqrt(numeratorGrad/denominatorGrad);
346 ervell2 = sqrt(numeratorFlux/denominatorFlux);
347 sumflux = flux0 + flux1 + fluy0 + fluy1 - sumf;
348 errflx0 = abs((flux0 + exactflux0)/exactflux0);
349 errflx1 = abs((flux1 + exactflux1)/exactflux1);
350 errfly0 = abs((fluy0 + exactfluy0)/exactfluy0);
351 errfly1 = abs((fluy1 + exactfluy1)/exactfluy1);
352 eren = abs(ener1 - ener2)/max(ener1, ener2);
353
354 // generate a VTK file of exact solution on the maximal grid level
355 if (exactOutput)
356 {
357 Dune::VTKWriter<GV> vtkwriter(gridview);
358 char fname[128];
359 sprintf(fname,"%d","exactSol-numRefine", grid.maxLevel());
360 vtkwriter.addCellData(exactSol,"exact pressure solution~");
361 vtkwriter.write(fname,Dune::VTK::ascii);
362 }
363
364 return;
365 }
366};
367
371struct ResultEvaluation
372{
373private:
374 template<int dim>
375 struct FaceLayout
376 {
377 bool contains (Dune::GeometryType gt)
378 {
379 return gt.dim() == dim-1;
380 }
381 };
382
383public:
384 double relativeL2Error;
385 double relativeL2ErrorIn;
386 double relativeL2ErrorOut;
387 double ergrad;
388 double ervell2;
389 double ervell2In;
390 double uMinExact;
391 double uMaxExact;
392 double uMin;
393 double uMax;
394 double flux0;
395 double flux1;
396 double fluy0;
397 double fluy1;
398 double fluz0;
399 double fluz1;
400 double sumf;
401 double sumflux;
402 double exactflux0;
403 double exactflux1;
404 double exactfluy0;
405 double exactfluy1;
406 double exactfluz0;
407 double exactfluz1;
408 double errflx0;
409 double errflx1;
410 double errfly0;
411 double errfly1;
412 double erflm;
413 double ener1;
414 double hMax;
415
416 template<class GridView, class Problem>
417 void evaluate(const GridView& gridView,
418 Problem& problem, bool consecutiveNumbering = false)
419 {
420 using Grid = typename GridView::Grid;
421 using Scalar = typename Grid::ctype;
422 enum {dim=Grid::dimension};
423 using Element = typename Grid::template Codim<0>::Entity;
424 using Geometry = typename Element::Geometry;
425 using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
426 using SolVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
427 using JacobianInverseTransposed = typename Geometry::JacobianInverseTransposed;
428 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
429
430 ElementMapper elementMapper(gridView, Dune::mcmgElementLayout());
431 SolVector exactSol(gridView.size(0));
432
433
434 uMinExact = 1e100;
435 uMaxExact = -1e100;
436 uMin = 1e100;
437 uMax = -1e100;
438 flux0 = 0;
439 flux1 = 0;
440 fluy0 = 0;
441 fluy1 = 0;
442 fluz0 = 0;
443 fluz1 = 0;
444 sumf = 0;
445 sumflux = 0;
446 exactflux0 = 0;
447 exactflux1 = 0;
448 exactfluy0 = 0;
449 exactfluy1 = 0;
450 exactfluz0 = 0;
451 exactfluz1 = 0;
452 erflm = 0;
453 ener1 = 0;
454 hMax = 0;
455 Scalar numerator = 0;
456 Scalar denominator = 0;
457 Scalar numeratorIn = 0;
458 Scalar denominatorIn = 0;
459 Scalar numeratorOut = 0;
460 Scalar denominatorOut = 0;
461 double numeratorGrad = 0;
462 double denominatorGrad = 0;
463 double numeratorFlux = 0;
464 double denominatorFlux = 0;
465 double numeratorFluxIn = 0;
466 double denominatorFluxIn = 0;
467 bool exactOutput = true;
468 for (const auto& element : elements(gridView))
469 {
470 // element geometry
471 const Geometry& geometry = element.geometry();
472
473 Dune::GeometryType geomType = geometry.type();
474
475 const Dune::FieldVector<Scalar,dim>& local = ReferenceElements::general(geomType).position(0, 0);
476 Dune::FieldVector<Scalar,dim> globalPos = geometry.global(local);
477
478 Scalar volume = geometry.volume();
479
480 int eIdx = problem.variables().index(element);
481
482 Scalar approxPressure = problem.variables().cellData(eIdx).globalPressure();
483 Scalar exactPressure = problem.exact(globalPos);
484
485 // evaluate exact solution vector
486 exactSol[eIdx] = exactPressure;
487
488 // output local relative error for each cell
489 //std::cout << "local relative error for cell "<< eIdx << " is: "
490 // << (approxPressure - exactPressure)/exactPressure << std::endl;
491
492 numerator += volume*(approxPressure - exactPressure)*(approxPressure - exactPressure);
493 denominator += volume*exactPressure*exactPressure;
494
495 // calculate the relative error for inner cells
496 bool bndCell = false;
497 for (const auto& intersection : intersections(gridView, element))
498 {
499 if (intersection.boundary())
500 {
501 bndCell = true;
502 break;
503 }
504 }
505
506 if (bndCell)
507 {
508 numeratorOut += volume*(approxPressure - exactPressure)*(approxPressure - exactPressure);
509 denominatorOut += volume*exactPressure*exactPressure;
510 }
511 else
512 {
513 numeratorIn += volume*(approxPressure - exactPressure)*(approxPressure - exactPressure);
514 denominatorIn += volume*exactPressure*exactPressure;
515 }
516
517 // update uMinExact and uMaxExact
518 using std::min;
519 using std::max;
520 uMinExact = min(uMinExact, exactPressure);
521 uMaxExact = max(uMaxExact, exactPressure);
522
523 // update uMin and uMax
524 uMin = min(uMin, approxPressure);
525 uMax = max(uMax, approxPressure);
526
527 typename Problem::PrimaryVariables sourceVec;
528 problem.source(sourceVec, element);
529 sumf += volume*sourceVec[0];
530
531 // get the absolute permeability
532 Dune::FieldMatrix<double,dim,dim> K = problem.spatialParams().intrinsicPermeability(element);
533
534 int isIdx = -1;
535 Dune::FieldVector<Scalar,2*dim> fluxVector;
536 Dune::FieldVector<Scalar,dim> exactGradient;
537 for (const auto& intersection : intersections(gridView, element))
538 {
539 // local number of facet
540 int fIdx = intersection.indexInInside();
541
542 if (consecutiveNumbering)
543 isIdx++;
544 else
545 isIdx = fIdx;
546
547 Dune::FieldVector<double,dim> faceGlobal = intersection.geometry().center();
548 double faceVol = intersection.geometry().volume();
549
550 // get normal vector
551 Dune::FieldVector<double,dim> unitOuterNormal = intersection.centerUnitOuterNormal();
552
553 // get the exact gradient
554 exactGradient = problem.exactGrad(faceGlobal);
555
556 // get the negative exact velocity
557 Dune::FieldVector<double,dim> KGrad(0);
558 K.umv(exactGradient, KGrad);
559
560 // calculate the exact normal velocity
561 double exactFlux = KGrad*unitOuterNormal;
562
563 // get the approximate normalvelocity
564 double approximateFlux = problem.variables().cellData(eIdx).fluxData().velocityTotal(fIdx)*unitOuterNormal;
565
566 // calculate the difference in the normal velocity
567 double fluxDiff = exactFlux + approximateFlux;
568
569 // update mean value error
570 using std::abs;
571 using std::max;
572 erflm = max(erflm, abs(fluxDiff));
573
574 numeratorFlux += volume*fluxDiff*fluxDiff;
575 denominatorFlux += volume*exactFlux*exactFlux;
576
577 if (!bndCell)
578 {
579 numeratorFluxIn += volume*fluxDiff*fluxDiff;
580 denominatorFluxIn += volume*exactFlux*exactFlux;
581 }
582
583 // calculate the fluxes through the element faces
584 exactFlux *= faceVol;
585 approximateFlux *= faceVol;
586 fluxVector[fIdx] = approximateFlux;
587
588 if (!intersection.neighbor()) {
589 if (abs(faceGlobal[2]) < 1e-6) {
590 fluz0 += approximateFlux;
591 exactfluz0 += exactFlux;
592 }
593 else if (abs(faceGlobal[2] - 1.0) < 1e-6) {
594 fluz1 += approximateFlux;
595 exactfluz1 += exactFlux;
596 }
597 if (abs(faceGlobal[1]) < 1e-6) {
598 fluy0 += approximateFlux;
599 exactfluy0 += exactFlux;
600 }
601 else if (abs(faceGlobal[1] - 1.0) < 1e-6) {
602 fluy1 += approximateFlux;
603 exactfluy1 += exactFlux;
604 }
605 else if (faceGlobal[0] < 1e-6) {
606 flux0 += approximateFlux;
607 exactflux0 += exactFlux;
608 }
609 else if (abs(faceGlobal[0] - 1.0) < 1e-6) {
610 flux1 += approximateFlux;
611 exactflux1 += exactFlux;
612 }
613 }
614 }
615
616 // calculate velocity on reference element
617 Dune::FieldVector<Scalar,dim> refVelocity;
618 if (geometry.corners() == 8) {
619 refVelocity[0] = 0.5*(fluxVector[1] - fluxVector[0]);
620 refVelocity[1] = 0.5*(fluxVector[3] - fluxVector[2]);
621 refVelocity[2] = 0.5*(fluxVector[5] - fluxVector[4]);
622 }
623 else {
624
625 }
626
627 // get the transposed Jacobian of the element mapping
628 const JacobianInverseTransposed& jacobianInv = geometry.jacobianInverseTransposed(local);
629 JacobianInverseTransposed jacobianT(jacobianInv);
630 jacobianT.invert();
631
632 // calculate the element velocity by the Piola transformation
633 Dune::FieldVector<Scalar,dim> elementVelocity(0);
634 jacobianT.umtv(refVelocity, elementVelocity);
635 elementVelocity /= geometry.integrationElement(local);
636
637 // get the approximate gradient
638 Dune::FieldVector<Scalar,dim> approximateGradient;
639 K.solve(approximateGradient, elementVelocity);
640
641 // get the exact gradient
642 exactGradient = problem.exactGrad(globalPos);
643
644 // the difference between exact and approximate gradient
645 Dune::FieldVector<Scalar,dim> gradDiff(exactGradient);
646 gradDiff += approximateGradient;
647
648 // add to energy
649 ener1 += volume*(approximateGradient*elementVelocity);
650
651 numeratorGrad += volume*(gradDiff*gradDiff);
652 denominatorGrad += volume*(exactGradient*exactGradient);
653
654 // calculate the maximum of the diagonal length of all elements on leaf grid
655 for (int i = 0; i < element.geometry().corners(); ++i)
656 {
657 Dune::FieldVector<Scalar,dim> corner1 = element.geometry().corner(i);
658
659 for (int j = 0; j < element.geometry().corners(); ++j)
660 {
661 // get all corners of current element and compare the distances between them
662 Dune::FieldVector<Scalar,dim> corner2 = element.geometry().corner(j);
663
664 // distance vector between corners
665 Dune::FieldVector<Scalar,dim> distVec = corner1 - corner2;
666 Scalar dist = distVec.two_norm();
667
668 if (hMax < dist)
669 hMax = dist;
670 }
671 }
672 }
673
674 using std::sqrt;
675 using std::abs;
676 relativeL2Error = sqrt(numerator/denominator);
677 relativeL2ErrorIn = sqrt(numeratorIn/denominatorIn);
678 relativeL2ErrorOut = sqrt(numeratorOut/denominatorOut);
679 ergrad = sqrt(numeratorGrad/denominatorGrad);
680 ervell2 = sqrt(numeratorFlux/denominatorFlux);
681 ervell2In = sqrt(numeratorFluxIn/denominatorFluxIn);
682 sumflux = flux0 + flux1 + fluy0 + fluy1 - sumf;
683 errflx0 = abs((flux0 + exactflux0)/exactflux0);
684 errflx1 = abs((flux1 + exactflux1)/exactflux1);
685 errfly0 = abs((fluy0 + exactfluy0)/exactfluy0);
686 errfly1 = abs((fluy1 + exactfluy1)/exactfluy1);
687
688 // generate a VTK file of exact solution
689 if (exactOutput)
690 {
691 Dune::VTKWriter<GridView> vtkwriter(gridView);
692 char fname[128];
693 sprintf(fname, "exactSol-numRefine%d", gridView.grid().maxLevel());
694 vtkwriter.addCellData(exactSol,"exact pressure solution~");
695 vtkwriter.write(fname,Dune::VTK::ascii);
696 }
697
698 return;
699 }
700
701
702 template<class GridView, class Problem, class SolVector, class VelVector>
703 void evaluateCP(const GridView& gridView, Problem& problem,
704 const SolVector& solution, const VelVector& velocity, bool switchNormals = false)
705 {
706 using Grid = typename GridView::Grid;
707 using Scalar = typename Grid::ctype;
708 enum {dim=Grid::dimension, maxIntersections = 12};
709 using Element = typename Grid::template Codim<0>::Entity;
710 using Geometry = typename Element::Geometry;
711 using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
712 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
713 using ReferenceFaces = Dune::ReferenceElements<Scalar, dim-1>;
714
715 ElementMapper elementMapper(gridView, Dune::mcmgElementLayout());
716
717 uMinExact = 1e100;
718 uMaxExact = -1e100;
719 uMin = 1e100;
720 uMax = -1e100;
721 sumf = 0;
722 sumflux = 0;
723 erflm = 0;
724 Scalar maxFlux = -1;
725 Scalar numerator = 0;
726 Scalar denominator = 0;
727 for (const auto& element : elements(gridView))
728 {
729 // element geometry
730 const Geometry& geometry = element.geometry();
731
732 Dune::GeometryType geomType = geometry.type();
733
734 const Dune::FieldVector<Scalar,dim>& local = ReferenceElements::general(geomType).position(0, 0);
735 Dune::FieldVector<Scalar,dim> globalPos = geometry.global(local);
736
737 Scalar volume = geometry.integrationElement(local)
738 *ReferenceElements::general(geomType).volume();
739
740 int eIdx = elementMapper.index(element);
741
742 Scalar approxPressure = solution[eIdx];
743 Scalar exactPressure = problem.exact(globalPos);
744 numerator += volume*(approxPressure - exactPressure)*(approxPressure - exactPressure);
745 denominator += volume*exactPressure*exactPressure;
746
747 // update uMinExact and uMaxExact
748 using std::min;
749 using std::max;
750 uMinExact = min(uMinExact, exactPressure);
751 uMaxExact = max(uMaxExact, exactPressure);
752
753 // update uMin and uMax
754 uMin = min(uMin, approxPressure);
755 uMax = max(uMax, approxPressure);
756
757 sumf += volume*problem.source(globalPos, element, local)[0];
758
759 // get the absolute permeability
760 Dune::FieldMatrix<Scalar,dim,dim> K = problem.spatialParams().K(globalPos, element, local);
761
762 int i = -1;
763 Dune::FieldVector<Scalar,dim> exactGradient;
764 for (const auto& intersection : intersections(gridView, element))
765 {
766 // get geometry type of face
767 Dune::GeometryType gtf = intersection.geometryInInside().type();
768
769 // local number of facet
770 i++;
771
772 // center in face's reference element
773 const Dune::FieldVector<Scalar,dim-1>& faceLocalNm1 = ReferenceFaces::general(gtf).position(0,0);
774
775 // center of face in global coordinates
776 Dune::FieldVector<Scalar,dim> faceGlobal = intersection.geometry().global(faceLocalNm1);
777
778 // get normal vector
779 Dune::FieldVector<Scalar,dim> unitOuterNormal = intersection.unitOuterNormal(faceLocalNm1);
780
781 if (switchNormals)
782 unitOuterNormal *= -1.0;
783
784 // get the exact gradient
785 exactGradient = problem.exactGrad(faceGlobal);
786 exactGradient.axpy(-problem.wettingPhase().density(), problem.gravity());
787
788 // get the negative exact velocity
789 Dune::FieldVector<Scalar,dim> KGrad(0);
790 K.umv(exactGradient, KGrad);
791
792 // calculate the exact normal velocity
793 Scalar exactFlux = KGrad*unitOuterNormal;
794
795 // get the approximate normalvelocity
796 Scalar approximateFlux = problem.cellData(eIdx).fluxData().velocityTotal(i)*unitOuterNormal;
797
798 // calculate the difference in the normal velocity
799 Scalar fluxDiff = exactFlux + approximateFlux;
800
801 using std::abs;
802 // if (abs(fluxDiff) > 1e-16)
803 // {
804 // std::cout << "faceGlobal " << faceGlobal << ": exact "
805 // << exactFlux << ", approximate " << approximateFlux << std::endl;
806 // }
807
808 // update mean value error
809 erflm = max(erflm, abs(fluxDiff));
810
811 maxFlux = max(maxFlux, abs(exactFlux));
812
813 sumflux += approximateFlux;
814 }
815 }
816
817 using std::sqrt;
818 using std::max;
819 relativeL2Error = sqrt(numerator/denominator);
820 sumflux -= sumf;
821 erflm /= max(1e-6, maxFlux);
822
823 return;
824 }
825};
826
827}
828
829#endif
onemodelproblem.hh
Base class for definition of an sequential diffusion (pressure) or transport problem.
Dumux::intersections
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
Dumux
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Dumux::ResultEvaluation::ener1
double ener1
Definition: resultevaluation.hh:60
Dumux::ResultEvaluation::relativeL2Error
double relativeL2Error
Definition: resultevaluation.hh:40
Dumux::ResultEvaluation::exactfluz0
double exactfluz0
Definition: resultevaluation3d.hh:406
Dumux::ResultEvaluation::relativeL2ErrorOut
double relativeL2ErrorOut
Definition: resultevaluation3d.hh:386
Dumux::ResultEvaluation::fluy0
double fluy0
Definition: resultevaluation.hh:47
Dumux::ResultEvaluation::exactfluz1
double exactfluz1
Definition: resultevaluation3d.hh:407
Dumux::ResultEvaluation::uMinExact
double uMinExact
Definition: resultevaluation3d.hh:390
Dumux::ResultEvaluation::uMin
double uMin
Definition: resultevaluation.hh:43
Dumux::ResultEvaluation::exactflux1
double exactflux1
Definition: resultevaluation.hh:52
Dumux::ResultEvaluation::errfly0
double errfly0
Definition: resultevaluation.hh:57
Dumux::ResultEvaluation::ervell2
double ervell2
Definition: resultevaluation.hh:42
Dumux::ResultEvaluation::uMaxExact
double uMaxExact
Definition: resultevaluation3d.hh:391
Dumux::ResultEvaluation::sumf
double sumf
Definition: resultevaluation.hh:49
Dumux::ResultEvaluation::ervell2In
double ervell2In
Definition: resultevaluation3d.hh:389
Dumux::ResultEvaluation::evaluate
void evaluate(const GridView &gridView, Problem &problem, bool consecutiveNumbering=false)
Definition: resultevaluation3d.hh:417
Dumux::ResultEvaluation::flux0
double flux0
Definition: resultevaluation.hh:45
Dumux::ResultEvaluation::erflm
double erflm
Definition: resultevaluation.hh:59
Dumux::ResultEvaluation::fluy1
double fluy1
Definition: resultevaluation.hh:48
Dumux::ResultEvaluation::ergrad
double ergrad
Definition: resultevaluation.hh:41
Dumux::ResultEvaluation::relativeL2ErrorIn
double relativeL2ErrorIn
Definition: resultevaluation3d.hh:385
Dumux::ResultEvaluation::errfly1
double errfly1
Definition: resultevaluation.hh:58
Dumux::ResultEvaluation::errflx0
double errflx0
Definition: resultevaluation.hh:55
Dumux::ResultEvaluation::exactflux0
double exactflux0
Definition: resultevaluation.hh:51
Dumux::ResultEvaluation::sumflux
double sumflux
Definition: resultevaluation.hh:50
Dumux::ResultEvaluation::fluz1
double fluz1
Definition: resultevaluation3d.hh:399
Dumux::ResultEvaluation::errflx1
double errflx1
Definition: resultevaluation.hh:56
Dumux::ResultEvaluation::exactfluy1
double exactfluy1
Definition: resultevaluation.hh:54
Dumux::ResultEvaluation::uMax
double uMax
Definition: resultevaluation.hh:44
Dumux::ResultEvaluation::evaluateCP
void evaluateCP(const GridView &gridView, Problem &problem, const SolVector &solution, const VelVector &velocity, bool switchNormals=false)
Definition: resultevaluation3d.hh:703
Dumux::ResultEvaluation::hMax
double hMax
Definition: resultevaluation3d.hh:414
Dumux::ResultEvaluation::fluz0
double fluz0
Definition: resultevaluation3d.hh:398
Dumux::ResultEvaluation::flux1
double flux1
Definition: resultevaluation.hh:46
Dumux::ResultEvaluation::exactfluy0
double exactfluy0
Definition: resultevaluation.hh:53
Dumux::BenchmarkResult
calculate errors for a FVCA5 benchmark problem
Definition: resultevaluation3d.hh:38
Dumux::BenchmarkResult::ener1
double ener1
Definition: resultevaluation3d.hh:76
Dumux::BenchmarkResult::ergrad
double ergrad
Definition: resultevaluation3d.hh:51
Dumux::BenchmarkResult::errflx1
double errflx1
Definition: resultevaluation3d.hh:72
Dumux::BenchmarkResult::fluy0
double fluy0
Definition: resultevaluation3d.hh:59
Dumux::BenchmarkResult::fluy1
double fluy1
Definition: resultevaluation3d.hh:60
Dumux::BenchmarkResult::uMinExact
double uMinExact
Definition: resultevaluation3d.hh:53
Dumux::BenchmarkResult::exactfluz0
double exactfluz0
Definition: resultevaluation3d.hh:69
Dumux::BenchmarkResult::flux0
double flux0
Definition: resultevaluation3d.hh:57
Dumux::BenchmarkResult::uMaxExact
double uMaxExact
Definition: resultevaluation3d.hh:54
Dumux::BenchmarkResult::sumflux
double sumflux
Definition: resultevaluation3d.hh:64
Dumux::BenchmarkResult::exactflux1
double exactflux1
Definition: resultevaluation3d.hh:66
Dumux::BenchmarkResult::flux1
double flux1
Definition: resultevaluation3d.hh:58
Dumux::BenchmarkResult::eren
double eren
Definition: resultevaluation3d.hh:78
Dumux::BenchmarkResult::evaluate
void evaluate(const Grid &grid, ProblemType &problem, SolutionType &solution, bool pureNeumann=false)
Definition: resultevaluation3d.hh:82
Dumux::BenchmarkResult::errflx0
double errflx0
Definition: resultevaluation3d.hh:71
Dumux::BenchmarkResult::errfly0
double errfly0
Definition: resultevaluation3d.hh:73
Dumux::BenchmarkResult::uMax
double uMax
Definition: resultevaluation3d.hh:56
Dumux::BenchmarkResult::fluz0
double fluz0
Definition: resultevaluation3d.hh:61
Dumux::BenchmarkResult::ervell2
double ervell2
Definition: resultevaluation3d.hh:52
Dumux::BenchmarkResult::sumf
double sumf
Definition: resultevaluation3d.hh:63
Dumux::BenchmarkResult::exactfluz1
double exactfluz1
Definition: resultevaluation3d.hh:70
Dumux::BenchmarkResult::exactfluy1
double exactfluy1
Definition: resultevaluation3d.hh:68
Dumux::BenchmarkResult::relativeL2Error
double relativeL2Error
Definition: resultevaluation3d.hh:50
Dumux::BenchmarkResult::uMin
double uMin
Definition: resultevaluation3d.hh:55
Dumux::BenchmarkResult::fluz1
double fluz1
Definition: resultevaluation3d.hh:62
Dumux::BenchmarkResult::erflm
double erflm
Definition: resultevaluation3d.hh:75
Dumux::BenchmarkResult::exactfluy0
double exactfluy0
Definition: resultevaluation3d.hh:67
Dumux::BenchmarkResult::errfly1
double errfly1
Definition: resultevaluation3d.hh:74
Dumux::BenchmarkResult::exactflux0
double exactflux0
Definition: resultevaluation3d.hh:65
Dumux::BenchmarkResult::uMean
double uMean
Definition: resultevaluation3d.hh:79
Dumux::BenchmarkResult::ener2
double ener2
Definition: resultevaluation3d.hh:77
Generated by doxygen 1.9.3