My Project
Loading...
Searching...
No Matches
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
35#include <opm/models/discretization/common/restrictprolong.hh>
36
39
40#include <opm/models/utils/simulatorutils.hpp>
41
42#include <functional>
43#include <iostream>
44#include <limits>
45#include <string>
46
47namespace Opm::Properties {
48
49template <class TypeTag, class MyTypeTag>
50struct NewtonMethod;
51
52} // namespace Opm::Properties
53
54namespace Opm {
55
65template<class TypeTag>
67{
68private:
69 using Implementation = GetPropType<TypeTag, Properties::Problem>;
71
72 static const int vtkOutputFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
74
80
83
88
89 enum {
90 dim = GridView::dimension,
91 dimWorld = GridView::dimensionworld
92 };
93
94 using Element = typename GridView::template Codim<0>::Entity;
95 using Vertex = typename GridView::template Codim<dim>::Entity;
96 using VertexIterator = typename GridView::template Codim<dim>::Iterator;
97
98 using CoordScalar = typename GridView::Grid::ctype;
99 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
100
101public:
102 // the default restriction and prolongation for adaptation is simply an empty one
104
105private:
106 // copying a problem is not a good idea
107 FvBaseProblem(const FvBaseProblem& ) = delete;
108
109public:
117 explicit FvBaseProblem(Simulator& simulator)
118 : nextTimeStepSize_(0.0)
119 , gridView_(simulator.gridView())
120 , elementMapper_(gridView_, Dune::mcmgElementLayout())
121 , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
122 , boundingBoxMin_(std::numeric_limits<double>::max())
123 , boundingBoxMax_(-std::numeric_limits<double>::max())
124 , simulator_(simulator)
125 , defaultVtkWriter_(0)
126 {
127 // calculate the bounding box of the local partition of the grid view
128 VertexIterator vIt = gridView_.template begin<dim>();
129 const VertexIterator vEndIt = gridView_.template end<dim>();
130 for (; vIt!=vEndIt; ++vIt) {
131 for (unsigned i=0; i<dim; i++) {
132 boundingBoxMin_[i] = std::min(boundingBoxMin_[i], vIt->geometry().corner(0)[i]);
133 boundingBoxMax_[i] = std::max(boundingBoxMax_[i], vIt->geometry().corner(0)[i]);
134 }
135 }
136
137 // communicate to get the bounding box of the whole domain
138 for (unsigned i = 0; i < dim; ++i) {
139 boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
140 boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
141 }
142
143 if (enableVtkOutput_()) {
144 bool asyncVtkOutput =
145 simulator_.gridView().comm().size() == 1 &&
146 Parameters::Get<Parameters::EnableAsyncVtkOutput>();
147
148 // asynchonous VTK output currently does not work in conjunction with grid
149 // adaptivity because the async-IO code assumes that the grid stays
150 // constant. complain about that case.
151 bool enableGridAdaptation = Parameters::Get<Parameters::EnableGridAdaptation>();
152 if (asyncVtkOutput && enableGridAdaptation)
153 throw std::runtime_error("Asynchronous VTK output currently cannot be used "
154 "at the same time as grid adaptivity");
155
156 std::string outputDir = asImp_().outputDir();
157
158 defaultVtkWriter_ =
159 new VtkMultiWriter(asyncVtkOutput, gridView_, outputDir, asImp_().name());
160 }
161 }
162
164 { delete defaultVtkWriter_; }
165
170 static void registerParameters()
171 {
172 Model::registerParameters();
173 Parameters::Register<Parameters::MaxTimeStepSize<Scalar>>
174 ("The maximum size to which all time steps are limited to [s]");
175 Parameters::Register<Parameters::MinTimeStepSize<Scalar>>
176 ("The minimum size to which all time steps are limited to [s]");
177 Parameters::Register<Parameters::MaxTimeStepDivisions>
178 ("The maximum number of divisions by two of the timestep size "
179 "before the simulation bails out");
180 Parameters::Register<Parameters::EnableAsyncVtkOutput>
181 ("Dispatch a separate thread to write the VTK output");
182 Parameters::Register<Parameters::ContinueOnConvergenceError>
183 ("Continue with a non-converged solution instead of giving up "
184 "if we encounter a time step size smaller than the minimum time "
185 "step size.");
186 }
187
196 { return true; }
197
206 std::string outputDir() const
207 {
208 return simulatorOutputDir();
209 }
210
217 static std::string helpPreamble(int,
218 const char **argv)
219 {
220 std::string desc = Implementation::briefDescription();
221 if (!desc.empty())
222 desc = desc + "\n";
223
224 return
225 "Usage: "+std::string(argv[0]) + " [OPTIONS]\n"
226 + desc;
227 }
228
235 static std::string briefDescription()
236 { return ""; }
237
238 // TODO (?): detailedDescription()
239
256 static int handlePositionalParameter(std::function<void(const std::string&,
257 const std::string&)>,
258 std::set<std::string>&,
259 std::string& errorMsg,
260 int,
261 const char** argv,
262 int paramIdx,
263 int)
264 {
265 errorMsg = std::string("Illegal parameter \"")+argv[paramIdx]+"\".";
266 return 0;
267 }
268
275 { }
276
281 void prefetch(const Element&) const
282 {
283 // do nothing by default
284 }
285
290 {
291 elementMapper_.update(gridView_);
292 vertexMapper_.update(gridView_);
293
294 if (enableVtkOutput_())
295 defaultVtkWriter_->gridChanged();
296 }
297
307 template <class Context>
308 void boundary(BoundaryRateVector&,
309 const Context&,
310 unsigned,
311 unsigned) const
312 { throw std::logic_error("Problem does not provide a boundary() method"); }
313
324 template <class Context>
325 void constraints(Constraints&,
326 const Context&,
327 unsigned,
328 unsigned) const
329 { throw std::logic_error("Problem does not provide a constraints() method"); }
330
343 template <class Context>
344 void source(RateVector&,
345 const Context&,
346 unsigned,
347 unsigned) const
348 { throw std::logic_error("Problem does not provide a source() method"); }
349
360 template <class Context>
361 void initial(PrimaryVariables&,
362 const Context&,
363 unsigned,
364 unsigned) const
365 { throw std::logic_error("Problem does not provide a initial() method"); }
366
382 template <class Context>
383 Scalar extrusionFactor(const Context&,
384 unsigned,
385 unsigned) const
386 { return asImp_().extrusionFactor(); }
387
388 Scalar extrusionFactor() const
389 { return 1.0; }
390
397
402 { }
403
408 { }
409
414 { }
415
420 { }
421
429 { }
430
437 {
438 std::cerr << "The end of episode " << simulator().episodeIndex() + 1 << " has been "
439 << "reached, but the problem does not override the endEpisode() method. "
440 << "Doing nothing!\n";
441 }
442
446 void finalize()
447 {
448 const auto& executionTimer = simulator().executionTimer();
449
450 Scalar executionTime = executionTimer.realTimeElapsed();
451 Scalar setupTime = simulator().setupTimer().realTimeElapsed();
452 Scalar prePostProcessTime = simulator().prePostProcessTimer().realTimeElapsed();
453 Scalar localCpuTime = executionTimer.cpuTimeElapsed();
454 Scalar globalCpuTime = executionTimer.globalCpuTimeElapsed();
455 Scalar writeTime = simulator().writeTimer().realTimeElapsed();
456 Scalar linearizeTime = simulator().linearizeTimer().realTimeElapsed();
457 Scalar solveTime = simulator().solveTimer().realTimeElapsed();
458 Scalar updateTime = simulator().updateTimer().realTimeElapsed();
459 unsigned numProcesses = static_cast<unsigned>(this->gridView().comm().size());
461 if (gridView().comm().rank() == 0) {
462 std::cout << std::setprecision(3)
463 << "Simulation of problem '" << asImp_().name() << "' finished.\n"
464 << "\n"
465 << "------------------------ Timing ------------------------\n"
466 << "Setup time: " << setupTime << " seconds"
467 << humanReadableTime(setupTime)
468 << ", " << setupTime/(executionTime + setupTime)*100 << "%\n"
469 << "Simulation time: " << executionTime << " seconds"
471 << ", " << executionTime/(executionTime + setupTime)*100 << "%\n"
472 << " Linearization time: " << linearizeTime << " seconds"
474 << ", " << linearizeTime/executionTime*100 << "%\n"
475 << " Linear solve time: " << solveTime << " seconds"
477 << ", " << solveTime/executionTime*100 << "%\n"
478 << " Newton update time: " << updateTime << " seconds"
480 << ", " << updateTime/executionTime*100 << "%\n"
481 << " Pre/postprocess time: " << prePostProcessTime << " seconds"
483 << ", " << prePostProcessTime/executionTime*100 << "%\n"
484 << " Output write time: " << writeTime << " seconds"
486 << ", " << writeTime/executionTime*100 << "%\n"
487 << "First process' simulation CPU time: " << localCpuTime << " seconds"
489 << "Number of processes: " << numProcesses << "\n"
490 << "Threads per processes: " << threadsPerProcess << "\n"
491 << "Total CPU time: " << globalCpuTime << " seconds"
493 << "\n"
494 << "----------------------------------------------------------------\n"
495 << std::endl;
496 }
497 }
498
504 {
505 unsigned maxFails = asImp_().maxTimeIntegrationFailures();
506 Scalar minTimeStepSize = asImp_().minTimeStepSize();
507
508 std::string errorMessage;
509 for (unsigned i = 0; i < maxFails; ++i) {
510 bool converged = model().update();
511 if (converged)
512 return;
513
514 Scalar dt = simulator().timeStepSize();
515 Scalar nextDt = dt / 2.0;
516 if (dt < minTimeStepSize*(1 + 1e-9)) {
517 if (asImp_().continueOnConvergenceError()) {
518 if (gridView().comm().rank() == 0)
519 std::cout << "Newton solver did not converge with minimum time step of "
520 << dt << " seconds. Continuing with unconverged solution!\n"
521 << std::flush;
522 return;
523 }
524 else {
526 "Time integration did not succeed with the minumum time step size of "
527 + std::to_string(double(minTimeStepSize)) + " seconds. Giving up!";
528 break; // give up: we can't make the time step smaller anymore!
529 }
530 }
531 else if (nextDt < minTimeStepSize)
533 simulator().setTimeStepSize(nextDt);
534
535 // update failed
536 if (gridView().comm().rank() == 0)
537 std::cout << "Newton solver did not converge with "
538 << "dt=" << dt << " seconds. Retrying with time step of "
539 << nextDt << " seconds\n" << std::flush;
540 }
541
542 if (errorMessage.empty())
544 "Newton solver didn't converge after "
545 +std::to_string(maxFails)+" time-step divisions. dt="
546 +std::to_string(double(simulator().timeStepSize()));
547 throw std::runtime_error(errorMessage);
548 }
549
553 Scalar minTimeStepSize() const
554 { return Parameters::Get<Parameters::MinTimeStepSize<Scalar>>(); }
555
561 { return Parameters::Get<Parameters::MaxTimeStepDivisions>(); }
562
569 { return Parameters::Get<Parameters::ContinueOnConvergenceError>(); }
570
575 { nextTimeStepSize_ = dt; }
576
582 Scalar nextTimeStepSize() const
583 {
584 if (nextTimeStepSize_ > 0.0)
585 return nextTimeStepSize_;
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
719
728 {
729 return 0;
730 }
731
745 template <class Restarter>
747 {
748 if (enableVtkOutput_())
749 defaultVtkWriter_->serialize(res);
750 }
751
762 template <class Restarter>
764 {
765 if (enableVtkOutput_())
766 defaultVtkWriter_->deserialize(res);
767 }
768
775 void writeOutput(bool verbose = true)
776 {
777 if (!enableVtkOutput_())
778 return;
779
780 if (verbose && gridView().comm().rank() == 0)
781 std::cout << "Writing visualization results for the current time step.\n"
782 << std::flush;
783
784 // calculate the time _after_ the time was updated
785 Scalar t = simulator().time() + simulator().timeStepSize();
786
787 defaultVtkWriter_->beginWrite(t);
788 model().prepareOutputFields();
789 model().appendOutputFields(*defaultVtkWriter_);
790 defaultVtkWriter_->endWrite();
791
792 }
793
799 { return defaultVtkWriter_; }
800
801protected:
802 Scalar nextTimeStepSize_;
803
804 bool enableVtkOutput_() const
805 { return Parameters::Get<Parameters::EnableVtkOutput>(); }
806
807private:
809 Implementation& asImp_()
810 { return *static_cast<Implementation *>(this); }
811
813 const Implementation& asImp_() const
814 { return *static_cast<const Implementation *>(this); }
815
816 // Grid management stuff
817 const GridView gridView_;
818 ElementMapper elementMapper_;
819 VertexMapper vertexMapper_;
820 GlobalPosition boundingBoxMin_;
821 GlobalPosition boundingBoxMax_;
822
823 // Attributes required for the actual simulation
824 Simulator& simulator_;
825 mutable VtkMultiWriter *defaultVtkWriter_;
826};
827
828} // namespace Opm
829
830#endif
Definition restrictprolong.hh:142
Base class for all problems which use a finite volume spatial discretization.
Definition fvbaseproblem.hh:67
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:256
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:775
void boundary(BoundaryRateVector &, const Context &, unsigned, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition fvbaseproblem.hh:308
void finalize()
Called after the simulation has been run sucessfully.
Definition fvbaseproblem.hh:446
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:170
void source(RateVector &, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition fvbaseproblem.hh:344
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:195
unsigned maxTimeIntegrationFailures() const
Returns the maximum number of subsequent failures for the time integration before giving up.
Definition fvbaseproblem.hh:560
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:413
void endTimeStep()
Called by the simulator after each time integration.
Definition fvbaseproblem.hh:428
FvBaseProblem(Simulator &simulator)
Definition fvbaseproblem.hh:117
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:503
Model & model()
Returns numerical model used for the problem.
Definition fvbaseproblem.hh:689
void endEpisode()
Called when the end of an simulation episode is reached.
Definition fvbaseproblem.hh:436
void beginTimeStep()
Called by the simulator before each time integration.
Definition fvbaseproblem.hh:407
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:217
void deserialize(Restarter &res)
This method restores the complete state of the problem from disk.
Definition fvbaseproblem.hh:763
void prefetch(const Element &) const
Allows to improve the performance by prefetching all data which is associated with a given element.
Definition fvbaseproblem.hh:281
Scalar extrusionFactor(const Context &, unsigned, unsigned) const
Return how much the domain is extruded at a given sub-control volume.
Definition fvbaseproblem.hh:383
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:582
Scalar minTimeStepSize() const
Returns the minimum allowable size of a time step.
Definition fvbaseproblem.hh:553
void beginEpisode()
Called at the beginning of an simulation episode.
Definition fvbaseproblem.hh:401
void initial(PrimaryVariables &, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition fvbaseproblem.hh:361
std::string outputDir() const
Determine the directory for simulation output.
Definition fvbaseproblem.hh:206
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:289
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:235
void endIteration()
Called by the simulator after each Newton-Raphson update.
Definition fvbaseproblem.hh:419
void setNextTimeStepSize(Scalar dt)
Impose the next time step size to be used externally.
Definition fvbaseproblem.hh:574
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:568
void constraints(Constraints &, const Context &, unsigned, unsigned) const
Evaluate the constraints for a control volume.
Definition fvbaseproblem.hh:325
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition fvbaseproblem.hh:274
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
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.
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:395
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:798
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:66
void endWrite(bool onlyDiscard=false) override
Finalizes the current writer.
Definition vtkmultiwriter.hh:384
void deserialize(Restarter &res)
Read the multi-writer's state from a restart file.
Definition vtkmultiwriter.hh:437
void gridChanged()
Updates the internal data structures after mesh refinement.
Definition vtkmultiwriter.hh:165
void beginWrite(double t) override
Called whenever a new time step must be written.
Definition vtkmultiwriter.hh:174
void serialize(Restarter &res)
Write the multi-writer's state to a restart file.
Definition vtkmultiwriter.hh:403
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.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
std::string humanReadableTime(double timeInSeconds, bool isAmendment)
Given a time step size in seconds, return it in a format which is more easily parsable by humans.
Definition simulatorutils.cpp:45
std::string simulatorOutputDir()
Determine and check the configured directory for simulation output.
Definition simulatorutils.cpp:119
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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:235
Load or save a state of a problem to/from the harddisk.
Specify the maximum size of a time integration [s].
Definition fvbaseparameters.hh:106
Simplifies writing multi-file VTK datasets.