OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1512 from andlaus/ewoms_parameters

Browse Source 
Switch flow to the eWoms parameter system
This commit is contained in:
Atgeirr Flø Rasmussen 2018-08-16 10:36:11 +02:00 committed by GitHub
commit 31280d579b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
21 changed files with 1581 additions and 688 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -275,6 +275,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/BlackoilModelBase_impl.hpp
opm/autodiff/BlackoilModelEnums.hpp
opm/autodiff/BlackoilModelParameters.hpp
opm/autodiff/BlackoilModelParametersEbos.hpp
opm/autodiff/BlackoilPressureModel.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/Compat.hpp
@ -297,6 +298,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/GridInit.hpp
opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/ISTLSolver.hpp
opm/autodiff/ISTLSolverEbos.hpp
opm/autodiff/IterationReport.hpp
opm/autodiff/moduleVersion.hpp
opm/autodiff/multiPhaseUpwind.hpp

29
compareECLFiles.cmake
View File

@ -218,14 +218,19 @@ add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater
REL_TOL ${rel_tol}
DIR aquifer-oilwater)
foreach(SIM flow flow_legacy)
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1
SIMULATOR ${SIM}
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS tolerance_wells=1e-6 max_iter=20)
endforeach()
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS --tolerance-wells=1e-6 --flow-newton-max-iterations=20)
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1
SIMULATOR flow_legacy
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS tolerance_wells=1e-6 max_iter=20)
foreach(SIM flow flow_legacy)
add_test_compareECLFiles(CASENAME spe9
@ -246,14 +251,14 @@ add_test_compareECLFiles(CASENAME msw_2d_h
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS use_multisegment_well=true)
TEST_ARGS --use-multisegment-well=true)
add_test_compareECLFiles(CASENAME msw_3d_hfa
FILENAME 3D_MSW
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${rel_tol}
TEST_ARGS use_multisegment_well=true)
TEST_ARGS --use-multisegment-well=true)
add_test_compareECLFiles(CASENAME polymer_oilwater
FILENAME 2D_OILWATER_POLYMER
@ -267,14 +272,14 @@ add_test_compareECLFiles(CASENAME polymer_simple2D
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS max_iter=20)
TEST_ARGS --flow-newton-max-iterations=20)
add_test_compareECLFiles(CASENAME spe5
FILENAME SPE5CASE1
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${coarse_rel_tol}
TEST_ARGS max_iter=20)
TEST_ARGS --flow-newton-max-iterations=20)
add_test_compareECLFiles(CASENAME wecon_wtest
FILENAME 3D_WECON

59
examples/flow.cpp
View File

@ -28,6 +28,10 @@
#include <opm/simulators/flow_ebos_energy.hpp>
#include <opm/simulators/flow_ebos_oilwater_polymer.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/autodiff/FlowMainEbos.hpp>
#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/parametersystem.hh>
#include <opm/autodiff/MissingFeatures.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/material/common/ResetLocale.hpp>
@ -43,6 +47,14 @@
#include <dune/common/parallel/mpihelper.hh>
#endif
BEGIN_PROPERTIES
// this is a dummy type tag that is used to setup the parameters before the actual
// simulator.
NEW_TYPE_TAG(FlowEarlyBird, INHERITS_FROM(EclFlowProblem));
END_PROPERTIES
namespace detail
{
boost::filesystem::path simulationCaseName( const std::string& casename ) {
@ -90,38 +102,35 @@ int main(int argc, char** argv)
int mpiRank = mpiHelper.rank();
#endif
const bool outputCout = (mpiRank == 0);
// we always want to use the default locale, and thus spare us the trouble
// with incorrect locale settings.
Opm::resetLocale();
Opm::ParameterGroup param(argc, argv, false, outputCout);
// this is a work-around for a catch 22: we do not know what code path to use without
// parsing the deck, but we don't know the deck without having access to the
// parameters and this requires to know the type tag to be used. To solve this, we
// use a type tag just for parsing the parameters before we instantiate the actual
// simulator object. (Which parses the parameters again, but since this is done in an
// identical manner it does not matter.)
typedef TTAG(FlowEarlyBird) PreTypeTag;
typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem;
// See if a deck was specified on the command line.
if (!param.unhandledArguments().empty()) {
if (param.unhandledArguments().size() != 1) {
if (outputCout)
std::cerr << "You can only specify a single input deck on the command line.\n";
return EXIT_FAILURE;
} else {
const auto casename = detail::simulationCaseName( param.unhandledArguments()[ 0 ] );
param.insertParameter("deck_filename", casename.string() );
}
}
PreProblem::setBriefDescription("Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project.");
// We must have an input deck. Grid and props will be read from that.
if (!param.has("deck_filename")) {
if (outputCout)
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return EXIT_FAILURE;
}
int status = Opm::FlowMainEbos<PreTypeTag>::setupParameters_(argc, argv);
if (status != 0)
// if setupParameters_ returns a value smaller than 0, there was no error, but
// the program should abort. This is the case e.g. for the --help and the
// --print-properties parameters.
return (status >= 0)?status:0;
std::string deckFilename = param.get<std::string>("deck_filename");
bool outputCout = false;
if (mpiRank == 0)
outputCout = EWOMS_GET_PARAM(PreTypeTag, bool, EnableTerminalOutput);
std::string deckFilename = EWOMS_GET_PARAM(PreTypeTag, std::string, EclDeckFileName);
typedef typename GET_PROP_TYPE(PreTypeTag, Vanguard) PreVanguard;
deckFilename = PreVanguard::canonicalDeckPath(deckFilename).string();
// Create Deck and EclipseState.
try {

19
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -27,7 +27,10 @@
#include <ebos/eclproblem.hh>
#include <ewoms/common/start.hh>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
#include <opm/autodiff/NonlinearSolverEbos.hpp>
#include <opm/autodiff/BlackoilModelParametersEbos.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp>
#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
@ -47,7 +50,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/autodiff/ISTLSolver.hpp>
#include <opm/autodiff/ISTLSolverEbos.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <dune/istl/owneroverlapcopy.hh>
@ -65,10 +68,9 @@
//#include <fstream>
BEGIN_PROPERTIES
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem, FlowNonLinearSolver, FlowIstlSolver, FlowModelParameters, FlowTimeSteppingParameters));
SET_STRING_PROP(EclFlowProblem, OutputDir, "");
SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
@ -82,7 +84,8 @@ SET_BOOL_PROP(EclFlowProblem, EnablePolymer, false);
SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true);
SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false);
}}
END_PROPERTIES
namespace Opm {
/// A model implementation for three-phase black oil.
@ -98,7 +101,7 @@ namespace Opm {
// --------- Types and enums ---------
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilModelParametersEbos<TypeTag> ModelParameters;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
@ -124,7 +127,7 @@ namespace Opm {
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx > ISTLSolverType;
typedef ISTLSolverEbos<TypeTag> ISTLSolverType;
//typedef typename SolutionVector :: value_type PrimaryVariables ;
// --------- Public methods ---------

229
opm/autodiff/BlackoilModelParametersEbos.hpp Normal file
View File

@ -0,0 +1,229 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
#define OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED
#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/parametersystem.hh>
#include <string>
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowModelParameters);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(EclDeckFileName);
NEW_PROP_TAG(DpMaxRel);
NEW_PROP_TAG(DsMax);
NEW_PROP_TAG(DrMaxRel);
NEW_PROP_TAG(DbhpMaxRel);
NEW_PROP_TAG(DwellFractionMax);
NEW_PROP_TAG(MaxResidualAllowed);
NEW_PROP_TAG(ToleranceMb);
NEW_PROP_TAG(ToleranceCnv);
NEW_PROP_TAG(ToleranceCnvRelaxed);
NEW_PROP_TAG(ToleranceWells);
NEW_PROP_TAG(ToleranceWellControl);
NEW_PROP_TAG(MaxWelleqIter);
NEW_PROP_TAG(UseMultisegmentWell);
NEW_PROP_TAG(MaxSinglePrecisionDays);
NEW_PROP_TAG(MaxStrictIter);
NEW_PROP_TAG(SolveWelleqInitially);
NEW_PROP_TAG(UpdateEquationsScaling);
NEW_PROP_TAG(UseUpdateStabilization);
NEW_PROP_TAG(MatrixAddWellContributions);
NEW_PROP_TAG(PreconditionerAddWellContributions);
// parameters for multisegment wells
NEW_PROP_TAG(TolerancePressureMsWells);
NEW_PROP_TAG(MaxPressureChangeMsWells);
NEW_PROP_TAG(UseInnerIterationsMsWells);
NEW_PROP_TAG(MaxInnerIterMsWells);
SET_SCALAR_PROP(FlowModelParameters, DpMaxRel, 0.3);
SET_SCALAR_PROP(FlowModelParameters, DsMax, 0.2);
SET_SCALAR_PROP(FlowModelParameters, DrMaxRel, 1e9);
SET_SCALAR_PROP(FlowModelParameters, DbhpMaxRel, 1.0);
SET_SCALAR_PROP(FlowModelParameters, DwellFractionMax, 0.2);
SET_SCALAR_PROP(FlowModelParameters, MaxResidualAllowed, 1e7);
SET_SCALAR_PROP(FlowModelParameters, ToleranceMb, 1e-5);
SET_SCALAR_PROP(FlowModelParameters, ToleranceCnv,1e-2);
SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvRelaxed, 1e9);
SET_SCALAR_PROP(FlowModelParameters, ToleranceWells, 1e-4);
SET_SCALAR_PROP(FlowModelParameters, ToleranceWellControl, 1e-7);
SET_INT_PROP(FlowModelParameters, MaxWelleqIter, 15);
SET_BOOL_PROP(FlowModelParameters, UseMultisegmentWell, false);
SET_SCALAR_PROP(FlowModelParameters, MaxSinglePrecisionDays, 20.0);
SET_INT_PROP(FlowModelParameters, MaxStrictIter, 8);
SET_BOOL_PROP(FlowModelParameters, SolveWelleqInitially, true);
SET_BOOL_PROP(FlowModelParameters, UpdateEquationsScaling, false);
SET_BOOL_PROP(FlowModelParameters, UseUpdateStabilization, true);
SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false);
SET_BOOL_PROP(FlowModelParameters, PreconditionerAddWellContributions, false);
SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5);
SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5);
SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true);
SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 10);
END_PROPERTIES
namespace Opm
{
/// Solver parameters for the BlackoilModel.
template <class TypeTag>
struct BlackoilModelParametersEbos
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
/// Max relative change in pressure in single iteration.
double dp_max_rel_;
/// Max absolute change in saturation in single iteration.
double ds_max_;
/// Max relative change in gas-oil or oil-gas ratio in single iteration.
double dr_max_rel_;
/// Max relative change in bhp in single iteration.
double dbhp_max_rel_;
/// Max absolute change in well volume fraction in single iteration.
double dwell_fraction_max_;
/// Absolute max limit for residuals.
double max_residual_allowed_;
/// Relative mass balance tolerance (total mass balance error).
double tolerance_mb_;
/// Local convergence tolerance (max of local saturation errors).
double tolerance_cnv_;
/// Relaxed local convergence tolerance (used when iter >= max_strict_iter_).
double tolerance_cnv_relaxed_;
/// Well convergence tolerance.
double tolerance_wells_;
/// Tolerance for the well control equations
// TODO: it might need to distinguish between rate control and pressure control later
double tolerance_well_control_;
/// Tolerance for the pressure equations for multisegment wells
double tolerance_pressure_ms_wells_;
/// Maximum pressure change over an iteratio for ms wells
double max_pressure_change_ms_wells_;
/// Whether to use inner iterations for ms wells
bool use_inner_iterations_ms_wells_;
/// Maximum inner iteration number for ms wells
int max_inner_iter_ms_wells_;
/// Maximum iteration number of the well equation solution
int max_welleq_iter_;
/// Tolerance for time step in seconds where single precision can be used
/// for solving for the Jacobian
double maxSinglePrecisionTimeStep_;
/// Maximum number of Newton iterations before we give up on the CNV convergence criterion
int max_strict_iter_;
/// Solve well equation initially
bool solve_welleq_initially_;
/// Update scaling factors for mass balance equations
bool update_equations_scaling_;
/// Try to detect oscillation or stagnation.
bool use_update_stabilization_;
/// Whether to use MultisegmentWell to handle multisegment wells
/// it is something temporary before the multisegment well model is considered to be
/// well developed and tested.
/// if it is false, we will handle multisegment wells as standard wells, which will be
/// the default behavoir for the moment. Later, we might set it to be true by default if necessary
bool use_multisegment_well_;
/// The file name of the deck
std::string deck_file_name_;
// Whether to add influences of wells between cells to the matrix and preconditioner matrix
bool matrix_add_well_contributions_;
// Whether to add influences of wells between cells to the preconditioner matrix only
bool preconditioner_add_well_contributions_;
/// Construct from user parameters or defaults.
BlackoilModelParametersEbos()
{
dp_max_rel_ = EWOMS_GET_PARAM(TypeTag, Scalar, DpMaxRel);
ds_max_ = EWOMS_GET_PARAM(TypeTag, Scalar, DsMax);
dr_max_rel_ = EWOMS_GET_PARAM(TypeTag, Scalar, DrMaxRel);
dbhp_max_rel_= EWOMS_GET_PARAM(TypeTag, Scalar, DbhpMaxRel);
dwell_fraction_max_ = EWOMS_GET_PARAM(TypeTag, Scalar, DwellFractionMax);
max_residual_allowed_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxResidualAllowed);
tolerance_mb_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceMb);
tolerance_cnv_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnv);
tolerance_cnv_relaxed_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnvRelaxed);
tolerance_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceWells);
tolerance_well_control_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceWellControl);
max_welleq_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxWelleqIter);
use_multisegment_well_ = EWOMS_GET_PARAM(TypeTag, bool, UseMultisegmentWell);
tolerance_pressure_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, TolerancePressureMsWells);
max_pressure_change_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells);
use_inner_iterations_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseInnerIterationsMsWells);
max_inner_iter_ms_wells_ = EWOMS_GET_PARAM(TypeTag, int, MaxInnerIterMsWells);
maxSinglePrecisionTimeStep_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays) *24*60*60;
max_strict_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxStrictIter);
solve_welleq_initially_ = EWOMS_GET_PARAM(TypeTag, bool, SolveWelleqInitially);
update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling);
use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization);
matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
preconditioner_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, PreconditionerAddWellContributions);
deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
}
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DpMaxRel, "Maximum relative change of pressure in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DsMax, "Maximum absolute change of any saturation in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DrMaxRel, "Maximum relative change of the gas-in-oil or oil-in-gas ratio in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DbhpMaxRel, "Maximum relative change of the bottom-hole pressure in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DwellFractionMax, "Maximum absolute change of a well's volume fraction in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxResidualAllowed, "Absolute maximum tolerated for residuals without cutting the time step size");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceMb, "Tolerated mass balance error relative to total mass present");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnv, "Local convergence tolerance (Maximum of local saturation errors)");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnvRelaxed, "Relaxed local convergence tolerance that applies for iterations after the iterations with the strict tolerance");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceWells, "Well convergence tolerance");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceWellControl, "Tolerance for the well control equations");
EWOMS_REGISTER_PARAM(TypeTag, int, MaxWelleqIter, "Maximum number of iterations to determine solution the well equations");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseMultisegmentWell, "Use the well model for multi-segment wells instead of the one for single-segment wells");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, TolerancePressureMsWells, "Tolerance for the pressure equations for multi-segment wells");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells, "Maximum relative pressure change for a single iteration of the multi-segment well model");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseInnerIterationsMsWells, "Use nested iterations for multi-segment wells");
EWOMS_REGISTER_PARAM(TypeTag, int, MaxInnerIterMsWells, "Maximum number of inner iterations for multi-segment wells");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays, "Maximum time step size where single precision floating point arithmetic can be used solving for the linear systems of equations");
EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIter, "Maximum number of Newton iterations before relaxed tolerances are used for the CNV convergence criterion");
EWOMS_REGISTER_PARAM(TypeTag, bool, SolveWelleqInitially, "Fully solve the well equations before each iteration of the reservoir model");
EWOMS_REGISTER_PARAM(TypeTag, bool, UpdateEquationsScaling, "Update scaling factors for mass balance equations during the run");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method");
EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices");
EWOMS_REGISTER_PARAM(TypeTag, bool, PreconditionerAddWellContributions, "Explicitly specify the influences of wells between cells for the preconditioner matrix only");
}
};
} // namespace Opm
#endif // OPM_BLACKOILMODELPARAMETERS_EBOS_HEADER_INCLUDED

2
opm/autodiff/BlackoilWellModel.hpp
View File

@ -69,7 +69,7 @@ namespace Opm {
public:
// --------- Types ---------
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilModelParametersEbos<TypeTag> ModelParameters;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;

3
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -2,8 +2,6 @@
namespace Opm {
template<typename TypeTag>
BlackoilWellModel<TypeTag>::
BlackoilWellModel(Simulator& ebosSimulator,
@ -33,7 +31,6 @@ namespace Opm {
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::

456
opm/autodiff/FlowMainEbos.hpp
View File

@ -56,13 +56,30 @@
#include <dune/common/parallel/mpihelper.hh>
#endif
BEGIN_PROPERTIES;
NEW_PROP_TAG(OutputMode);
NEW_PROP_TAG(EnableDryRun);
NEW_PROP_TAG(OutputInterval);
NEW_PROP_TAG(UseAmg);
SET_STRING_PROP(EclFlowProblem, OutputMode, "all");
// TODO: enumeration parameters. we use strings for now.
SET_STRING_PROP(EclFlowProblem, EnableDryRun, "auto");
SET_INT_PROP(EclFlowProblem, OutputInterval, 1);
END_PROPERTIES;
namespace Opm
{
// The FlowMain class is the ebos based black-oil simulator.
template <class TypeTag>
class FlowMainEbos
{
enum FileOutputValue{
enum FileOutputMode
{
//! \brief No output to files.
OUTPUT_NONE = 0,
//! \brief Output only to log files, no eclipse output.
@ -84,42 +101,134 @@ namespace Opm
typedef Opm::SimulatorFullyImplicitBlackoilEbos<TypeTag> Simulator;
typedef typename Simulator::ReservoirState ReservoirState;
/// This is the main function of Flow.
/// It runs a complete simulation, with the given grid and
/// simulator classes, based on user command-line input. The
/// content of this function used to be in the main() function of
/// flow.cpp.
// Read the command line parameters. Throws an exception if something goes wrong.
static int setupParameters_(int argc, char** argv)
{
// register the flow specific parameters
EWOMS_REGISTER_PARAM(TypeTag, std::string, OutputMode,
"Specify which messages are going to be printed. Valid values are: none, log, all (default)");
EWOMS_REGISTER_PARAM(TypeTag, std::string, EnableDryRun,
"Specify if the simulation ought to be actually run, or just pretended to be");
EWOMS_REGISTER_PARAM(TypeTag, int, OutputInterval,
"Specify the number of report steps between two consecutive writes of restart data");
Simulator::registerParameters();
typedef typename BlackoilModelEbos<TypeTag>::ISTLSolverType ISTLSolverType;
ISTLSolverType::registerParameters();
// register the parameters inherited from ebos
Ewoms::registerAllParameters_<TypeTag>(/*finalizeRegistration=*/false);
// hide the parameters unused by flow. TODO: this is a pain to maintain
EWOMS_HIDE_PARAM(TypeTag, EnableGravity);
EWOMS_HIDE_PARAM(TypeTag, EnableGridAdaptation);
// this parameter is actually used in eWoms, but the flow well model
// hard-codes the assumption that the intensive quantities cache is enabled,
// so flow crashes. Let's hide the parameter for that reason.
EWOMS_HIDE_PARAM(TypeTag, EnableIntensiveQuantityCache);
// thermodynamic hints are not implemented/required by the eWoms blackoil
// model
EWOMS_HIDE_PARAM(TypeTag, EnableThermodynamicHints);
// in flow only the deck file determines the end time of the simulation
EWOMS_HIDE_PARAM(TypeTag, EndTime);
// time stepping is not (yet) done by the eWoms code in flow
EWOMS_HIDE_PARAM(TypeTag, InitialTimeStepSize);
EWOMS_HIDE_PARAM(TypeTag, MaxTimeStepDivisions);
EWOMS_HIDE_PARAM(TypeTag, MaxTimeStepSize);
EWOMS_HIDE_PARAM(TypeTag, MinTimeStepSize);
EWOMS_HIDE_PARAM(TypeTag, PredeterminedTimeStepsFile);
// flow currently uses its own linear solver
EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxError);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxIterations);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverOverlapSize);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverTolerance);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverVerbosity);
EWOMS_HIDE_PARAM(TypeTag, PreconditionerRelaxation);
// flow also does not use the eWoms Newton method
EWOMS_HIDE_PARAM(TypeTag, NewtonMaxError);
EWOMS_HIDE_PARAM(TypeTag, NewtonMaxIterations);
EWOMS_HIDE_PARAM(TypeTag, NewtonRawTolerance);
EWOMS_HIDE_PARAM(TypeTag, NewtonTargetIterations);
EWOMS_HIDE_PARAM(TypeTag, NewtonVerbose);
EWOMS_HIDE_PARAM(TypeTag, NewtonWriteConvergence);
// the default eWoms checkpoint/restart mechanism does not work with flow
EWOMS_HIDE_PARAM(TypeTag, RestartTime);
EWOMS_HIDE_PARAM(TypeTag, RestartWritingInterval);
EWOMS_END_PARAM_REGISTRATION(TypeTag);
// read in the command line parameters
int status = Ewoms::setupParameters_<TypeTag>(argc, const_cast<const char**>(argv), /*doRegistration=*/false);
if (status == 0) {
// deal with --print-properties and --print-parameters
bool doExit = false;
int mpiRank = 0;
#if HAVE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
#endif
if (EWOMS_GET_PARAM(TypeTag, int, PrintProperties) == 1) {
doExit = true;
if (mpiRank == 0)
Ewoms::Properties::printValues<TypeTag>();
}
if (EWOMS_GET_PARAM(TypeTag, int, PrintParameters) == 1) {
doExit = true;
if (mpiRank == 0)
Ewoms::Parameters::printValues<TypeTag>();
}
if (doExit)
return -1;
}
return status;
}
/// This is the main function of Flow. It runs a complete simulation with the
/// given grid and simulator classes, based on the user-specified command-line
/// input.
int execute(int argc, char** argv)
{
try {
setupParallelism();
printStartupMessage();
const bool ok = setupParameters(argc, argv);
if (!ok) {
return EXIT_FAILURE;
}
// deal with some administrative boilerplate
setupEbosSimulator();
int status = setupParameters_(argc, argv);
if (status)
return status;
setupParallelism();
setupOutput();
printStartupMessage();
setupEbosSimulator();
setupLogging();
printPRTHeader();
runDiagnostics();
setupLinearSolver();
createSimulator();
// Run.
auto ret = runSimulator();
// do the actual work
runSimulator();
// clean up
mergeParallelLogFiles();
return ret;
return EXIT_SUCCESS;
}
catch (const std::exception &e) {
catch (const std::exception& e) {
std::ostringstream message;
message << "Program threw an exception: " << e.what();
if( output_cout_ )
{
if (output_cout_) {
// in some cases exceptions are thrown before the logging system is set
// up.
if (OpmLog::hasBackend("STREAMLOG")) {
@ -138,34 +247,14 @@ namespace Opm
void setupParallelism()
{
// determine the rank of the current process and the number of processes
// involved in the simulation. MPI must have already been initialized here.
// involved in the simulation. MPI must have already been initialized
// here. (yes, the name of this method is misleading.)
#if HAVE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank_);
int mpi_size;
MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
MPI_Comm_size(MPI_COMM_WORLD, &mpi_size_);
#else
mpi_rank_ = 0;
const int mpi_size = 1;
#endif
output_cout_ = ( mpi_rank_ == 0 );
must_distribute_ = ( mpi_size > 1 );
#ifdef _OPENMP
// OpenMP setup.
if (!getenv("OMP_NUM_THREADS")) {
// Default to at most 2 threads, regardless of
// number of cores (unless ENV(OMP_NUM_THREADS) is defined)
int num_cores = omp_get_num_procs();
int num_threads = std::min(2, num_cores);
omp_set_num_threads(num_threads);
}
// omp_get_num_threads() only works as expected within a parallel region.
const int num_omp_threads = omp_get_max_threads();
if (mpi_size == 1) {
std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl;
} else {
std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl;
}
mpi_size_ = 1;
#endif
}
@ -191,71 +280,42 @@ namespace Opm
}
}
// Read parameters, see if a deck was specified on the command line, and if
// it was, insert it into parameters.
// Writes to:
// param_
// Returns true if ok, false if not.
bool setupParameters(int argc, char** argv)
{
param_ = ParameterGroup(argc, argv, false, output_cout_);
// See if a deck was specified on the command line.
if (!param_.unhandledArguments().empty()) {
if (param_.unhandledArguments().size() != 1) {
std::cerr << "You can only specify a single input deck on the command line.\n";
return false;
} else {
const auto casename = this->simulationCaseName( param_.unhandledArguments()[ 0 ] );
param_.insertParameter("deck_filename", casename.string() );
}
}
// We must have an input deck. Grid and props will be read from that.
if (!param_.has("deck_filename")) {
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return false;
}
return true;
}
// Set output_to_files_ and set/create output dir. Write parameter file.
// Extract the minimum priority and determines if log files ought to be created.
// Writes to:
// output_to_files_
// output_dir_
// Throws std::runtime_error if failed to create (if requested) output dir.
// output_
void setupOutput()
{
const std::string output = param_.getDefault("output", std::string("all"));
static std::map<std::string, FileOutputValue> string2OutputEnum =
const std::string outputModeString =
EWOMS_GET_PARAM(TypeTag, std::string, OutputMode);
static std::map<std::string, FileOutputMode> stringToOutputMode =
{ {"none", OUTPUT_NONE },
{"false", OUTPUT_LOG_ONLY },
{"log", OUTPUT_LOG_ONLY },
{"all" , OUTPUT_ALL },
{"true" , OUTPUT_ALL }};
auto converted = string2OutputEnum.find(output);
if ( converted != string2OutputEnum.end() )
{
output_ = string2OutputEnum[output];
auto outputModeIt = stringToOutputMode.find(outputModeString);
if (outputModeIt != stringToOutputMode.end()) {
output_ = outputModeIt->second;
}
else
{
std::cerr << "Value " << output <<
" passed to option output was invalid. Using \"all\" instead."
else {
output_ = OUTPUT_ALL;
std::cerr << "Value " << outputModeString <<
" is not a recognized output mode. Using \"all\" instead."
<< std::endl;
}
output_to_files_ = output_cout_ && (output_ != OUTPUT_NONE);
output_cout_ = false;
if (mpi_rank_ == 0) {
output_cout_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
output_to_files_ = (output_ != OUTPUT_NONE);
}
}
// Setup OpmLog backend with output_dir.
// Setup the OpmLog backends
void setupLogging()
{
std::string deck_filename = param_.get<std::string>("deck_filename");
std::string deck_filename = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
// create logFile
using boost::filesystem::path;
path fpath(deck_filename);
@ -273,8 +333,7 @@ namespace Opm
logFileStream << output_dir << "/" << baseName;
debugFileStream << output_dir << "/" << baseName;
if ( must_distribute_ && mpi_rank_ != 0 )
{
if (mpi_rank_ != 0) {
// Added rank to log file for non-zero ranks.
// This prevents message loss.
debugFileStream << "."<< mpi_rank_;
@ -286,19 +345,17 @@ namespace Opm
logFile_ = logFileStream.str();
if( output_ > OUTPUT_NONE)
{
if (output_ > OUTPUT_NONE) {
std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::NoDebugMessageTypes, false, output_cout_);
OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
prtLog->setMessageLimiter(std::make_shared<MessageLimiter>());
prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
}
if( output_ >= OUTPUT_LOG_ONLY && !param_.getDefault("no_debug_log", false) )
{
if (output_ >= OUTPUT_LOG_ONLY) {
std::string debugFile = debugFileStream.str();
std::shared_ptr<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
OpmLog::addBackend( "DEBUGLOG" , debugLog);
OpmLog::addBackend("DEBUGLOG", debugLog);
}
std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(std::cout, Log::StdoutMessageTypes);
@ -312,18 +369,12 @@ namespace Opm
{Log::MessageType::Bug, msgLimits.getBugPrintLimit(0)}};
streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true));
if ( output_cout_ )
{
// Read Parameters.
OpmLog::debug("\n--------------- Reading parameters ---------------\n");
}
}
// Print an ASCII-art header to the PRT and DEBUG files.
void printPRTHeader()
{
// Print header for PRT file.
if ( output_cout_ ) {
if (output_cout_) {
const std::string version = moduleVersionName();
const double megabyte = 1024 * 1024;
unsigned num_cpu = std::thread::hardware_concurrency();
@ -346,7 +397,7 @@ namespace Opm
ss << " and is part of OPM.\nFor more information visit: https://opm-project.org \n\n";
ss << "Flow Version = " + version + "\n";
if (uname(&arch) == 0) {
ss << "System = " << arch.nodename << " (Number of cores: " << num_cpu;
ss << "System = " << arch.nodename << " (Number of logical cores: " << num_cpu;
ss << ", RAM: " << std::fixed << std::setprecision (2) << mem_size << " MB) \n";
ss << "Architecture = " << arch.sysname << " " << arch.machine << " (Release: " << arch.release;
ss << ", Version: " << arch.version << " )\n";
@ -355,6 +406,10 @@ namespace Opm
ss << "User = " << user << std::endl;
}
ss << "Simulation started on " << tmstr << " hrs\n";
ss << "Parameters used by Flow:\n";
Ewoms::Parameters::printValues<TypeTag>(ss);
OpmLog::note(ss.str());
}
}
@ -364,20 +419,15 @@ namespace Opm
// force closing of all log files.
OpmLog::removeAllBackends();
if( mpi_rank_ != 0 || !must_distribute_ || !output_to_files_ )
{
if (mpi_rank_ != 0 || mpi_size_ < 2 || !output_to_files_) {
return;
}
namespace fs = boost::filesystem;
fs::path output_path(".");
const std::string& output_dir = eclState().getIOConfig().getOutputDir();
if ( param_.has("output_dir") )
{
output_path = fs::path(output_dir);
}
fs::path deck_filename(param_.get<std::string>("deck_filename"));
fs::path deck_filename(EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName));
std::for_each(fs::directory_iterator(output_path),
fs::directory_iterator(),
detail::ParallelFileMerger(output_path, deck_filename.stem().string()));
@ -385,47 +435,6 @@ namespace Opm
void setupEbosSimulator()
{
std::vector<const char*> argv;
argv.push_back("flow_ebos");
std::string deckFileParam("--ecl-deck-file-name=");
const std::string& deckFileName = param_.get<std::string>("deck_filename");
deckFileParam += deckFileName;
argv.push_back(deckFileParam.c_str());
std::string outputDirParam("--output-dir=");
if (param_.has("output_dir")) {
const std::string& output_dir = param_.get<std::string>("output_dir");
outputDirParam += output_dir;
argv.push_back(outputDirParam.c_str());
}
std::string asyncOutputParam("--enable-async-ecl-output=");
if (param_.has("async_output")) {
const std::string& value = param_.get<std::string>("async_output");
asyncOutputParam += value;
argv.push_back(asyncOutputParam.c_str());
}
std::string outputDoublePrecisionParam("--ecl-output-double-precision=");
if (param_.has("restart_double_si")) {
const std::string& value = param_.get<std::string>("restart_double_si");
outputDoublePrecisionParam += value;
argv.push_back(outputDoublePrecisionParam.c_str());
}
#if defined(_OPENMP)
std::string numThreadsParam("--threads-per-process=");
int numThreads = omp_get_max_threads();
numThreadsParam += std::to_string(numThreads);
argv.push_back(numThreadsParam.c_str());
#endif // defined(_OPENMP)
EbosSimulator::registerParameters();
Ewoms::setupParameters_<TypeTag>(argv.size(), &argv[0]);
EbosThreadManager::init();
ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false));
ebosSimulator_->model().applyInitialSolution();
@ -435,10 +444,22 @@ namespace Opm
}
// Possible to force initialization only behavior (NOSIM).
if (param_.has("nosim")) {
const bool nosim = param_.get<bool>("nosim");
const std::string& dryRunString = EWOMS_GET_PARAM(TypeTag, std::string, EnableDryRun);
if (dryRunString != "" && dryRunString != "auto") {
bool yesno;
if (dryRunString == "true"
|| dryRunString == "t"
|| dryRunString == "1")
yesno = true;
else if (dryRunString == "false"
|| dryRunString == "f"
|| dryRunString == "0")
yesno = false;
else
throw std::invalid_argument("Invalid value for parameter EnableDryRun: '"
+dryRunString+"'");
auto& ioConfig = eclState().getIOConfig();
ioConfig.overrideNOSIM( nosim );
ioConfig.overrideNOSIM(yesno);
}
}
catch (const std::invalid_argument& e) {
@ -446,12 +467,6 @@ namespace Opm
std::cerr << "Exception caught: " << e.what() << std::endl;
throw;
}
// Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
if (param_.has("output_interval")) {
const int output_interval = param_.get<int>("output_interval");
eclState().getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) );
}
}
const Deck& deck() const
@ -475,8 +490,7 @@ namespace Opm
// OpmLog singleton.
void runDiagnostics()
{
if( ! output_cout_ )
{
if (!output_cout_) {
return;
}
@ -486,8 +500,7 @@ namespace Opm
}
// Run the simulator.
// Returns EXIT_SUCCESS if it does not throw.
int runSimulator()
void runSimulator()
{
const auto& schedule = this->schedule();
const auto& timeMap = schedule.getTimeMap();
@ -498,6 +511,18 @@ namespace Opm
const auto& initConfig = eclState().getInitConfig();
simtimer.init(timeMap, (size_t)initConfig.getRestartStep());
if (output_cout_) {
std::ostringstream oss;
// This allows a user to catch typos and misunderstandings in the
// use of simulator parameters.
if (Ewoms::Parameters::printUnused<TypeTag>(oss)) {
std::cout << "----------------- Unrecognized parameters: -----------------\n";
std::cout << oss.str();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
if (!ioConfig.initOnly()) {
if (output_cout_) {
std::string msg;
@ -513,13 +538,6 @@ namespace Opm
ss << "\n\n================ End of simulation ===============\n\n";
successReport.reportFullyImplicit(ss, &failureReport);
OpmLog::info(ss.str());
if (param_.anyUnused()) {
// This allows a user to catch typos and misunderstandings in the
// use of simulator parameters.
std::cout << "-------------------- Unused parameters: --------------------\n";
param_.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} else {
@ -528,39 +546,35 @@ namespace Opm
}
}
return EXIT_SUCCESS;
}
// Setup linear solver.
// Writes to:
// fis_solver_
// linearSolver_
void setupLinearSolver()
{
typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
const std::string cprSolver = "cpr";
if (!param_.has("solver_approach") )
{
if ( eclState().getSimulationConfig().useCPR() )
{
/* Deactivate selection of CPR via eclipse keyword
as this preconditioner is still considered experimental
and fails miserably for some models.
param_.insertParameter("solver_approach", cprSolver);
*/
if ( output_cout_ )
{
std::ostringstream message;
message << "Ignoring request for CPRPreconditioner "
<< "via Eclipse keyword as it is considered "
<<" experimental. To activate use "
<<"\"solver_approach=cprSolver\" command "
<<"line parameter.";
OpmLog::info(message.str());
}
}
}
typedef typename BlackoilModelEbos<TypeTag>::ISTLSolverType ISTLSolverType;
extractParallelGridInformationToISTL(grid(), parallel_information_);
fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );
auto *tmp = new ISTLSolverType(parallel_information_);
linearSolver_.reset(tmp);
// Deactivate selection of CPR via eclipse keyword
// as this preconditioner is still considered experimental
// and fails miserably for some models.
if (output_cout_
&& eclState().getSimulationConfig().useCPR()
&& !tmp->parameters().use_cpr_)
{
std::ostringstream message;
message << "Ignoring request for CPRPreconditioner "
<< "via Eclipse keyword as it is considered "
<<" experimental. To activate use "
<<"\"--flow-use-cpr=true\" command "
<<"line parameter.";
OpmLog::info(message.str());
}
}
/// This is the main function of Flow.
@ -570,37 +584,7 @@ namespace Opm
void createSimulator()
{
// Create the simulator instance.
simulator_.reset(new Simulator(*ebosSimulator_,
param_,
*fis_solver_));
}
private:
boost::filesystem::path simulationCaseName( const std::string& casename ) {
namespace fs = boost::filesystem;
const auto exists = []( const fs::path& f ) -> bool {
if( !fs::exists( f ) ) return false;
if( fs::is_regular_file( f ) ) return true;
return fs::is_symlink( f )
&& fs::is_regular_file( fs::read_symlink( f ) );
};
auto simcase = fs::path( casename );
if( exists( simcase ) ) {
return simcase;
}
for( const auto& ext : { std::string("data"), std::string("DATA") } ) {
if( exists( simcase.replace_extension( ext ) ) ) {
return simcase;
}
}
throw std::invalid_argument( "Cannot find input case " + casename );
simulator_.reset(new Simulator(*ebosSimulator_, *linearSolver_));
}
unsigned long long getTotalSystemMemory()
@ -614,15 +598,15 @@ namespace Opm
Grid& grid()
{ return ebosSimulator_->vanguard().grid(); }
private:
std::unique_ptr<EbosSimulator> ebosSimulator_;
int mpi_rank_ = 0;
int mpi_size_ = 1;
bool output_cout_ = false;
FileOutputValue output_ = OUTPUT_ALL;
bool must_distribute_ = false;
ParameterGroup param_;
FileOutputMode output_ = OUTPUT_ALL;
bool output_to_files_ = false;
boost::any parallel_information_;
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
std::unique_ptr<NewtonIterationBlackoilInterface> linearSolver_;
std::unique_ptr<Simulator> simulator_;
std::string logFile_;
};

290
opm/autodiff/ISTLSolver.hpp
View File

@ -20,6 +20,7 @@
#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
#define OPM_ISTLSOLVER_HEADER_INCLUDED
#include <opm/autodiff/ISTLSolverEbos.hpp>
#include <opm/autodiff/BlackoilAmg.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
@ -41,297 +42,8 @@
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
namespace Dune
{
namespace FMatrixHelp {
//! invert 4x4 Matrix without changing the original matrix
template <typename K>
static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
{
inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
matrix[1][1] * matrix[2][3] * matrix[3][2] -
matrix[2][1] * matrix[1][2] * matrix[3][3] +
matrix[2][1] * matrix[1][3] * matrix[3][2] +
matrix[3][1] * matrix[1][2] * matrix[2][3] -
matrix[3][1] * matrix[1][3] * matrix[2][2];
inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
matrix[1][0] * matrix[2][3] * matrix[3][2] +
matrix[2][0] * matrix[1][2] * matrix[3][3] -
matrix[2][0] * matrix[1][3] * matrix[3][2] -
matrix[3][0] * matrix[1][2] * matrix[2][3] +
matrix[3][0] * matrix[1][3] * matrix[2][2];
inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
matrix[1][0] * matrix[2][3] * matrix[3][1] -
matrix[2][0] * matrix[1][1] * matrix[3][3] +
matrix[2][0] * matrix[1][3] * matrix[3][1] +
matrix[3][0] * matrix[1][1] * matrix[2][3] -
matrix[3][0] * matrix[1][3] * matrix[2][1];
inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
matrix[1][0] * matrix[2][2] * matrix[3][1] +
matrix[2][0] * matrix[1][1] * matrix[3][2] -
matrix[2][0] * matrix[1][2] * matrix[3][1] -
matrix[3][0] * matrix[1][1] * matrix[2][2] +
matrix[3][0] * matrix[1][2] * matrix[2][1];
inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
matrix[0][1] * matrix[2][3] * matrix[3][2] +
matrix[2][1] * matrix[0][2] * matrix[3][3] -
matrix[2][1] * matrix[0][3] * matrix[3][2] -
matrix[3][1] * matrix[0][2] * matrix[2][3] +
matrix[3][1] * matrix[0][3] * matrix[2][2];
inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
matrix[0][0] * matrix[2][3] * matrix[3][2] -
matrix[2][0] * matrix[0][2] * matrix[3][3] +
matrix[2][0] * matrix[0][3] * matrix[3][2] +
matrix[3][0] * matrix[0][2] * matrix[2][3] -
matrix[3][0] * matrix[0][3] * matrix[2][2];
inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
matrix[0][0] * matrix[2][3] * matrix[3][1] +
matrix[2][0] * matrix[0][1] * matrix[3][3] -
matrix[2][0] * matrix[0][3] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[2][3] +
matrix[3][0] * matrix[0][3] * matrix[2][1];
inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
matrix[0][0] * matrix[2][2] * matrix[3][1] -
matrix[2][0] * matrix[0][1] * matrix[3][2] +
matrix[2][0] * matrix[0][2] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[2][2] -
matrix[3][0] * matrix[0][2] * matrix[2][1];
inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
matrix[0][1] * matrix[1][3] * matrix[3][2] -
matrix[1][1] * matrix[0][2] * matrix[3][3] +
matrix[1][1] * matrix[0][3] * matrix[3][2] +
matrix[3][1] * matrix[0][2] * matrix[1][3] -
matrix[3][1] * matrix[0][3] * matrix[1][2];
inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
matrix[0][0] * matrix[1][3] * matrix[3][2] +
matrix[1][0] * matrix[0][2] * matrix[3][3] -
matrix[1][0] * matrix[0][3] * matrix[3][2] -
matrix[3][0] * matrix[0][2] * matrix[1][3] +
matrix[3][0] * matrix[0][3] * matrix[1][2];
inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
matrix[0][0] * matrix[1][3] * matrix[3][1] -
matrix[1][0] * matrix[0][1] * matrix[3][3] +
matrix[1][0] * matrix[0][3] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[1][3] -
matrix[3][0] * matrix[0][3] * matrix[1][1];
inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
matrix[0][0] * matrix[1][2] * matrix[3][1] +
matrix[1][0] * matrix[0][1] * matrix[3][2] -
matrix[1][0] * matrix[0][2] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[1][2] +
matrix[3][0] * matrix[0][2] * matrix[1][1];
inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
matrix[0][1] * matrix[1][3] * matrix[2][2] +
matrix[1][1] * matrix[0][2] * matrix[2][3] -
matrix[1][1] * matrix[0][3] * matrix[2][2] -
matrix[2][1] * matrix[0][2] * matrix[1][3] +
matrix[2][1] * matrix[0][3] * matrix[1][2];
inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
matrix[0][0] * matrix[1][3] * matrix[2][2] -
matrix[1][0] * matrix[0][2] * matrix[2][3] +
matrix[1][0] * matrix[0][3] * matrix[2][2] +
matrix[2][0] * matrix[0][2] * matrix[1][3] -
matrix[2][0] * matrix[0][3] * matrix[1][2];
inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
matrix[0][0] * matrix[1][3] * matrix[2][1] +
matrix[1][0] * matrix[0][1] * matrix[2][3] -
matrix[1][0] * matrix[0][3] * matrix[2][1] -
matrix[2][0] * matrix[0][1] * matrix[1][3] +
matrix[2][0] * matrix[0][3] * matrix[1][1];
inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
matrix[0][0] * matrix[1][2] * matrix[2][1] -
matrix[1][0] * matrix[0][1] * matrix[2][2] +
matrix[1][0] * matrix[0][2] * matrix[2][1] +
matrix[2][0] * matrix[0][1] * matrix[1][2] -
matrix[2][0] * matrix[0][2] * matrix[1][1];
K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
// return identity for singular or nearly singular matrices.
if (std::abs(det) < 1e-40) {
for (int i = 0; i < 4; ++i){
inverse[i][i] = 1.0;
}
return 1.0;
}
K inv_det = 1.0 / det;
inverse *= inv_det;
return det;
}
} // end FMatrixHelp
namespace ISTLUtility {
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
{
FieldMatrix<K,1,1> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
{
FieldMatrix<K,2,2> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
{
FieldMatrix<K,3,3> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
{
FieldMatrix<K,4,4> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling matrix.invert
template <typename K, int n>
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
{
Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
matrix.invert();
}
} // end ISTLUtility
template <class Scalar, int n, int m>
class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
{
public:
typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
using BaseType :: operator= ;
using BaseType :: rows;
using BaseType :: cols;
explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
void invert()
{
ISTLUtility::invertMatrix( *this );
}
const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
BaseType& asBase() { return static_cast< BaseType& > (*this); }
};
template<class K, int n, int m>
void
print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
typename FieldMatrix<K,n,m>::size_type I,
typename FieldMatrix<K,n,m>::size_type J,
typename FieldMatrix<K,n,m>::size_type therow, int width,
int precision)
{
print_row(s, A.asBase(), I, J, therow, width, precision);
}
template<class K, int n, int m>
K& firstmatrixelement (MatrixBlock<K,n,m>& A)
{
return firstmatrixelement( A.asBase() );
}
template<typename Scalar, int n, int m>
struct MatrixDimension< MatrixBlock< Scalar, n, m > >
: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
{
};
#if HAVE_UMFPACK
/// \brief UMFPack specialization for MatrixBlock to make AMG happy
///
/// Without this the empty default implementation would be used.
template<typename T, typename A, int n, int m>
class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
: public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
{
typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
public:
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
UMFPack(const RealMatrix& matrix, int verbose, bool)
: Base(reinterpret_cast<const Matrix&>(matrix), verbose)
{}
};
#endif
#if HAVE_SUPERLU
/// \brief SuperLU specialization for MatrixBlock to make AMG happy
///
/// Without this the empty default implementation would be used.
template<typename T, typename A, int n, int m>
class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
: public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
{
typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
public:
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
: Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
{}
};
#endif
} // end namespace Dune
namespace Opm
{
namespace Detail
{
//! calculates ret = A^T * B
template< class K, int m, int n, int p >
static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, n, m > &A,
const Dune::FieldMatrix< K, n, p > &B,
Dune::FieldMatrix< K, m, p > &ret )
{
typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
for( size_type i = 0; i < m; ++i )
{
for( size_type j = 0; j < p; ++j )
{
ret[ i ][ j ] = K( 0 );
for( size_type k = 0; k < n; ++k )
ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
}
}
}
}
/// This class solves the fully implicit black-oil system by
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed

705
opm/autodiff/ISTLSolverEbos.hpp Normal file
View File

@ -0,0 +1,705 @@
/*
Copyright 2016 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
#define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
#include <opm/autodiff/BlackoilAmg.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <ewoms/common/parametersystem.hh>
#include <ewoms/common/propertysystem.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowIstlSolver, INHERITS_FROM(FlowIstlSolverParams));
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(GlobalEqVector);
NEW_PROP_TAG(JacobianMatrix);
NEW_PROP_TAG(Indices);
END_PROPERTIES
namespace Dune
{
namespace FMatrixHelp {
//! invert 4x4 Matrix without changing the original matrix
template <typename K>
static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
{
inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
matrix[1][1] * matrix[2][3] * matrix[3][2] -
matrix[2][1] * matrix[1][2] * matrix[3][3] +
matrix[2][1] * matrix[1][3] * matrix[3][2] +
matrix[3][1] * matrix[1][2] * matrix[2][3] -
matrix[3][1] * matrix[1][3] * matrix[2][2];
inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
matrix[1][0] * matrix[2][3] * matrix[3][2] +
matrix[2][0] * matrix[1][2] * matrix[3][3] -
matrix[2][0] * matrix[1][3] * matrix[3][2] -
matrix[3][0] * matrix[1][2] * matrix[2][3] +
matrix[3][0] * matrix[1][3] * matrix[2][2];
inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
matrix[1][0] * matrix[2][3] * matrix[3][1] -
matrix[2][0] * matrix[1][1] * matrix[3][3] +
matrix[2][0] * matrix[1][3] * matrix[3][1] +
matrix[3][0] * matrix[1][1] * matrix[2][3] -
matrix[3][0] * matrix[1][3] * matrix[2][1];
inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
matrix[1][0] * matrix[2][2] * matrix[3][1] +
matrix[2][0] * matrix[1][1] * matrix[3][2] -
matrix[2][0] * matrix[1][2] * matrix[3][1] -
matrix[3][0] * matrix[1][1] * matrix[2][2] +
matrix[3][0] * matrix[1][2] * matrix[2][1];
inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
matrix[0][1] * matrix[2][3] * matrix[3][2] +
matrix[2][1] * matrix[0][2] * matrix[3][3] -
matrix[2][1] * matrix[0][3] * matrix[3][2] -
matrix[3][1] * matrix[0][2] * matrix[2][3] +
matrix[3][1] * matrix[0][3] * matrix[2][2];
inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
matrix[0][0] * matrix[2][3] * matrix[3][2] -
matrix[2][0] * matrix[0][2] * matrix[3][3] +
matrix[2][0] * matrix[0][3] * matrix[3][2] +
matrix[3][0] * matrix[0][2] * matrix[2][3] -
matrix[3][0] * matrix[0][3] * matrix[2][2];
inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
matrix[0][0] * matrix[2][3] * matrix[3][1] +
matrix[2][0] * matrix[0][1] * matrix[3][3] -
matrix[2][0] * matrix[0][3] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[2][3] +
matrix[3][0] * matrix[0][3] * matrix[2][1];
inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
matrix[0][0] * matrix[2][2] * matrix[3][1] -
matrix[2][0] * matrix[0][1] * matrix[3][2] +
matrix[2][0] * matrix[0][2] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[2][2] -
matrix[3][0] * matrix[0][2] * matrix[2][1];
inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
matrix[0][1] * matrix[1][3] * matrix[3][2] -
matrix[1][1] * matrix[0][2] * matrix[3][3] +
matrix[1][1] * matrix[0][3] * matrix[3][2] +
matrix[3][1] * matrix[0][2] * matrix[1][3] -
matrix[3][1] * matrix[0][3] * matrix[1][2];
inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
matrix[0][0] * matrix[1][3] * matrix[3][2] +
matrix[1][0] * matrix[0][2] * matrix[3][3] -
matrix[1][0] * matrix[0][3] * matrix[3][2] -
matrix[3][0] * matrix[0][2] * matrix[1][3] +
matrix[3][0] * matrix[0][3] * matrix[1][2];
inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
matrix[0][0] * matrix[1][3] * matrix[3][1] -
matrix[1][0] * matrix[0][1] * matrix[3][3] +
matrix[1][0] * matrix[0][3] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[1][3] -
matrix[3][0] * matrix[0][3] * matrix[1][1];
inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
matrix[0][0] * matrix[1][2] * matrix[3][1] +
matrix[1][0] * matrix[0][1] * matrix[3][2] -
matrix[1][0] * matrix[0][2] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[1][2] +
matrix[3][0] * matrix[0][2] * matrix[1][1];
inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
matrix[0][1] * matrix[1][3] * matrix[2][2] +
matrix[1][1] * matrix[0][2] * matrix[2][3] -
matrix[1][1] * matrix[0][3] * matrix[2][2] -
matrix[2][1] * matrix[0][2] * matrix[1][3] +
matrix[2][1] * matrix[0][3] * matrix[1][2];
inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
matrix[0][0] * matrix[1][3] * matrix[2][2] -
matrix[1][0] * matrix[0][2] * matrix[2][3] +
matrix[1][0] * matrix[0][3] * matrix[2][2] +
matrix[2][0] * matrix[0][2] * matrix[1][3] -
matrix[2][0] * matrix[0][3] * matrix[1][2];
inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
matrix[0][0] * matrix[1][3] * matrix[2][1] +
matrix[1][0] * matrix[0][1] * matrix[2][3] -
matrix[1][0] * matrix[0][3] * matrix[2][1] -
matrix[2][0] * matrix[0][1] * matrix[1][3] +
matrix[2][0] * matrix[0][3] * matrix[1][1];
inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
matrix[0][0] * matrix[1][2] * matrix[2][1] -
matrix[1][0] * matrix[0][1] * matrix[2][2] +
matrix[1][0] * matrix[0][2] * matrix[2][1] +
matrix[2][0] * matrix[0][1] * matrix[1][2] -
matrix[2][0] * matrix[0][2] * matrix[1][1];
K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
// return identity for singular or nearly singular matrices.
if (std::abs(det) < 1e-40) {
for (int i = 0; i < 4; ++i){
inverse[i][i] = 1.0;
}
return 1.0;
}
K inv_det = 1.0 / det;
inverse *= inv_det;
return det;
}
} // end FMatrixHelp
namespace ISTLUtility {
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
{
FieldMatrix<K,1,1> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
{
FieldMatrix<K,2,2> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
{
FieldMatrix<K,3,3> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
{
FieldMatrix<K,4,4> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling matrix.invert
template <typename K, int n>
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
{
Dune::FMatrixPrecision<K>::set_singular_limit(1.e-20);
matrix.invert();
}
} // end ISTLUtility
template <class Scalar, int n, int m>
class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
{
public:
typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
using BaseType :: operator= ;
using BaseType :: rows;
using BaseType :: cols;
explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
void invert()
{
ISTLUtility::invertMatrix( *this );
}
const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
BaseType& asBase() { return static_cast< BaseType& > (*this); }
};
template<class K, int n, int m>
void
print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
typename FieldMatrix<K,n,m>::size_type I,
typename FieldMatrix<K,n,m>::size_type J,
typename FieldMatrix<K,n,m>::size_type therow, int width,
int precision)
{
print_row(s, A.asBase(), I, J, therow, width, precision);
}
template<class K, int n, int m>
K& firstmatrixelement (MatrixBlock<K,n,m>& A)
{
return firstmatrixelement( A.asBase() );
}
template<typename Scalar, int n, int m>
struct MatrixDimension< MatrixBlock< Scalar, n, m > >
: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
{
};
#if HAVE_UMFPACK
/// \brief UMFPack specialization for MatrixBlock to make AMG happy
///
/// Without this the empty default implementation would be used.
template<typename T, typename A, int n, int m>
class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
: public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
{
typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
public:
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
UMFPack(const RealMatrix& matrix, int verbose, bool)
: Base(reinterpret_cast<const Matrix&>(matrix), verbose)
{}
};
#endif
#if HAVE_SUPERLU
/// \brief SuperLU specialization for MatrixBlock to make AMG happy
///
/// Without this the empty default implementation would be used.
template<typename T, typename A, int n, int m>
class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
: public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
{
typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
public:
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
: Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
{}
};
#endif
} // end namespace Dune
namespace Opm
{
namespace Detail
{
//! calculates ret = A^T * B
template< class K, int m, int n, int p >
static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, n, m > &A,
const Dune::FieldMatrix< K, n, p > &B,
Dune::FieldMatrix< K, m, p > &ret )
{
typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
for( size_type i = 0; i < m; ++i )
{
for( size_type j = 0; j < p; ++j )
{
ret[ i ][ j ] = K( 0 );
for( size_type k = 0; k < n; ++k )
ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
}
}
}
}
/// This class solves the fully implicit black-oil system by
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed
/// number of cell variables np .
/// \tparam MatrixBlockType The type of the matrix block used.
/// \tparam VectorBlockType The type of the vector block used.
/// \tparam pressureIndex The index of the pressure component in the vector
/// vector block. It is used to guide the AMG coarsening.
/// Default is zero.
template <class TypeTag>
class ISTLSolverEbos : public NewtonIterationBlackoilInterface
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix;
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { pressureIndex = Indices::pressureSwitchIdx };
public:
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector;
static void registerParameters()
{
NewtonIterationBlackoilInterleavedParameters::registerParameters<TypeTag>();
}
/// Construct a system solver.
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
ISTLSolverEbos(const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 )
, parallelInformation_(parallelInformation_arg)
, isIORank_(isIORank(parallelInformation_arg))
{
parameters_.template init<TypeTag>();
}
const NewtonIterationBlackoilInterleavedParameters& parameters() const
{ return parameters_; }
// dummy method that is not implemented for this class
SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual&) const
{
OPM_THROW(std::logic_error,"This method is not implemented");
return SolutionVector();
}
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] residual residual object containing A and b.
/// \return the solution x
/// \copydoc NewtonIterationBlackoilInterface::iterations
int iterations () const { return iterations_; }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
const boost::any& parallelInformation() const { return parallelInformation_; }
public:
/// \brief construct the CPR preconditioner and the solver.
/// \tparam P The type of the parallel information.
/// \param parallelInformation the information about the parallelization.
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
template<Dune::SolverCategory::Category category=Dune::SolverCategory::sequential,
class LinearOperator, class POrComm>
#else
template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
#endif
void constructPreconditionerAndSolve(LinearOperator& linearOperator,
Vector& x, Vector& istlb,
const POrComm& parallelInformation_arg,
Dune::InverseOperatorResult& result) const
{
// Construct scalar product.
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
#else
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
#endif
// Communicate if parallel.
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
#if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_)
{
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
typedef typename CPRSelectorType::Operator MatrixOperator;
std::unique_ptr< MatrixOperator > opA;
if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
{
// create new operator in case linear operator and matrix operator differ
opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
}
const double relax = parameters_.ilu_relaxation_;
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
if ( parameters_.use_cpr_ )
{
using Matrix = typename MatrixOperator::matrix_type;
using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
using AMG = typename ISTLUtility
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
Criterion crit(15, 2000);
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
}
else
{
typedef typename CPRSelectorType::AMG AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
}
}
else
#endif
{
// Construct preconditioner.
auto precond = constructPrecond(linearOperator, parallelInformation_arg);
// Solve.
solve(linearOperator, x, istlb, *sp, *precond, result);
}
}
// 3x3 matrix block inversion was unstable at least 2.3 until and including
// 2.5.0. There may still be some issue with the 4x4 matrix block inversion
// we therefore still use the block inversion in OPM
typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
Matrix::block_type::rows,
Matrix::block_type::cols> >,
Vector, Vector> SeqPreconditioner;
template <class Operator>
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
{
const double relax = parameters_.ilu_relaxation_;
const int ilu_fillin = parameters_.ilu_fillin_level_;
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
const bool ilu_redblack = parameters_.ilu_redblack_;
const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
return precond;
}
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
#else
typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
Matrix::block_type::rows,
Matrix::block_type::cols> >,
Vector, Vector, Comm> ParPreconditioner;
#endif
template <class Operator>
std::unique_ptr<ParPreconditioner>
constructPrecond(Operator& opA, const Comm& comm) const
{
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = parameters_.ilu_relaxation_;
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
const bool ilu_redblack = parameters_.ilu_redblack_;
const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
}
#endif
template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, comm, amg );
}
template <class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax,
const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax,
milu, comm, amg );
}
template <class C, class LinearOperator, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax,
const MILU_VARIANT milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax,
comm, amg, parameters_ );
}
template <class C, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const MILU_VARIANT milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, milu,
comm, amg, parameters_ );
}
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
{
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.
// GMRes solver
int verbosity = ( isIORank_ ) ? parameters_.linear_solver_verbosity_ : 0;
if ( parameters_.newton_use_gmres_ ) {
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_restart_,
parameters_.linear_solver_maxiter_,
verbosity);
// Solve system.
linsolve.apply(x, istlb, result);
}
else { // BiCGstab solver
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_maxiter_,
verbosity);
// Solve system.
linsolve.apply(x, istlb, result);
}
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
void solve(Matrix& A, Vector& x, Vector& b ) const
{
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
Operator opA(A, istlComm);
solve( opA, x, b, istlComm );
}
else
#endif
{
// Construct operator, scalar product and vectors needed.
Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
solve( opA, x, b );
}
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
template <class Operator, class Comm >
void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
{
Dune::InverseOperatorResult result;
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
const size_t size = opA.getmat().N();
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
// As we use a dune-istl with block size np the number of components
// per parallel is only one.
info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
size, 1);
// Construct operator, scalar product and vectors needed.
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
}
else
#endif
{
OPM_THROW(std::logic_error,"this method if for parallel solve only");
}
checkConvergence( result );
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
template <class Operator>
void solve(Operator& opA, Vector& x, Vector& b ) const
{
Dune::InverseOperatorResult result;
// Construct operator, scalar product and vectors needed.
Dune::Amg::SequentialInformation info;
constructPreconditionerAndSolve(opA, x, b, info, result);
checkConvergence( result );
}
void checkConvergence( const Dune::InverseOperatorResult& result ) const
{
// store number of iterations
iterations_ = result.iterations;
// Check for failure of linear solver.
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
const std::string msg("Convergence failure for linear solver.");
OPM_THROW_NOLOG(LinearSolverProblem, msg);
}
}
protected:
mutable int iterations_;
boost::any parallelInformation_;
bool isIORank_;
NewtonIterationBlackoilInterleavedParameters parameters_;
}; // end ISTLSolver
} // namespace Opm
#endif

78
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@ -30,9 +30,48 @@
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <ewoms/common/parametersystem.hh>
#include <array>
#include <memory>
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowIstlSolverParams);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(LinearSolverReduction);
NEW_PROP_TAG(IluRelaxation);
NEW_PROP_TAG(LinearSolverMaxIter);
NEW_PROP_TAG(LinearSolverRestart);
NEW_PROP_TAG(FlowLinearSolverVerbosity);
NEW_PROP_TAG(IluFillinLevel);
NEW_PROP_TAG(MIluVariant);
NEW_PROP_TAG(IluRedblack);
NEW_PROP_TAG(IluReorderSpheres);
NEW_PROP_TAG(UseGmres);
NEW_PROP_TAG(LinearSolverRequireFullSparsityPattern);
NEW_PROP_TAG(LinearSolverIgnoreConvergenceFailure);
NEW_PROP_TAG(UseAmg);
NEW_PROP_TAG(UseCpr);
SET_SCALAR_PROP(FlowIstlSolverParams, LinearSolverReduction, 1e-2);
SET_SCALAR_PROP(FlowIstlSolverParams, IluRelaxation, 0.9);
SET_INT_PROP(FlowIstlSolverParams, LinearSolverMaxIter, 150);
SET_INT_PROP(FlowIstlSolverParams, LinearSolverRestart, 40);
SET_INT_PROP(FlowIstlSolverParams, FlowLinearSolverVerbosity, 0);
SET_INT_PROP(FlowIstlSolverParams, IluFillinLevel, 0);
SET_STRING_PROP(FlowIstlSolverParams, MIluVariant, "ILU");
SET_BOOL_PROP(FlowIstlSolverParams, IluRedblack, false);
SET_BOOL_PROP(FlowIstlSolverParams, IluReorderSpheres, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseGmres, false);
SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverRequireFullSparsityPattern, false);
SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverIgnoreConvergenceFailure, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false);
END_PROPERTIES
namespace Opm
{
/// This class carries all parameters for the NewtonIterationBlackoilInterleaved class
@ -54,6 +93,45 @@ namespace Opm
bool linear_solver_use_amg_;
bool use_cpr_;
template <class TypeTag>
void init()
{
// TODO: these parameters have undocumented non-trivial dependencies
linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction);
ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation);
linear_solver_maxiter_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter);
linear_solver_restart_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverRestart);
linear_solver_verbosity_ = EWOMS_GET_PARAM(TypeTag, int, FlowLinearSolverVerbosity);
ilu_fillin_level_ = EWOMS_GET_PARAM(TypeTag, int, IluFillinLevel);
ilu_milu_ = convertString2Milu(EWOMS_GET_PARAM(TypeTag, std::string, MIluVariant));
ilu_redblack_ = EWOMS_GET_PARAM(TypeTag, bool, IluRedblack);
ilu_reorder_sphere_ = EWOMS_GET_PARAM(TypeTag, bool, IluReorderSpheres);
newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres);
require_full_sparsity_pattern_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern);
ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure);
linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg);
use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr);
}
template <class TypeTag>
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve");
EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIter, "The maximum number of iterations of the linear solver");
EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverRestart, "The number of iterations after which GMRES is restarted");
EWOMS_REGISTER_PARAM(TypeTag, int, FlowLinearSolverVerbosity, "The verbosity level of the linear solver (0: off, 2: all)");
EWOMS_REGISTER_PARAM(TypeTag, int, IluFillinLevel, "The fill-in level of the linear solver's ILU preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, std::string, MIluVariant, "Specify which variant of the modified-ILU preconditioner ought to be used");
EWOMS_REGISTER_PARAM(TypeTag, bool, IluRedblack, "Use red-black partioning for the ILU preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, bool, IluReorderSpheres, "Reorder the entries of the matrix in the ILU preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver");
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern, "Produce the full sparsity pattern for the linear solver");
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner");
}
NewtonIterationBlackoilInterleavedParameters() { reset(); }
// read values from parameter class
NewtonIterationBlackoilInterleavedParameters( const ParameterGroup& param )

58
opm/autodiff/NonlinearSolverEbos.hpp
View File

@ -26,18 +26,41 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
#include <ewoms/common/parametersystem.hh>
#include <ewoms/common/propertysystem.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <memory>
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowNonLinearSolver);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(NewtonMaxRelax);
NEW_PROP_TAG(FlowNewtonMaxIterations);
NEW_PROP_TAG(FlowNewtonMinIterations);
NEW_PROP_TAG(NewtonRelaxationType);
SET_SCALAR_PROP(FlowNonLinearSolver, NewtonMaxRelax, 0.5);
SET_INT_PROP(FlowNonLinearSolver, FlowNewtonMaxIterations, 10);
SET_INT_PROP(FlowNonLinearSolver, FlowNewtonMinIterations, 1);
SET_STRING_PROP(FlowNonLinearSolver, NewtonRelaxationType, "dampen");
END_PROPERTIES
namespace Opm {
/// A nonlinear solver class suitable for general fully-implicit models,
/// as well as pressure, transport and sequential models.
template <class PhysicalModel>
template <class TypeTag, class PhysicalModel>
class NonlinearSolverEbos
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
// Available relaxation scheme types.
enum RelaxType {
@ -46,7 +69,7 @@ namespace Opm {
};
// Solver parameters controlling nonlinear process.
struct SolverParametersEbos
struct SolverParameters
{
RelaxType relaxType_;
double relaxMax_;
@ -55,28 +78,33 @@ namespace Opm {
int maxIter_; // max nonlinear iterations
int minIter_; // min nonlinear iterations
explicit SolverParametersEbos(const ParameterGroup& param)
SolverParameters()
{
// set default values
reset();
// overload with given parameters
relaxMax_ = param.getDefault("relax_max", relaxMax_);
maxIter_ = param.getDefault("max_iter", maxIter_);
minIter_ = param.getDefault("min_iter", minIter_);
relaxMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxRelax);
maxIter_ = EWOMS_GET_PARAM(TypeTag, int, FlowNewtonMaxIterations);
minIter_ = EWOMS_GET_PARAM(TypeTag, int, FlowNewtonMinIterations);
std::string relaxationType = param.getDefault("relax_type", std::string("dampen"));
if (relaxationType == "dampen") {
const auto& relaxationTypeString = EWOMS_GET_PARAM(TypeTag, std::string, NewtonRelaxationType);
if (relaxationTypeString == "dampen") {
relaxType_ = Dampen;
} else if (relaxationType == "sor") {
} else if (relaxationTypeString == "sor") {
relaxType_ = SOR;
} else {
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationType);
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationTypeString);
}
}
SolverParametersEbos()
{ reset(); }
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonMaxRelax, "The maximum relaxation factor of a Newton iteration used by flow");
EWOMS_REGISTER_PARAM(TypeTag, int, FlowNewtonMaxIterations, "The maximum number of Newton iterations per time step used by flow");
EWOMS_REGISTER_PARAM(TypeTag, int, FlowNewtonMinIterations, "The minimum number of Newton iterations per time step used by flow");
EWOMS_REGISTER_PARAM(TypeTag, std::string, NewtonRelaxationType, "The type of relaxation used by flow's Newton method");
}
void reset()
{
@ -104,7 +132,7 @@ namespace Opm {
///
/// \param[in] param parameters controlling nonlinear process
/// \param[in, out] model physical simulation model.
NonlinearSolverEbos(const SolverParametersEbos& param,
NonlinearSolverEbos(const SolverParameters& param,
std::unique_ptr<PhysicalModel> model)
: param_(param)
, model_(std::move(model))
@ -329,13 +357,13 @@ namespace Opm {
{ return param_.minIter_; }
/// Set parameters to override those given at construction time.
void setParameters(const SolverParametersEbos& param)
void setParameters(const SolverParameters& param)
{ param_ = param; }
private:
// --------- Data members ---------
SimulatorReport failureReport_;
SolverParametersEbos param_;
SolverParameters param_;
std::unique_ptr<PhysicalModel> model_;
int linearizations_;
int nonlinearIterations_;

89
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -25,7 +25,7 @@
#include <opm/autodiff/IterationReport.hpp>
#include <opm/autodiff/NonlinearSolverEbos.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilModelParametersEbos.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp>
@ -36,7 +36,17 @@
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/unused.hh>
BEGIN_PROPERTIES;
NEW_PROP_TAG(EnableTerminalOutput);
NEW_PROP_TAG(EnableAdaptiveTimeStepping);
NEW_PROP_TAG(EnableTuning);
SET_BOOL_PROP(EclFlowProblem, EnableTerminalOutput, true);
SET_BOOL_PROP(EclFlowProblem, EnableAdaptiveTimeStepping, true);
SET_BOOL_PROP(EclFlowProblem, EnableTuning, false);
END_PROPERTIES;
namespace Opm {
@ -56,13 +66,15 @@ public:
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef AdaptiveTimeSteppingEbos<TypeTag> TimeStepper;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilState ReservoirState;
typedef BlackoilModelEbos<TypeTag> Model;
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolverEbos<Model> Solver;
typedef NonlinearSolverEbos<TypeTag, Model> Solver;
typedef typename Model::ModelParameters ModelParameters;
typedef typename Solver::SolverParameters SolverParameters;
typedef BlackoilWellModel<TypeTag> WellModel;
typedef BlackoilAquiferModel<TypeTag> AquiferModel;
@ -85,30 +97,34 @@ public:
///
/// \param[in] props fluid and rock properties
/// \param[in] linsolver linear solver
/// \param[in] has_disgas true for dissolved gas option
/// \param[in] has_vapoil true for vaporized oil option
/// \param[in] eclipse_state the object which represents an internalized ECL deck
/// \param[in] output_writer
/// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow
SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
const ParameterGroup& param,
NewtonIterationBlackoilInterface& linsolver)
NewtonIterationBlackoilInterface& linearSolver)
: ebosSimulator_(ebosSimulator)
, param_(param)
, modelParam_(param)
, solverParam_(param)
, solver_(linsolver)
, phaseUsage_(phaseUsageFromDeck(eclState()))
, terminalOutput_(param.getDefault("output_terminal", true))
, linearSolver_(linearSolver)
{
#if HAVE_MPI
if (solver_.parallelInformation().type() == typeid(ParallelISTLInformation)) {
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
// Only rank 0 does print to std::cout
terminalOutput_ = terminalOutput_ && (info.communicator().rank() == 0);
}
#endif
phaseUsage_ = phaseUsageFromDeck(eclState());
// Only rank 0 does print to std::cout
const auto& comm = grid().comm();
terminalOutput_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
terminalOutput_ = terminalOutput_ && (comm.rank() == 0);
}
static void registerParameters()
{
ModelParameters::registerParameters();
SolverParameters::registerParameters();
TimeStepper::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableTerminalOutput,
"Print high-level information about the simulation's progress to the terminal");
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAdaptiveTimeStepping,
"Use adaptive time stepping between report steps");
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableTuning,
"Honor some aspects of the TUNING keyword.");
}
/// Run the simulation.
@ -139,14 +155,15 @@ public:
// adaptive time stepping
const auto& events = schedule().getEvents();
std::unique_ptr< AdaptiveTimeSteppingEbos > adaptiveTimeStepping;
const bool useTUNING = param_.getDefault("use_TUNING", false);
if (param_.getDefault("timestep.adaptive", true)) {
if (useTUNING) {
adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(schedule().getTuning(), timer.currentStepNum(), param_, terminalOutput_));
std::unique_ptr<TimeStepper > adaptiveTimeStepping;
bool enableAdaptive = EWOMS_GET_PARAM(TypeTag, bool, EnableAdaptiveTimeStepping);
bool enableTUNING = EWOMS_GET_PARAM(TypeTag, bool, EnableTuning);
if (enableAdaptive) {
if (enableTUNING) {
adaptiveTimeStepping.reset(new TimeStepper(schedule().getTuning(), timer.currentStepNum(), terminalOutput_));
}
else {
adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(param_, terminalOutput_));
adaptiveTimeStepping.reset(new TimeStepper(terminalOutput_));
}
double suggestedStepSize = -1.0;
@ -239,7 +256,7 @@ public:
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if (adaptiveTimeStepping) {
if (useTUNING) {
if (enableTUNING) {
if (events.hasEvent(ScheduleEvents::TUNING_CHANGE,timer.currentStepNum())) {
adaptiveTimeStepping->updateTUNING(schedule().getTuning(), timer.currentStepNum());
}
@ -332,7 +349,7 @@ protected:
modelParam_,
wellModel,
aquifer_model,
solver_,
linearSolver_,
terminalOutput_));
return std::unique_ptr<Solver>(new Solver(solverParam_, std::move(model)));
@ -367,23 +384,19 @@ protected:
// Data.
Simulator& ebosSimulator_;
std::unique_ptr<WellConnectionAuxiliaryModule<TypeTag>> wellAuxMod_;
typedef typename Solver::SolverParametersEbos SolverParametersEbos;
SimulatorReport failureReport_;
const ParameterGroup param_;
ModelParameters modelParam_;
SolverParametersEbos solverParam_;
SolverParameters solverParam_;
// Observed objects.
NewtonIterationBlackoilInterface& solver_;
NewtonIterationBlackoilInterface& linearSolver_;
PhaseUsage phaseUsage_;
// Misc. data
bool terminalOutput_;
bool terminalOutput_;
};
} // namespace Opm
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#endif // OPM_SIMULATOR_FULLY_IMPLICIT_BLACKOIL_EBOS_HPP

4
opm/autodiff/WellInterface.hpp
View File

@ -41,7 +41,7 @@
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilModelParametersEbos.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
@ -69,7 +69,7 @@ namespace Opm
using WellState = WellStateFullyImplicitBlackoil;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilModelParametersEbos<TypeTag> ModelParameters;
static const int Water = BlackoilPhases::Aqua;
static const int Oil = BlackoilPhases::Liquid;

159
opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp
View File

@ -6,6 +6,7 @@
#include <iostream>
#include <utility>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
@ -14,11 +15,57 @@
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowTimeSteppingParameters);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(SolverRestartFactor);
NEW_PROP_TAG(SolverGrowthFactor);
NEW_PROP_TAG(SolverMaxGrowth);
NEW_PROP_TAG(SolverMaxTimeStepInDays);
NEW_PROP_TAG(SolverMaxRestarts);
NEW_PROP_TAG(SolverVerbosity);
NEW_PROP_TAG(TimeStepVerbosity);
NEW_PROP_TAG(InitialTimeStepInDays);
NEW_PROP_TAG(FullTimeStepInitially);
NEW_PROP_TAG(TimeStepAfterEventInDays);
NEW_PROP_TAG(TimeStepControl);
NEW_PROP_TAG(TimeStepControlTolerance);
NEW_PROP_TAG(TimeStepControlTargetIterations);
NEW_PROP_TAG(TimeStepControlTargetNewtonIterations);
NEW_PROP_TAG(TimeStepControlDecayRate);
NEW_PROP_TAG(TimeStepControlGrowthRate);
NEW_PROP_TAG(TimeStepControlFileName);
SET_SCALAR_PROP(FlowTimeSteppingParameters, SolverRestartFactor, 0.33);
SET_SCALAR_PROP(FlowTimeSteppingParameters, SolverGrowthFactor, 2.0);
SET_SCALAR_PROP(FlowTimeSteppingParameters, SolverMaxGrowth, 3.0);
SET_SCALAR_PROP(FlowTimeSteppingParameters, SolverMaxTimeStepInDays, 365.0);
SET_INT_PROP(FlowTimeSteppingParameters, SolverMaxRestarts, 10);
SET_INT_PROP(FlowTimeSteppingParameters, SolverVerbosity, 1);
SET_INT_PROP(FlowTimeSteppingParameters, TimeStepVerbosity, 1);
SET_SCALAR_PROP(FlowTimeSteppingParameters, InitialTimeStepInDays, 1.0);
SET_BOOL_PROP(FlowTimeSteppingParameters, FullTimeStepInitially, false);
SET_SCALAR_PROP(FlowTimeSteppingParameters, TimeStepAfterEventInDays, -1.0);
SET_STRING_PROP(FlowTimeSteppingParameters, TimeStepControl, "pid");
SET_SCALAR_PROP(FlowTimeSteppingParameters, TimeStepControlTolerance, 1e-1);
SET_INT_PROP(FlowTimeSteppingParameters, TimeStepControlTargetIterations, 30);
SET_INT_PROP(FlowTimeSteppingParameters, TimeStepControlTargetNewtonIterations, 8);
SET_SCALAR_PROP(FlowTimeSteppingParameters, TimeStepControlDecayRate, 0.75);
SET_SCALAR_PROP(FlowTimeSteppingParameters, TimeStepControlGrowthRate, 1.25);
SET_STRING_PROP(FlowTimeSteppingParameters, TimeStepControlFileName, "timesteps");
END_PROPERTIES
namespace Opm {
// AdaptiveTimeStepping
//---------------------
template<class TypeTag>
class AdaptiveTimeSteppingEbos
{
template <class Solver>
@ -48,23 +95,21 @@ namespace Opm {
public:
//! \brief contructor taking parameter object
AdaptiveTimeSteppingEbos(const ParameterGroup& param,
const bool terminalOutput = true)
AdaptiveTimeSteppingEbos(const bool terminalOutput = true)
: timeStepControl_()
, restartFactor_( param.getDefault("solver.restartfactor", double(0.33)))
, growthFactor_( param.getDefault("solver.growthfactor", double(2)))
, maxGrowth_( param.getDefault("timestep.control.maxgrowth", double(3.0)))
// default is 1 year, convert to seconds
, maxTimeStep_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0), unit::day))
, solverRestartMax_( param.getDefault("solver.restart", int(10)))
, solverVerbose_( param.getDefault("solver.verbose", bool(true)) && terminalOutput)
, timestepVerbose_( param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
, suggestedNextTimestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", 1.0), unit::day))
, fullTimestepInitially_( param.getDefault("full_timestep_initially", bool(false)))
, timestepAfterEvent_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0), unit::day))
, restartFactor_(EWOMS_GET_PARAM(TypeTag, double, SolverRestartFactor)) // 0.33
, growthFactor_(EWOMS_GET_PARAM(TypeTag, double, SolverGrowthFactor)) // 2.0
, maxGrowth_(EWOMS_GET_PARAM(TypeTag, double, SolverMaxGrowth)) // 3.0
, maxTimeStep_(EWOMS_GET_PARAM(TypeTag, double, SolverMaxTimeStepInDays)*24*60*60) // 365.25
, solverRestartMax_(EWOMS_GET_PARAM(TypeTag, int, SolverMaxRestarts)) // 10
, solverVerbose_(EWOMS_GET_PARAM(TypeTag, int, SolverVerbosity) > 0 && terminalOutput) // 2
, timestepVerbose_(EWOMS_GET_PARAM(TypeTag, int, TimeStepVerbosity) > 0 && terminalOutput) // 2
, suggestedNextTimestep_(EWOMS_GET_PARAM(TypeTag, double, InitialTimeStepInDays)*24*60*60) // 1.0
, fullTimestepInitially_(EWOMS_GET_PARAM(TypeTag, bool, FullTimeStepInitially)) // false
, timestepAfterEvent_(EWOMS_GET_PARAM(TypeTag, double, TimeStepAfterEventInDays)*24*60*60) // 1e30
, useNewtonIteration_(false)
{
init_(param);
init_();
}
@ -72,26 +117,62 @@ namespace Opm {
//! \brief contructor taking parameter object
//! \param tuning Pointer to ecl TUNING keyword
//! \param timeStep current report step
//! \param param The parameter object
AdaptiveTimeSteppingEbos(const Tuning& tuning,
size_t timeStep,
const ParameterGroup& param,
const bool terminalOutput = true)
: timeStepControl_()
, restartFactor_(tuning.getTSFCNV(timeStep))
, growthFactor_(tuning.getTFDIFF(timeStep))
, maxGrowth_(tuning.getTSFMAX(timeStep))
// default is 1 year, convert to seconds
, maxTimeStep_(tuning.getTSMAXZ(timeStep))
, solverRestartMax_(param.getDefault("solver.restart", int(10)))
, solverVerbose_(param.getDefault("solver.verbose", bool(true)) && terminalOutput)
, timestepVerbose_(param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
, suggestedNextTimestep_(tuning.getTSINIT(timeStep))
, fullTimestepInitially_(param.getDefault("full_timestep_initially", bool(false)))
, timestepAfterEvent_(tuning.getTMAXWC(timeStep))
, maxTimeStep_(EWOMS_GET_PARAM(TypeTag, double, SolverMaxTimeStepInDays)*24*60*60) // 365.25
, solverRestartMax_(EWOMS_GET_PARAM(TypeTag, int, SolverMaxRestarts)) // 10
, solverVerbose_(EWOMS_GET_PARAM(TypeTag, int, SolverVerbosity) > 0 && terminalOutput) // 2
, timestepVerbose_(EWOMS_GET_PARAM(TypeTag, int, TimeStepVerbosity) > 0 && terminalOutput) // 2
, suggestedNextTimestep_(EWOMS_GET_PARAM(TypeTag, double, InitialTimeStepInDays)*24*60*60) // 1.0
, fullTimestepInitially_(EWOMS_GET_PARAM(TypeTag, bool, FullTimeStepInitially)) // false
, timestepAfterEvent_(EWOMS_GET_PARAM(TypeTag, double, TimeStepAfterEventInDays)*24*60*60) // 1e30
, useNewtonIteration_(false)
{
init_(param);
init_();
}
static void registerParameters()
{
// TODO: make sure the help messages are correct (and useful)
EWOMS_REGISTER_PARAM(TypeTag, double, SolverRestartFactor,
"The factor time steps are elongated after restarts");
EWOMS_REGISTER_PARAM(TypeTag, double, SolverGrowthFactor,
"The factor time steps are elongated after a successful substep");
EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxGrowth,
"The maximum factor time steps are elongated after a report step");
EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxTimeStepInDays,
"The maximum size of a time step in days");
EWOMS_REGISTER_PARAM(TypeTag, int, SolverMaxRestarts,
"The maximum number of breakdowns before a substep is given up and the simulator is terminated");
EWOMS_REGISTER_PARAM(TypeTag, int, SolverVerbosity,
"Specify the \"chattiness\" of the non-linear solver itself");
EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepVerbosity,
"Specify the \"chattiness\" during the time integration");
EWOMS_REGISTER_PARAM(TypeTag, double, InitialTimeStepInDays,
"The size of the initial time step in days");
EWOMS_REGISTER_PARAM(TypeTag, bool, FullTimeStepInitially,
"Always attempt to finish a report step using a single substep");
EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepAfterEventInDays,
"Time step size of the first time step after an event occurs during the simulation in days");
EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControl,
"The algorithm used to determine time-step sizes. valid options are: 'pid' (default), 'pid+iteration', 'pid+newtoniteration', 'iterationcount' and 'hardcoded'");
EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlTolerance,
"The tolerance used by the time step size control algorithm");
EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetIterations,
"The number of linear iterations which the time step control scheme should aim for (if applicable)");
EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations,
"The number of Newton iterations which the time step control scheme should aim for (if applicable)");
EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayRate,
"The decay rate of the time step size of the number of target iterations is exceeded");
EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthRate,
"The growth rate of the time step size of the number of target iterations is undercut");
EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControlFileName,
"The name of the file which contains the hardcoded time steps sizes");
}
/** \brief step method that acts like the solver::step method
@ -322,35 +403,32 @@ namespace Opm {
protected:
void init_(const ParameterGroup& param)
void init_()
{
// valid are "pid" and "pid+iteration"
std::string control = param.getDefault("timestep.control", std::string("pid"));
// iterations is the accumulation of all linear iterations over all newton steops per time step
const int defaultTargetIterations = 30;
const int defaultTargetNewtonIterations = 8;
std::string control = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControl); // "pid"
const double tol = param.getDefault("timestep.control.tol", double(1e-1));
const double tol = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlTolerance); // 1e-1
if (control == "pid") {
timeStepControl_ = TimeStepControlType(new PIDTimeStepControl(tol));
}
else if (control == "pid+iteration") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
}
else if (control == "pid+newtoniteration") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations);
const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
useNewtonIteration_ = true;
}
else if (control == "iterationcount") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
const double decayrate = param.getDefault("timestep.control.decayrate", double(0.75));
const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25));
const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
}
else if (control == "hardcoded") {
const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps"));
const std::string filename = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControlFileName); // "timesteps"
timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));
}
@ -361,7 +439,6 @@ namespace Opm {
assert(growthFactor_ >= 1.0);
}
typedef std::unique_ptr<TimeStepControlInterface> TimeStepControlType;
SimulatorReport failureReport_; //!< statistics for the failed substeps of the last timestep
@ -370,9 +447,9 @@ namespace Opm {
double growthFactor_; //!< factor to multiply time step when solver recovered from failed convergence
double maxGrowth_; //!< factor that limits the maximum growth of a time step
double maxTimeStep_; //!< maximal allowed time step size
const int solverRestartMax_; //!< how many restart of solver are allowed
const bool solverVerbose_; //!< solver verbosity
const bool timestepVerbose_; //!< timestep verbosity
int solverRestartMax_; //!< how many restart of solver are allowed
bool solverVerbose_; //!< solver verbosity
bool timestepVerbose_; //!< timestep verbosity
double suggestedNextTimestep_; //!< suggested size of next timestep
bool fullTimestepInitially_; //!< beginning with the size of the time step from data file
double timestepAfterEvent_; //!< suggested size of timestep after an event

6
tests/run-init-regressionTest.sh
View File

@ -21,7 +21,11 @@ TEST_ARGS="$@"
rm -Rf ${RESULT_PATH}
mkdir -p ${RESULT_PATH}
cd ${RESULT_PATH}
${BINPATH}/${EXE_NAME} ${TEST_ARGS} nosim=true output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${BINPATH}/${EXE_NAME} ${TEST_ARGS} --enable-dry-run=true --output-dir=${RESULT_PATH}
else
${BINPATH}/${EXE_NAME} ${TEST_ARGS} nosim=true output_dir=${RESULT_PATH}
fi
cd ..
ecode=0

13
tests/run-parallel-regressionTest.sh
View File

@ -19,11 +19,20 @@ TEST_ARGS="$@"
rm -Rf ${RESULT_PATH}
mkdir -p ${RESULT_PATH}
cd ${RESULT_PATH}
${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA linear_solver_reduction=1e-7 tolerance_cnv=5e-6 tolerance_mb=1e-8 output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8 --output-dir=${RESULT_PATH}
else
${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA linear_solver_reduction=1e-7 tolerance_cnv=5e-6 tolerance_mb=1e-8 output_dir=${RESULT_PATH}
fi
test $? -eq 0 || exit 1
mkdir mpi
cd mpi
mpirun -np 4 ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA linear_solver_reduction=1e-7 tolerance_cnv=5e-6 tolerance_mb=1e-8 output_dir=${RESULT_PATH}/mpi
if test "${EXE_NAME}" = "flow"; then
mpirun -np 4 ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8 --output-dir=${RESULT_PATH}/mpi
else
mpirun -np 4 ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA linear_solver_reduction=1e-7 tolerance_cnv=5e-6 tolerance_mb=1e-8 output_dir=${RESULT_PATH}/mpi
fi
test $? -eq 0 || exit 1
cd ..

7
tests/run-porv-acceptanceTest.sh
View File

@ -21,9 +21,14 @@ TEST_ARGS="$@"
rm -Rf ${RESULT_PATH}
mkdir -p ${RESULT_PATH}
cd ${RESULT_PATH}
${BINPATH}/${EXE_NAME} ${TEST_ARGS} nosim=true output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${BINPATH}/${EXE_NAME} ${TEST_ARGS} --enable-dry-run=true --output-dir=${RESULT_PATH}
else
${BINPATH}/${EXE_NAME} ${TEST_ARGS} nosim=true output_dir=${RESULT_PATH}
fi
cd ..
ecode=0
${COMPARE_ECL_COMMAND} -t INIT -k PORV ${RESULT_PATH}/${FILENAME} ${INPUT_DATA_PATH}/opm-porevolume-reference/${EXE_NAME}/${FILENAME} ${ABS_TOL} ${REL_TOL}
if [ $? -ne 0 ]

6
tests/run-regressionTest.sh
View File

@ -17,7 +17,11 @@ TEST_ARGS="$@"
mkdir -p ${RESULT_PATH}
cd ${RESULT_PATH}
${BINPATH}/${EXE_NAME} ${TEST_ARGS} output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${BINPATH}/${EXE_NAME} ${TEST_ARGS} --output-dir=${RESULT_PATH}
else
${BINPATH}/${EXE_NAME} ${TEST_ARGS} output_dir=${RESULT_PATH}
fi
test $? -eq 0 || exit 1
cd ..

14
tests/run-restart-regressionTest.sh
View File

@ -29,11 +29,21 @@ then
else
CMD_PREFIX=""
fi
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA timestep.adaptive=false output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA --enable-adaptive-time-stepping=false --output-dir=${RESULT_PATH}
else
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${TEST_ARGS}.DATA timestep.adaptive=false output_dir=${RESULT_PATH}
fi
test $? -eq 0 || exit 1
${OPM_PACK_COMMAND} -o ${BASE_NAME} ${TEST_ARGS}_RESTART.DATA
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${BASE_NAME} timestep.adaptive=false output_dir=${RESULT_PATH}
if test "${EXE_NAME}" = "flow"; then
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${BASE_NAME} --enable-adaptive-time-stepping=false --output-dir=${RESULT_PATH}
else
${CMD_PREFIX} ${BINPATH}/${EXE_NAME} ${BASE_NAME} timestep.adaptive=false output_dir=${RESULT_PATH}
fi
test $? -eq 0 || exit 1
ecode=0

41
tests/test_wellmodel.cpp
View File

@ -44,7 +44,8 @@
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/FlowMainEbos.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
@ -56,12 +57,11 @@
#include <opm/autodiff/StandardWell.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
// maybe should just include BlackoilModelEbos.hpp
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
}
}
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
using StandardWell = Opm::StandardWell<TTAG(EclFlowProblem)>;
@ -120,6 +120,25 @@ struct SetupTest {
int current_timestep;
};
struct GlobalFixture {
GlobalFixture()
{
int argcDummy = 1;
const char *tmp[] = {"test_wellmodel"};
char **argvDummy = const_cast<char**>(tmp);
// MPI setup.
#if HAVE_DUNE_FEM
Dune::Fem::MPIManager::initialize(argcDummy, argvDummy);
#else
Dune::MPIHelper::instance(argcDummy, argvDummy);
#endif
Opm::FlowMainEbos<TTAG(EclFlowProblem)>::setupParameters_(argcDummy, argvDummy);
}
};
BOOST_GLOBAL_FIXTURE(GlobalFixture);
BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const SetupTest setup_test;
@ -127,7 +146,7 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
BOOST_CHECK_EQUAL( wells_ecl.size(), 2);
const Opm::Well* well = wells_ecl[1];
const Opm::BlackoilModelParameters param;
const Opm::BlackoilModelParametersEbos<TTAG(EclFlowProblem) > param;
// For the conversion between the surface volume rate and resrevoir voidage rate
typedef Opm::BlackOilFluidSystem<double> FluidSystem;
@ -144,8 +163,8 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const int num_comp = wells->number_of_phases;
BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param , *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param , *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells , param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr , param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
}
@ -158,7 +177,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{
const int nw = wells_struct ? (wells_struct->number_of_wells) : 0;
const Opm::BlackoilModelParameters param;
const Opm::BlackoilModelParametersEbos<TTAG(EclFlowProblem)> param;
for (int w = 0; w < nw; ++w) {
const std::string well_name(wells_struct->name[w]);