From 98c5b8fba749a840110035815d3d6ec41b24a6f8 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 29 Nov 2013 16:33:46 +0100 Subject: [PATCH] reformat the source code this was done semi-automatically. The line length has been set to 80 characters, but with a quite low penalty for crossing this limit. --- examples/tutorial1.cpp | 15 +- examples/tutorial1problem.hh | 128 +++++---- tests/models/co2injection_flash.cc | 13 +- tests/models/co2injection_flash_ni.cc | 17 +- tests/models/co2injection_immiscible.cc | 8 +- tests/models/co2injection_immiscible_ni.cc | 5 +- tests/models/co2injection_ncp.cc | 5 +- tests/models/co2injection_ncp_ni.cc | 5 +- tests/models/co2injection_pvs.cc | 11 +- tests/models/co2injection_pvs_ni.cc | 11 +- tests/models/cuvette_pvs.cc | 2 +- tests/models/diffusion_flash.cc | 2 +- tests/models/diffusion_ncp.cc | 2 +- tests/models/diffusion_pvs.cc | 2 +- tests/models/finger_immiscible.cc | 8 +- tests/models/fracture_discretefracture.cc | 4 +- tests/models/groundwater_immiscible.cc | 8 +- tests/models/infiltration_pvs.cc | 9 +- tests/models/lens_immiscible.cc | 5 +- tests/models/lens_richards.cc | 2 +- tests/models/obstacle_immiscible.cc | 5 +- tests/models/obstacle_ncp.cc | 5 +- tests/models/obstacle_pvs.cc | 8 +- tests/models/outflow_pvs.cc | 5 +- tests/models/powerinjection_darcy.cc | 11 +- tests/models/powerinjection_forchheimer.cc | 11 +- tests/models/problems/co2injectionflash.hh | 10 +- tests/models/problems/co2injectionproblem.hh | 156 ++++++----- tests/models/problems/cuvetteproblem.hh | 172 ++++++------ tests/models/problems/diffusionproblem.hh | 46 ++-- tests/models/problems/fingergridcreator.hh | 136 ++++++---- tests/models/problems/fingerproblem.hh | 85 +++--- tests/models/problems/fractureproblem.hh | 118 +++++---- tests/models/problems/groundwaterproblem.hh | 99 ++++--- tests/models/problems/infiltrationproblem.hh | 108 ++++---- tests/models/problems/lensgridcreator.hh | 135 ++++++---- tests/models/problems/lensproblem.hh | 106 ++++---- .../problems/navierstokestestproblem.hh | 60 +++-- tests/models/problems/obstacleproblem.hh | 105 ++++---- tests/models/problems/outflowproblem.hh | 52 ++-- .../models/problems/powerinjectionproblem.hh | 49 ++-- tests/models/problems/reservoirproblem.hh | 245 +++++++++--------- tests/models/problems/richardslensproblem.hh | 72 ++--- tests/models/problems/stokes2ctestproblem.hh | 54 ++-- tests/models/problems/stokesnitestproblem.hh | 74 +++--- tests/models/problems/stokestestproblem.hh | 51 ++-- tests/models/problems/waterairproblem.hh | 114 ++++---- tests/models/reservoir_blackoil.cc | 5 +- tests/models/test_navierstokes.cc | 5 +- tests/models/test_quadrature.cc | 136 +++++----- tests/models/test_stokes.cc | 2 +- tests/models/test_stokes2c.cc | 2 +- tests/models/test_stokesni.cc | 2 +- tests/models/waterair_pvs_ni.cc | 5 +- 54 files changed, 1363 insertions(+), 1148 deletions(-) diff --git a/examples/tutorial1.cpp b/examples/tutorial1.cpp index 23a67c289..55a4f3e52 100644 --- a/examples/tutorial1.cpp +++ b/examples/tutorial1.cpp @@ -19,14 +19,17 @@ /*! * \file * - * \brief Main file of the tutorial for a fully coupled twophase VCVF discretization. + * \brief Main file of the tutorial for a fully coupled twophase VCVF + *discretization. */ -#include "config.h" /*@\label{tutorial-coupled:include-begin}@*/ +#include "config.h" /*@\label{tutorial-coupled:include-begin}@*/ #include /*@\label{tutorial-coupled:include-end}@*/ -#include "tutorial1problem.hh" /*@\label{tutorial-coupled:include-problem-header}@*/ +#include "tutorial1problem.hh" /*@\label{tutorial-coupled:include-problem-header}@*/ -int main(int argc, char** argv) +int main(int argc, char **argv) { - typedef TTAG(TutorialProblemCoupled) TypeTag; /*@\label{tutorial-coupled:set-type-tag}@*/ - return Ewoms::start(argc, argv); /*@\label{tutorial-coupled:call-start}@*/ + typedef TTAG( + TutorialProblemCoupled) TypeTag; /*@\label{tutorial-coupled:set-type-tag}@*/ + return Ewoms::start(argc, + argv); /*@\label{tutorial-coupled:call-start}@*/ } diff --git a/examples/tutorial1problem.hh b/examples/tutorial1problem.hh index f06f80dc1..aad6ab6fc 100644 --- a/examples/tutorial1problem.hh +++ b/examples/tutorial1problem.hh @@ -21,8 +21,10 @@ * * \copydoc Ewoms::TutorialProblemCoupled */ -#ifndef EWOMS_TUTORIAL1_PROBLEM_HH // guardian macro /*@\label{tutorial-coupled:guardian1}@*/ -#define EWOMS_TUTORIAL1_PROBLEM_HH // guardian macro /*@\label{tutorial-coupled:guardian2}@*/ +#ifndef EWOMS_TUTORIAL1_PROBLEM_HH // guardian macro + // /*@\label{tutorial-coupled:guardian1}@*/ +#define EWOMS_TUTORIAL1_PROBLEM_HH // guardian macro + // /*@\label{tutorial-coupled:guardian2}@*/ // The numerical model #include @@ -32,13 +34,16 @@ #include // Headers required for the capillary pressure law -#include /*@\label{tutorial-coupled:rawLawInclude}@*/ +#include \ + /*@\label{tutorial-coupled:rawLawInclude}@*/ #include #include // For the DUNE grid -#include /*@\label{tutorial-coupled:include-grid-manager}@*/ -#include /*@\label{tutorial-coupled:include-grid-creator}@*/ +#include \ + /*@\label{tutorial-coupled:include-grid-manager}@*/ +#include \ + /*@\label{tutorial-coupled:include-grid-creator}@*/ // For Dune::FieldMatrix #include @@ -53,25 +58,37 @@ class TutorialProblemCoupled; namespace Opm { namespace Properties { // Create a new type tag for the problem -NEW_TYPE_TAG(TutorialProblemCoupled, INHERITS_FROM(VcfvImmiscibleTwoPhase)); /*@\label{tutorial-coupled:create-type-tag}@*/ +NEW_TYPE_TAG(TutorialProblemCoupled, + INHERITS_FROM(VcfvImmiscibleTwoPhase)); /*@\label + {tutorial-coupled:create-type-tag}@*/ // Set the "Problem" property -SET_PROP(TutorialProblemCoupled, Problem) /*@\label{tutorial-coupled:set-problem}@*/ -{ typedef Ewoms::TutorialProblemCoupled type;}; +SET_PROP(TutorialProblemCoupled, + Problem) /*@\label{tutorial-coupled:set-problem}@*/ +{ + typedef Ewoms::TutorialProblemCoupled type; +}; // Set grid and the grid creator to be used -SET_TYPE_PROP(TutorialProblemCoupled, Grid, Dune::YaspGrid); /*@\label{tutorial-coupled:set-grid}@*/ -SET_TYPE_PROP(TutorialProblemCoupled, GridCreator, Ewoms::CubeGridCreator); /*@\label{tutorial-coupled:set-gridcreator}@*/ +SET_TYPE_PROP(TutorialProblemCoupled, Grid, + Dune::YaspGrid); /*@\label{tutorial-coupled:set-grid}@*/ +SET_TYPE_PROP(TutorialProblemCoupled, GridCreator, + Ewoms::CubeGridCreator); /*@\label{tutorial-coupled:set-gridcreator}@*/ // Set the wetting phase /*@\label{tutorial-coupled:2p-system-start}@*/ -SET_TYPE_PROP(TutorialProblemCoupled, WettingPhase, /*@\label{tutorial-coupled:wettingPhase}@*/ - Opm::LiquidPhase >); +SET_TYPE_PROP( + TutorialProblemCoupled, + WettingPhase, /*@\label{tutorial-coupled:wettingPhase}@*/ + Opm::LiquidPhase >); // Set the non-wetting phase -SET_TYPE_PROP(TutorialProblemCoupled, NonwettingPhase, /*@\label{tutorial-coupled:nonwettingPhase}@*/ - Opm::LiquidPhase >); /*@\label{tutorial-coupled:2p-system-end}@*/ +SET_TYPE_PROP( + TutorialProblemCoupled, + NonwettingPhase, /*@\label{tutorial-coupled:nonwettingPhase}@*/ + Opm::LiquidPhase >); /*@\label + {tutorial-coupled:2p-system-end}@*/ // Set the material law SET_PROP(TutorialProblemCoupled, MaterialLaw) @@ -83,11 +100,13 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations - typedef Opm::RegularizedBrooksCorey RawMaterialLaw; /*@\label{tutorial-coupled:rawlaw}@*/ + typedef Opm::RegularizedBrooksCorey + RawMaterialLaw; /*@\label{tutorial-coupled:rawlaw}@*/ public: // Convert absolute saturations into effective ones before passing @@ -96,23 +115,27 @@ public: }; // Disable gravity -SET_BOOL_PROP(TutorialProblemCoupled, EnableGravity, false); /*@\label{tutorial-coupled:gravity}@*/ +SET_BOOL_PROP(TutorialProblemCoupled, EnableGravity, + false); /*@\label{tutorial-coupled:gravity}@*/ -// define how long the simulation should run [s] /*@\label{tutorial-coupled:default-params-begin}@*/ +// define how long the simulation should run [s] +// /*@\label{tutorial-coupled:default-params-begin}@*/ SET_SCALAR_PROP(TutorialProblemCoupled, EndTime, 100e3); // define the size of the initial time step [s] SET_SCALAR_PROP(TutorialProblemCoupled, InitialTimeStepSize, 125.0); // define the physical size of the problem's domain [m] -SET_SCALAR_PROP(TutorialProblemCoupled, DomainSizeX, 300.0); /*@\label{tutorial-coupled:grid-default-params-begin}@*/ +SET_SCALAR_PROP(TutorialProblemCoupled, DomainSizeX, + 300.0); /*@\label{tutorial-coupled:grid-default-params-begin}@*/ SET_SCALAR_PROP(TutorialProblemCoupled, DomainSizeY, 60.0); SET_SCALAR_PROP(TutorialProblemCoupled, DomainSizeZ, 0.0); // // define the number of cells used for discretizing the physical domain SET_INT_PROP(TutorialProblemCoupled, CellsX, 100); SET_INT_PROP(TutorialProblemCoupled, CellsY, 1); -SET_INT_PROP(TutorialProblemCoupled, CellsZ, 1); /*@\label{tutorial-coupled:default-params-end}@*/ +SET_INT_PROP(TutorialProblemCoupled, CellsZ, + 1); /*@\label{tutorial-coupled:default-params-end}@*/ } // namespace Properties } // namespace Opm @@ -120,7 +143,8 @@ namespace Ewoms { //! Tutorial problem using the fully-implicit immiscible model. template class TutorialProblemCoupled - : public GET_PROP_TYPE(TypeTag, BaseProblem) /*@\label{tutorial-coupled:def-problem}@*/ + : public GET_PROP_TYPE(TypeTag, + BaseProblem) /*@\label{tutorial-coupled:def-problem}@*/ { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; @@ -136,11 +160,13 @@ class TutorialProblemCoupled typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; /*@\label{tutorial-coupled:matLawObjectType}@*/ + typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) + MaterialLawParams; /*@\label{tutorial-coupled:matLawObjectType}@*/ // phase indices enum { numPhases = FluidSystem::numPhases }; @@ -154,20 +180,21 @@ class TutorialProblemCoupled public: //! The constructor of the problem TutorialProblemCoupled(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()), #else : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()), #endif - , eps_(3e-6) + eps_(3e-6) { // Use an isotropic and homogeneous intrinsic permeability K_ = this->toDimMatrix_(1e-7); // Parameters of the Brooks-Corey law - materialParams_.setEntryPressure(500.0); // entry pressure [Pa] /*@\label{tutorial-coupled:setLawParams}@*/ + materialParams_.setEntryPressure(500.0); // entry pressure [Pa] + // /*@\label{tutorial-coupled:setLawParams}@*/ materialParams_.setLambda(2); // shape parameter // Set the residual saturations @@ -178,7 +205,8 @@ public: materialParams_.finalize(); } - //! Specifies the problem name. This is used for files generated by the simulation. + //! Specifies the problem name. This is used for files generated by the + // simulation. const char *name() const { return "tutorial_coupled"; } @@ -189,30 +217,34 @@ public: //! Returns the intrinsic permeability tensor [m^2] at a position. template - const DimMatrix &intrinsicPermeability(const Context &context, /*@\label{tutorial-coupled:permeability}@*/ - int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability( + const Context &context, /*@\label{tutorial-coupled:permeability}@*/ + int spaceIdx, int timeIdx) const { return K_; } //! Defines the porosity [-] of the medium at a given position template - Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const /*@\label{tutorial-coupled:porosity}@*/ + Scalar porosity(const Context &context, int spaceIdx, + int timeIdx) const /*@\label{tutorial-coupled:porosity}@*/ { return 0.2; } //! Returns the parameter object for the material law at a given position template - const MaterialLawParams& materialLawParams(const Context &context, /*@\label{tutorial-coupled:matLawParams}@*/ - int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams( + const Context &context, /*@\label{tutorial-coupled:matLawParams}@*/ + int spaceIdx, int timeIdx) const { return materialParams_; } //! Evaluates the boundary conditions. template - void boundary(BoundaryRateVector &values, - const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); if (pos[0] < eps_) { // Free-flow conditions on left boundary - const auto &materialParams = this->materialLawParams(context, spaceIdx, timeIdx); + const auto &materialParams + = this->materialLawParams(context, spaceIdx, timeIdx); Opm::ImmiscibleFluidState fs; Scalar Sw = 1.0; @@ -231,7 +263,7 @@ public: // forced outflow at the right boundary RateVector massRate(0.0); - massRate[contiWEqIdx] = 0.0; // [kg / (s m^2)] + massRate[contiWEqIdx] = 0.0; // [kg / (s m^2)] massRate[contiNEqIdx] = 3e-2; // [kg / (s m^2)] values.setMassRate(massRate); @@ -240,10 +272,12 @@ public: values.setNoFlow(); } - //! Evaluates the source term for all conserved quantities at a given position + //! Evaluates the source term for all conserved quantities at a given + // position //! of the domain [kg/(m^3 * s)]. Positive values mean that mass is created. template - void source(RateVector &source, const Context &context, int spaceIdx, int timeIdx) const + void source(RateVector &source, const Context &context, int spaceIdx, + int timeIdx) const { source[contiWEqIdx] = 0.0; source[contiNEqIdx] = 0.0; @@ -251,8 +285,8 @@ public: //! Evaluates the initial value at a given position in the domain. template - void initial(PrimaryVariables &values, - const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { Opm::ImmiscibleFluidState fs; @@ -266,7 +300,9 @@ public: // set pressure of the wetting phase to 200 kPa = 2 bar Scalar pC[numPhases]; - MaterialLaw::capillaryPressures(pC, materialLawParams(context, spaceIdx, timeIdx), fs); + MaterialLaw::capillaryPressures(pC, materialLawParams(context, spaceIdx, + timeIdx), + fs); fs.setPressure(wPhaseIdx, 200e3); fs.setPressure(nPhaseIdx, 200e3 + pC[nPhaseIdx] - pC[nPhaseIdx]); @@ -276,7 +312,7 @@ public: private: DimMatrix K_; // Object that holds the values/parameters of the selected material law. - MaterialLawParams materialParams_; /*@\label{tutorial-coupled:matParamsObject}@*/ + MaterialLawParams materialParams_; /*@\label{tutorial-coupled:matParamsObject}@*/ // small epsilon value Scalar eps_; diff --git a/tests/models/co2injection_flash.cc b/tests/models/co2injection_flash.cc index 0a21b49d0..84da85e0d 100644 --- a/tests/models/co2injection_flash.cc +++ b/tests/models/co2injection_flash.cc @@ -35,17 +35,18 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionFlashProblem, INHERITS_FROM(VcfvFlash, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionFlashProblem, + INHERITS_FROM(VcfvFlash, Co2InjectionBaseProblem)); // for the flash model we want to use thermodynamic hints or it will // get _very_ slow. SET_BOOL_PROP(Co2InjectionFlashProblem, EnableHints, true); // use the flash solver adapted to the CO2 injection problem -SET_TYPE_PROP(Co2InjectionFlashProblem, - FlashSolver, - Ewoms::Co2InjectionFlash ); +SET_TYPE_PROP( + Co2InjectionFlashProblem, FlashSolver, + Ewoms::Co2InjectionFlash); // the flash model has serious problems with the numerical // precision. if quadruple precision math is available, we use it, @@ -58,7 +59,7 @@ SET_SCALAR_PROP(Co2InjectionFlashProblem, NewtonRelativeTolerance, 1e-5); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionFlashProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_flash_ni.cc b/tests/models/co2injection_flash_ni.cc index e330e7858..1c882f263 100644 --- a/tests/models/co2injection_flash_ni.cc +++ b/tests/models/co2injection_flash_ni.cc @@ -32,7 +32,8 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionFlashNIProblem, INHERITS_FROM(VcfvFlash, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionFlashNIProblem, + INHERITS_FROM(VcfvFlash, Co2InjectionBaseProblem)); SET_BOOL_PROP(Co2InjectionFlashNIProblem, EnableEnergy, true); @@ -41,10 +42,10 @@ SET_BOOL_PROP(Co2InjectionFlashNIProblem, EnableEnergy, true); SET_BOOL_PROP(Co2InjectionFlashNIProblem, EnableHints, true); // use the CO2 injection problem adapted flash solver -SET_TYPE_PROP(Co2InjectionFlashNIProblem, - FlashSolver, - Ewoms::Co2InjectionFlash ); +SET_TYPE_PROP( + Co2InjectionFlashNIProblem, FlashSolver, + Ewoms::Co2InjectionFlash); // the flash model has serious problems with the numerical // precision. if quadruple precision math is available, we use it, @@ -54,10 +55,10 @@ SET_TYPE_PROP(Co2InjectionFlashNIProblem, Scalar, quad); #else SET_SCALAR_PROP(Co2InjectionFlashNIProblem, NewtonRelativeTolerance, 1e-5); #endif +} +} -} } - -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionFlashNIProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_immiscible.cc b/tests/models/co2injection_immiscible.cc index 6e0ab0fb7..d3ff5572f 100644 --- a/tests/models/co2injection_immiscible.cc +++ b/tests/models/co2injection_immiscible.cc @@ -29,13 +29,15 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionImmiscibleProblem, INHERITS_FROM(VcfvImmiscible, Co2InjectionBaseProblem)); -} } +NEW_TYPE_TAG(Co2InjectionImmiscibleProblem, + INHERITS_FROM(VcfvImmiscible, Co2InjectionBaseProblem)); +} +} //////////////////////// // the main function //////////////////////// -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionImmiscibleProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_immiscible_ni.cc b/tests/models/co2injection_immiscible_ni.cc index 1028ec314..4cbfc5053 100644 --- a/tests/models/co2injection_immiscible_ni.cc +++ b/tests/models/co2injection_immiscible_ni.cc @@ -30,7 +30,8 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionImmiscibleNIProblem, INHERITS_FROM(VcfvImmiscible, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionImmiscibleNIProblem, + INHERITS_FROM(VcfvImmiscible, Co2InjectionBaseProblem)); SET_BOOL_PROP(Co2InjectionImmiscibleNIProblem, EnableEnergy, true); } @@ -39,7 +40,7 @@ SET_BOOL_PROP(Co2InjectionImmiscibleNIProblem, EnableEnergy, true); //////////////////////// // the main function //////////////////////// -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionImmiscibleNIProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_ncp.cc b/tests/models/co2injection_ncp.cc index 60e9e573b..3fa9a09c0 100644 --- a/tests/models/co2injection_ncp.cc +++ b/tests/models/co2injection_ncp.cc @@ -30,11 +30,12 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionNcpProblem, INHERITS_FROM(VcfvNcp, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionNcpProblem, + INHERITS_FROM(VcfvNcp, Co2InjectionBaseProblem)); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionNcpProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_ncp_ni.cc b/tests/models/co2injection_ncp_ni.cc index 7797b3760..3c430f194 100644 --- a/tests/models/co2injection_ncp_ni.cc +++ b/tests/models/co2injection_ncp_ni.cc @@ -30,12 +30,13 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionNcpNIProblem, INHERITS_FROM(VcfvNcp, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionNcpNIProblem, + INHERITS_FROM(VcfvNcp, Co2InjectionBaseProblem)); SET_BOOL_PROP(Co2InjectionNcpNIProblem, EnableEnergy, true); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionNcpNIProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_pvs.cc b/tests/models/co2injection_pvs.cc index 9e2b41f82..3759dfd1b 100644 --- a/tests/models/co2injection_pvs.cc +++ b/tests/models/co2injection_pvs.cc @@ -19,7 +19,8 @@ /*! * \file * - * \brief Test for the isothermal primary variable switching VCVF discretization. + * \brief Test for the isothermal primary variable switching VCVF + *discretization. */ #include "config.h" @@ -29,10 +30,12 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionPvsProblem, INHERITS_FROM(VcfvPvs, Co2InjectionBaseProblem)); -} } +NEW_TYPE_TAG(Co2InjectionPvsProblem, + INHERITS_FROM(VcfvPvs, Co2InjectionBaseProblem)); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionPvsProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/co2injection_pvs_ni.cc b/tests/models/co2injection_pvs_ni.cc index 689bdc113..8fbe45ce8 100644 --- a/tests/models/co2injection_pvs_ni.cc +++ b/tests/models/co2injection_pvs_ni.cc @@ -19,7 +19,8 @@ /*! * \file * - * \brief Test for the non-isothermal primary variable switching VCVF discretization. + * \brief Test for the non-isothermal primary variable switching VCVF + *discretization. */ #include "config.h" @@ -29,12 +30,14 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(Co2InjectionPvsNIProblem, INHERITS_FROM(VcfvPvs, Co2InjectionBaseProblem)); +NEW_TYPE_TAG(Co2InjectionPvsNIProblem, + INHERITS_FROM(VcfvPvs, Co2InjectionBaseProblem)); SET_BOOL_PROP(Co2InjectionPvsNIProblem, EnableEnergy, true); -} } +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Co2InjectionPvsNIProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/cuvette_pvs.cc b/tests/models/cuvette_pvs.cc index 1049d8ddb..e61690dcd 100644 --- a/tests/models/cuvette_pvs.cc +++ b/tests/models/cuvette_pvs.cc @@ -33,7 +33,7 @@ NEW_TYPE_TAG(CuvetteProblem, INHERITS_FROM(VcfvPvs, CuvetteBaseProblem)); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(CuvetteProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/diffusion_flash.cc b/tests/models/diffusion_flash.cc index 9a3944c9f..43218ce0e 100644 --- a/tests/models/diffusion_flash.cc +++ b/tests/models/diffusion_flash.cc @@ -33,7 +33,7 @@ NEW_TYPE_TAG(DiffusionProblem, INHERITS_FROM(VcfvFlash, DiffusionBaseProblem)); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(DiffusionProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/diffusion_ncp.cc b/tests/models/diffusion_ncp.cc index 0cb58fde1..a2e5ef1ce 100644 --- a/tests/models/diffusion_ncp.cc +++ b/tests/models/diffusion_ncp.cc @@ -33,7 +33,7 @@ NEW_TYPE_TAG(DiffusionProblem, INHERITS_FROM(VcfvNcp, DiffusionBaseProblem)); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(DiffusionProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/diffusion_pvs.cc b/tests/models/diffusion_pvs.cc index e50339433..fbfe528ba 100644 --- a/tests/models/diffusion_pvs.cc +++ b/tests/models/diffusion_pvs.cc @@ -33,7 +33,7 @@ NEW_TYPE_TAG(DiffusionProblem, INHERITS_FROM(VcfvPvs, DiffusionBaseProblem)); } } -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(DiffusionProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/finger_immiscible.cc b/tests/models/finger_immiscible.cc index 683bea690..21b1a0c15 100644 --- a/tests/models/finger_immiscible.cc +++ b/tests/models/finger_immiscible.cc @@ -29,10 +29,12 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(FingerProblem, INHERITS_FROM(VcfvImmiscibleTwoPhase, FingerBaseProblem)); -}} +NEW_TYPE_TAG(FingerProblem, + INHERITS_FROM(VcfvImmiscibleTwoPhase, FingerBaseProblem)); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(FingerProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/fracture_discretefracture.cc b/tests/models/fracture_discretefracture.cc index da18d0707..dea6aee35 100644 --- a/tests/models/fracture_discretefracture.cc +++ b/tests/models/fracture_discretefracture.cc @@ -25,13 +25,13 @@ #if !HAVE_ALUGRID #warning "The ALUGrid Dune grid manager is required for this test." -int main(int argc, char** argv) +int main(int argc, char **argv) {} #else #include #include "problems/fractureproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(FractureProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/groundwater_immiscible.cc b/tests/models/groundwater_immiscible.cc index 527226f3b..63260c81e 100644 --- a/tests/models/groundwater_immiscible.cc +++ b/tests/models/groundwater_immiscible.cc @@ -29,10 +29,12 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(GroundWaterProblem, INHERITS_FROM(VcfvImmiscibleOnePhase, GroundWaterBaseProblem)); -}} +NEW_TYPE_TAG(GroundWaterProblem, + INHERITS_FROM(VcfvImmiscibleOnePhase, GroundWaterBaseProblem)); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(GroundWaterProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/infiltration_pvs.cc b/tests/models/infiltration_pvs.cc index b6f1261ea..e486b09ab 100644 --- a/tests/models/infiltration_pvs.cc +++ b/tests/models/infiltration_pvs.cc @@ -29,12 +29,13 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(VcfvPvs, InfiltrationBaseProblem)); -}} +NEW_TYPE_TAG(InfiltrationProblem, + INHERITS_FROM(VcfvPvs, InfiltrationBaseProblem)); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(InfiltrationProblem) ProblemTypeTag; return Ewoms::start(argc, argv); } - diff --git a/tests/models/lens_immiscible.cc b/tests/models/lens_immiscible.cc index ace5490ca..5ff15f574 100644 --- a/tests/models/lens_immiscible.cc +++ b/tests/models/lens_immiscible.cc @@ -31,9 +31,10 @@ namespace Opm { namespace Properties { NEW_TYPE_TAG(LensProblem, INHERITS_FROM(VcfvImmiscibleTwoPhase, LensBaseProblem)); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(LensProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/lens_richards.cc b/tests/models/lens_richards.cc index 7df9c6fb3..91a5671bf 100644 --- a/tests/models/lens_richards.cc +++ b/tests/models/lens_richards.cc @@ -26,7 +26,7 @@ #include #include "problems/richardslensproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(RichardsLensProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/obstacle_immiscible.cc b/tests/models/obstacle_immiscible.cc index 2b146c640..e5eaf1f8a 100644 --- a/tests/models/obstacle_immiscible.cc +++ b/tests/models/obstacle_immiscible.cc @@ -31,9 +31,10 @@ namespace Opm { namespace Properties { NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(VcfvImmiscible, ObstacleBaseProblem)); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(ObstacleProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/obstacle_ncp.cc b/tests/models/obstacle_ncp.cc index 07ff4294c..a3a51f77c 100644 --- a/tests/models/obstacle_ncp.cc +++ b/tests/models/obstacle_ncp.cc @@ -31,9 +31,10 @@ namespace Opm { namespace Properties { NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(VcfvNcp, ObstacleBaseProblem)); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(ObstacleProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/obstacle_pvs.cc b/tests/models/obstacle_pvs.cc index e43249b80..895a22761 100644 --- a/tests/models/obstacle_pvs.cc +++ b/tests/models/obstacle_pvs.cc @@ -20,7 +20,8 @@ /** * \file * - * \brief Test for the isothermal primary variable switching VCVF discretization. + * \brief Test for the isothermal primary variable switching VCVF + *discretization. */ #include "config.h" @@ -34,9 +35,10 @@ NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(VcfvPvs, ObstacleBaseProblem)); // Verbosity of the PVS model (0=silent, 1=medium, 2=chatty) SET_INT_PROP(ObstacleProblem, PvsVerbosity, 1); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(ObstacleProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/outflow_pvs.cc b/tests/models/outflow_pvs.cc index 3d7a32051..32f4c823c 100644 --- a/tests/models/outflow_pvs.cc +++ b/tests/models/outflow_pvs.cc @@ -33,9 +33,10 @@ NEW_TYPE_TAG(OutflowProblem, INHERITS_FROM(VcfvPvs, OutflowBaseProblem)); // Verbosity of the PVS model (0=silent, 1=medium, 2=chatty) SET_INT_PROP(OutflowProblem, PvsVerbosity, 1); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(OutflowProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/powerinjection_darcy.cc b/tests/models/powerinjection_darcy.cc index 4111ba975..5a6798389 100644 --- a/tests/models/powerinjection_darcy.cc +++ b/tests/models/powerinjection_darcy.cc @@ -29,12 +29,15 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(PowerInjectionProblem, INHERITS_FROM(VcfvImmiscibleTwoPhase, PowerInjectionBaseProblem)); +NEW_TYPE_TAG(PowerInjectionProblem, + INHERITS_FROM(VcfvImmiscibleTwoPhase, PowerInjectionBaseProblem)); -SET_TYPE_PROP(PowerInjectionProblem, VelocityModule, Ewoms::VcfvDarcyVelocityModule); -} } +SET_TYPE_PROP(PowerInjectionProblem, VelocityModule, + Ewoms::VcfvDarcyVelocityModule); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(PowerInjectionProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/powerinjection_forchheimer.cc b/tests/models/powerinjection_forchheimer.cc index c1e7c0ff7..9e5e8ef83 100644 --- a/tests/models/powerinjection_forchheimer.cc +++ b/tests/models/powerinjection_forchheimer.cc @@ -29,12 +29,15 @@ namespace Opm { namespace Properties { -NEW_TYPE_TAG(PowerInjectionProblem, INHERITS_FROM(VcfvImmiscibleTwoPhase, PowerInjectionBaseProblem)); +NEW_TYPE_TAG(PowerInjectionProblem, + INHERITS_FROM(VcfvImmiscibleTwoPhase, PowerInjectionBaseProblem)); -SET_TYPE_PROP(PowerInjectionProblem, VelocityModule, Ewoms::VcfvForchheimerVelocityModule); -} } +SET_TYPE_PROP(PowerInjectionProblem, VelocityModule, + Ewoms::VcfvForchheimerVelocityModule); +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(PowerInjectionProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/problems/co2injectionflash.hh b/tests/models/problems/co2injectionflash.hh index d234b7133..6db7685d6 100644 --- a/tests/models/problems/co2injectionflash.hh +++ b/tests/models/problems/co2injectionflash.hh @@ -43,25 +43,25 @@ class Co2InjectionFlash : public Opm::NcpFlash typedef typename ParentType::ComponentVector ComponentVector; enum { numPhases = FluidSystem::numPhases }; + public: /*! * \brief Guess initial values for all quantities. */ template - static void guessInitial(FluidState &fluidState, - ParameterCache ¶mCache, + static void guessInitial(FluidState &fluidState, ParameterCache ¶mCache, const ComponentVector &globalMolarities) { ParentType::guessInitial(fluidState, paramCache, globalMolarities); - for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - // pressure. something close to the reservoid pressure as initial guess + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // pressure. something close to the reservoid pressure as initial + // guess fluidState.setPressure(phaseIdx, 1.0135e6); } paramCache.updateAllPressures(fluidState); } - }; } // namespace Ewoms diff --git a/tests/models/problems/co2injectionproblem.hh b/tests/models/problems/co2injectionproblem.hh index 84a95ccc0..5ad7391a7 100644 --- a/tests/models/problems/co2injectionproblem.hh +++ b/tests/models/problems/co2injectionproblem.hh @@ -55,7 +55,8 @@ class Co2InjectionProblem; namespace Co2Injection { #include -}} // namespace Ewoms +} +} // namespace Ewoms namespace Opm { namespace Properties { @@ -77,18 +78,21 @@ NEW_PROP_TAG(SimulationName); SET_TYPE_PROP(Co2InjectionBaseProblem, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(Co2InjectionBaseProblem, Problem, Ewoms::Co2InjectionProblem); +SET_TYPE_PROP(Co2InjectionBaseProblem, Problem, + Ewoms::Co2InjectionProblem); // Set fluid configuration SET_PROP(Co2InjectionBaseProblem, FluidSystem) -{ private: +{ +private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef Ewoms::Co2Injection::CO2Tables CO2Tables; static const bool useComplexRelations = false; + public: typedef Opm::FluidSystems::BrineCO2 type; - //typedef Opm::FluidSystems::H2ON2 type; + // typedef Opm::FluidSystems::H2ON2 type; }; // Set the material Law @@ -99,7 +103,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::gPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations @@ -181,8 +186,7 @@ namespace Ewoms { * hydrostatic pressure is assumed. */ template -class Co2InjectionProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class Co2InjectionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -190,30 +194,23 @@ class Co2InjectionProblem typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - enum { - // Grid and world dimension - dim = GridView::dimension, - dimWorld = GridView::dimensionworld - }; + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; // copy some indices for convenience typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { - numPhases = FluidSystem::numPhases, - - gPhaseIdx = FluidSystem::gPhaseIdx, - lPhaseIdx = FluidSystem::lPhaseIdx, - - CO2Idx = FluidSystem::CO2Idx, - BrineIdx = FluidSystem::BrineIdx, - - conti0EqIdx = Indices::conti0EqIdx, - contiCO2EqIdx = conti0EqIdx + CO2Idx - }; + enum { numPhases = FluidSystem::numPhases }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; + enum { lPhaseIdx = FluidSystem::lPhaseIdx }; + enum { CO2Idx = FluidSystem::CO2Idx }; + enum { BrineIdx = FluidSystem::BrineIdx }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { contiCO2EqIdx = conti0EqIdx + CO2Idx }; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -229,7 +226,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ Co2InjectionProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -239,20 +236,23 @@ public: { eps_ = 1e-6; - temperatureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow); - temperatureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh); + temperatureLow_ + = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow); + temperatureHigh_ + = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh); nTemperature_ = EWOMS_GET_PARAM(TypeTag, int, FluidSystemNumTemperature); nPressure_ = EWOMS_GET_PARAM(TypeTag, int, FluidSystemNumPressure); pressureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureLow); - pressureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh); + pressureHigh_ + = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh); maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth); temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature); name_ = EWOMS_GET_PARAM(TypeTag, std::string, SimulationName); // initialize the tables of the fluid system - //FluidSystem::init(); + // FluidSystem::init(); FluidSystem::init(/*Tmin=*/temperatureLow_, /*Tmax=*/temperatureHigh_, /*nT=*/nTemperature_, @@ -297,17 +297,33 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow, "The lower temperature [K] for tabulation of the fluid system"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh, "The upper temperature [K] for tabulation of the fluid system"); - EWOMS_REGISTER_PARAM(TypeTag, int, FluidSystemNumTemperature, "The number of intervals between the lower and upper temperature"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow, + "The lower temperature [K] for tabulation of the " + "fluid system"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh, + "The upper temperature [K] for tabulation of the " + "fluid system"); + EWOMS_REGISTER_PARAM(TypeTag, int, FluidSystemNumTemperature, + "The number of intervals between the lower and " + "upper temperature"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureLow, "The lower pressure [Pa] for tabulation of the fluid system"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureHigh, "The upper pressure [Pa] for tabulation of the fluid system"); - EWOMS_REGISTER_PARAM(TypeTag, int, FluidSystemNumPressure, "The number of intervals between the lower and upper pressure"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureLow, + "The lower pressure [Pa] for tabulation of the " + "fluid system"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureHigh, + "The upper pressure [Pa] for tabulation of the " + "fluid system"); + EWOMS_REGISTER_PARAM(TypeTag, int, FluidSystemNumPressure, + "The number of intervals between the lower and " + "upper pressure"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature, "The temperature [K] in the reservoir"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth, "The maximum depth [m] of the reservoir"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName, "The name of the simulation used for the output files"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature, + "The temperature [K] in the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth, + "The maximum depth [m] of the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName, + "The name of the simulation used for the output " + "files"); } /*! @@ -339,8 +355,8 @@ public: // Write mass balance information for rank 0 if (this->gridView().comm().rank() == 0) { - std::cout<<"Storage: liquid=[" << storageL << "]" - << " gas=[" << storageG << "]\n"; + std::cout << "Storage: liquid=[" << storageL << "]" + << " gas=[" << storageG << "]\n"; } } @@ -360,7 +376,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -384,7 +401,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -398,18 +416,18 @@ public: * In this case, we assume the rock-matrix to be granite. */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 790 // specific heat capacity of granite [J / (kg K)] - * 2700; // density of granite [kg/m^3] + return 790 // specific heat capacity of granite [J / (kg K)] + * 2700; // density of granite [kg/m^3] } /*! * \copydoc VcfvMultiPhaseProblem::heatConductionParams */ template - const HeatConductionLawParams& + const HeatConductionLawParams & heatConductionParams(const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -429,8 +447,7 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -449,7 +466,9 @@ public: Opm::ImmiscibleFluidState fs; fs.setSaturation(gPhaseIdx, 1.0); - fs.setPressure(gPhaseIdx,context. volVars(spaceIdx, timeIdx).fluidState().pressure(gPhaseIdx)); + fs.setPressure(gPhaseIdx, + context.volVars(spaceIdx, timeIdx).fluidState().pressure( + gPhaseIdx)); fs.setTemperature(temperature(context, spaceIdx, timeIdx)); typename FluidSystem::ParameterCache paramCache; paramCache.updatePhase(fs, gPhaseIdx); @@ -475,13 +494,15 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); - //const auto &matParams = this->materialLawParams(context, spaceIdx, timeIdx); - //values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true); + // const auto &matParams = this->materialLawParams(context, spaceIdx, + // timeIdx); + // values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true); values.assignNaive(fs); } @@ -492,16 +513,16 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} private: template - void initialFluidState_(FluidState &fs, const Context &context, int spaceIdx, int timeIdx) const + void initialFluidState_(FluidState &fs, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -520,11 +541,12 @@ private: // set pressures ////// Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, 1e5); - Scalar depth = maxDepth_ - pos[dim -1]; - Scalar pl = 1e5 - densityL*this->gravity()[dim - 1]*depth; + Scalar depth = maxDepth_ - pos[dim - 1]; + Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth; Scalar pC[numPhases]; - const auto &matParams = this->materialLawParams(context, spaceIdx, timeIdx); + const auto &matParams + = this->materialLawParams(context, spaceIdx, timeIdx); MaterialLaw::capillaryPressures(pC, matParams, fs); fs.setPressure(lPhaseIdx, pl + (pC[lPhaseIdx] - pC[lPhaseIdx])); @@ -539,8 +561,7 @@ private: typename FluidSystem::ParameterCache paramCache; typedef Opm::ComputeFromReferencePhase CFRP; - CFRP::solve(fs, - paramCache, + CFRP::solve(fs, paramCache, /*refPhaseIdx=*/lPhaseIdx, /*setViscosity=*/true, /*setEnthalpy=*/true); @@ -563,8 +584,9 @@ private: Scalar lambdaWater = 0.6; Scalar lambdaGranite = 2.8; - Scalar lambdaWet = std::pow(lambdaGranite, (1-poro)) * std::pow(lambdaWater, poro); - Scalar lambdaDry = std::pow(lambdaGranite, (1-poro)); + Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro)) + * std::pow(lambdaWater, poro); + Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro)); params.setFullySaturatedLambda(gPhaseIdx, lambdaDry); params.setFullySaturatedLambda(lPhaseIdx, lambdaWet); @@ -572,7 +594,7 @@ private: } bool isFineMaterial_(const GlobalPosition &pos) const - { return pos[dim-1] > fineLayerBottom_; } + { return pos[dim - 1] > fineLayerBottom_; } DimMatrix fineK_; DimMatrix coarseK_; @@ -594,7 +616,7 @@ private: int nTemperature_; int nPressure_; - std::string name_ ; + std::string name_; Scalar pressureLow_, pressureHigh_; Scalar temperatureLow_, temperatureHigh_; diff --git a/tests/models/problems/cuvetteproblem.hh b/tests/models/problems/cuvetteproblem.hh index a45d91a87..6c9bb6ef6 100644 --- a/tests/models/problems/cuvetteproblem.hh +++ b/tests/models/problems/cuvetteproblem.hh @@ -65,9 +65,9 @@ SET_TYPE_PROP(CuvetteBaseProblem, Grid, Dune::YaspGrid<2>); SET_TYPE_PROP(CuvetteBaseProblem, Problem, Ewoms::CuvetteProblem); // Set the fluid system -SET_TYPE_PROP(CuvetteBaseProblem, - FluidSystem, - Opm::FluidSystems::H2OAirMesitylene); +SET_TYPE_PROP( + CuvetteBaseProblem, FluidSystem, + Opm::FluidSystems::H2OAirMesitylene); // Enable gravity SET_BOOL_PROP(CuvetteBaseProblem, EnableGravity, true); @@ -94,7 +94,7 @@ public: // material law API typedef Opm::ThreePAdapter type; - //typedef Opm::MpLinearMaterial type; + // typedef Opm::MpLinearMaterial type; }; // Set the heat conduction law @@ -148,7 +148,7 @@ namespace Ewoms { * problem to about 2-3 hours simulation time. Complete remediation * of the domain requires much longer (about 10 days simulated time). */ -template +template class CuvetteProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -158,33 +158,30 @@ class CuvetteProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw; - typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams; + typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) + HeatConductionLawParams; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; // copy some indices for convenience typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { - numPhases = FluidSystem::numPhases, - numComponents = FluidSystem::numComponents, + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { wPhaseIdx = FluidSystem::wPhaseIdx }; + enum { nPhaseIdx = FluidSystem::nPhaseIdx }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; + enum { H2OIdx = FluidSystem::H2OIdx }; + enum { airIdx = FluidSystem::airIdx }; + enum { NAPLIdx = FluidSystem::NAPLIdx }; + enum { conti0EqIdx = Indices::conti0EqIdx }; - wPhaseIdx = FluidSystem::wPhaseIdx, - nPhaseIdx = FluidSystem::nPhaseIdx, - gPhaseIdx = FluidSystem::gPhaseIdx, - - H2OIdx = FluidSystem::H2OIdx, - airIdx = FluidSystem::airIdx, - NAPLIdx = FluidSystem::NAPLIdx, - - conti0EqIdx = Indices::conti0EqIdx, - - // Grid and world dimension - dim = GridView::dimension, - dimWorld = GridView::dimensionworld - }; + // Grid and world dimension + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; typedef typename GridView::ctype CoordScalar; typedef Dune::FieldVector GlobalPosition; @@ -195,14 +192,14 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ CuvetteProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()), #else : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()), #endif - , eps_(1e-6) + eps_(1e-6) { if (Valgrind::IsRunning()) FluidSystem::init(/*minT=*/283.15, /*maxT=*/500.0, /*nT=*/20, @@ -219,7 +216,7 @@ public: finePorosity_ = 0.42; coarsePorosity_ = 0.42; - // parameters for the capillary pressure law +// parameters for the capillary pressure law #if 1 // three-phase van Genuchten law fineMaterialParams_.setVgAlpha(0.0005); @@ -294,7 +291,7 @@ public: */ const char *name() const { - static std::string tmp = std::string("cuvette_")+this->model().name(); + static std::string tmp = std::string("cuvette_") + this->model().name(); return tmp.c_str(); } @@ -316,7 +313,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -341,7 +339,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -354,7 +353,7 @@ public: * \copydoc VcfvMultiPhaseProblem::heatConductionParams */ template - const HeatConductionLawParams& + const HeatConductionLawParams & heatConductionParams(const Context &context, int spaceIdx, int timeIdx) const { return heatCondParams_; } @@ -362,11 +361,11 @@ public: * \copydoc VcfvMultiPhaseProblem::heatCapacitySolid */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 850 // specific heat capacity [J / (kg K)] - * 2650; // density of sand [kg/m^3] + return 850 // specific heat capacity [J / (kg K)] + * 2650; // density of sand [kg/m^3] } //! \} @@ -380,7 +379,8 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -392,27 +392,26 @@ public: values.setFreeFlow(context, spaceIdx, timeIdx, fs); values.setNoFlow(); } - else if (onLeftBoundary_(pos)) - { + else if (onLeftBoundary_(pos)) { // injection RateVector molarRate; // inject with the same composition as the gas phase of // the injection fluid state Scalar molarInjectionRate = 0.3435; // [mol/(m^2 s)] - for (int compIdx = 0; compIdx < numComponents; ++ compIdx) - molarRate[conti0EqIdx + compIdx] = - - molarInjectionRate - * injectFluidState_.moleFraction(gPhaseIdx, compIdx); + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + molarRate[conti0EqIdx + compIdx] + = -molarInjectionRate + * injectFluidState_.moleFraction(gPhaseIdx, compIdx); // calculate the total mass injection rate [kg / (m^2 s) - Scalar massInjectionRate = - molarInjectionRate - * injectFluidState_.averageMolarMass(gPhaseIdx); + Scalar massInjectionRate + = molarInjectionRate + * injectFluidState_.averageMolarMass(gPhaseIdx); // set the boundary rate vector values.setMolarRate(molarRate); - values.setEnthalpyRate(- injectFluidState_.enthalpy(gPhaseIdx) + values.setEnthalpyRate(-injectFluidState_.enthalpy(gPhaseIdx) * massInjectionRate); // [J / (m^2 s)] } else @@ -430,7 +429,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { Opm::CompositionalFluidState fs; @@ -447,7 +447,8 @@ public: * everywhere. */ template - void source(RateVector &rate, const Context &context, int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -467,9 +468,8 @@ private: bool isContaminated_(const GlobalPosition &pos) const { - return - (0.20 <= pos[0]) && (pos[0] <= 0.80) - && (0.4 <= pos[1]) && (pos[1] <= 0.65); + return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1]) + && (pos[1] <= 0.65); } bool isFineMaterial_(const GlobalPosition &pos) const @@ -478,14 +478,13 @@ private: return true; else if (pos[1] <= 0.15 && 1.20 <= pos[0]) return true; - else return false; + else + return false; } template - void initialFluidState_(FluidState &fs, - const Context &context, - int spaceIdx, - int timeIdx) const + void initialFluidState_(FluidState &fs, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -493,23 +492,24 @@ private: Scalar pw = 1e5; - if(isContaminated_(pos)) { + if (isContaminated_(pos)) { fs.setSaturation(wPhaseIdx, 0.12); fs.setSaturation(nPhaseIdx, 0.07); fs.setSaturation(gPhaseIdx, 1 - 0.12 - 0.07); // set the capillary pressures - const auto &matParams = materialLawParams(context, spaceIdx, timeIdx); + const auto &matParams + = materialLawParams(context, spaceIdx, timeIdx); Scalar pc[numPhases]; MaterialLaw::capillaryPressures(pc, matParams, fs); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - fs.setPressure(phaseIdx, - pw + (pc[phaseIdx] - pc[wPhaseIdx])); + fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[wPhaseIdx])); // compute the phase compositions typedef Opm::MiscibleMultiPhaseComposition MMPC; typename FluidSystem::ParameterCache paramCache; - MMPC::solve(fs, paramCache, /*setViscosity=*/true, /*setEnthalpy=*/true); + MMPC::solve(fs, paramCache, /*setViscosity=*/true, + /*setEnthalpy=*/true); } else { fs.setSaturation(wPhaseIdx, 0.12); @@ -517,33 +517,33 @@ private: fs.setSaturation(nPhaseIdx, 0); // set the capillary pressures - const auto &matParams = materialLawParams(context, spaceIdx, timeIdx); + const auto &matParams + = materialLawParams(context, spaceIdx, timeIdx); Scalar pc[numPhases]; MaterialLaw::capillaryPressures(pc, matParams, fs); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - fs.setPressure(phaseIdx, - pw + (pc[phaseIdx] - pc[wPhaseIdx])); + fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[wPhaseIdx])); // compute the phase compositions typedef Opm::MiscibleMultiPhaseComposition MMPC; typename FluidSystem::ParameterCache paramCache; - MMPC::solve(fs, paramCache, /*setViscosity=*/true, /*setEnthalpy=*/true); + MMPC::solve(fs, paramCache, /*setViscosity=*/true, + /*setEnthalpy=*/true); // set the contaminant mole fractions to zero. this is a // little bit hacky... - for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0); if (phaseIdx == nPhaseIdx) continue; Scalar sumx = 0; - for (int compIdx = 0; compIdx < numComponents; ++ compIdx) + for (int compIdx = 0; compIdx < numComponents; ++compIdx) sumx += fs.moleFraction(phaseIdx, compIdx); - for (int compIdx = 0; compIdx < numComponents; ++ compIdx) - fs.setMoleFraction(phaseIdx, - compIdx, + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + fs.setMoleFraction(phaseIdx, compIdx, fs.moleFraction(phaseIdx, compIdx) / sumx); } } @@ -570,12 +570,13 @@ private: for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { Scalar lambdaSaturated; if (FluidSystem::isLiquid(phaseIdx)) { - Scalar lambdaFluid = - FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); - lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro); + Scalar lambdaFluid + = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)) + + std::pow(lambdaFluid, poro); } else - lambdaSaturated = std::pow(lambdaGranite, (1-poro)); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)); params.setFullySaturatedLambda(phaseIdx, lambdaSaturated); if (!FluidSystem::isLiquid(phaseIdx)) @@ -585,26 +586,21 @@ private: void initInjectFluidState_() { - injectFluidState_.setTemperature(383.0); // [K] - injectFluidState_.setPressure(gPhaseIdx, 1e5); // [Pa] + injectFluidState_.setTemperature(383.0); // [K] + injectFluidState_.setPressure(gPhaseIdx, 1e5); // [Pa] injectFluidState_.setSaturation(gPhaseIdx, 1.0); // [-] Scalar xgH2O = 0.417; - injectFluidState_.setMoleFraction(gPhaseIdx, - H2OIdx, - xgH2O); // [-] - injectFluidState_.setMoleFraction(gPhaseIdx, - airIdx, - 1 - xgH2O); // [-] - injectFluidState_.setMoleFraction(gPhaseIdx, - NAPLIdx, - 0.0); // [-] + injectFluidState_.setMoleFraction(gPhaseIdx, H2OIdx, xgH2O); // [-] + injectFluidState_.setMoleFraction(gPhaseIdx, airIdx, 1 - xgH2O); // [-] + injectFluidState_.setMoleFraction(gPhaseIdx, NAPLIdx, 0.0); // [-] // set the specific enthalpy of the gas phase typename FluidSystem::ParameterCache paramCache; paramCache.updatePhase(injectFluidState_, gPhaseIdx); - Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gPhaseIdx); + Scalar h + = FluidSystem::enthalpy(injectFluidState_, paramCache, gPhaseIdx); injectFluidState_.setEnthalpy(gPhaseIdx, h); } diff --git a/tests/models/problems/diffusionproblem.hh b/tests/models/problems/diffusionproblem.hh index eceec2bdf..f5e88a279 100644 --- a/tests/models/problems/diffusionproblem.hh +++ b/tests/models/problems/diffusionproblem.hh @@ -82,7 +82,8 @@ private: typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::gPhaseIdx> + Traits; public: typedef Opm::LinearMaterial type; @@ -108,8 +109,8 @@ SET_SCALAR_PROP(DiffusionBaseProblem, EndTime, 1e6); // The default for the initial time step size of the simulation SET_SCALAR_PROP(DiffusionBaseProblem, InitialTimeStepSize, 1000); -}} // namespace Opm, Properties - +} +} // namespace Opm, Properties namespace Ewoms { /*! @@ -123,8 +124,7 @@ namespace Ewoms { * diffusion. */ template -class DiffusionProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class DiffusionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -152,7 +152,8 @@ class DiffusionProblem }; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -167,7 +168,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ DiffusionProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -212,7 +213,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { return K_; } /*! @@ -226,15 +228,15 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { return materialParams_; } /*! * \copydoc VcfvMultiPhaseProblem::temperature */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return temperature_; } //! \} @@ -250,8 +252,7 @@ public: * This problem sets no-flow boundaries everywhere. */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { values.setNoFlow(); } @@ -266,9 +267,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); if (onLeftSide_(pos)) @@ -284,16 +284,15 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} private: bool onLeftSide_(const GlobalPosition &pos) const - { return pos[0] < (this->bboxMin()[0] + this->bboxMax()[0])/2; } + { return pos[0] < (this->bboxMin()[0] + this->bboxMax()[0]) / 2; } void setupInitialFluidStates_() { @@ -314,15 +313,16 @@ private: typedef Opm::ComputeFromReferencePhase CFRP; typename FluidSystem::ParameterCache paramCache; - CFRP::solve(leftInitialFluidState_, paramCache, gPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/false); + CFRP::solve(leftInitialFluidState_, paramCache, gPhaseIdx, + /*setViscosity=*/false, /*setEnthalpy=*/false); // create the initial fluid state for the right half of the domain rightInitialFluidState_.assign(leftInitialFluidState_); xH2O = 0.0; rightInitialFluidState_.setMoleFraction(gPhaseIdx, H2OIdx, xH2O); rightInitialFluidState_.setMoleFraction(gPhaseIdx, N2Idx, 1 - xH2O); - CFRP::solve(rightInitialFluidState_, paramCache, gPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/false); - + CFRP::solve(rightInitialFluidState_, paramCache, gPhaseIdx, + /*setViscosity=*/false, /*setEnthalpy=*/false); } DimMatrix K_; diff --git a/tests/models/problems/fingergridcreator.hh b/tests/models/problems/fingergridcreator.hh index 07a2ccbb6..b27d2c519 100644 --- a/tests/models/problems/fingergridcreator.hh +++ b/tests/models/problems/fingergridcreator.hh @@ -50,8 +50,7 @@ namespace Opm { ////////// // Specify the properties for the finger problem ////////// -namespace Properties -{ +namespace Properties { // declare the properties required by the for the finger grid creator NEW_PROP_TAG(Grid); NEW_PROP_TAG(Scalar); @@ -88,16 +87,24 @@ public: */ static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, "The number of global refinements of the grid executed after it was loaded"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, "The size of the domain in x direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, "The number of intervalls in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, + "The number of intervalls in x direction"); if (dim > 1) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, "The size of the domain in y direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, "The number of intervalls in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, + "The number of intervalls in y direction"); } if (dim > 2) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, "The size of the domain in z direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, "The number of intervalls in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, + "The number of intervalls in z direction"); } } @@ -123,20 +130,21 @@ public: cellRes[2] = EWOMS_GET_PARAM(TypeTag, int, CellsZ); } - unsigned numRefinments = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + unsigned numRefinments + = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); Dune::GridFactory factory(grid_); if (dim == 3) { - Dune::FieldVector pos; + Dune::FieldVector pos; for (int k = 0; k <= cellRes[0]; k++) { - pos[2] = upperRight[2]*double(k)/cellRes[2]; + pos[2] = upperRight[2] * double(k) / cellRes[2]; for (int j = 0; j <= cellRes[1]; j++) { - pos[1] = upperRight[1]*double(j)/cellRes[1]; + pos[1] = upperRight[1] * double(j) / cellRes[1]; for (int i = 0; i <= cellRes[0]; i++) { - pos[0] = upperRight[0]*double(i)/cellRes[0]; + pos[0] = upperRight[0] * double(i) / cellRes[0]; factory.insertVertex(pos); } } @@ -144,12 +152,12 @@ public: } else { assert(dim == 2); - Dune::FieldVector pos; + Dune::FieldVector pos; for (int j = 0; j <= cellRes[1]; j++) { - pos[1] = upperRight[1]*double(j)/cellRes[1]; + pos[1] = upperRight[1] * double(j) / cellRes[1]; for (int i = 0; i <= cellRes[0]; i++) { - pos[0] = upperRight[0]*double(i)/cellRes[0]; + pos[0] = upperRight[0] * double(i) / cellRes[0]; factory.insertVertex(pos); } } @@ -166,14 +174,14 @@ public: int m = cellRes[0] + 1; int n = cellRes[1] + 1; for (int k = 0; k < cellRes[2]; ++k) { - int i0 = k*m*n + j*m + i; - int i1 = k*m*n + j*m + (i+1); - int i2 = k*m*n + (j+1)*m + i; - int i3 = k*m*n + (j+1)*m + (i+1); - int i4 = (k+1)*m*n + j*m + i; - int i5 = (k+1)*m*n + j*m + (i+1); - int i6 = (k+1)*m*n + (j+1)*m + i; - int i7 = (k+1)*m*n + (j+1)*m + (i+1); + int i0 = k * m * n + j * m + i; + int i1 = k * m * n + j * m + (i + 1); + int i2 = k * m * n + (j + 1) * m + i; + int i3 = k * m * n + (j + 1) * m + (i + 1); + int i4 = (k + 1) * m * n + j * m + i; + int i5 = (k + 1) * m * n + j * m + (i + 1); + int i6 = (k + 1) * m * n + (j + 1) * m + i; + int i7 = (k + 1) * m * n + (j + 1) * m + (i + 1); #if FINGER_CUBES v[0] = i0; @@ -184,44 +192,57 @@ public: v[5] = i5; v[6] = i6; v[7] = i7; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::cube, 3), v); #else v[0] = i0; v[1] = i1; v[2] = i2; v[3] = i4; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i4; v[1] = i5; v[2] = i6; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i2; v[1] = i5; v[2] = i4; v[3] = i1; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i2; v[1] = i3; v[2] = i7; v[3] = i5; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i5; v[1] = i7; v[2] = i6; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i1; v[1] = i3; v[2] = i5; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); #endif } } @@ -229,26 +250,29 @@ public: assert(dim == 2); int m = cellRes[0] + 1; - int i0 = j*m + i; - int i1 = j*m + (i+1); - int i2 = (j+1)*m + i; - int i3 = (j+1)*m + (i+1); + int i0 = j * m + i; + int i1 = j * m + (i + 1); + int i2 = (j + 1) * m + i; + int i3 = (j + 1) * m + (i + 1); #if FINGER_CUBES v[0] = i0; v[1] = i1; v[2] = i2; v[3] = i3; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::cube, 2), v); #else v[0] = i0; v[1] = i1; v[2] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 2), v); v[0] = i1; v[1] = i3; v[2] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 2), v); #endif } } @@ -304,16 +328,24 @@ public: */ static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, "The number of global refinements of the grid executed after it was loaded"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, "The size of the domain in x direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, "The number of intervalls in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, + "The number of intervalls in x direction"); if (dim > 1) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, "The size of the domain in y direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, "The number of intervalls in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, + "The number of intervalls in y direction"); } if (dim > 2) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, "The size of the domain in z direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, "The number of intervalls in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, + "The number of intervalls in z direction"); } } @@ -322,7 +354,7 @@ public: */ static void makeGrid() { -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) std::bitset isPeriodic(false); std::array cellRes; #else @@ -346,16 +378,16 @@ public: cellRes[2] = EWOMS_GET_PARAM(TypeTag, int, CellsZ); } - unsigned numRefinments = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + unsigned numRefinments + = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); grid_ = new Dune::YaspGrid( #ifdef HAVE_MPI Dune::MPIHelper::getCommunicator(), #endif - upperRight, // upper right - cellRes, // number of cells - isPeriodic, - 0); // overlap size + upperRight, // upper right + cellRes, // number of cells + isPeriodic, 0); // overlap size grid_->globalRefine(numRefinments); } diff --git a/tests/models/problems/fingerproblem.hh b/tests/models/problems/fingerproblem.hh index 7b5c667c5..eed2acc7d 100644 --- a/tests/models/problems/fingerproblem.hh +++ b/tests/models/problems/fingerproblem.hh @@ -59,7 +59,8 @@ NEW_TYPE_TAG(FingerBaseProblem); SET_TYPE_PROP(FingerBaseProblem, GridCreator, Ewoms::FingerGridCreator); // Retrieve the grid type from the grid creator -SET_TYPE_PROP(FingerBaseProblem, Grid, typename GET_PROP_TYPE(TypeTag, GridCreator)::Grid); +SET_TYPE_PROP(FingerBaseProblem, Grid, + typename GET_PROP_TYPE(TypeTag, GridCreator)::Grid); // declare the properties specific for the finger problem NEW_PROP_TAG(InitialWaterSaturation); @@ -72,6 +73,7 @@ SET_PROP(FingerBaseProblem, WettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; @@ -81,6 +83,7 @@ SET_PROP(FingerBaseProblem, NonwettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::GasPhase > type; }; @@ -92,7 +95,8 @@ SET_PROP(FingerBaseProblem, MaterialLaw) typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // use the parker-lenhard hysteresis law typedef Opm::ParkerLenhard ParkerLenhard; @@ -100,10 +104,10 @@ SET_PROP(FingerBaseProblem, MaterialLaw) }; // Enable partial reassembly of the jacobian matrix? -//SET_BOOL_PROP(FingerBaseProblem, EnablePartialReassemble, true); +// SET_BOOL_PROP(FingerBaseProblem, EnablePartialReassemble, true); // Enable reuse of jacobian matrices? -//SET_BOOL_PROP(FingerBaseProblem, EnableJacobianRecycling, true); +// SET_BOOL_PROP(FingerBaseProblem, EnableJacobianRecycling, true); // Write the solutions of individual newton iterations? SET_BOOL_PROP(FingerBaseProblem, NewtonWriteConvergence, false); @@ -156,10 +160,9 @@ namespace Ewoms { * discretization is fine enough. */ template -class FingerProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class FingerProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { -//!\cond SKIP_THIS + //!\cond SKIP_THIS typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; @@ -189,7 +192,8 @@ class FingerProblem typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP(TypeTag, MaterialLaw)::ParkerLenhard ParkerLenhard; @@ -200,14 +204,14 @@ class FingerProblem typedef Dune::FieldVector GlobalPosition; typedef Dune::FieldMatrix DimMatrix; -//!\endcond + //!\endcond public: /*! * \copydoc Doxygen::defaultProblemConstructor */ FingerProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -230,7 +234,7 @@ public: * \copydoc VcfvProblem::name */ std::string name() const - { return std::string("finger_")+this->model().name(); } + { return std::string("finger_") + this->model().name(); } /*! * \copydoc VcfvMultiPhaseProblem::registerParameters @@ -239,7 +243,9 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation, "The initial saturation in the domain [] of the wetting phase"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation, + "The initial saturation in the domain [] of the " + "wetting phase"); } /*! @@ -260,11 +266,7 @@ public: // initialize the material parameter objects of the individual // finite volumes -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) - int n = GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView().size(dimWorld); -#else - int n = GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView().size(dimWorld); -#endif + int n = this->model().numDofs(); materialParams_.resize(n); for (int i = 0; i < n; ++i) { materialParams_[i].setMicParams(&micParams_); @@ -292,12 +294,12 @@ public: auto elemIt = this->gridView().template begin<0>(); const auto &elemEndIt = this->gridView().template end<0>(); - for (; elemIt != elemEndIt; ++elemIt) - { + for (; elemIt != elemEndIt; ++elemIt) { elemCtx.updateAll(*elemIt); for (int scvIdx = 0; scvIdx < elemCtx.numScv(); ++scvIdx) { int globalIdx = elemCtx.globalSpaceIndex(scvIdx, /*timeIdx=*/0); - const auto &fs = elemCtx.volVars(scvIdx, /*timeIdx=*/0).fluidState(); + const auto &fs + = elemCtx.volVars(scvIdx, /*timeIdx=*/0).fluidState(); ParkerLenhard::update(materialParams_[globalIdx], fs); } } @@ -314,29 +316,30 @@ public: * \copydoc VcfvMultiPhaseProblem::temperature */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return temperature_; } /*! * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { return K_; } /*! * \copydoc VcfvMultiPhaseProblem::porosity */ - template - Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const - { return 0.4; } + template + Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const + { return 0.4; } /*! * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); return materialParams_[globalSpaceIdx]; @@ -353,13 +356,13 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.cvCenter(spaceIdx, timeIdx); - if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos)) { + if (onLeftBoundary_(pos) || onRightBoundary_(pos) + || onLowerBoundary_(pos)) { values.setNoFlow(); } else { @@ -371,7 +374,7 @@ public: // override the value for the liquid phase by forced // imbibition of water on inlet boundary segments if (onInlet_(pos)) { - values[contiWEqIdx] = - 0.001; // [kg/(m^2 s)] + values[contiWEqIdx] = -0.001; // [kg/(m^2 s)] } } @@ -386,9 +389,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { // assign the primary variables values.assignNaive(initialFluidState_); @@ -398,8 +400,7 @@ public: * \copydoc VcfvProblem::constraints */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -421,9 +422,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -443,15 +443,16 @@ private: bool onInlet_(const GlobalPosition &pos) const { Scalar width = this->bboxMax()[0] - this->bboxMin()[0]; - Scalar lambda = (this->bboxMax()[0] - pos[0])/width; + Scalar lambda = (this->bboxMax()[0] - pos[0]) / width; if (!onUpperBoundary_(pos)) return false; Scalar xInject[] = { 0.25, 0.75 }; Scalar injectLen[] = { 0.1, 0.1 }; - for (unsigned i = 0; i < sizeof(xInject)/sizeof(Scalar); ++ i) { - if (xInject[i] - injectLen[i]/2 < lambda && lambda < xInject[i] + injectLen[i]/2) + for (unsigned i = 0; i < sizeof(xInject) / sizeof(Scalar); ++i) { + if (xInject[i] - injectLen[i] / 2 < lambda + && lambda < xInject[i] + injectLen[i] / 2) return true; } return false; diff --git a/tests/models/problems/fractureproblem.hh b/tests/models/problems/fractureproblem.hh index e4f406481..f9af6a82c 100644 --- a/tests/models/problems/fractureproblem.hh +++ b/tests/models/problems/fractureproblem.hh @@ -63,9 +63,9 @@ NEW_TYPE_TAG(FractureProblem, INHERITS_FROM(VcfvDiscreteFracture)); SET_TYPE_PROP(FractureProblem, GridCreator, Ewoms::ArtGridCreator); // Set the grid type -SET_TYPE_PROP(FractureProblem, - Grid, - Dune::ALUGrid); +SET_TYPE_PROP( + FractureProblem, Grid, + Dune::ALUGrid); // Set the problem property SET_TYPE_PROP(FractureProblem, Problem, Ewoms::FractureProblem); @@ -98,13 +98,14 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations typedef Opm::RegularizedBrooksCorey EffectiveLaw; - //typedef RegularizedVanGenuchten EffectiveLaw; - //typedef LinearMaterial EffectiveLaw; + // typedef RegularizedVanGenuchten EffectiveLaw; + // typedef LinearMaterial EffectiveLaw; public: typedef Opm::EffToAbsLaw type; @@ -155,9 +156,8 @@ namespace Ewoms { * the fractures and gradually pushes the oil out on the right side, * where the pressure is kept constant. */ -template -class FractureProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +template +class FractureProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; @@ -167,13 +167,15 @@ class FractureProblem typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; - typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams; + typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) + HeatConductionLawParams; enum { // phase indices @@ -193,10 +195,10 @@ class FractureProblem typedef Dune::FieldVector GlobalPosition; typedef Dune::FieldMatrix DimMatrix; - template + template struct FaceLayout { - bool contains (Dune::GeometryType gt) + bool contains(Dune::GeometryType gt) { return gt.dim() == dim - 1; } }; typedef Dune::MultipleCodimMultipleGeomTypeMapper FaceMapper; @@ -208,7 +210,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ FractureProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -250,8 +252,8 @@ public: matrixMaterialParams_.finalize(); fractureMaterialParams_.finalize(); - matrixK_ = this->toDimMatrix_(1e-15); //m^2 - fractureK_ = this->toDimMatrix_(1e5*1e-15); //m^2 + matrixK_ = this->toDimMatrix_(1e-15); // m^2 + fractureK_ = this->toDimMatrix_(1e5 * 1e-15); // m^2 matrixPorosity_ = 0.10; fracturePorosity_ = 0.25; @@ -296,8 +298,7 @@ public: * \copydoc VcfvMultiPhaseProblem::temperature */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return temperature_; } // \} @@ -311,7 +312,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { return matrixK_; } /*! @@ -320,7 +322,9 @@ public: * \copydoc Doxygen::contextParams */ template - const DimMatrix &fractureIntrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &fractureIntrinsicPermeability(const Context &context, + int spaceIdx, + int timeIdx) const { return fractureK_; } /*! @@ -336,14 +340,16 @@ public: * \copydoc Doxygen::contextParams */ template - Scalar fracturePorosity(const Context &context, int spaceIdx, int timeIdx) const + Scalar fracturePorosity(const Context &context, int spaceIdx, + int timeIdx) const { return fracturePorosity_; } /*! * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { return matrixMaterialParams_; } /*! @@ -352,7 +358,9 @@ public: * \copydoc Doxygen::contextParams */ template - const MaterialLawParams& fractureMaterialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &fractureMaterialLawParams(const Context &context, + int spaceIdx, + int timeIdx) const { return fractureMaterialParams_; } /*! @@ -374,14 +382,15 @@ public: * \param timeIdx The index used by the time discretization. */ template - Scalar fractureWidth(const Context &context, int spaceIdx1, int spaceIdx2, int timeIdx) const + Scalar fractureWidth(const Context &context, int spaceIdx1, int spaceIdx2, + int timeIdx) const { return fractureWidth_; } /*! * \copydoc VcfvMultiPhaseProblem::heatConductionParams */ template - const HeatConductionLawParams& + const HeatConductionLawParams & heatConductionParams(const Context &context, int spaceIdx, int timeIdx) const { return heatCondParams_; } @@ -391,11 +400,11 @@ public: * In this case, we assume the rock-matrix to be granite. */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 790 // specific heat capacity of granite [J / (kg K)] - * 2700; // density of granite [kg/m^3] + return 790 // specific heat capacity of granite [J / (kg K)] + * 2700; // density of granite [kg/m^3] } // \} @@ -409,21 +418,20 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); - if (onRightBoundary_(pos)) - { + if (onRightBoundary_(pos)) { // on the right boundary, we impose a free-flow // (i.e. Dirichlet) condition FluidState fluidState; fluidState.setTemperature(temperature_); fluidState.setSaturation(wPhaseIdx, 0.0); - fluidState.setSaturation(nPhaseIdx, 1.0 - fluidState.saturation(wPhaseIdx)); + fluidState.setSaturation(nPhaseIdx, + 1.0 - fluidState.saturation(wPhaseIdx)); fluidState.setPressure(wPhaseIdx, 1e5); fluidState.setPressure(nPhaseIdx, fluidState.pressure(wPhaseIdx)); @@ -448,8 +456,7 @@ public: * \copydoc VcfvProblem::constraints */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -472,25 +479,31 @@ public: fractureFluidState.setTemperature(temperature_ + 10); fractureFluidState.setSaturation(wPhaseIdx, 1.0); - fractureFluidState.setSaturation(nPhaseIdx, 1.0 - fractureFluidState.saturation(wPhaseIdx)); + fractureFluidState.setSaturation(nPhaseIdx, + 1.0 - fractureFluidState.saturation( + wPhaseIdx)); Scalar pCFracture[numPhases]; - MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_, fractureFluidState); + MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_, + fractureFluidState); fractureFluidState.setPressure(wPhaseIdx, /*pressure=*/1e5); - fractureFluidState.setPressure(nPhaseIdx, fractureFluidState.pressure(wPhaseIdx) + (pCFracture[nPhaseIdx] - pCFracture[wPhaseIdx])); + fractureFluidState.setPressure(nPhaseIdx, + fractureFluidState.pressure(wPhaseIdx) + + (pCFracture[nPhaseIdx] + - pCFracture[wPhaseIdx])); constraints.setAllConstraint(); - constraints.assignNaiveFromFracture(fractureFluidState, matrixMaterialParams_); + constraints.assignNaiveFromFracture(fractureFluidState, + matrixMaterialParams_); } /*! * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { FluidState fluidState; fluidState.setTemperature(temperature_); @@ -498,7 +511,8 @@ public: fluidState.setPressure(nPhaseIdx, fluidState.pressure(wPhaseIdx)); fluidState.setSaturation(wPhaseIdx, 0.0); - fluidState.setSaturation(nPhaseIdx, 1.0 - fluidState.saturation(wPhaseIdx)); + fluidState.setSaturation(nPhaseIdx, + 1.0 - fluidState.saturation(wPhaseIdx)); values.assignNaive(fluidState); } @@ -510,9 +524,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } // \} @@ -551,17 +564,18 @@ private: for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { Scalar lambdaSaturated; if (FluidSystem::isLiquid(phaseIdx)) { - Scalar lambdaFluid = - FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); - lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro); + Scalar lambdaFluid + = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)) + + std::pow(lambdaFluid, poro); } else - lambdaSaturated = std::pow(lambdaGranite, (1-poro)); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)); params.setFullySaturatedLambda(phaseIdx, lambdaSaturated); } - Scalar lambdaVac = std::pow(lambdaGranite, (1-poro)); + Scalar lambdaVac = std::pow(lambdaGranite, (1 - poro)); params.setVacuumLambda(lambdaVac); } diff --git a/tests/models/problems/groundwaterproblem.hh b/tests/models/problems/groundwaterproblem.hh index 9ab943c9e..9a76dc814 100644 --- a/tests/models/problems/groundwaterproblem.hh +++ b/tests/models/problems/groundwaterproblem.hh @@ -62,15 +62,17 @@ SET_PROP(GroundWaterBaseProblem, Fluid) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; // Set the grid type SET_TYPE_PROP(GroundWaterBaseProblem, Grid, Dune::YaspGrid<2>); -//SET_TYPE_PROP(GroundWaterBaseProblem, Grid, Dune::SGrid<2, 2>); +// SET_TYPE_PROP(GroundWaterBaseProblem, Grid, Dune::SGrid<2, 2>); -SET_TYPE_PROP(GroundWaterBaseProblem, Problem, Ewoms::GroundWaterProblem); +SET_TYPE_PROP(GroundWaterBaseProblem, Problem, + Ewoms::GroundWaterProblem); SET_SCALAR_PROP(GroundWaterBaseProblem, LensLowerLeftX, 0.25); SET_SCALAR_PROP(GroundWaterBaseProblem, LensLowerLeftY, 0.25); @@ -81,8 +83,10 @@ SET_SCALAR_PROP(GroundWaterBaseProblem, LensUpperRightZ, 0.75); SET_SCALAR_PROP(GroundWaterBaseProblem, Permeability, 1e-10); SET_SCALAR_PROP(GroundWaterBaseProblem, PermeabilityLens, 1e-12); // Linear solver settings -SET_TYPE_PROP(GroundWaterBaseProblem, LinearSolverWrapper, Ewoms::Linear::SolverWrapperConjugatedGradients ); -SET_TYPE_PROP(GroundWaterBaseProblem, PreconditionerWrapper, Ewoms::Linear::PreconditionerWrapperILU ); +SET_TYPE_PROP(GroundWaterBaseProblem, LinearSolverWrapper, + Ewoms::Linear::SolverWrapperConjugatedGradients); +SET_TYPE_PROP(GroundWaterBaseProblem, PreconditionerWrapper, + Ewoms::Linear::PreconditionerWrapperILU); SET_INT_PROP(GroundWaterBaseProblem, LinearSolverVerbosity, 0); // Enable gravity @@ -113,8 +117,7 @@ namespace Ewoms { * occupied by a rectangular lens of lower permeability. */ template -class GroundWaterProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class GroundWaterProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -135,7 +138,8 @@ class GroundWaterProblem typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GridView::ctype CoordScalar; @@ -148,7 +152,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ GroundWaterProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -166,12 +170,16 @@ public: lensUpperRight_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX); if (dim > 1) - lensUpperRight_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY); + lensUpperRight_[1] + = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY); if (dim > 2) - lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY); + lensUpperRight_[2] + = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY); - intrinsicPerm_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, Permeability)); - intrinsicPermLens_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, PermeabilityLens)); + intrinsicPerm_ + = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, Permeability)); + intrinsicPermLens_ = this->toDimMatrix_( + EWOMS_GET_PARAM(TypeTag, Scalar, PermeabilityLens)); } /*! @@ -181,21 +189,36 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX, "The x-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX, "The x-coordinate of the lens' upper-right corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX, + "The x-coordinate of the lens' lower-left corner " + "[m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX, + "The x-coordinate of the lens' upper-right corner " + "[m]."); if (dimWorld > 1) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY, "The y-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY, "The y-coordinate of the lens' upper-right corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY, + "The y-coordinate of the lens' lower-left " + "corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY, + "The y-coordinate of the lens' upper-right " + "corner [m]."); } if (dimWorld > 2) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ, "The z-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ, "The z-coordinate of the lens' upper-right corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ, + "The z-coordinate of the lens' lower-left " + "corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ, + "The z-coordinate of the lens' upper-right " + "corner [m]."); } - EWOMS_REGISTER_PARAM(TypeTag, Scalar, Permeability, "The intrinsic permeability [m^2] of the ambient material."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, PermeabilityLens, "The intrinsic permeability [m^2] of the lens."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, Permeability, + "The intrinsic permeability [m^2] of the ambient " + "material."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, PermeabilityLens, + "The intrinsic permeability [m^2] of the lens."); } /*! @@ -231,8 +254,12 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const - { return isInLens_(context.pos(spaceIdx, timeIdx))?intrinsicPermLens_:intrinsicPerm_; } + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const + { + return isInLens_(context.pos(spaceIdx, timeIdx)) ? intrinsicPermLens_ + : intrinsicPerm_; + } //! \} /*! @@ -244,8 +271,7 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); @@ -258,7 +284,8 @@ public: else // on upper boundary pressure = 1e5; - Opm::ImmiscibleFluidState fs; + Opm::ImmiscibleFluidState fs; fs.setSaturation(/*phaseIdx=*/0, 1.0); fs.setPressure(/*phaseIdx=*/0, pressure); fs.setTemperature(T); @@ -270,7 +297,6 @@ public: // no flow boundary values.setNoFlow(); } - } //! \} @@ -284,35 +310,34 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { - //const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); - values[pressure0Idx] = 1.0e+5;// + 9.81*1.23*(20-globalPos[dim-1]); + // const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); + values[pressure0Idx] = 1.0e+5; // + 9.81*1.23*(20-globalPos[dim-1]); } /*! * \copydoc VcfvProblem::source */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} private: bool onLowerBoundary_(const GlobalPosition &pos) const - { return pos[dim-1] < eps_; } + { return pos[dim - 1] < eps_; } bool onUpperBoundary_(const GlobalPosition &pos) const - { return pos[dim-1] > this->bboxMax()[dim-1] - eps_; } + { return pos[dim - 1] > this->bboxMax()[dim - 1] - eps_; } bool isInLens_(const GlobalPosition &pos) const { - return - lensLowerLeft_[0] <= pos[0] && pos[0] <= lensUpperRight_[0] && - lensLowerLeft_[1] <= pos[1] && pos[1] <= lensUpperRight_[1]; + return lensLowerLeft_[0] <= pos[0] && pos[0] <= lensUpperRight_[0] + && lensLowerLeft_[1] <= pos[1] && pos[1] <= lensUpperRight_[1]; } GlobalPosition lensLowerLeft_; diff --git a/tests/models/problems/infiltrationproblem.hh b/tests/models/problems/infiltrationproblem.hh index 48b7ef7b8..071cfaf5e 100644 --- a/tests/models/problems/infiltrationproblem.hh +++ b/tests/models/problems/infiltrationproblem.hh @@ -55,12 +55,13 @@ NEW_TYPE_TAG(InfiltrationBaseProblem); SET_TYPE_PROP(InfiltrationBaseProblem, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(InfiltrationBaseProblem, Problem, Ewoms::InfiltrationProblem); +SET_TYPE_PROP(InfiltrationBaseProblem, Problem, + Ewoms::InfiltrationProblem); // Set the fluid system -SET_TYPE_PROP(InfiltrationBaseProblem, - FluidSystem, - Opm::FluidSystems::H2OAirMesitylene); +SET_TYPE_PROP( + InfiltrationBaseProblem, FluidSystem, + Opm::FluidSystems::H2OAirMesitylene); // Enable gravity? SET_BOOL_PROP(InfiltrationBaseProblem, EnableGravity, true); @@ -113,7 +114,8 @@ SET_SCALAR_PROP(InfiltrationBaseProblem, EndTime, 6e3); SET_SCALAR_PROP(InfiltrationBaseProblem, InitialTimeStepSize, 60); // The default DGF file to load -SET_STRING_PROP(InfiltrationBaseProblem, GridFile, "./grids/infiltration_50x3.dgf"); +SET_STRING_PROP(InfiltrationBaseProblem, GridFile, + "./grids/infiltration_50x3.dgf"); } // namespace Properties } // namespace Opm @@ -140,7 +142,7 @@ namespace Ewoms { * (Dirichlet), Top and bottom are Neumann boundaries, all no-flow * except for the small infiltration zone in the upper left part. */ -template +template class InfiltrationProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -151,7 +153,8 @@ class InfiltrationProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; @@ -188,21 +191,21 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ InfiltrationProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()), #else : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()), #endif - , eps_(1e-6) + eps_(1e-6) { temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius FluidSystem::init(/*tempMin=*/temperature_ - 1, /*tempMax=*/temperature_ + 1, /*nTemp=*/3, - /*pressMin=*/0.8*1e5, - /*pressMax=*/3*1e5, + /*pressMin=*/0.8 * 1e5, + /*pressMax=*/3 * 1e5, /*nPress=*/200); // intrinsic permeabilities @@ -264,7 +267,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -278,7 +282,7 @@ public: template Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const { - //const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); + // const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); // if (isFineMaterial_(pos)) // return finePorosity_; // else @@ -290,7 +294,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { return materialParams_; } /*! @@ -299,11 +304,11 @@ public: * In this case, we assume the rock-matrix to be quartz. */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 850. // specific heat capacity [J / (kg K)] - * 2650.; // density of sand [kg/m^3] + return 850. // specific heat capacity [J / (kg K)] + * 2650.; // density of sand [kg/m^3] } //! \} @@ -317,7 +322,8 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -350,7 +356,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { Opm::CompositionalFluidState fs; @@ -368,7 +375,8 @@ public: * everywhere. */ template - void source(RateVector &rate, const Context &context, int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -390,15 +398,17 @@ private: { return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; } template - void initialFluidState_(FluidState &fs, const Context &context, int spaceIdx, int timeIdx) const + void initialFluidState_(FluidState &fs, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition pos = context.pos(spaceIdx, timeIdx); Scalar y = pos[1]; Scalar x = pos[0]; Scalar densityW = 1000.0; - Scalar pc = 9.81 * densityW * (y - (5 - 5e-4*x)); - if (pc < 0.0) pc = 0.0; + Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x)); + if (pc < 0.0) + pc = 0.0; // set pressures const auto &matParams = materialLawParams(context, spaceIdx, timeIdx); @@ -427,57 +437,57 @@ private: if (onLeftBoundary_(pos)) pg += 10e3; MaterialLaw::capillaryPressures(pcAll, matParams, fs); - for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gPhaseIdx])); // set composition of gas phase fs.setMoleFraction(gPhaseIdx, H2OIdx, 1e-6); - fs.setMoleFraction(gPhaseIdx, airIdx, 1 - fs.moleFraction(gPhaseIdx, H2OIdx)); + fs.setMoleFraction(gPhaseIdx, airIdx, + 1 - fs.moleFraction(gPhaseIdx, H2OIdx)); fs.setMoleFraction(gPhaseIdx, NAPLIdx, 0); typedef Opm::ComputeFromReferencePhase CFRP; typename FluidSystem::ParameterCache paramCache; - CFRP::solve(fs, - paramCache, - gPhaseIdx, + CFRP::solve(fs, paramCache, gPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/false); - fs.setMoleFraction(wPhaseIdx, H2OIdx, 1 - fs.moleFraction(wPhaseIdx, H2OIdx)); + fs.setMoleFraction(wPhaseIdx, H2OIdx, + 1 - fs.moleFraction(wPhaseIdx, H2OIdx)); } static Scalar invertPCGW_(Scalar pcIn, const MaterialLawParams &pcParams) { - Scalar lower,upper; + Scalar lower, upper; int k; int maxIt = 50; Scalar bisLimit = 1.; Scalar Sw, pcGW; - lower=0.0; upper=1.0; - for (k=1; k<=25; k++) - { - Sw = 0.5*(upper+lower); + lower = 0.0; + upper = 1.0; + for (k = 1; k <= 25; k++) { + Sw = 0.5 * (upper + lower); pcGW = MaterialLaw::pCGW(pcParams, Sw); - Scalar delta = pcGW-pcIn; - if (delta<0.) delta*=-1.; - if (deltapcIn) lower=Sw; - else upper=Sw; + if (pcGW > pcIn) + lower = Sw; + else + upper = Sw; } return 0; } bool isFineMaterial_(const GlobalPosition &pos) const { - return - 70. <= pos[0] && pos[0] <= 85. && - 7.0 <= pos[1] && pos[1] <= 7.50; + return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; } DimMatrix fineK_; diff --git a/tests/models/problems/lensgridcreator.hh b/tests/models/problems/lensgridcreator.hh index 26730c899..5ef4ae243 100644 --- a/tests/models/problems/lensgridcreator.hh +++ b/tests/models/problems/lensgridcreator.hh @@ -84,10 +84,12 @@ class LensGridCreator public: #if LENS_CUBES typedef Dune::UGGrid Grid; - //typedef Dune::ALUGrid Grid; +// typedef Dune::ALUGrid +// Grid; #else typedef Dune::UGGrid Grid; - //typedef Dune::ALUGrid Grid; +// typedef Dune::ALUGrid +// Grid; #endif /*! @@ -95,16 +97,24 @@ public: */ static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, "The number of global refinements of the grid executed after it was loaded"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, "The size of the domain in x direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, "The number of intervalls in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, + "The number of intervalls in x direction"); if (dim > 1) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, "The size of the domain in y direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, "The number of intervalls in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, + "The number of intervalls in y direction"); } if (dim > 2) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, "The size of the domain in z direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, "The number of intervalls in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, + "The number of intervalls in z direction"); } } @@ -131,15 +141,15 @@ public: Dune::GridFactory factory; if (dim == 3) { - Dune::FieldVector pos; + Dune::FieldVector pos; for (int k = 0; k <= cellRes[0]; k++) { - pos[2] = upperRight[2]*double(k)/cellRes[2]; + pos[2] = upperRight[2] * double(k) / cellRes[2]; for (int j = 0; j <= cellRes[1]; j++) { - pos[1] = upperRight[1]*double(j)/cellRes[1]; + pos[1] = upperRight[1] * double(j) / cellRes[1]; for (int i = 0; i <= cellRes[0]; i++) { - pos[0] = upperRight[0]*double(i)/cellRes[0]; + pos[0] = upperRight[0] * double(i) / cellRes[0]; factory.insertVertex(pos); } } @@ -147,12 +157,12 @@ public: } else { assert(dim == 2); - Dune::FieldVector pos; + Dune::FieldVector pos; for (int j = 0; j <= cellRes[1]; j++) { - pos[1] = upperRight[1]*double(j)/cellRes[1]; + pos[1] = upperRight[1] * double(j) / cellRes[1]; for (int i = 0; i <= cellRes[0]; i++) { - pos[0] = upperRight[0]*double(i)/cellRes[0]; + pos[0] = upperRight[0] * double(i) / cellRes[0]; factory.insertVertex(pos); } } @@ -169,14 +179,14 @@ public: int m = cellRes[0] + 1; int n = cellRes[1] + 1; for (int k = 0; k < cellRes[2]; ++k) { - int i0 = k*m*n + j*m + i; - int i1 = k*m*n + j*m + (i+1); - int i2 = k*m*n + (j+1)*m + i; - int i3 = k*m*n + (j+1)*m + (i+1); - int i4 = (k+1)*m*n + j*m + i; - int i5 = (k+1)*m*n + j*m + (i+1); - int i6 = (k+1)*m*n + (j+1)*m + i; - int i7 = (k+1)*m*n + (j+1)*m + (i+1); + int i0 = k * m * n + j * m + i; + int i1 = k * m * n + j * m + (i + 1); + int i2 = k * m * n + (j + 1) * m + i; + int i3 = k * m * n + (j + 1) * m + (i + 1); + int i4 = (k + 1) * m * n + j * m + i; + int i5 = (k + 1) * m * n + j * m + (i + 1); + int i6 = (k + 1) * m * n + (j + 1) * m + i; + int i7 = (k + 1) * m * n + (j + 1) * m + (i + 1); #if LENS_CUBES v[0] = i0; @@ -187,44 +197,57 @@ public: v[5] = i5; v[6] = i6; v[7] = i7; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::cube, 3), v); #else v[0] = i0; v[1] = i1; v[2] = i2; v[3] = i4; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i4; v[1] = i5; v[2] = i6; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i2; v[1] = i5; v[2] = i4; v[3] = i1; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i2; v[1] = i3; v[2] = i7; v[3] = i5; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i5; v[1] = i7; v[2] = i6; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); v[0] = i1; v[1] = i3; v[2] = i5; v[3] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,3), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 3), + v); #endif } } @@ -232,26 +255,29 @@ public: assert(dim == 2); int m = cellRes[0] + 1; - int i0 = j*m + i; - int i1 = j*m + (i+1); - int i2 = (j+1)*m + i; - int i3 = (j+1)*m + (i+1); + int i0 = j * m + i; + int i1 = j * m + (i + 1); + int i2 = (j + 1) * m + i; + int i3 = (j + 1) * m + (i + 1); #if LENS_CUBES v[0] = i0; v[1] = i1; v[2] = i2; v[3] = i3; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::cube, 2), v); #else v[0] = i0; v[1] = i1; v[2] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 2), v); v[0] = i1; v[1] = i3; v[2] = i2; - factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v); + factory.insertElement( + Dune::GeometryType(Dune::GeometryType::simplex, 2), v); #endif } } @@ -259,7 +285,8 @@ public: grid_ = factory.createGrid(); - unsigned numRefinements = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + unsigned numRefinements + = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); grid_->globalRefine(numRefinements); } @@ -309,16 +336,24 @@ public: */ static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, "The number of global refinements of the grid executed after it was loaded"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, "The size of the domain in x direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, "The number of intervalls in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, GridGlobalRefinements, + "The number of global refinements of the grid " + "executed after it was loaded"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeX, + "The size of the domain in x direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsX, + "The number of intervalls in x direction"); if (dim > 1) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, "The size of the domain in y direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, "The number of intervalls in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeY, + "The size of the domain in y direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsY, + "The number of intervalls in y direction"); } if (dim > 2) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, "The size of the domain in z direction"); - EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, "The number of intervalls in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, DomainSizeZ, + "The size of the domain in z direction"); + EWOMS_REGISTER_PARAM(TypeTag, int, CellsZ, + "The number of intervalls in z direction"); } } @@ -327,7 +362,7 @@ public: */ static void makeGrid() { -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) std::bitset isPeriodic(false); std::array cellRes; #else @@ -351,15 +386,15 @@ public: cellRes[2] = EWOMS_GET_PARAM(TypeTag, int, CellsZ); } - unsigned numRefinements = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); + unsigned numRefinements + = EWOMS_GET_PARAM(TypeTag, unsigned, GridGlobalRefinements); grid_ = new Dune::YaspGrid( #ifdef HAVE_MPI /*mpiCommunicator=*/Dune::MPIHelper::getCommunicator(), #endif /*upperRightCorner=*/upperRight, - /*numCells=*/cellRes, - isPeriodic, + /*numCells=*/cellRes, isPeriodic, /*overlap=*/1); grid_->globalRefine(numRefinements); } diff --git a/tests/models/problems/lensproblem.hh b/tests/models/problems/lensproblem.hh index 726809607..338e0ce1e 100644 --- a/tests/models/problems/lensproblem.hh +++ b/tests/models/problems/lensproblem.hh @@ -68,7 +68,8 @@ NEW_PROP_TAG(LensUpperRightZ); SET_TYPE_PROP(LensBaseProblem, GridCreator, Ewoms::LensGridCreator); // Retrieve the grid type from the grid creator -SET_TYPE_PROP(LensBaseProblem, Grid, typename GET_PROP_TYPE(TypeTag, GridCreator)::Grid); +SET_TYPE_PROP(LensBaseProblem, Grid, + typename GET_PROP_TYPE(TypeTag, GridCreator)::Grid); // Set the problem property SET_TYPE_PROP(LensBaseProblem, Problem, Ewoms::LensProblem); @@ -78,6 +79,7 @@ SET_PROP(LensBaseProblem, WettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; @@ -87,6 +89,7 @@ SET_PROP(LensBaseProblem, NonwettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; @@ -99,7 +102,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations @@ -114,10 +118,10 @@ public: SET_TAG_PROP(LensBaseProblem, LinearSolver, ParallelAmgBackend); // Enable partial reassembly of the jacobian matrix? -//SET_BOOL_PROP(LensBaseProblem, EnablePartialReassemble, true); +// SET_BOOL_PROP(LensBaseProblem, EnablePartialReassemble, true); // Enable reuse of jacobian matrices? -//SET_BOOL_PROP(LensBaseProblem, EnableJacobianRecycling, true); +// SET_BOOL_PROP(LensBaseProblem, EnableJacobianRecycling, true); // Write the solutions of individual newton iterations? SET_BOOL_PROP(LensBaseProblem, NewtonWriteConvergence, false); @@ -177,8 +181,7 @@ namespace Ewoms { * saturation on both sides is zero. */ template -class LensProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class LensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -208,7 +211,8 @@ class LensProblem }; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -223,7 +227,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ LensProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -242,7 +246,8 @@ public: if (dimWorld == 3) { lensLowerLeft_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftZ); - lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightZ); + lensUpperRight_[2] + = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightZ); } // parameters for the Van Genuchten law @@ -271,14 +276,26 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX, "The x-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY, "The y-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX, "The x-coordinate of the lens' upper-right corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY, "The y-coordinate of the lens' upper-right corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX, + "The x-coordinate of the lens' lower-left corner " + "[m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY, + "The y-coordinate of the lens' lower-left corner " + "[m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX, + "The x-coordinate of the lens' upper-right corner " + "[m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY, + "The y-coordinate of the lens' upper-right corner " + "[m]."); if (dimWorld == 3) { - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ, "The z-coordinate of the lens' lower-left corner [m]."); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ, "The z-coordinate of the lens' upper-right corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ, + "The z-coordinate of the lens' lower-left " + "corner [m]."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ, + "The z-coordinate of the lens' upper-right " + "corner [m]."); } }; @@ -291,7 +308,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); @@ -311,7 +329,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); @@ -324,8 +343,7 @@ public: * \copydoc VcfvMultiPhaseProblem::temperature */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return temperature_; } //! \} @@ -356,7 +374,7 @@ public: // Write mass balance information for rank 0 if (this->gridView().comm().rank() == 0) { - std::cout<<"Storage: " << storage << std::endl; + std::cout << "Storage: " << storage << std::endl; } } @@ -371,43 +389,43 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (onLeftBoundary_(pos) || onRightBoundary_(pos)) { // free flow boundary - Scalar densityW = WettingPhase::density(temperature_, /*pressure=*/1e5); + Scalar densityW + = WettingPhase::density(temperature_, /*pressure=*/1e5); Scalar T = temperature(context, spaceIdx, timeIdx); Scalar pw, Sw; // set wetting phase pressure and saturation - if (onLeftBoundary_(pos)) - { + if (onLeftBoundary_(pos)) { Scalar height = this->bboxMax()[1] - this->bboxMin()[1]; Scalar depth = this->bboxMax()[1] - pos[1]; - Scalar alpha = (1 + 1.5/height); + Scalar alpha = (1 + 1.5 / height); // hydrostatic pressure scaled by alpha - pw = 1e5 - alpha*densityW*this->gravity()[1]*depth; + pw = 1e5 - alpha * densityW * this->gravity()[1] * depth; Sw = 1.0; } else { Scalar depth = this->bboxMax()[1] - pos[1]; // hydrostatic pressure - pw = 1e5 - densityW*this->gravity()[1]*depth; + pw = 1e5 - densityW * this->gravity()[1] * depth; Sw = 1.0; } // specify a full fluid state using pw and Sw - const MaterialLawParams &matParams = - this->materialLawParams(context, spaceIdx, timeIdx); + const MaterialLawParams &matParams + = this->materialLawParams(context, spaceIdx, timeIdx); - Opm::ImmiscibleFluidState fs; + Opm::ImmiscibleFluidState fs; fs.setSaturation(wPhaseIdx, Sw); fs.setSaturation(nPhaseIdx, 1 - Sw); fs.setTemperature(T); @@ -415,7 +433,7 @@ public: Scalar pC[numPhases]; MaterialLaw::capillaryPressures(pC, matParams, fs); fs.setPressure(wPhaseIdx, pw); - fs.setPressure(nPhaseIdx, pw + pC[nPhaseIdx] - pC[wPhaseIdx]); + fs.setPressure(nPhaseIdx, pw + pC[nPhaseIdx] - pC[wPhaseIdx]); // impose an freeflow boundary condition values.setFreeFlow(context, spaceIdx, timeIdx, fs); @@ -432,7 +450,6 @@ public: // no flow boundary values.setNoFlow(); } - } //! \} @@ -446,9 +463,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); Scalar depth = this->bboxMax()[1] - pos[1]; @@ -467,11 +483,11 @@ public: Scalar densityW = FluidSystem::density(fs, paramCache, wPhaseIdx); // hydrostatic pressure (assuming incompressibility) - Scalar pw = 1e5 - densityW*this->gravity()[1]*depth; + Scalar pw = 1e5 - densityW * this->gravity()[1] * depth; // calculate the capillary pressure - const MaterialLawParams &matParams = - this->materialLawParams(context, spaceIdx, timeIdx); + const MaterialLawParams &matParams + = this->materialLawParams(context, spaceIdx, timeIdx); Scalar pC[numPhases]; MaterialLaw::capillaryPressures(pC, matParams, fs); @@ -490,9 +506,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -501,7 +516,8 @@ private: bool isInLens_(const GlobalPosition &pos) const { for (int i = 0; i < dim; ++i) { - if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i] + eps_) + if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i] + + eps_) return false; } return true; @@ -522,8 +538,8 @@ private: bool onInlet_(const GlobalPosition &pos) const { Scalar width = this->bboxMax()[0] - this->bboxMin()[0]; - Scalar lambda = (this->bboxMax()[0] - pos[0])/width; - return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0/3.0; + Scalar lambda = (this->bboxMax()[0] - pos[0]) / width; + return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0; } GlobalPosition lensLowerLeft_; diff --git a/tests/models/problems/navierstokestestproblem.hh b/tests/models/problems/navierstokestestproblem.hh index cdd52da60..01a8600c2 100644 --- a/tests/models/problems/navierstokestestproblem.hh +++ b/tests/models/problems/navierstokestestproblem.hh @@ -45,15 +45,18 @@ namespace Properties { NEW_TYPE_TAG(NavierStokesTestProblem, INHERITS_FROM(VcfvNavierStokes)); // Set the grid type -SET_TYPE_PROP(NavierStokesTestProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>); +SET_TYPE_PROP(NavierStokesTestProblem, Grid, + Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>); // Set the property which defines the type of the physical problem -SET_TYPE_PROP(NavierStokesTestProblem, Problem, Ewoms::NavierStokesTestProblem); +SET_TYPE_PROP(NavierStokesTestProblem, Problem, + Ewoms::NavierStokesTestProblem); SET_PROP(NavierStokesTestProblem, Fluid) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::GasPhase > type; }; @@ -71,8 +74,10 @@ SET_SCALAR_PROP(NavierStokesTestProblem, EndTime, 1e-3); SET_SCALAR_PROP(NavierStokesTestProblem, InitialTimeStepSize, 1e-3); // Default grid file to load -SET_STRING_PROP(NavierStokesTestProblem, GridFile, "grids/test_navierstokes.dgf"); -}} +SET_STRING_PROP(NavierStokesTestProblem, GridFile, + "grids/test_navierstokes.dgf"); +} +} namespace Ewoms { /*! @@ -102,7 +107,8 @@ class NavierStokesTestProblem : public StokesProblem typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -125,7 +131,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ NavierStokesTestProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -151,8 +157,7 @@ public: * This problem assumes a constant temperature of 10 degrees Celsius. */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return 273.15 + 10; } //! \} @@ -166,18 +171,19 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { -/* const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); + /* const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); - values.setOutflow(massBalanceIdx); - values.setDirichlet(momentumXIdx); - values.setDirichlet(momentumYIdx); - // set pressure for all vertices at the bottom - if (onLowerBoundary_(pos)) { - values.setDirichlet(massBalanceIdx); - } -*/ + values.setOutflow(massBalanceIdx); + values.setDirichlet(momentumXIdx); + values.setDirichlet(momentumYIdx); + // set pressure for all vertices at the bottom + if (onLowerBoundary_(pos)) { + values.setDirichlet(massBalanceIdx); + } + */ values.setNoFlow(context, spaceIdx, timeIdx); } @@ -192,9 +198,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { initial_(values); } /*! @@ -203,14 +208,16 @@ public: * For this problem, we fix the velocity of upper boundary. */ template - void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const + void constraints(Constraints &constraints, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); if (onUpperBoundary_(pos)) { // lid moves from left to right const Scalar lidVelocity = 1.0; - constraints.setConstraint(momentum0EqIdx, velocity0Idx + 0, lidVelocity); + constraints.setConstraint(momentum0EqIdx, velocity0Idx + 0, + lidVelocity); constraints.setConstraint(momentum0EqIdx + 1, velocity0Idx + 1, 0); constraints.setConstraint(conti0EqIdx, pressureIdx, 1e5); } @@ -220,9 +227,8 @@ public: * \copydoc VcfvProblem::source */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -246,7 +252,7 @@ private: { return globalPos[1] < this->bboxMin()[1] + eps_; } bool onUpperBoundary_(const GlobalPosition &globalPos) const - { return globalPos[1] > this->bboxMax()[1] - eps_; } + { return globalPos[1] > this->bboxMax()[1] - eps_; } Scalar eps_; }; diff --git a/tests/models/problems/obstacleproblem.hh b/tests/models/problems/obstacleproblem.hh index f3fa9ac20..ab99ad5dd 100644 --- a/tests/models/problems/obstacleproblem.hh +++ b/tests/models/problems/obstacleproblem.hh @@ -58,13 +58,10 @@ NEW_TYPE_TAG(ObstacleBaseProblem); SET_TYPE_PROP(ObstacleBaseProblem, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(ObstacleBaseProblem, - Problem, - Ewoms::ObstacleProblem); +SET_TYPE_PROP(ObstacleBaseProblem, Problem, Ewoms::ObstacleProblem); // Set fluid configuration -SET_TYPE_PROP(ObstacleBaseProblem, - FluidSystem, +SET_TYPE_PROP(ObstacleBaseProblem, FluidSystem, Opm::FluidSystems::H2ON2); // Set the material Law @@ -77,7 +74,7 @@ private: typedef Opm::TwoPhaseMaterialTraits - MaterialTraits; + MaterialTraits; typedef Opm::LinearMaterial EffMaterialLaw; @@ -139,8 +136,7 @@ namespace Ewoms { * and the right boundary where a free flow condition is assumed. */ template -class ObstacleProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class ObstacleProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -149,7 +145,8 @@ class ObstacleProblem typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -160,12 +157,9 @@ class ObstacleProblem // Grid and world dimension dim = GridView::dimension, dimWorld = GridView::dimensionworld, - numPhases = GET_PROP_VALUE(TypeTag, NumPhases), - gPhaseIdx = FluidSystem::gPhaseIdx, lPhaseIdx = FluidSystem::lPhaseIdx, - H2OIdx = FluidSystem::H2OIdx, N2Idx = FluidSystem::N2Idx }; @@ -180,7 +174,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ ObstacleProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -218,10 +212,10 @@ public: // parameters for the linear law, i.e. minimum and maximum // pressures - fineMaterialParams_.setPcMinSat(lPhaseIdx,0.0); - fineMaterialParams_.setPcMaxSat(lPhaseIdx,0.0); - coarseMaterialParams_.setPcMinSat(lPhaseIdx,0.0); - coarseMaterialParams_.setPcMaxSat(lPhaseIdx,0.0); + fineMaterialParams_.setPcMinSat(lPhaseIdx, 0.0); + fineMaterialParams_.setPcMaxSat(lPhaseIdx, 0.0); + coarseMaterialParams_.setPcMinSat(lPhaseIdx, 0.0); + coarseMaterialParams_.setPcMaxSat(lPhaseIdx, 0.0); /* // entry pressures for Brooks-Corey @@ -254,13 +248,9 @@ public: this->model().globalPhaseStorage(phaseStorage, phaseIdx); if (this->gridView().comm().rank() == 0) { - std::cout - <<"Storage in " - << FluidSystem::phaseName(phaseIdx) - << "Phase: [" - << phaseStorage - << "]" - << "\n"; + std::cout << "Storage in " << FluidSystem::phaseName(phaseIdx) + << "Phase: [" << phaseStorage << "]" + << "\n"; } } @@ -270,9 +260,8 @@ public: // Write mass balance information for rank 0 if (this->gridView().comm().rank() == 0) { - std::cout - <<"Storage total: [" << storage << "]" - << "\n"; + std::cout << "Storage total: [" << storage << "]" + << "\n"; } } @@ -287,7 +276,8 @@ public: const std::string name() const { std::ostringstream oss; - oss << "obstacle" << "_" << this->model().name(); + oss << "obstacle" + << "_" << this->model().name(); return oss.str(); } @@ -304,7 +294,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { if (isFineMaterial_(context.pos(spaceIdx, timeIdx))) return fineK_; @@ -328,7 +319,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams &materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -344,18 +336,18 @@ public: * medium is granite. */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 790 // specific heat capacity of granite [J / (kg K)] - * 2700; // density of granite [kg/m^3] + return 790 // specific heat capacity of granite [J / (kg K)] + * 2700; // density of granite [kg/m^3] } /*! * \copydoc VcfvMultiPhaseProblem::heatConductionParams */ template - const HeatConductionLawParams& + const HeatConductionLawParams & heatConductionParams(const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -375,7 +367,8 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -398,7 +391,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const auto &matParams = materialLawParams(context, spaceIdx, timeIdx); values.assignMassConservative(outletFluidState_, matParams); @@ -411,9 +405,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = 0.0; } //! \} @@ -424,11 +417,7 @@ private: * fine-permeability region or not. */ bool isFineMaterial_(const GlobalPosition &pos) const - { - return - 10 <= pos[0] && pos[0] <= 20 && - 0 <= pos[1] && pos[1] <= 35; - } + { return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; } bool onInlet_(const GlobalPosition &globalPos) const { @@ -446,12 +435,15 @@ private: void initFluidStates_() { - initFluidState_(inletFluidState_, coarseMaterialParams_, /*isInlet=*/true); - initFluidState_(outletFluidState_, coarseMaterialParams_, /*isInlet=*/false); + initFluidState_(inletFluidState_, coarseMaterialParams_, + /*isInlet=*/true); + initFluidState_(outletFluidState_, coarseMaterialParams_, + /*isInlet=*/false); } template - void initFluidState_(FluidState &fs, const MaterialLawParams &matParams, bool isInlet) + void initFluidState_(FluidState &fs, const MaterialLawParams &matParams, + bool isInlet) { int refPhaseIdx; int otherPhaseIdx; @@ -496,18 +488,16 @@ private: // calulate the capillary pressure PhaseVector pC; MaterialLaw::capillaryPressures(pC, matParams, fs); - fs.setPressure(otherPhaseIdx, - fs.pressure(refPhaseIdx) - + (pC[otherPhaseIdx] - pC[refPhaseIdx])); + fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx) + + (pC[otherPhaseIdx] - pC[refPhaseIdx])); // make the fluid state consistent with local thermodynamic // equilibrium - typedef Opm::ComputeFromReferencePhase ComputeFromReferencePhase; + typedef Opm::ComputeFromReferencePhase + ComputeFromReferencePhase; typename FluidSystem::ParameterCache paramCache; - ComputeFromReferencePhase::solve(fs, - paramCache, - refPhaseIdx, + ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/false); } @@ -517,8 +507,9 @@ private: Scalar lambdaWater = 0.6; Scalar lambdaGranite = 2.8; - Scalar lambdaWet = std::pow(lambdaGranite, (1-poro)) * std::pow(lambdaWater, poro); - Scalar lambdaDry = std::pow(lambdaGranite, (1-poro)); + Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro)) + * std::pow(lambdaWater, poro); + Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro)); params.setFullySaturatedLambda(gPhaseIdx, lambdaDry); params.setFullySaturatedLambda(lPhaseIdx, lambdaWet); diff --git a/tests/models/problems/outflowproblem.hh b/tests/models/problems/outflowproblem.hh index 8f5700bc8..e7694ecfd 100644 --- a/tests/models/problems/outflowproblem.hh +++ b/tests/models/problems/outflowproblem.hh @@ -51,8 +51,10 @@ SET_TYPE_PROP(OutflowBaseProblem, Problem, Ewoms::OutflowProblem); // Set fluid system SET_PROP(OutflowBaseProblem, FluidSystem) -{ private: +{ +private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: // Two-component single phase fluid system typedef Opm::FluidSystems::H2ON2LiquidPhase type; @@ -94,8 +96,7 @@ namespace Ewoms { * used. */ template -class OutflowProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class OutflowProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -103,7 +104,8 @@ class OutflowProblem typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -129,17 +131,18 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ OutflowProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()), #else : ParentType(timeManager, - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()), #endif - , eps_(1e-6) + eps_(1e-6) { temperature_ = 273.15 + 20; - FluidSystem::init(/*minT=*/temperature_ - 1, /*maxT=*/temperature_ + 2, /*numT=*/3, + FluidSystem::init(/*minT=*/temperature_ - 1, /*maxT=*/temperature_ + 2, + /*numT=*/3, /*minp=*/0.8e5, /*maxp=*/2.5e5, /*nump=*/500); // set parameters of porous medium @@ -174,7 +177,8 @@ public: * This problem uses a constant intrinsic permeability. */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { return perm_; } /*! @@ -216,14 +220,14 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); if (onLeftBoundary_(globalPos)) { - Opm::CompositionalFluidState fs; + Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); fs.setPressure(/*phaseIdx=*/0, fs.pressure(/*phaseIdx=*/0) + 1e5); @@ -235,7 +239,8 @@ public: values.setFreeFlow(context, spaceIdx, timeIdx, fs); } else if (onRightBoundary_(globalPos)) { - Opm::CompositionalFluidState fs; + Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); // impose an outflow boundary condition @@ -257,7 +262,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); @@ -272,9 +278,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -287,14 +292,15 @@ private: { return pos[0] > this->bboxMax()[0] - eps_; } template - void initialFluidState_(FluidState &fs, - const Context &context, + void initialFluidState_(FluidState &fs, const Context &context, int spaceIdx, int timeIdx) const { Scalar T = temperature(context, spaceIdx, timeIdx); - //Scalar rho = FluidSystem::H2O::liquidDensity(T, /*pressure=*/1.5e5); - //Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] - this->bboxMax()[dim - 1]; - //Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] - this->bboxMax()[dim - 1]; + // Scalar rho = FluidSystem::H2O::liquidDensity(T, /*pressure=*/1.5e5); + // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] - + // this->bboxMax()[dim - 1]; + // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] - + // this->bboxMax()[dim - 1]; fs.setSaturation(/*phaseIdx=*/0, 1.0); fs.setPressure(/*phaseIdx=*/0, 1e5 /* + rho*z */); diff --git a/tests/models/problems/powerinjectionproblem.hh b/tests/models/problems/powerinjectionproblem.hh index 097c7e276..4224574ea 100644 --- a/tests/models/problems/powerinjectionproblem.hh +++ b/tests/models/problems/powerinjectionproblem.hh @@ -57,16 +57,19 @@ NEW_TYPE_TAG(PowerInjectionBaseProblem); SET_TYPE_PROP(PowerInjectionBaseProblem, Grid, Dune::YaspGrid); // set the GridCreator property -SET_TYPE_PROP(PowerInjectionBaseProblem, GridCreator, Ewoms::CubeGridCreator); +SET_TYPE_PROP(PowerInjectionBaseProblem, GridCreator, + Ewoms::CubeGridCreator); // Set the problem property -SET_TYPE_PROP(PowerInjectionBaseProblem, Problem, Ewoms::PowerInjectionProblem); +SET_TYPE_PROP(PowerInjectionBaseProblem, Problem, + Ewoms::PowerInjectionProblem); // Set the wetting phase SET_PROP(PowerInjectionBaseProblem, WettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; @@ -76,6 +79,7 @@ SET_PROP(PowerInjectionBaseProblem, NonwettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::GasPhase > type; }; @@ -88,7 +92,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations @@ -136,8 +141,7 @@ namespace Ewoms { * Systems, University of Stuttgart, 2011 */ template -class PowerInjectionProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class PowerInjectionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -166,7 +170,8 @@ class PowerInjectionProblem }; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -181,7 +186,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ PowerInjectionProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -236,7 +241,7 @@ public: // Write mass balance information for rank 0 if (this->gridView().comm().rank() == 0) { - std::cout<<"Storage: " << storage << std::endl; + std::cout << "Storage: " << storage << std::endl; } } //! \} @@ -250,14 +255,16 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { return K_; } /*! * \copydoc VcfvForchheimerBaseProblem::ergunCoefficient */ template - Scalar ergunCoefficient(const Context &context, int spaceIdx, int timeIdx) const + Scalar ergunCoefficient(const Context &context, int spaceIdx, + int timeIdx) const { return 0.3866; } /*! @@ -271,15 +278,15 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { return materialParams_; } /*! * \copydoc VcfvMultiPhaseProblem::temperature */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return temperature_; } //! \} @@ -296,8 +303,7 @@ public: * left and a free-flow boundary on the right. */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -310,11 +316,10 @@ public: // impose a forced flow boundary values.setMassRate(massRate); } - else { + else { // free flow boundary with initial condition on the right values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_); } - } //! \} @@ -328,9 +333,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { // assign the primary variables values.assignNaive(initialFluidState_); @@ -343,9 +347,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} diff --git a/tests/models/problems/reservoirproblem.hh b/tests/models/problems/reservoirproblem.hh index 7aab922f3..cc43bb582 100644 --- a/tests/models/problems/reservoirproblem.hh +++ b/tests/models/problems/reservoirproblem.hh @@ -71,10 +71,11 @@ private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef Opm::ThreePhaseMaterialTraits Traits; + typedef Opm:: + ThreePhaseMaterialTraits Traits; public: typedef Opm::LinearMaterial type; @@ -132,8 +133,7 @@ namespace Ewoms { * which is 50% above the reservoir pressure. */ template -class ReservoirProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class ReservoirProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -141,32 +141,28 @@ class ReservoirProblem typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - enum { - // Grid and world dimension - dim = GridView::dimension, - dimWorld = GridView::dimensionworld - }; + // Grid and world dimension + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; // copy some indices for convenience - enum { - numPhases = FluidSystem::numPhases, - numComponents = FluidSystem::numComponents, - - gPhaseIdx = FluidSystem::gPhaseIdx, - oPhaseIdx = FluidSystem::oPhaseIdx, - wPhaseIdx = FluidSystem::wPhaseIdx, - - gCompIdx = FluidSystem::gCompIdx, - oCompIdx = FluidSystem::oCompIdx, - wCompIdx = FluidSystem::wCompIdx - }; + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; + enum { oPhaseIdx = FluidSystem::oPhaseIdx }; + enum { wPhaseIdx = FluidSystem::wPhaseIdx }; + enum { gCompIdx = FluidSystem::gCompIdx }; + enum { oCompIdx = FluidSystem::oCompIdx }; + enum { wCompIdx = FluidSystem::wCompIdx }; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GET_PROP_TYPE(TypeTag, BlackOilFluidState) BlackOilFluidState; + typedef typename GET_PROP_TYPE(TypeTag, + BlackOilFluidState) BlackOilFluidState; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -180,7 +176,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ ReservoirProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -195,71 +191,66 @@ public: name_ = EWOMS_GET_PARAM(TypeTag, std::string, SimulationName); FluidSystem::initBegin(); - std::vector > Bg = { - { 1.013529e+05, 9.998450e-01 }, - { 2.757903e+06, 3.075500e-02 }, - { 5.515806e+06, 1.537947e-02 }, - { 8.273709e+06, 1.021742e-02 }, - { 1.103161e+07, 7.662783e-03 }, - { 1.378951e+07, 6.151899e-03 }, - { 1.654742e+07, 5.108709e-03 }, - { 1.930532e+07, 4.378814e-03 }, - { 2.206322e+07, 3.857780e-03 }, - { 2.482113e+07, 3.388401e-03 }, - { 2.757903e+07, 3.049842e-03 } - }; - std::vector > Bo = { - { 1.013529e+05, 1.000000e+00 }, - { 2.757903e+06, 1.012000e+00 }, - { 5.515806e+06, 1.025500e+00 }, - { 8.273709e+06, 1.038000e+00 }, - { 1.103161e+07, 1.051000e+00 }, - { 1.378951e+07, 1.063000e+00 }, - { 1.654742e+07, 1.075000e+00 }, - { 1.930532e+07, 1.087000e+00 }, - { 2.206322e+07, 1.098500e+00 }, - { 2.482113e+07, 1.110000e+00 }, - { 2.757903e+07, 1.120000e+00 } - }; - std::vector > Rs = { - { 1.013529e+05, 0.000000e+00 }, - { 2.757903e+06, 2.938776e+01 }, - { 5.515806e+06, 5.966605e+01 }, - { 8.273709e+06, 8.905380e+01 }, - { 1.103161e+07, 1.184416e+02 }, - { 1.378951e+07, 1.474731e+02 }, - { 1.654742e+07, 1.754360e+02 }, - { 1.930532e+07, 2.012616e+02 }, - { 2.206322e+07, 2.261967e+02 }, - { 2.482113e+07, 2.475696e+02 }, - { 2.757903e+07, 2.671614e+02 } - }; - std::vector > muo = { - { 1.013529e+05, 1.200000e-03 }, - { 2.757903e+06, 1.170000e-03 }, - { 5.515806e+06, 1.140000e-03 }, - { 8.273709e+06, 1.110000e-03 }, - { 1.103161e+07, 1.080000e-03 }, - { 1.378951e+07, 1.060000e-03 }, - { 1.654742e+07, 1.030000e-03 }, - { 1.930532e+07, 1.000000e-03 }, - { 2.206322e+07, 9.800000e-04 }, - { 2.482113e+07, 9.500000e-04 }, - { 2.757903e+07, 9.400000e-04 } - }; - std::vector > mug = { - { 1.013529e+05, 1.250000e-05 }, - { 2.757903e+06, 1.300000e-05 }, - { 5.515806e+06, 1.350000e-05 }, - { 8.273709e+06, 1.400000e-05 }, - { 1.103161e+07, 1.450000e-05 }, - { 1.378951e+07, 1.500000e-05 }, - { 1.654742e+07, 1.550000e-05 }, - { 1.930532e+07, 1.600000e-05 }, - { 2.206322e+07, 1.650000e-05 }, - { 2.482113e+07, 1.700000e-05 }, - { 2.757903e+07, 1.750000e-05 }, - }; + std::vector > Bg + = { { 1.013529e+05, 9.998450e-01 }, + { 2.757903e+06, 3.075500e-02 }, + { 5.515806e+06, 1.537947e-02 }, + { 8.273709e+06, 1.021742e-02 }, + { 1.103161e+07, 7.662783e-03 }, + { 1.378951e+07, 6.151899e-03 }, + { 1.654742e+07, 5.108709e-03 }, + { 1.930532e+07, 4.378814e-03 }, + { 2.206322e+07, 3.857780e-03 }, + { 2.482113e+07, 3.388401e-03 }, + { 2.757903e+07, 3.049842e-03 } }; + std::vector > Bo + = { { 1.013529e+05, 1.000000e+00 }, + { 2.757903e+06, 1.012000e+00 }, + { 5.515806e+06, 1.025500e+00 }, + { 8.273709e+06, 1.038000e+00 }, + { 1.103161e+07, 1.051000e+00 }, + { 1.378951e+07, 1.063000e+00 }, + { 1.654742e+07, 1.075000e+00 }, + { 1.930532e+07, 1.087000e+00 }, + { 2.206322e+07, 1.098500e+00 }, + { 2.482113e+07, 1.110000e+00 }, + { 2.757903e+07, 1.120000e+00 } }; + std::vector > Rs + = { { 1.013529e+05, 0.000000e+00 }, + { 2.757903e+06, 2.938776e+01 }, + { 5.515806e+06, 5.966605e+01 }, + { 8.273709e+06, 8.905380e+01 }, + { 1.103161e+07, 1.184416e+02 }, + { 1.378951e+07, 1.474731e+02 }, + { 1.654742e+07, 1.754360e+02 }, + { 1.930532e+07, 2.012616e+02 }, + { 2.206322e+07, 2.261967e+02 }, + { 2.482113e+07, 2.475696e+02 }, + { 2.757903e+07, 2.671614e+02 } }; + std::vector > muo + = { { 1.013529e+05, 1.200000e-03 }, + { 2.757903e+06, 1.170000e-03 }, + { 5.515806e+06, 1.140000e-03 }, + { 8.273709e+06, 1.110000e-03 }, + { 1.103161e+07, 1.080000e-03 }, + { 1.378951e+07, 1.060000e-03 }, + { 1.654742e+07, 1.030000e-03 }, + { 1.930532e+07, 1.000000e-03 }, + { 2.206322e+07, 9.800000e-04 }, + { 2.482113e+07, 9.500000e-04 }, + { 2.757903e+07, 9.400000e-04 } }; + std::vector > mug + = { { 1.013529e+05, 1.250000e-05 }, + { 2.757903e+06, 1.300000e-05 }, + { 5.515806e+06, 1.350000e-05 }, + { 8.273709e+06, 1.400000e-05 }, + { 1.103161e+07, 1.450000e-05 }, + { 1.378951e+07, 1.500000e-05 }, + { 1.654742e+07, 1.550000e-05 }, + { 1.930532e+07, 1.600000e-05 }, + { 2.206322e+07, 1.650000e-05 }, + { 2.482113e+07, 1.700000e-05 }, + { 2.757903e+07, 1.750000e-05 }, }; FluidSystem::setGasFormationVolumeFactor(Bg); FluidSystem::setOilFormationVolumeFactor(Bo); FluidSystem::setGasDissolutionFactor(Rs); @@ -270,7 +261,7 @@ public: FluidSystem::setSurfaceDensities(/*oil=*/720.51, /*water=*/1009.32, /*gas=*/1.1245); - for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) FluidSystem::setReferenceVolumeFactor(phaseIdx, 1.0); FluidSystem::initEnd(); @@ -285,7 +276,7 @@ public: finePorosity_ = 0.2; coarsePorosity_ = 0.3; - for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { fineMaterialParams_.setPcMinSat(phaseIdx, 0.0); fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0); @@ -307,9 +298,13 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature, "The temperature [K] in the reservoir"); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth, "The maximum depth [m] of the reservoir"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName, "The name of the simulation used for the output files"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature, + "The temperature [K] in the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth, + "The maximum depth [m] of the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName, + "The name of the simulation used for the output " + "files"); } /*! @@ -319,7 +314,8 @@ public: * above one with low permeability. */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -343,7 +339,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -388,8 +385,7 @@ public: * extraction and production wells, so all boundaries are no-flow. */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { // no flow on top and bottom @@ -410,7 +406,8 @@ public: * the whole domain. */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { ////// // set the primary variables @@ -429,23 +426,19 @@ public: * reservoir. */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); Scalar x = pos[0] - this->bboxMin()[0]; - Scalar y = pos[dim-1] - this->bboxMin()[dim-1]; - Scalar height = this->bboxMax()[dim-1] - this->bboxMin()[dim-1]; + Scalar y = pos[dim - 1] - this->bboxMin()[dim - 1]; + Scalar height = this->bboxMax()[dim - 1] - this->bboxMin()[dim - 1]; Scalar width = this->bboxMax()[0] - this->bboxMin()[0]; - if ((onLeftBoundary_(pos) - || onRightBoundary_(pos)) - && y < height/2) - { + if ((onLeftBoundary_(pos) || onRightBoundary_(pos)) && y < height / 2) { // injectors auto fs = initialFluidState_; - Scalar pInj = pReservoir_*1.5; + Scalar pInj = pReservoir_ * 1.5; fs.setPressure(wPhaseIdx, pInj); fs.setPressure(oPhaseIdx, pInj); fs.setPressure(gPhaseIdx, pInj); @@ -461,21 +454,20 @@ public: // set the composition of the oil phase to the initial // composition for (int compIdx = 0; compIdx < numComponents; ++compIdx) - fs.setMoleFraction(oPhaseIdx, - compIdx, - initialFluidState_.moleFraction(oPhaseIdx, compIdx)); + fs.setMoleFraction(oPhaseIdx, compIdx, + initialFluidState_.moleFraction(oPhaseIdx, + compIdx)); fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0); constraints.setAllConstraint(); constraints.assignNaive(fs); } - else if (width/2 - 1 < x && x < width/2 + 1 && y > height/2) - { + else if (width / 2 - 1 < x && x < width / 2 + 1 && y > height / 2) { // producer auto fs = initialFluidState_; - Scalar pProd = pReservoir_/1.5; + Scalar pProd = pReservoir_ / 1.5; fs.setPressure(wPhaseIdx, pProd); fs.setPressure(oPhaseIdx, pProd); fs.setPressure(gPhaseIdx, pProd); @@ -486,9 +478,9 @@ public: // set the compositions to the initial composition for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) for (int compIdx = 0; compIdx < numComponents; ++compIdx) - fs.setMoleFraction(phaseIdx, - compIdx, - initialFluidState_.moleFraction(phaseIdx, compIdx)); + fs.setMoleFraction(phaseIdx, compIdx, + initialFluidState_.moleFraction(phaseIdx, + compIdx)); constraints.setAllConstraint(); constraints.assignNaive(fs); @@ -501,9 +493,8 @@ public: * For this problem, the source term of all components is 0 everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -558,18 +549,18 @@ private: Scalar pSat = pReservoir_; // the saturation pressure of the oil Scalar Bo = FluidSystem::oilFormationVolumeFactor(pSat); Scalar Rs = FluidSystem::gasDissolutionFactor(pSat); - Scalar rhoo = FluidSystem::surfaceDensity(oPhaseIdx)/Bo; + Scalar rhoo = FluidSystem::surfaceDensity(oPhaseIdx) / Bo; Scalar rhogref = FluidSystem::surfaceDensity(gPhaseIdx); // calculate composition of oil phase in terms of mass // fractions. - Scalar XoG = Rs*rhogref / rhoo; + Scalar XoG = Rs * rhogref / rhoo; // convert mass to mole fractions Scalar MG = FluidSystem::molarMass(gCompIdx); Scalar MO = FluidSystem::molarMass(oCompIdx); - Scalar xoG = XoG*MO/((MO - MG)*XoG + MG); + Scalar xoG = XoG * MO / ((MO - MG) * XoG + MG); Scalar xoO = 1 - xoG; // finally set the oil-phase composition @@ -587,7 +578,7 @@ private: { return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); } bool isFineMaterial_(const GlobalPosition &pos) const - { return pos[dim-1] > layerBottom_; } + { return pos[dim - 1] > layerBottom_; } DimMatrix fineK_; DimMatrix coarseK_; @@ -606,7 +597,7 @@ private: Scalar maxDepth_; Scalar eps_; - std::string name_ ; + std::string name_; }; } // namespace Ewoms diff --git a/tests/models/problems/richardslensproblem.hh b/tests/models/problems/richardslensproblem.hh index 4da9070d1..f5c7dd8fa 100644 --- a/tests/models/problems/richardslensproblem.hh +++ b/tests/models/problems/richardslensproblem.hh @@ -53,13 +53,16 @@ SET_TYPE_PROP(RichardsLensProblem, Grid, Dune::YaspGrid<2>); // Set the physical problem to be solved SET_PROP(RichardsLensProblem, Problem) -{ typedef Ewoms::RichardsLensProblem type; }; +{ + typedef Ewoms::RichardsLensProblem type; +}; // Set the wetting phase SET_PROP(RichardsLensProblem, WettingPhase) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::LiquidPhase > type; }; @@ -72,7 +75,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::nPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations @@ -112,7 +116,8 @@ SET_SCALAR_PROP(RichardsLensProblem, EndTime, 3000); SET_SCALAR_PROP(RichardsLensProblem, InitialTimeStepSize, 100); // The default DGF file to load -SET_STRING_PROP(RichardsLensProblem, GridFile, "./grids/richardslens_24x16.dgf"); +SET_STRING_PROP(RichardsLensProblem, GridFile, + "./grids/richardslens_24x16.dgf"); } // namespace Properties } // namespace Opm @@ -134,14 +139,14 @@ namespace Ewoms { * instead of a \c DNAPL infiltrates from the top. */ template -class RichardsLensProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class RichardsLensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; @@ -152,17 +157,15 @@ class RichardsLensProblem // copy some indices for convenience pressureWIdx = Indices::pressureWIdx, contiWEqIdx = Indices::contiWEqIdx, - wPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex), nPhaseIdx = 1 - wPhaseIdx, - numPhases = FluidSystem::numPhases, // Grid and world dimension dimWorld = GridView::dimensionworld }; - //get the material law from the property system + // get the material law from the property system typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; //! The parameters of the material law to be used typedef typename MaterialLaw::Params MaterialLawParams; @@ -177,13 +180,14 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ RichardsLensProblem(TimeManager &timeManager) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()), #else - GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()) + : ParentType(timeManager, + GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafView()), #endif - , pnRef_(1e5) + pnRef_(1e5) { eps_ = 3e-6; pnRef_ = 1e5; @@ -206,10 +210,10 @@ public: // parameters for the linear law // minimum and maximum pressures -// lensMaterialParams_.setEntryPC(0); -// outerMaterialParams_.setEntryPC(0); -// lensMaterialParams_.setMaxPC(0); -// outerMaterialParams_.setMaxPC(0); + // lensMaterialParams_.setEntryPC(0); + // outerMaterialParams_.setEntryPC(0); + // lensMaterialParams_.setMaxPC(0); + // outerMaterialParams_.setMaxPC(0); lensK_ = this->toDimMatrix_(1e-12); outerK_ = this->toDimMatrix_(5e-12); @@ -237,7 +241,8 @@ public: * \copydoc VcfvMultiPhaseProblem::intrinsicPermeability */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isInLens_(pos)) @@ -256,7 +261,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); if (isInLens_(pos)) @@ -270,7 +276,8 @@ public: * \copydetails Doxygen::contextParams */ template - Scalar referencePressure(const Context &context, int spaceIdx, int timeIdx) const + Scalar referencePressure(const Context &context, int spaceIdx, + int timeIdx) const { return pnRef_; } //! \} @@ -284,14 +291,14 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); - if (onLeftBoundary_(pos) || - onRightBoundary_(pos)) - { - const auto &materialParams = this->materialLawParams(context, spaceIdx, timeIdx); + if (onLeftBoundary_(pos) || onRightBoundary_(pos)) { + const auto &materialParams + = this->materialLawParams(context, spaceIdx, timeIdx); Scalar Sw = 0.0; Opm::ImmiscibleFluidState fs; @@ -328,11 +335,11 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { - const auto &materialParams = this->materialLawParams(context, spaceIdx, timeIdx); + const auto &materialParams + = this->materialLawParams(context, spaceIdx, timeIdx); Scalar Sw = 0.0; Opm::ImmiscibleFluidState fs; @@ -351,7 +358,8 @@ public: * everywhere. */ template - void source(RateVector &rate, const Context &context, int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } //! \} @@ -372,8 +380,8 @@ private: bool onInlet_(const GlobalPosition &pos) const { Scalar width = this->bboxMax()[0] - this->bboxMin()[0]; - Scalar lambda = (this->bboxMax()[0] - pos[0])/width; - return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0/3.0; + Scalar lambda = (this->bboxMax()[0] - pos[0]) / width; + return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0; } bool isInLens_(const GlobalPosition &pos) const diff --git a/tests/models/problems/stokes2ctestproblem.hh b/tests/models/problems/stokes2ctestproblem.hh index e0afed1a7..a55022635 100644 --- a/tests/models/problems/stokes2ctestproblem.hh +++ b/tests/models/problems/stokes2ctestproblem.hh @@ -51,13 +51,11 @@ SET_TYPE_PROP(Stokes2cTestProblem, Grid, Dune::YaspGrid<2>); SET_TYPE_PROP(Stokes2cTestProblem, Problem, Ewoms::Stokes2cTestProblem); //! Select the fluid system -SET_TYPE_PROP(Stokes2cTestProblem, - FluidSystem, +SET_TYPE_PROP(Stokes2cTestProblem, FluidSystem, Opm::FluidSystems::H2OAir); //! Select the phase to be considered -SET_INT_PROP(Stokes2cTestProblem, - StokesPhaseIndex, +SET_INT_PROP(Stokes2cTestProblem, StokesPhaseIndex, GET_PROP_TYPE(TypeTag, FluidSystem)::gPhaseIdx); // Disable gravity @@ -90,15 +88,15 @@ namespace Ewoms { * exhibiting slightly higher humitiy than the ones on the right. */ template -class Stokes2cTestProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class Stokes2cTestProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; @@ -110,11 +108,9 @@ class Stokes2cTestProblem // copy some indices for convenience conti0EqIdx = Indices::conti0EqIdx, momentum0EqIdx = Indices::momentum0EqIdx, - velocity0Idx = Indices::velocity0Idx, moleFrac1Idx = Indices::moleFrac1Idx, pressureIdx = Indices::pressureIdx, - H2OIdx = FluidSystem::H2OIdx, AirIdx = FluidSystem::AirIdx }; @@ -127,7 +123,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ Stokes2cTestProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -158,8 +154,7 @@ public: * This problem assumes a temperature of 10 degrees Celsius. */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return 273.15 + 10; /* -> 10 deg C */ } // \} @@ -177,13 +172,14 @@ public: * upper edge. */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (onLowerBoundary_(pos)) values.setOutFlow(context, spaceIdx, timeIdx); - else if(onUpperBoundary_(pos)) { + else if (onUpperBoundary_(pos)) { // upper boundary is constraint! values = 0.0; } @@ -208,18 +204,19 @@ public: * 0.5% is set. */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &globalPos = context.pos(spaceIdx, timeIdx); values = 0.0; - //parabolic profile + // parabolic profile const Scalar v1 = 1.0; - values[velocity0Idx + 1] = - - v1*(globalPos[0] - this->bboxMin()[0])*(this->bboxMax()[0] - globalPos[0]) - / (0.25*(this->bboxMax()[0] - this->bboxMin()[0])*(this->bboxMax()[0] - this->bboxMin()[0])); + values[velocity0Idx + 1] + = -v1 * (globalPos[0] - this->bboxMin()[0]) + * (this->bboxMax()[0] - globalPos[0]) + / (0.25 * (this->bboxMax()[0] - this->bboxMin()[0]) + * (this->bboxMax()[0] - this->bboxMin()[0])); Scalar moleFrac[numComponents]; if (onUpperBoundary_(globalPos)) @@ -240,9 +237,8 @@ public: * is 0 everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } /*! @@ -252,8 +248,7 @@ public: * initial conditions. */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -262,8 +257,11 @@ public: PrimaryVariables initCond; initial(initCond, context, spaceIdx, timeIdx); - constraints.setConstraint(pressureIdx, conti0EqIdx, initCond[pressureIdx]);; - constraints.setConstraint(moleFrac1Idx, conti0EqIdx + 1, initCond[moleFrac1Idx]); + constraints.setConstraint(pressureIdx, conti0EqIdx, + initCond[pressureIdx]); + ; + constraints.setConstraint(moleFrac1Idx, conti0EqIdx + 1, + initCond[moleFrac1Idx]); for (int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) constraints.setConstraint(velocity0Idx + axisIdx, momentum0EqIdx + axisIdx, diff --git a/tests/models/problems/stokesnitestproblem.hh b/tests/models/problems/stokesnitestproblem.hh index 9eef9f45e..d6acab4ab 100644 --- a/tests/models/problems/stokesnitestproblem.hh +++ b/tests/models/problems/stokesnitestproblem.hh @@ -48,13 +48,11 @@ SET_TYPE_PROP(StokesNITestProblem, Grid, Dune::YaspGrid<2>); SET_TYPE_PROP(StokesNITestProblem, Problem, Ewoms::StokesNITestProblem); //! Select the fluid system -SET_TYPE_PROP(StokesNITestProblem, - FluidSystem, +SET_TYPE_PROP(StokesNITestProblem, FluidSystem, Opm::FluidSystems::H2OAir); //! Select the phase to be considered -SET_INT_PROP(StokesNITestProblem, - StokesPhaseIndex, +SET_INT_PROP(StokesNITestProblem, StokesPhaseIndex, GET_PROP_TYPE(TypeTag, FluidSystem)::gPhaseIdx); // Enable gravity @@ -90,8 +88,7 @@ namespace Ewoms { * conditions. */ template -class StokesNITestProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class StokesNITestProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; @@ -100,14 +97,14 @@ class StokesNITestProblem typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; enum { // Number of equations and grid dimension numEq = GET_PROP_VALUE(TypeTag, NumEq), - dimWorld = GridView::dimensionworld - }; + dimWorld = GridView::dimensionworld }; enum { // primary variable indices pressureIdx = Indices::pressureIdx, @@ -133,7 +130,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ StokesNITestProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -170,16 +167,14 @@ public: * \copydoc VcfvProblem::boundary */ template - void boundary(BoundaryRateVector &values, - const Context &context, - int spaceIdx, - int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (onUpperBoundary_(pos)) values.setOutFlow(context, spaceIdx, timeIdx); - else if(onLowerBoundary_(pos)) { + else if (onLowerBoundary_(pos)) { // lower boundary is constraint! values = 0.0; } @@ -200,9 +195,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -224,18 +218,18 @@ public: // parabolic velocity profile const Scalar maxVelocity = 1.0; - Scalar a = - 4*maxVelocity/(width*width); - Scalar b = - a*width; + Scalar a = -4 * maxVelocity / (width * width); + Scalar b = -a * width; Scalar c = 0; DimVector velocity(0.0); - velocity[1] = a * x*x + b * x + c; + velocity[1] = a * x * x + b * x + c; // hydrostatic pressure Scalar rho = 1.189; - Scalar pressure = 1e5 - rho*this->gravity()[1]*y; + Scalar pressure = 1e5 - rho * this->gravity()[1] * y; - for (int axisIdx = 0; axisIdx < dimWorld; ++ axisIdx) + for (int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) values[velocity0Idx + axisIdx] = velocity[axisIdx]; values[pressureIdx] = pressure; @@ -250,9 +244,8 @@ public: * is 0 everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } /*! @@ -262,20 +255,23 @@ public: * adjacent to the inlet. */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); - if (onLowerBoundary_(pos) || onUpperBoundary_(pos)) - { + if (onLowerBoundary_(pos) || onUpperBoundary_(pos)) { PrimaryVariables initCond; initial(initCond, context, spaceIdx, timeIdx); - constraints.setConstraint(temperatureIdx, energyEqIdx, initCond[temperatureIdx]);; - constraints.setConstraint(pressureIdx, conti0EqIdx, initCond[pressureIdx]); - constraints.setConstraint(moleFrac1Idx, conti0EqIdx+1, initCond[moleFrac1Idx]);; + constraints.setConstraint(temperatureIdx, energyEqIdx, + initCond[temperatureIdx]); + ; + constraints.setConstraint(pressureIdx, conti0EqIdx, + initCond[pressureIdx]); + constraints.setConstraint(moleFrac1Idx, conti0EqIdx + 1, + initCond[moleFrac1Idx]); + ; for (int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) constraints.setConstraint(velocity0Idx + axisIdx, momentum0EqIdx + axisIdx, @@ -300,18 +296,12 @@ private: bool onBoundary_(const GlobalPosition &pos) const { - return onLeftBoundary_(pos) - || onRightBoundary_(pos) - || onLowerBoundary_(pos) - || onUpperBoundary_(pos); + return onLeftBoundary_(pos) || onRightBoundary_(pos) + || onLowerBoundary_(pos) || onUpperBoundary_(pos); } bool inLens_(const GlobalPosition &pos) const - { - return - pos[0]<0.75 && pos[0]>0.25 && - pos[1]<0.75 && pos[1]>0.25; - } + { return pos[0] < 0.75 && pos[0] > 0.25 && pos[1] < 0.75 && pos[1] > 0.25; } Scalar eps_; }; diff --git a/tests/models/problems/stokestestproblem.hh b/tests/models/problems/stokestestproblem.hh index a4e05317c..9661db5ab 100644 --- a/tests/models/problems/stokestestproblem.hh +++ b/tests/models/problems/stokestestproblem.hh @@ -56,6 +56,7 @@ SET_PROP(StokesTestProblem, Fluid) { private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: typedef Opm::GasPhase > type; }; @@ -92,15 +93,15 @@ namespace Ewoms { * free-flow on the left and no-flow at the top and bottom boundaries. */ template -class StokesTestProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +class StokesTestProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; @@ -128,7 +129,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ StokesTestProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -154,8 +155,7 @@ public: * This problem assumes a constant temperature of 10 degrees Celsius. */ template - Scalar temperature(const Context &context, - int spaceIdx, int timeIdx) const + Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const { return 273.15 + 10; } // -> 10 deg C //! \} @@ -173,7 +173,8 @@ public: * a parabolic velocity profile via constraints. */ template - void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const + void boundary(BoundaryRateVector &values, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -183,16 +184,16 @@ public: // parabolic velocity profile const Scalar maxVelocity = 1.0; - Scalar a = - 4*maxVelocity/(height*height); - Scalar b = - a*height; + Scalar a = -4 * maxVelocity / (height * height); + Scalar b = -a * height; Scalar c = 0; DimVector velocity(0.0); - velocity[0] = a * y*y + b * y + c; + velocity[0] = a * y * y + b * y + c; if (onRightBoundary_(pos)) values.setOutFlow(context, spaceIdx, timeIdx); - else if(onLeftBoundary_(pos)) { + else if (onLeftBoundary_(pos)) { // left boundary is constraint! values = 0.0; } @@ -213,9 +214,8 @@ public: * \copydoc VcfvProblem::initial */ template - void initial(PrimaryVariables &values, - const Context &context, - int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -225,12 +225,12 @@ public: // parabolic velocity profile on boundaries const Scalar maxVelocity = 1.0; - Scalar a = - 4*maxVelocity/(height*height); - Scalar b = - a*height; + Scalar a = -4 * maxVelocity / (height * height); + Scalar b = -a * height; Scalar c = 0; DimVector velocity(0.0); - velocity[0] = a * y*y + b * y + c; + velocity[0] = a * y * y + b * y + c; for (int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) values[velocity0Idx + axisIdx] = velocity[axisIdx]; @@ -244,9 +244,8 @@ public: * is 0 everywhere. */ template - void source(RateVector &rate, - const Context &context, - int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = Scalar(0.0); } /*! @@ -256,8 +255,7 @@ public: * velocity profile using constraints. */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); @@ -266,7 +264,9 @@ public: PrimaryVariables initCond; initial(initCond, context, spaceIdx, timeIdx); - constraints.setConstraint(pressureIdx, conti0EqIdx, initCond[pressureIdx]);; + constraints.setConstraint(pressureIdx, conti0EqIdx, + initCond[pressureIdx]); + ; for (int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) constraints.setConstraint(velocity0Idx + axisIdx, momentum0EqIdx + axisIdx, @@ -291,9 +291,8 @@ private: bool onBoundary_(const GlobalPosition &pos) const { - return - onLeftBoundary_(pos) || onRightBoundary_(pos) || - onLowerBoundary_(pos) || onUpperBoundary_(pos); + return onLeftBoundary_(pos) || onRightBoundary_(pos) + || onLowerBoundary_(pos) || onUpperBoundary_(pos); } Scalar eps_; diff --git a/tests/models/problems/waterairproblem.hh b/tests/models/problems/waterairproblem.hh index 61d3a568e..30fa186f7 100644 --- a/tests/models/problems/waterairproblem.hh +++ b/tests/models/problems/waterairproblem.hh @@ -68,7 +68,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::TwoPhaseMaterialTraits Traits; + /*nonWettingPhaseIdx=*/FluidSystem::gPhaseIdx> + Traits; // define the material law which is parameterized by effective // saturations @@ -149,9 +150,8 @@ namespace Ewoms { * saturation of zero and a geothermal temperature gradient of 0.03 * K/m. */ -template -class WaterAirProblem - : public GET_PROP_TYPE(TypeTag, BaseProblem) +template +class WaterAirProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) { typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType; @@ -162,31 +162,30 @@ class WaterAirProblem // copy some indices for convenience typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { - numPhases = FluidSystem::numPhases, + enum { numPhases = FluidSystem::numPhases }; - // energy related indices - temperatureIdx = Indices::temperatureIdx, - energyEqIdx = Indices::energyEqIdx, + // energy related indices + enum { temperatureIdx = Indices::temperatureIdx }; + enum { energyEqIdx = Indices::energyEqIdx }; - // component indices - H2OIdx = FluidSystem::H2OIdx, - AirIdx = FluidSystem::AirIdx, + // component indices + enum { H2OIdx = FluidSystem::H2OIdx }; + enum { AirIdx = FluidSystem::AirIdx }; - // phase indices - lPhaseIdx = FluidSystem::lPhaseIdx, - gPhaseIdx = FluidSystem::gPhaseIdx, + // phase indices + enum { lPhaseIdx = FluidSystem::lPhaseIdx }; + enum { gPhaseIdx = FluidSystem::gPhaseIdx }; - // equation indices - conti0EqIdx = Indices::conti0EqIdx, + // equation indices + enum { conti0EqIdx = Indices::conti0EqIdx }; - // Grid and world dimension - dim = GridView::dimension, - dimWorld = GridView::dimensionworld - }; + // Grid and world dimension + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; + typedef typename GET_PROP_TYPE(TypeTag, + BoundaryRateVector) BoundaryRateVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints; typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; @@ -195,7 +194,8 @@ class WaterAirProblem typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw; - typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams; + typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) + HeatConductionLawParams; typedef typename GridView::ctype CoordScalar; typedef Dune::FieldVector GlobalPosition; @@ -207,7 +207,7 @@ public: * \copydoc Doxygen::defaultProblemConstructor */ WaterAirProblem(TimeManager &timeManager) -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) : ParentType(timeManager, GET_PROP_TYPE(TypeTag, GridCreator)::grid().leafGridView()) #else @@ -276,7 +276,8 @@ public: * permeable than the lower one. */ template - const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, int timeIdx) const + const DimMatrix &intrinsicPermeability(const Context &context, int spaceIdx, + int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -301,7 +302,8 @@ public: * \copydoc VcfvMultiPhaseProblem::materialLawParams */ template - const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const + const MaterialLawParams &materialLawParams(const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); if (isFineMaterial_(pos)) @@ -316,18 +318,18 @@ public: * In this case, we assume the rock-matrix to be granite. */ template - Scalar heatCapacitySolid(const Context &context, int spaceIdx, int timeIdx) const + Scalar heatCapacitySolid(const Context &context, int spaceIdx, + int timeIdx) const { - return - 790 // specific heat capacity of granite [J / (kg K)] - * 2700; // density of granite [kg/m^3] + return 790 // specific heat capacity of granite [J / (kg K)] + * 2700; // density of granite [kg/m^3] } /*! * \copydoc VcfvMultiPhaseProblem::heatConductionParams */ template - const HeatConductionLawParams& + const HeatConductionLawParams & heatConductionParams(const Context &context, int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); @@ -352,15 +354,12 @@ public: * right boundaries of the domain. */ template - void boundary(BoundaryRateVector &values, - const Context &context, + void boundary(BoundaryRateVector &values, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.cvCenter(spaceIdx, timeIdx); - assert(onLeftBoundary_(pos) || - onLowerBoundary_(pos) || - onRightBoundary_(pos) || - onUpperBoundary_(pos)); + assert(onLeftBoundary_(pos) || onLowerBoundary_(pos) + || onRightBoundary_(pos) || onUpperBoundary_(pos)); if (onInlet_(pos)) { RateVector massRate(0.0); @@ -370,8 +369,6 @@ public: values.setMassRate(massRate); } else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) { - //int globalIdx = context.elemContext().globalSpaceIndex(context.insideScvIndex(spaceIdx,timeIdx), timeIdx); - Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); @@ -397,9 +394,10 @@ public: * liquid water and assume hydrostatic pressure. */ template - void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const + void initial(PrimaryVariables &values, const Context &context, int spaceIdx, + int timeIdx) const { - //int globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + // int globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx); Opm::CompositionalFluidState fs; initialFluidState_(fs, context, spaceIdx, timeIdx); @@ -415,14 +413,14 @@ public: * of the finite-volumes which are closest to the inlet constant. */ template - void constraints(Constraints &constraints, - const Context &context, + void constraints(Constraints &constraints, const Context &context, int spaceIdx, int timeIdx) const { const auto &pos = context.pos(spaceIdx, timeIdx); if (onInlet_(pos)) { - constraints.setConstraint(temperatureIdx, energyEqIdx, 380);; + constraints.setConstraint(temperatureIdx, energyEqIdx, 380); + ; } } @@ -433,8 +431,8 @@ public: * everywhere. */ template - void source(RateVector &rate, - const Context &context, int spaceIdx, int timeIdx) const + void source(RateVector &rate, const Context &context, int spaceIdx, + int timeIdx) const { rate = 0; } //! \} @@ -459,12 +457,13 @@ private: { return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); } template - void initialFluidState_(FluidState &fs, const Context &context, int spaceIdx, int timeIdx) const + void initialFluidState_(FluidState &fs, const Context &context, + int spaceIdx, int timeIdx) const { const GlobalPosition &pos = context.pos(spaceIdx, timeIdx); Scalar densityW = 1000.0; - fs.setPressure(lPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81); + fs.setPressure(lPhaseIdx, 1e5 + (maxDepth_ - pos[1]) * densityW * 9.81); fs.setSaturation(lPhaseIdx, 1.0); fs.setMoleFraction(lPhaseIdx, H2OIdx, 1.0); fs.setMoleFraction(lPhaseIdx, AirIdx, 0.0); @@ -472,18 +471,20 @@ private: if (inHighTemperatureRegion_(pos)) fs.setTemperature(380); else - fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03); + fs.setTemperature(283.0 + (maxDepth_ - pos[1]) * 0.03); // set the gas saturation and pressure fs.setSaturation(gPhaseIdx, 0); Scalar pc[numPhases]; const auto &matParams = materialLawParams(context, spaceIdx, timeIdx); MaterialLaw::capillaryPressures(pc, matParams, fs); - fs.setPressure(gPhaseIdx, fs.pressure(lPhaseIdx) + (pc[gPhaseIdx] - pc[lPhaseIdx])); + fs.setPressure(gPhaseIdx, fs.pressure(lPhaseIdx) + + (pc[gPhaseIdx] - pc[lPhaseIdx])); typename FluidSystem::ParameterCache paramCache; typedef Opm::ComputeFromReferencePhase CFRP; - CFRP::solve(fs, paramCache, lPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/true); + CFRP::solve(fs, paramCache, lPhaseIdx, /*setViscosity=*/false, + /*setEnthalpy=*/true); } void computeHeatCondParams_(HeatConductionLawParams ¶ms, Scalar poro) @@ -507,12 +508,13 @@ private: for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { Scalar lambdaSaturated; if (FluidSystem::isLiquid(phaseIdx)) { - Scalar lambdaFluid = - FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); - lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro); + Scalar lambdaFluid + = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)) + + std::pow(lambdaFluid, poro); } else - lambdaSaturated = std::pow(lambdaGranite, (1-poro)); + lambdaSaturated = std::pow(lambdaGranite, (1 - poro)); params.setFullySaturatedLambda(phaseIdx, lambdaSaturated); if (!FluidSystem::isLiquid(phaseIdx)) @@ -521,7 +523,7 @@ private: } bool isFineMaterial_(const GlobalPosition &pos) const - { return pos[dim-1] > layerBottom_; } + { return pos[dim - 1] > layerBottom_; } DimMatrix fineK_; DimMatrix coarseK_; diff --git a/tests/models/reservoir_blackoil.cc b/tests/models/reservoir_blackoil.cc index eeb924edb..f0685da23 100644 --- a/tests/models/reservoir_blackoil.cc +++ b/tests/models/reservoir_blackoil.cc @@ -30,9 +30,10 @@ namespace Opm { namespace Properties { NEW_TYPE_TAG(ReservoirProblem, INHERITS_FROM(VcfvBlackOil, ReservoirBaseProblem)); -}} +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(ReservoirProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/test_navierstokes.cc b/tests/models/test_navierstokes.cc index 2543e552a..fded9bc31 100644 --- a/tests/models/test_navierstokes.cc +++ b/tests/models/test_navierstokes.cc @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief Test for the isothermal Navier-Stokes VCVF discretization; this test case is + * \brief Test for the isothermal Navier-Stokes VCVF discretization; this test + * case is * known as lid-driven cavity-flow in literature. */ #include "config.h" @@ -26,7 +27,7 @@ #include #include "problems/navierstokestestproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(NavierStokesTestProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/test_quadrature.cc b/tests/models/test_quadrature.cc index f49a78fcf..f0b9ede5f 100644 --- a/tests/models/test_quadrature.cc +++ b/tests/models/test_quadrature.cc @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief A test for numerical integration using the vertex-centered finite volume geometries. + * \brief A test for numerical integration using the vertex-centered finite + * volume geometries. */ #include "config.h" @@ -53,16 +54,14 @@ void testIdenityMapping() { QuadratureGeom foo; - Scalar corners[][3] = { - { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 1, 1, 0 }, - { 0, 0, 1 }, - { 1, 0, 1 }, - { 0, 1, 1 }, - { 1, 1, 1 } - }; + Scalar corners[][3] = { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 1, 0 }, + { 1, 1, 0 }, + { 0, 0, 1 }, + { 1, 0, 1 }, + { 0, 1, 1 }, + { 1, 1, 1 } }; foo.setCorners(corners, 8); std::cout << "testing identity mapping...\n"; @@ -72,9 +71,9 @@ void testIdenityMapping() for (int k = 0; k < n; ++k) { LocalPosition localPos; - localPos[0] = Scalar(i)/(n - 1); - localPos[1] = Scalar(j)/(n - 1); - localPos[2] = Scalar(k)/(n - 1); + localPos[0] = Scalar(i) / (n - 1); + localPos[1] = Scalar(j) / (n - 1); + localPos[2] = Scalar(k) / (n - 1); GlobalPosition globalPos = foo.global(localPos); @@ -107,10 +106,12 @@ void writeTetrahedronSubControlVolumes(const Grid &grid) for (; eIt != eEndIt; ++eIt) { fvElemGeom.update(gridView, *eIt); for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; ++scvIdx) { - const auto &scvLocalGeom = *(fvElemGeom.subContVol[scvIdx].localGeometry); + const auto &scvLocalGeom + = *(fvElemGeom.subContVol[scvIdx].localGeometry); - for (int i = 0; i < scvLocalGeom.numCorners; ++ i) { - GlobalPosition pos(eIt->geometry().global(scvLocalGeom.corner(i))); + for (int i = 0; i < scvLocalGeom.numCorners; ++i) { + GlobalPosition pos( + eIt->geometry().global(scvLocalGeom.corner(i))); gf2.insertVertex(pos); } } @@ -121,15 +122,16 @@ void writeTetrahedronSubControlVolumes(const Grid &grid) for (; eIt != eEndIt; ++eIt) { fvElemGeom.update(gridView, *eIt); for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; ++scvIdx) { - const auto &scvLocalGeom = *fvElemGeom.subContVol[scvIdx].localGeometry; + const auto &scvLocalGeom + = *fvElemGeom.subContVol[scvIdx].localGeometry; std::vector vertexIndices; - for (int i = 0; i < scvLocalGeom.numCorners; ++ i) { + for (int i = 0; i < scvLocalGeom.numCorners; ++i) { vertexIndices.push_back(cornerOffset); - ++ cornerOffset; + ++cornerOffset; } - gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), + gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim), vertexIndices); } } @@ -147,22 +149,16 @@ void testTetrahedron() typedef Dune::ALUGrid Grid; typedef Dune::GridFactory GridFactory; GridFactory gf; - Scalar corners[][3] = { - { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 0, 0, 1 } - }; + Scalar corners[][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } }; - for (unsigned i = 0; i < sizeof(corners)/sizeof(corners[0]); ++i) { + for (unsigned i = 0; i < sizeof(corners) / sizeof(corners[0]); ++i) { GlobalPosition pos; for (unsigned j = 0; j < dim; ++j) pos[j] = corners[i][j]; gf.insertVertex(pos); } std::vector v = { 0, 1, 2, 3 }; - gf.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,dim), - v); + gf.insertElement(Dune::GeometryType(Dune::GeometryType::simplex, dim), v); auto *grid = gf.createGrid(); // write the sub-control volumes to a VTK file. @@ -193,10 +189,12 @@ void writeCubeSubControlVolumes(const Grid &grid) for (; eIt != eEndIt; ++eIt) { fvElemGeom.update(gridView, *eIt); for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; ++scvIdx) { - const auto &scvLocalGeom = *(fvElemGeom.subContVol[scvIdx].localGeometry); + const auto &scvLocalGeom + = *(fvElemGeom.subContVol[scvIdx].localGeometry); - for (int i = 0; i < scvLocalGeom.numCorners; ++ i) { - GlobalPosition pos(eIt->geometry().global(scvLocalGeom.corner(i))); + for (int i = 0; i < scvLocalGeom.numCorners; ++i) { + GlobalPosition pos( + eIt->geometry().global(scvLocalGeom.corner(i))); gf2.insertVertex(pos); } } @@ -207,15 +205,16 @@ void writeCubeSubControlVolumes(const Grid &grid) for (; eIt != eEndIt; ++eIt) { fvElemGeom.update(gridView, *eIt); for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; ++scvIdx) { - const auto &scvLocalGeom = *fvElemGeom.subContVol[scvIdx].localGeometry; + const auto &scvLocalGeom + = *fvElemGeom.subContVol[scvIdx].localGeometry; std::vector vertexIndices; - for (int i = 0; i < scvLocalGeom.numCorners; ++ i) { + for (int i = 0; i < scvLocalGeom.numCorners; ++i) { vertexIndices.push_back(cornerOffset); - ++ cornerOffset; + ++cornerOffset; } - gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), + gf2.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim), vertexIndices); } } @@ -233,26 +232,23 @@ void testCube() typedef Dune::ALUGrid Grid; typedef Dune::GridFactory GridFactory; GridFactory gf; - Scalar corners[][3] = { - { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 2, 0 }, - { 3, 3, 0 }, - { 0, 0, 4 }, - { 5, 0, 5 }, - { 0, 6, 6 }, - { 7, 7, 7 }, - }; + Scalar corners[][3] = { { 0, 0, 0 }, + { 1, 0, 0 }, + { 0, 2, 0 }, + { 3, 3, 0 }, + { 0, 0, 4 }, + { 5, 0, 5 }, + { 0, 6, 6 }, + { 7, 7, 7 }, }; - for (unsigned i = 0; i < sizeof(corners)/sizeof(corners[0]); ++i) { + for (unsigned i = 0; i < sizeof(corners) / sizeof(corners[0]); ++i) { GlobalPosition pos; for (unsigned j = 0; j < dim; ++j) pos[j] = corners[i][j]; gf.insertVertex(pos); } std::vector v = { 0, 1, 2, 3, 4, 5, 6, 7 }; - gf.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), - v); + gf.insertElement(Dune::GeometryType(Dune::GeometryType::cube, dim), v); auto *grid = gf.createGrid(); // write the sub-control volumes to a VTK file. @@ -266,7 +262,7 @@ void testQuadrature() { std::cout << "testing SCV quadrature...\n"; -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,3) +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) std::bitset isPeriodic(false); std::array cellRes; #else @@ -283,16 +279,15 @@ void testQuadrature() #ifdef HAVE_MPI Dune::MPIHelper::getCommunicator(), #endif - upperRight, // upper right - cellRes, // number of cells - isPeriodic, - 0); // overlap + upperRight, // upper right + cellRes, // number of cells + isPeriodic, 0); // overlap // compute approximate integral auto gridView = grid.leafView(); auto eIt = gridView.begin<0>(); const auto eEndIt = gridView.end<0>(); - Scalar result=0; + Scalar result = 0; // instanciate a FvElementGeometry typedef Ewoms::VcfvElementGeometry FvElementGeometry; FvElementGeometry fvElemGeom; @@ -303,43 +298,46 @@ void testQuadrature() // loop over all sub-control volumes for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; ++scvIdx) { - const auto &scvLocalGeom = *fvElemGeom.subContVol[scvIdx].localGeometry; + const auto &scvLocalGeom + = *fvElemGeom.subContVol[scvIdx].localGeometry; Dune::GeometryType geomType = scvLocalGeom.type(); static const int quadratureOrder = 2; - const auto &rule = Dune::QuadratureRules::rule(geomType, quadratureOrder); + const auto &rule + = Dune::QuadratureRules::rule(geomType, + quadratureOrder); // integrate f over the sub-control volume - for (auto it = rule.begin(); it != rule.end(); ++ it) - { + for (auto it = rule.begin(); it != rule.end(); ++it) { auto posScvLocal = it->position(); auto posElemLocal = scvLocalGeom.global(posScvLocal); auto posGlobal = elemGeom.global(posScvLocal); Scalar fval = f(posGlobal); Scalar weight = it->weight(); - Scalar detjac = - scvLocalGeom.integrationElement(posScvLocal) * - elemGeom.integrationElement(posElemLocal); + Scalar detjac = scvLocalGeom.integrationElement(posScvLocal) + * elemGeom.integrationElement(posElemLocal); result += fval * weight * detjac; } } } - std::cout << "result: " << result << " (expected value: " << 1.0/(1 << dim) << ")\n"; + std::cout << "result: " << result + << " (expected value: " << 1.0 / (1 << dim) << ")\n"; } -int main(int argc, char** argv) +int main(int argc, char **argv) { // initialize MPI, finalize is done automatically on exit Dune::MPIHelper::instance(argc, argv); testIdenityMapping(); - // test the quadrature in a tetrahedron. since the CLang compiler - // prior to version 3.2 generates incorrect code here, we do not - // do it if the compiler is clang 3.1 or older. -#if !__clang__ || __clang_major__ > 3 || (__clang_major__ == 3 && __clang_minor__ >= 3) +// test the quadrature in a tetrahedron. since the CLang compiler +// prior to version 3.2 generates incorrect code here, we do not +// do it if the compiler is clang 3.1 or older. +#if !__clang__ || __clang_major__ > 3 \ + || (__clang_major__ == 3 && __clang_minor__ >= 3) testTetrahedron(); testCube(); #endif diff --git a/tests/models/test_stokes.cc b/tests/models/test_stokes.cc index 7d816e259..fa6be4dee 100644 --- a/tests/models/test_stokes.cc +++ b/tests/models/test_stokes.cc @@ -24,7 +24,7 @@ #include #include "problems/stokestestproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(StokesTestProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/test_stokes2c.cc b/tests/models/test_stokes2c.cc index 46d86ddd8..d2a81d86a 100644 --- a/tests/models/test_stokes2c.cc +++ b/tests/models/test_stokes2c.cc @@ -25,7 +25,7 @@ #include #include "problems/stokes2ctestproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(Stokes2cTestProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/test_stokesni.cc b/tests/models/test_stokesni.cc index bfcdd9b02..f31ca76a1 100644 --- a/tests/models/test_stokesni.cc +++ b/tests/models/test_stokesni.cc @@ -25,7 +25,7 @@ #include #include "problems/stokesnitestproblem.hh" -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(StokesNITestProblem) ProblemTypeTag; return Ewoms::start(argc, argv); diff --git a/tests/models/waterair_pvs_ni.cc b/tests/models/waterair_pvs_ni.cc index f9f17dac9..ab1667eab 100644 --- a/tests/models/waterair_pvs_ni.cc +++ b/tests/models/waterair_pvs_ni.cc @@ -32,9 +32,10 @@ namespace Properties { NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(VcfvPvs, WaterAirBaseProblem)); SET_BOOL_PROP(WaterAirProblem, EnableEnergy, true); -} } +} +} -int main(int argc, char** argv) +int main(int argc, char **argv) { typedef TTAG(WaterAirProblem) ProblemTypeTag; return Ewoms::start(argc, argv);