fvbaseproblem.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 This file is part of the Open Porous Media project (OPM).
5
6 OPM 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 2 of the License, or
9 (at your option) any later version.
10
11 OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_FV_BASE_PROBLEM_HH
29#define EWOMS_FV_BASE_PROBLEM_HH
30
31#include <dune/common/fvector.hh>
32
36
39
41
42#include <functional>
43#include <iomanip>
44#include <iostream>
45#include <limits>
46#include <memory>
47#include <string>
48
49namespace Opm::Properties {
50
51template <class TypeTag, class MyTypeTag>
52struct NewtonMethod;
53
54} // namespace Opm::Properties
55
56namespace Opm {
57
67template<class TypeTag>
69{
70private:
71 using Implementation = GetPropType<TypeTag, Properties::Problem>;
73
74 static constexpr auto vtkOutputFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
76
82
85
90
91 enum {
92 dim = GridView::dimension,
93 dimWorld = GridView::dimensionworld
94 };
95
96 using Element = typename GridView::template Codim<0>::Entity;
97 using Vertex = typename GridView::template Codim<dim>::Entity;
98 using VertexIterator = typename GridView::template Codim<dim>::Iterator;
99
100 using CoordScalar = typename GridView::Grid::ctype;
101 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
102
103public:
104 // the default restriction and prolongation for adaptation is simply an empty one
106
107private:
108 // copying a problem is not a good idea
109 FvBaseProblem(const FvBaseProblem&) = delete;
110
111public:
119 explicit FvBaseProblem(Simulator& simulator)
120 : nextTimeStepSize_(0.0)
121 , gridView_(simulator.gridView())
122 , elementMapper_(gridView_, Dune::mcmgElementLayout())
123 , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
124 , boundingBoxMin_(std::numeric_limits<double>::max())
125 , boundingBoxMax_(-std::numeric_limits<double>::max())
126 , simulator_(simulator)
127 {
128 // calculate the bounding box of the local partition of the grid view
129 for (const auto& vertex : vertices(gridView_)) {
130 for (unsigned i = 0; i < dim; ++i) {
131 boundingBoxMin_[i] = std::min(boundingBoxMin_[i], vertex.geometry().corner(0)[i]);
132 boundingBoxMax_[i] = std::max(boundingBoxMax_[i], vertex.geometry().corner(0)[i]);
133 }
134 }
135
136 // communicate to get the bounding box of the whole domain
137 for (unsigned i = 0; i < dim; ++i) {
138 boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
139 boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
140 }
141
142 if (enableVtkOutput_()) {
143 const bool asyncVtkOutput =
144 simulator_.gridView().comm().size() == 1 &&
145 Parameters::Get<Parameters::EnableAsyncVtkOutput>();
146
147 // asynchonous VTK output currently does not work in conjunction with grid
148 // adaptivity because the async-IO code assumes that the grid stays
149 // constant. complain about that case.
150 const bool enableGridAdaptation = Parameters::Get<Parameters::EnableGridAdaptation>();
151 if (asyncVtkOutput && enableGridAdaptation) {
152 throw std::runtime_error("Asynchronous VTK output currently cannot be used "
153 "at the same time as grid adaptivity");
154 }
155
156 defaultVtkWriter_ = std::make_unique<VtkMultiWriter>(asyncVtkOutput,
157 gridView_,
158 asImp_().outputDir(),
159 asImp_().name());
160 }
161 }
162
167 static void registerParameters()
168 {
169 Model::registerParameters();
170 Parameters::Register<Parameters::MaxTimeStepSize<Scalar>>
171 ("The maximum size to which all time steps are limited to [s]");
172 Parameters::Register<Parameters::MinTimeStepSize<Scalar>>
173 ("The minimum size to which all time steps are limited to [s]");
174 Parameters::Register<Parameters::MaxTimeStepDivisions>
175 ("The maximum number of divisions by two of the timestep size "
176 "before the simulation bails out");
177 Parameters::Register<Parameters::EnableAsyncVtkOutput>
178 ("Dispatch a separate thread to write the VTK output");
179 Parameters::Register<Parameters::ContinueOnConvergenceError>
180 ("Continue with a non-converged solution instead of giving up "
181 "if we encounter a time step size smaller than the minimum time "
182 "step size.");
183 }
184
193 { return true; }
194
203 std::string outputDir() const
204 {
205 return simulatorOutputDir();
206 }
207
214 static std::string helpPreamble(int, const char **argv)
215 {
216 std::string desc = Implementation::briefDescription();
217 if (!desc.empty()) {
218 desc = desc + "\n";
219 }
220
221 return "Usage: " + std::string(argv[0]) + " [OPTIONS]\n" + desc;
222 }
223
230 static std::string briefDescription()
231 { return ""; }
232
233 // TODO (?): detailedDescription()
234
251 static int handlePositionalParameter(std::function<void(const std::string&,
252 const std::string&)>,
253 std::set<std::string>&,
254 std::string& errorMsg,
255 int,
256 const char** argv,
257 int paramIdx,
258 int)
259 {
260 errorMsg = std::string("Illegal parameter \"") + argv[paramIdx] + "\".";
261 return 0;
262 }
263
270 {}
271
276 void prefetch(const Element&) const
277 {
278 // do nothing by default
279 }
280
285 {
286 elementMapper_.update(gridView_);
287 vertexMapper_.update(gridView_);
288
289 if (enableVtkOutput_()) {
290 defaultVtkWriter_->gridChanged();
291 }
292 }
293
303 template <class Context>
304 void boundary(BoundaryRateVector&,
305 const Context&,
306 unsigned,
307 unsigned) const
308 { throw std::logic_error("Problem does not provide a boundary() method"); }
309
320 template <class Context>
321 void constraints(Constraints&,
322 const Context&,
323 unsigned,
324 unsigned) const
325 { throw std::logic_error("Problem does not provide a constraints() method"); }
326
339 template <class Context>
340 void source(RateVector&,
341 const Context&,
342 unsigned,
343 unsigned) const
344 { throw std::logic_error("Problem does not provide a source() method"); }
345
356 template <class Context>
357 void initial(PrimaryVariables&,
358 const Context&,
359 unsigned,
360 unsigned) const
361 { throw std::logic_error("Problem does not provide a initial() method"); }
362
378 template <class Context>
379 Scalar extrusionFactor(const Context&,
380 unsigned,
381 unsigned) const
382 { return asImp_().extrusionFactor(); }
383
384 Scalar extrusionFactor() const
385 { return 1.0; }
386
392 {}
393
398 {}
399
404 {}
405
410 {}
411
416 {}
417
425 {}
426
433 {
434 std::cerr << "The end of episode " << simulator().episodeIndex() + 1 << " has been "
435 << "reached, but the problem does not override the endEpisode() method. "
436 << "Doing nothing!\n";
437 }
438
442 void finalize()
443 {
444 const auto& executionTimer = simulator().executionTimer();
445
446 const Scalar executionTime = executionTimer.realTimeElapsed();
447 const Scalar setupTime = simulator().setupTimer().realTimeElapsed();
448 const Scalar prePostProcessTime = simulator().prePostProcessTimer().realTimeElapsed();
449 const Scalar localCpuTime = executionTimer.cpuTimeElapsed();
450 const Scalar globalCpuTime = executionTimer.globalCpuTimeElapsed();
451 const Scalar writeTime = simulator().writeTimer().realTimeElapsed();
452 const Scalar linearizeTime = simulator().linearizeTimer().realTimeElapsed();
453 const Scalar solveTime = simulator().solveTimer().realTimeElapsed();
454 const Scalar updateTime = simulator().updateTimer().realTimeElapsed();
455 const unsigned numProcesses = static_cast<unsigned>(this->gridView().comm().size());
456 const unsigned threadsPerProcess = ThreadManager::maxThreads();
457 if (gridView().comm().rank() == 0) {
458 std::cout << std::setprecision(3)
459 << "Simulation of problem '" << asImp_().name() << "' finished.\n"
460 << "\n"
461 << "------------------------ Timing ------------------------\n"
462 << "Setup time: " << setupTime << " seconds"
463 << humanReadableTime(setupTime)
464 << ", " << setupTime / (executionTime + setupTime) * 100 << "%\n"
465 << "Simulation time: " << executionTime << " seconds"
466 << humanReadableTime(executionTime)
467 << ", " << executionTime / (executionTime + setupTime) * 100 << "%\n"
468 << " Linearization time: " << linearizeTime << " seconds"
469 << humanReadableTime(linearizeTime)
470 << ", " << linearizeTime / executionTime * 100 << "%\n"
471 << " Linear solve time: " << solveTime << " seconds"
472 << humanReadableTime(solveTime)
473 << ", " << solveTime / executionTime * 100 << "%\n"
474 << " Newton update time: " << updateTime << " seconds"
475 << humanReadableTime(updateTime)
476 << ", " << updateTime / executionTime * 100 << "%\n"
477 << " Pre/postprocess time: " << prePostProcessTime << " seconds"
478 << humanReadableTime(prePostProcessTime)
479 << ", " << prePostProcessTime / executionTime * 100 << "%\n"
480 << " Output write time: " << writeTime << " seconds"
481 << humanReadableTime(writeTime)
482 << ", " << writeTime / executionTime * 100 << "%\n"
483 << "First process' simulation CPU time: " << localCpuTime << " seconds"
484 << humanReadableTime(localCpuTime) << "\n"
485 << "Number of processes: " << numProcesses << "\n"
486 << "Threads per processes: " << threadsPerProcess << "\n"
487 << "Total CPU time: " << globalCpuTime << " seconds"
488 << humanReadableTime(globalCpuTime) << "\n"
489 << "\n"
490 << "----------------------------------------------------------------\n"
491 << std::endl;
492 }
493 }
494
500 {
501 const unsigned maxFails = asImp_().maxTimeIntegrationFailures();
502 Scalar minTimeStep = asImp_().minTimeStepSize();
503
504 std::string errorMessage;
505 for (unsigned i = 0; i < maxFails; ++i) {
506 bool converged = model().update();
507 if (converged) {
508 return;
509 }
510
511 const Scalar dt = simulator().timeStepSize();
512 Scalar nextDt = dt / 2.0;
513 if (dt < minTimeStep * (1.0 + 1e-9)) {
514 if (asImp_().continueOnConvergenceError()) {
515 if (gridView().comm().rank() == 0) {
516 std::cout << "Newton solver did not converge with minimum time step of "
517 << dt << " seconds. Continuing with unconverged solution!\n"
518 << std::flush;
519 }
520 return;
521 }
522 else {
523 errorMessage = "Time integration did not succeed with the minumum time step size of " +
524 std::to_string(double(minTimeStep)) + " seconds. Giving up!";
525 break; // give up: we can't make the time step smaller anymore!
526 }
527 }
528 else if (nextDt < minTimeStep) {
529 nextDt = minTimeStep;
530 }
531 simulator().setTimeStepSize(nextDt);
532
533 // update failed
534 if (gridView().comm().rank() == 0) {
535 std::cout << "Newton solver did not converge with "
536 << "dt=" << dt << " seconds. Retrying with time step of "
537 << nextDt << " seconds\n" << std::flush;
538 }
539 }
540
541 if (errorMessage.empty()) {
542 errorMessage = "Newton solver didn't converge after " +
543 std::to_string(maxFails) + " time-step divisions. dt=" +
544 std::to_string(double(simulator().timeStepSize()));
545 }
546 throw std::runtime_error(errorMessage);
547 }
548
552 Scalar minTimeStepSize() const
553 { return Parameters::Get<Parameters::MinTimeStepSize<Scalar>>(); }
554
560 { return Parameters::Get<Parameters::MaxTimeStepDivisions>(); }
561
568 { return Parameters::Get<Parameters::ContinueOnConvergenceError>(); }
569
573 void setNextTimeStepSize(Scalar dt)
574 { nextTimeStepSize_ = dt; }
575
581 Scalar nextTimeStepSize() const
582 {
583 if (nextTimeStepSize_ > 0.0) {
584 return nextTimeStepSize_;
585 }
586
587 Scalar dtNext = std::min(Parameters::Get<Parameters::MaxTimeStepSize<Scalar>>(),
588 newtonMethod().suggestTimeStepSize(simulator().timeStepSize()));
589
590 if (dtNext < simulator().maxTimeStepSize() &&
591 simulator().maxTimeStepSize() < dtNext * 2)
592 {
593 dtNext = simulator().maxTimeStepSize() / 2 * 1.01;
594 }
595
596 return dtNext;
597 }
598
608 {
609 return simulator().timeStepIndex() > 0 &&
610 (simulator().timeStepIndex() % 10 == 0);
611 }
612
621 bool shouldWriteOutput() const
622 { return true; }
623
630 { model().advanceTimeLevel(); }
631
639 std::string name() const
640 { return "sim"; }
641
645 const GridView& gridView() const
646 { return gridView_; }
647
652 const GlobalPosition& boundingBoxMin() const
653 { return boundingBoxMin_; }
654
659 const GlobalPosition& boundingBoxMax() const
660 { return boundingBoxMax_; }
661
665 const VertexMapper& vertexMapper() const
666 { return vertexMapper_; }
667
671 const ElementMapper& elementMapper() const
672 { return elementMapper_; }
673
677 Simulator& simulator()
678 { return simulator_; }
679
683 const Simulator& simulator() const
684 { return simulator_; }
685
689 Model& model()
690 { return simulator_.model(); }
691
695 const Model& model() const
696 { return simulator_.model(); }
697
701 NewtonMethod& newtonMethod()
702 { return model().newtonMethod(); }
703
707 const NewtonMethod& newtonMethod() const
708 { return model().newtonMethod(); }
709 // \}
710
716 {
718 }
719
728 {
729 return 0;
730 }
731
745 template <class Restarter>
746 void serialize(Restarter& res)
747 {
748 if (enableVtkOutput_()) {
749 defaultVtkWriter_->serialize(res);
750 }
751 }
752
763 template <class Restarter>
764 void deserialize(Restarter& res)
765 {
766 if (enableVtkOutput_()) {
767 defaultVtkWriter_->deserialize(res);
768 }
769 }
770
777 void writeOutput(bool verbose = true)
778 {
779 if (!enableVtkOutput_()) {
780 return;
781 }
782
783 if (verbose && gridView().comm().rank() == 0) {
784 std::cout << "Writing visualization results for the current time step.\n"
785 << std::flush;
786 }
787
788 // calculate the time _after_ the time was updated
789 const Scalar t = simulator().time() + simulator().timeStepSize();
790
791 defaultVtkWriter_->beginWrite(t);
792 model().prepareOutputFields();
793 model().appendOutputFields(*defaultVtkWriter_);
794 defaultVtkWriter_->endWrite(false);
795 }
796
802 { return *defaultVtkWriter_; }
803
804protected:
806
807 bool enableVtkOutput_() const
808 { return Parameters::Get<Parameters::EnableVtkOutput>(); }
809
810private:
812 Implementation& asImp_()
813 { return *static_cast<Implementation*>(this); }
814
816 const Implementation& asImp_() const
817 { return *static_cast<const Implementation*>(this); }
818
819 // Grid management stuff
820 const GridView gridView_;
821 ElementMapper elementMapper_;
822 VertexMapper vertexMapper_;
823 GlobalPosition boundingBoxMin_;
824 GlobalPosition boundingBoxMax_;
825
826 // Attributes required for the actual simulation
827 Simulator& simulator_;
828 std::unique_ptr<VtkMultiWriter> defaultVtkWriter_{};
829};
830
831} // namespace Opm
832
833#endif
Definition: restrictprolong.hh:142
Base class for all problems which use a finite volume spatial discretization.
Definition: fvbaseproblem.hh:69
static int handlePositionalParameter(std::function< void(const std::string &, const std::string &)>, std::set< std::string > &, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Handles positional command line parameters.
Definition: fvbaseproblem.hh:251
std::string name() const
The problem name.
Definition: fvbaseproblem.hh:639
const Simulator & simulator() const
Returns Simulator object used by the simulation.
Definition: fvbaseproblem.hh:683
void writeOutput(bool verbose=true)
Write the relevant secondary variables of the current solution into an VTK output file.
Definition: fvbaseproblem.hh:777
void boundary(BoundaryRateVector &, const Context &, unsigned, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition: fvbaseproblem.hh:304
void finalize()
Called after the simulation has been run sucessfully.
Definition: fvbaseproblem.hh:442
Simulator & simulator()
Returns Simulator object used by the simulation.
Definition: fvbaseproblem.hh:677
const Model & model() const
Returns numerical model used for the problem.
Definition: fvbaseproblem.hh:695
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: fvbaseproblem.hh:167
void source(RateVector &, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: fvbaseproblem.hh:340
bool recycleFirstIterationStorage() const
Return if the storage term of the first iteration is identical to the storage term for the solution o...
Definition: fvbaseproblem.hh:192
unsigned maxTimeIntegrationFailures() const
Returns the maximum number of subsequent failures for the time integration before giving up.
Definition: fvbaseproblem.hh:559
unsigned markForGridAdaptation()
Mark grid cells for refinement or coarsening.
Definition: fvbaseproblem.hh:727
void serialize(Restarter &res)
This method writes the complete state of the problem to the harddisk.
Definition: fvbaseproblem.hh:746
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition: fvbaseproblem.hh:409
void endTimeStep()
Called by the simulator after each time integration.
Definition: fvbaseproblem.hh:424
FvBaseProblem(Simulator &simulator)
Definition: fvbaseproblem.hh:119
const GlobalPosition & boundingBoxMax() const
The coordinate of the corner of the GridView's bounding box with the largest values.
Definition: fvbaseproblem.hh:659
RestrictProlongOperator restrictProlongOperator()
return restriction and prolongation operator
Definition: fvbaseproblem.hh:715
void timeIntegration()
Called by Opm::Simulator in order to do a time integration on the model.
Definition: fvbaseproblem.hh:499
Model & model()
Returns numerical model used for the problem.
Definition: fvbaseproblem.hh:689
Scalar extrusionFactor() const
Definition: fvbaseproblem.hh:384
void endEpisode()
Called when the end of an simulation episode is reached.
Definition: fvbaseproblem.hh:432
void beginTimeStep()
Called by the simulator before each time integration.
Definition: fvbaseproblem.hh:403
static std::string helpPreamble(int, const char **argv)
Returns the string that is printed before the list of command line parameters in the help message.
Definition: fvbaseproblem.hh:214
void deserialize(Restarter &res)
This method restores the complete state of the problem from disk.
Definition: fvbaseproblem.hh:764
void prefetch(const Element &) const
Allows to improve the performance by prefetching all data which is associated with a given element.
Definition: fvbaseproblem.hh:276
Scalar extrusionFactor(const Context &, unsigned, unsigned) const
Return how much the domain is extruded at a given sub-control volume.
Definition: fvbaseproblem.hh:379
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: fvbaseproblem.hh:607
Scalar nextTimeStepSize() const
Called by Opm::Simulator whenever a solution for a time step has been computed and the simulation tim...
Definition: fvbaseproblem.hh:581
Scalar minTimeStepSize() const
Returns the minimum allowable size of a time step.
Definition: fvbaseproblem.hh:552
void beginEpisode()
Called at the beginning of an simulation episode.
Definition: fvbaseproblem.hh:397
void initial(PrimaryVariables &, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition: fvbaseproblem.hh:357
std::string outputDir() const
Determine the directory for simulation output.
Definition: fvbaseproblem.hh:203
Scalar nextTimeStepSize_
Definition: fvbaseproblem.hh:805
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition: fvbaseproblem.hh:665
void gridChanged()
Handle changes of the grid.
Definition: fvbaseproblem.hh:284
const GlobalPosition & boundingBoxMin() const
The coordinate of the corner of the GridView's bounding box with the smallest values.
Definition: fvbaseproblem.hh:652
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition: fvbaseproblem.hh:230
void endIteration()
Called by the simulator after each Newton-Raphson update.
Definition: fvbaseproblem.hh:415
void setNextTimeStepSize(Scalar dt)
Impose the next time step size to be used externally.
Definition: fvbaseproblem.hh:573
bool continueOnConvergenceError() const
Returns if we should continue with a non-converged solution instead of giving up if we encounter a ti...
Definition: fvbaseproblem.hh:567
void constraints(Constraints &, const Context &, unsigned, unsigned) const
Evaluate the constraints for a control volume.
Definition: fvbaseproblem.hh:321
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: fvbaseproblem.hh:269
const NewtonMethod & newtonMethod() const
Returns object which implements the Newton method.
Definition: fvbaseproblem.hh:707
const GridView & gridView() const
The GridView which used by the problem.
Definition: fvbaseproblem.hh:645
bool enableVtkOutput_() const
Definition: fvbaseproblem.hh:807
EmptyRestrictProlong RestrictProlongOperator
Definition: fvbaseproblem.hh:105
NewtonMethod & newtonMethod()
Returns object which implements the Newton method.
Definition: fvbaseproblem.hh:701
bool shouldWriteOutput() const
Returns true if the current solution should be written to disk (i.e. as a VTK file)
Definition: fvbaseproblem.hh:621
void initialSolutionApplied()
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition: fvbaseproblem.hh:391
void advanceTimeLevel()
Called by the simulator after everything which can be done about the current time step is finished an...
Definition: fvbaseproblem.hh:629
VtkMultiWriter & defaultVtkWriter() const
Method to retrieve the VTK writer which should be used to write the default ouput after each time ste...
Definition: fvbaseproblem.hh:801
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition: fvbaseproblem.hh:671
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hpp:66
Simplifies writing multi-file VTK datasets.
Definition: vtkmultiwriter.hh:65
Declare the properties used by the infrastructure code of the finite volume discretizations.
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: fvbaseprimaryvariables.hh:141
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:185
Definition: blackoilmodel.hh:79
Definition: blackoilboundaryratevector.hh:39
std::string humanReadableTime(double timeInSeconds, bool isAmendment=true)
Given a time step size in seconds, return it in a format which is more easily parsable by humans.
std::string simulatorOutputDir()
Determine and check the configured directory for simulation output.
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:233
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Specify the maximum size of a time integration [s].
Definition: fvbaseparameters.hh:106