3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
onemodelproblem.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
20#ifndef DUMUX_ONE_MODEL_PROBLEM_HH
21#define DUMUX_ONE_MODEL_PROBLEM_HH
22
23#include <dune/common/shared_ptr.hh>
27#include <dumux/io/restart.hh>
28
29
35namespace Dumux
36{
37
44template<class TypeTag>
46{
47private:
48 using Implementation = GetPropType<TypeTag, Properties::Problem>;
50 using Grid = typename GridView::Grid;
51
53
55
57
59
61
63 using VertexMapper = typename SolutionTypes::VertexMapper;
64 using ElementMapper = typename SolutionTypes::ElementMapper;
65
66 enum
67 {
68 dim = GridView::dimension,
69 dimWorld = GridView::dimensionworld
70 };
71 enum
72 {
73 wetting = 0, nonwetting = 1
74 };
75
76 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
77 using Element = typename GridView::Traits::template Codim<0>::Entity;
78 using Intersection = typename GridView::Intersection;
79
80 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
82
83 // private!! copy constructor
85 {}
86
87public:
88
90
94 OneModelProblem(Grid& grid, bool verbose = true)
95 : gridView_(grid.leafGridView()),
96 bBoxMin_(std::numeric_limits<double>::max()),
97 bBoxMax_(-std::numeric_limits<double>::max()),
98 variables_(grid.leafGridView()),
99 outputInterval_(1),
100 outputTimeInterval_(0)
101 {
102 // calculate the bounding box of the grid view
103 using std::max;
104 using std::min;
105 for (const auto& vertex : vertices(grid.leafGridView())) {
106 for (int i=0; i<dim; i++) {
107 bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
108 bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
109 }
110 }
111
112 timeManager_ = std::make_shared<TimeManager>(verbose);
113
114 model_ = std::make_shared<Model>(asImp_()) ;
115 maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
116 }
117
119
123 OneModelProblem(TimeManager& timeManager, Grid& grid)
124 : gridView_(grid.leafGridView()),
125 bBoxMin_(std::numeric_limits<double>::max()),
126 bBoxMax_(-std::numeric_limits<double>::max()),
127 variables_(grid.leafGridView()),
128 outputInterval_(1),
129 outputTimeInterval_(0)
130 {
131 // calculate the bounding box of the grid view
132 using std::max;
133 using std::min;
134 for (const auto& vertex : vertices(grid.leafGridView())) {
135 for (int i=0; i<dim; i++) {
136 bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
137 bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
138 }
139 }
140
141 timeManager_ = Dune::stackobject_to_shared_ptr<TimeManager>(timeManager);
142
143 model_ = std::make_shared<Model>(asImp_()) ;
144 maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
145 }
146
154 void boundaryTypes(BoundaryTypes &bcTypes,
155 const Intersection &intersection) const
156 {
157 // forward it to the method which only takes the global coordinate
158 asImp_().boundaryTypesAtPos(bcTypes, intersection.geometry().center());
159 }
160
168 void boundaryTypesAtPos(BoundaryTypes &bcTypes,
169 const GlobalPosition &globalPos) const
170 {
171 // Throw an exception (there is no reasonable default value
172 // for Dirichlet conditions)
173 DUNE_THROW(Dune::InvalidStateException,
174 "The problem does not provide "
175 "a boundaryTypesAtPos() method.");
176 }
177
187 void dirichlet(PrimaryVariables &values,
188 const Intersection &intersection) const
189 {
190 // forward it to the method which only takes the global coordinate
191 asImp_().dirichletAtPos(values, intersection.geometry().center());
192 }
193
203 void dirichletAtPos(PrimaryVariables &values,
204 const GlobalPosition &globalPos) const
205 {
206 // Throw an exception (there is no reasonable default value
207 // for Dirichlet conditions)
208 DUNE_THROW(Dune::InvalidStateException,
209 "The problem specifies that some boundary "
210 "segments are dirichlet, but does not provide "
211 "a dirichletAtPos() method.");
212 }
213
224 void neumann(PrimaryVariables &values,
225 const Intersection &intersection) const
226 {
227 // forward it to the interface with only the global position
228 asImp_().neumannAtPos(values, intersection.geometry().center());
229 }
230
241 void neumannAtPos(PrimaryVariables &values,
242 const GlobalPosition &globalPos) const
243 {
244 // Throw an exception (there is no reasonable default value
245 // for Neumann conditions)
246 DUNE_THROW(Dune::InvalidStateException,
247 "The problem specifies that some boundary "
248 "segments are neumann, but does not provide "
249 "a neumannAtPos() method.");
250 }
251
262 void source(PrimaryVariables &values,
263 const Element &element) const
264 {
265 // forward to generic interface
266 asImp_().sourceAtPos(values, element.geometry().center());
267 }
268
282 void sourceAtPos(PrimaryVariables &values,
283 const GlobalPosition &globalPos) const
284 { // Throw an exception (there is no initial condition)
285 DUNE_THROW(Dune::InvalidStateException,
286 "The problem does not provide "
287 "a sourceAtPos() method.");
288 }
289
299 void initial(PrimaryVariables &values,
300 const Element &element) const
301 {
302 // forward to generic interface
303 asImp_().initialAtPos(values, element.geometry().center());
304 }
305
316 void initialAtPos(PrimaryVariables &values,
317 const GlobalPosition &globalPos) const
318 {
319 // initialize with 0 by default
320 values = 0;
321 }
322
327 void init()
328 {
329 // set the initial condition of the model
330 variables_.initialize();
331 model().initialize();
332 }
333
339 {}
340
346 {}
347
357 {}
358
365 {}
366
372 Scalar maxTimeStepSize() const
373 { return maxTimeStepSize_; }
374
378 Scalar timeStepSize() const
379 { return timeManager().timeStepSize(); }
380
384 void setTimeStepSize(Scalar dt)
385 { timeManager().setTimeStepSize(dt); }
386
392 Scalar nextTimeStepSize(Scalar dt)
393 { return timeManager().timeStepSize();}
394
404 {
405 return
406 timeManager().timeStepIndex() > 0 &&
407 (timeManager().timeStepIndex() % 5 == 0);
408 }
409
415 void setOutputTimeInterval(const Scalar timeInterval)
416 {
417 outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
418 timeManager().startNextEpisode(outputTimeInterval_);
419 }
420
426 void setOutputInterval(int interval)
427 { outputInterval_ = interval; }
428
438 bool shouldWriteOutput() const
439 {
440 if (outputInterval_ > 0)
441 {
442 if (timeManager().timeStepIndex() % outputInterval_ == 0
443 || timeManager().willBeFinished()
444 || timeManager().episodeWillBeFinished())
445 {
446 return true;
447 }
448 }
449 else if (timeManager().willBeFinished()
450 || timeManager().episodeWillBeFinished() || timeManager().timeStepIndex() == 0)
451 {
452 return true;
453 }
454 return false;
455 }
456
458 {}
459
461 void writeOutput(bool verbose = true)
462 {
463 if (verbose && gridView().comm().rank() == 0)
464 std::cout << "Writing result file for current time step\n";
465 if (!resultWriter_)
466 resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
467 resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
468 model().addOutputVtkFields(*resultWriter_);
469 asImp_().addOutputVtkFields();
470 resultWriter_->endWrite();
471 }
472
477 {
478 if (outputTimeInterval_ > 0.0 && !timeManager().finished())
479 {
480 timeManager().startNextEpisode(outputTimeInterval_);
481 }
482 else if (!timeManager().finished())
483 {
484 std::cerr << "The end of an episode is reached, but the problem "
485 << "does not override the episodeEnd() method. "
486 << "Doing nothing!\n";
487 }
488 }
489
490 // \}
491
499 const std::string& name() const
500 {
501 return simname_;
502 }
503
511 void setName(const std::string& newName)
512 {
513 simname_ = newName;
514 }
515
519 const GridView &gridView() const
520 { return gridView_; }
521
525 const VertexMapper &vertexMapper() const
526 { return variables_.vertexMapper(); }
527
531 const ElementMapper &elementMapper() const
532 { return variables_.elementMapper(); }
533
538 const GlobalPosition &bBoxMin() const
539 { return bBoxMin_; }
540
545 const GlobalPosition &bBoxMax() const
546 { return bBoxMax_; }
547
551 TimeManager &timeManager()
552 { return *timeManager_; }
553
557 const TimeManager &timeManager() const
558 { return *timeManager_; }
559
563 Variables& variables ()
564 { return variables_; }
565
569 const Variables& variables () const
570 { return variables_; }
571
575 Model &model()
576 { return *model_; }
577
581 const Model &model() const
582 { return *model_; }
583 // \}
584
585
589 // \{
590
601 {
602 using Restarter = Restart;
603
604 Restarter res;
605 res.serializeBegin(asImp_());
606 std::cout << "Serialize to file " << res.fileName() << "\n";
607
608 timeManager().serialize(res);
609 resultWriter().serialize(res);
610 res.template deserializeEntities<0> (model(), gridView_);
611
612 res.serializeEnd();
613 }
614
621 void restart(double tRestart)
622 {
623 using Restarter = Restart;
624
625 Restarter res;
626 res.deserializeBegin(asImp_(), tRestart);
627 std::cout << "Deserialize from file " << res.fileName() << "\n";
628
629 timeManager().deserialize(res);
631 res.template deserializeEntities<0> (model(), gridView_);
632
633 res.deserializeEnd();
634 }
635
636 // \}
637
638protected:
640 {
641 if (!resultWriter_)
642 resultWriter_ = std::make_shared<VtkMultiWriter>(gridView_, asImp_().name());
643 return *resultWriter_;
644 }
645
647 {
648 if (!resultWriter_)
649 resultWriter_ = std::make_shared<VtkMultiWriter>(gridView_, asImp_().name());
650 return *resultWriter_;
651 }
652
653private:
655 Implementation &asImp_()
656 { return *static_cast<Implementation *>(this); }
657
659 const Implementation &asImp_() const
660 { return *static_cast<const Implementation *>(this); }
661
662 std::string simname_; // a string for the name of the current simulation,
663 // which could be set by means of an program argument,
664 // for example.
665 const GridView gridView_;
666
667 GlobalPosition bBoxMin_;
668 GlobalPosition bBoxMax_;
669
670 std::shared_ptr<TimeManager> timeManager_;
671 Scalar maxTimeStepSize_;
672
673 Variables variables_;
674
675 std::shared_ptr<Model> model_;
676
677 std::shared_ptr<VtkMultiWriter> resultWriter_;
678 int outputInterval_;
679 Scalar outputTimeInterval_;
680};
681
682}
683#endif
Generic class to read/write restart files.
Simplifies writing multi-file VTK datasets.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Load or save a state of a model to/from the harddisk.
Definition: restart.hh:44
void deserializeBegin(Problem &problem, double t)
Start reading a restart file at a certain simulated time.
Definition: restart.hh:168
void serializeBegin(Problem &problem)
Write the current state of the model to disk.
Definition: restart.hh:94
Simplifies writing multi-file VTK datasets.
Definition: vtkmultiwriter.hh:146
void serialize(Restarter &res)
Write the multi-writer's state to a restart file.
Definition: vtkmultiwriter.hh:374
void deserialize(Restarter &res)
Read the multi-writer's state from a restart file.
Definition: vtkmultiwriter.hh:408
Base class for definition of an sequential diffusion (pressure) or transport problem.
Definition: onemodelproblem.hh:46
Scalar timeStepSize() const
Returns the current time step size [seconds].
Definition: onemodelproblem.hh:378
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition: onemodelproblem.hh:525
void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
Evaluate the initial value for a control volume.
Definition: onemodelproblem.hh:316
void timeIntegration()
Called by TimeManager in order to do a time integration on the model.
Definition: onemodelproblem.hh:345
void setName(const std::string &newName)
Set the problem name.
Definition: onemodelproblem.hh:511
void setTimeStepSize(Scalar dt)
Sets the current time step size [seconds].
Definition: onemodelproblem.hh:384
void dirichlet(PrimaryVariables &values, const Intersection &intersection) const
Evaluate the boundary conditions for a dirichlet control volume.
Definition: onemodelproblem.hh:187
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition: onemodelproblem.hh:531
TimeManager & timeManager()
Returns TimeManager object used by the simulation.
Definition: onemodelproblem.hh:551
const TimeManager & timeManager() const
Returns TimeManager object used by the simulation.
Definition: onemodelproblem.hh:557
Scalar maxTimeStepSize() const
Returns the user specified maximum time step size.
Definition: onemodelproblem.hh:372
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: onemodelproblem.hh:403
VtkMultiWriter & resultWriter()
Definition: onemodelproblem.hh:639
void restart(double tRestart)
This method restores the complete state of the problem from disk.
Definition: onemodelproblem.hh:621
void postTimeStep()
Called by TimeManager whenever a solution for a timestep has been computed and the simulation time ha...
Definition: onemodelproblem.hh:356
void writeOutput(bool verbose=true)
Write the fields current solution into an VTK output file.
Definition: onemodelproblem.hh:461
Model & model()
Returns numerical model used for the problem.
Definition: onemodelproblem.hh:575
void advanceTimeLevel()
Called by the time manager after everything which can be done about the current time step is finished...
Definition: onemodelproblem.hh:364
void neumann(PrimaryVariables &values, const Intersection &intersection) const
Evaluate the boundary conditions for a neumann boundary segment.
Definition: onemodelproblem.hh:224
void setOutputInterval(int interval)
Sets the interval for Output.
Definition: onemodelproblem.hh:426
bool shouldWriteOutput() const
Returns true if the current solution should be written to disk (i.e. as a VTK file)
Definition: onemodelproblem.hh:438
VtkMultiWriter & resultWriter() const
Definition: onemodelproblem.hh:646
void episodeEnd()
Called when the end of an simulation episode is reached.
Definition: onemodelproblem.hh:476
OneModelProblem(TimeManager &timeManager, Grid &grid)
Constructs an object of type OneModelProblemProblem.
Definition: onemodelproblem.hh:123
const Variables & variables() const
Returns variables object.
Definition: onemodelproblem.hh:569
void setOutputTimeInterval(const Scalar timeInterval)
Sets a time interval for Output.
Definition: onemodelproblem.hh:415
void init()
Called by the TimeManager in order to initialize the problem.
Definition: onemodelproblem.hh:327
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
Evaluate the boundary conditions for a dirichlet control volume.
Definition: onemodelproblem.hh:203
Scalar nextTimeStepSize(Scalar dt)
Called by TimeManager whenever a solution for a timestep has been computed and the simulation time ha...
Definition: onemodelproblem.hh:392
void addOutputVtkFields()
Definition: onemodelproblem.hh:457
void source(PrimaryVariables &values, const Element &element) const
Evaluate the source term.
Definition: onemodelproblem.hh:262
void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition &globalPos) const
Specifies which kind of boundary condition should be used for which equation on a given boundary segm...
Definition: onemodelproblem.hh:168
const GlobalPosition & bBoxMin() const
The coordinate of the corner of the GridView's bounding box with the smallest values.
Definition: onemodelproblem.hh:538
const Model & model() const
Returns numerical model used for the problem.
Definition: onemodelproblem.hh:581
void boundaryTypes(BoundaryTypes &bcTypes, const Intersection &intersection) const
Specifies which kind of boundary condition should be used for which equation on a given boundary segm...
Definition: onemodelproblem.hh:154
void initial(PrimaryVariables &values, const Element &element) const
Evaluate the initial value for a control volume.
Definition: onemodelproblem.hh:299
const GlobalPosition & bBoxMax() const
The coordinate of the corner of the GridView's bounding box with the largest values.
Definition: onemodelproblem.hh:545
OneModelProblem(Grid &grid, bool verbose=true)
Constructs an object of type OneModelProblemProblem.
Definition: onemodelproblem.hh:94
const GridView & gridView() const
The GridView which used by the problem.
Definition: onemodelproblem.hh:519
const std::string & name() const
The problem name.
Definition: onemodelproblem.hh:499
void preTimeStep()
Called by TimeManager just before the time integration.
Definition: onemodelproblem.hh:338
void sourceAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: onemodelproblem.hh:282
void serialize()
This method writes the complete state of the problem to the harddisk.
Definition: onemodelproblem.hh:600
void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
Evaluate the boundary conditions for a neumann boundary segment.
Definition: onemodelproblem.hh:241
Variables & variables()
Returns variables object.
Definition: onemodelproblem.hh:563
Declares all properties used in Dumux.
Base file for properties related to sequential models.