28#ifndef EWOMS_FV_BASE_PROBLEM_HH
29#define EWOMS_FV_BASE_PROBLEM_HH
31#include <dune/common/fvector.hh>
51template <
class TypeTag,
class MyTypeTag>
67template<
class TypeTag>
74 static constexpr auto vtkOutputFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
92 dim = GridView::dimension,
93 dimWorld = GridView::dimensionworld
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;
100 using CoordScalar =
typename GridView::Grid::ctype;
101 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
122 , elementMapper_(gridView_,
Dune::mcmgElementLayout())
123 , vertexMapper_(gridView_,
Dune::mcmgVertexLayout())
124 , boundingBoxMin_(std::numeric_limits<double>::max())
125 , boundingBoxMax_(-std::numeric_limits<double>::max())
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]);
137 for (
unsigned i = 0; i < dim; ++i) {
138 boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
139 boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
143 const bool asyncVtkOutput =
144 simulator_.gridView().comm().size() == 1 &&
145 Parameters::Get<Parameters::EnableAsyncVtkOutput>();
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");
156 defaultVtkWriter_ = std::make_unique<VtkMultiWriter>(asyncVtkOutput,
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 "
216 std::string desc = Implementation::briefDescription();
221 return "Usage: " + std::string(argv[0]) +
" [OPTIONS]\n" + desc;
252 const std::string&)>,
253 std::set<std::string>&,
254 std::string& errorMsg,
260 errorMsg = std::string(
"Illegal parameter \"") + argv[paramIdx] +
"\".";
286 elementMapper_.update(gridView_);
287 vertexMapper_.update(gridView_);
290 defaultVtkWriter_->gridChanged();
303 template <
class Context>
308 {
throw std::logic_error(
"Problem does not provide a boundary() method"); }
320 template <
class Context>
325 {
throw std::logic_error(
"Problem does not provide a constraints() method"); }
339 template <
class Context>
344 {
throw std::logic_error(
"Problem does not provide a source() method"); }
356 template <
class Context>
361 {
throw std::logic_error(
"Problem does not provide a initial() method"); }
378 template <
class Context>
382 {
return asImp_().extrusionFactor(); }
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";
444 const auto& executionTimer =
simulator().executionTimer();
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());
457 if (
gridView().comm().rank() == 0) {
458 std::cout << std::setprecision(3)
459 <<
"Simulation of problem '" << asImp_().name() <<
"' finished.\n"
461 <<
"------------------------ Timing ------------------------\n"
462 <<
"Setup time: " << setupTime <<
" seconds"
464 <<
", " << setupTime / (executionTime + setupTime) * 100 <<
"%\n"
465 <<
"Simulation time: " << executionTime <<
" seconds"
467 <<
", " << executionTime / (executionTime + setupTime) * 100 <<
"%\n"
468 <<
" Linearization time: " << linearizeTime <<
" seconds"
470 <<
", " << linearizeTime / executionTime * 100 <<
"%\n"
471 <<
" Linear solve time: " << solveTime <<
" seconds"
473 <<
", " << solveTime / executionTime * 100 <<
"%\n"
474 <<
" Newton update time: " << updateTime <<
" seconds"
476 <<
", " << updateTime / executionTime * 100 <<
"%\n"
477 <<
" Pre/postprocess time: " << prePostProcessTime <<
" seconds"
479 <<
", " << prePostProcessTime / executionTime * 100 <<
"%\n"
480 <<
" Output write time: " << writeTime <<
" seconds"
482 <<
", " << writeTime / executionTime * 100 <<
"%\n"
483 <<
"First process' simulation CPU time: " << localCpuTime <<
" seconds"
485 <<
"Number of processes: " << numProcesses <<
"\n"
486 <<
"Threads per processes: " << threadsPerProcess <<
"\n"
487 <<
"Total CPU time: " << globalCpuTime <<
" seconds"
490 <<
"----------------------------------------------------------------\n"
501 const unsigned maxFails = asImp_().maxTimeIntegrationFailures();
502 Scalar minTimeStep = asImp_().minTimeStepSize();
504 std::string errorMessage;
505 for (
unsigned i = 0; i < maxFails; ++i) {
506 bool converged =
model().update();
511 const Scalar dt =
simulator().timeStepSize();
512 Scalar nextDt = dt / 2.0;
513 if (dt < minTimeStep * (1.0 + 1e-9)) {
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"
523 errorMessage =
"Time integration did not succeed with the minumum time step size of " +
528 else if (nextDt < minTimeStep) {
529 nextDt = minTimeStep;
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;
541 if (errorMessage.empty()) {
542 errorMessage =
"Newton solver didn't converge after " +
546 throw std::runtime_error(errorMessage);
553 {
return Parameters::Get<Parameters::MinTimeStepSize<Scalar>>(); }
560 {
return Parameters::Get<Parameters::MaxTimeStepDivisions>(); }
568 {
return Parameters::Get<Parameters::ContinueOnConvergenceError>(); }
590 if (dtNext <
simulator().maxTimeStepSize() &&
591 simulator().maxTimeStepSize() < dtNext * 2)
593 dtNext =
simulator().maxTimeStepSize() / 2 * 1.01;
609 return simulator().timeStepIndex() > 0 &&
630 {
model().advanceTimeLevel(); }
646 {
return gridView_; }
653 {
return boundingBoxMin_; }
660 {
return boundingBoxMax_; }
666 {
return vertexMapper_; }
672 {
return elementMapper_; }
678 {
return simulator_; }
684 {
return simulator_; }
690 {
return simulator_.model(); }
696 {
return simulator_.model(); }
702 {
return model().newtonMethod(); }
708 {
return model().newtonMethod(); }
745 template <
class Restarter>
749 defaultVtkWriter_->serialize(res);
763 template <
class Restarter>
767 defaultVtkWriter_->deserialize(res);
783 if (verbose &&
gridView().comm().rank() == 0) {
784 std::cout <<
"Writing visualization results for the current time step.\n"
791 defaultVtkWriter_->beginWrite(t);
792 model().prepareOutputFields();
793 model().appendOutputFields(*defaultVtkWriter_);
794 defaultVtkWriter_->endWrite(
false);
802 {
return *defaultVtkWriter_; }
808 {
return Parameters::Get<Parameters::EnableVtkOutput>(); }
812 Implementation& asImp_()
813 {
return *
static_cast<Implementation*
>(
this); }
816 const Implementation& asImp_()
const
817 {
return *
static_cast<const Implementation*
>(
this); }
820 const GridView gridView_;
821 ElementMapper elementMapper_;
822 VertexMapper vertexMapper_;
823 GlobalPosition boundingBoxMin_;
824 GlobalPosition boundingBoxMax_;
827 Simulator& simulator_;
828 std::unique_ptr<VtkMultiWriter> defaultVtkWriter_{};
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