28 #ifndef EWOMS_OBSTACLE_PROBLEM_HH
29 #define EWOMS_OBSTACLE_PROBLEM_HH
33 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
34 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
35 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
36 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <opm/material/heatconduction/Somerton.hpp>
41 #include <opm/common/Unused.hpp>
43 #include <dune/grid/yaspgrid.hh>
44 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
55 template <
class TypeTag>
60 namespace Properties {
71 Opm::FluidSystems::H2ON2<
typename GET_PROP_TYPE(TypeTag, Scalar)>);
74 SET_PROP(ObstacleBaseProblem, MaterialLaw)
79 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
80 typedef Opm::TwoPhaseMaterialTraits<Scalar,
81 FluidSystem::liquidPhaseIdx,
82 FluidSystem::gasPhaseIdx>
85 typedef Opm::LinearMaterial<MaterialTraits> EffMaterialLaw;
88 typedef Opm::EffToAbsLaw<EffMaterialLaw> type;
92 SET_PROP(ObstacleBaseProblem, HeatConductionLaw)
96 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
100 typedef Opm::Somerton<FluidSystem, Scalar> type;
113 SET_STRING_PROP(ObstacleBaseProblem, GridFile, "./data/obstacle_24x16.dgf");
144 template <
class TypeTag>
147 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
152 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
153 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
154 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
155 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
156 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
157 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
158 typedef typename
GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
159 typedef typename HeatConductionLaw::Params HeatConductionLawParams;
163 dim = GridView::dimension,
164 dimWorld = GridView::dimensionworld,
166 gasPhaseIdx = FluidSystem::gasPhaseIdx,
167 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
168 H2OIdx = FluidSystem::H2OIdx,
169 N2Idx = FluidSystem::N2Idx
172 typedef Dune::FieldVector<typename GridView::ctype, dimWorld> GlobalPosition;
173 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
174 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
175 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
183 : ParentType(simulator)
191 ParentType::finishInit();
194 temperature_ = 273.15 + 25;
197 Scalar Tmin = temperature_ - 1.0;
198 Scalar Tmax = temperature_ + 1.0;
201 Scalar pmin = 1.0e5 * 0.75;
202 Scalar pmax = 2.0e5 * 1.25;
205 FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
208 coarseK_ = this->toDimMatrix_(1e-12);
209 fineK_ = this->toDimMatrix_(1e-15);
213 coarsePorosity_ = 0.3;
216 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
217 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
218 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
219 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
223 fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
224 fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
225 coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
226 coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
238 fineMaterialParams_.finalize();
239 coarseMaterialParams_.finalize();
242 computeHeatCondParams_(fineHeatCondParams_, finePorosity_);
243 computeHeatCondParams_(coarseHeatCondParams_, coarsePorosity_);
254 this->model().checkConservativeness();
257 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
258 PrimaryVariables phaseStorage;
259 this->model().globalPhaseStorage(phaseStorage, phaseIdx);
261 if (this->gridView().comm().rank() == 0) {
262 std::cout <<
"Storage in " << FluidSystem::phaseName(phaseIdx)
263 <<
"Phase: [" << phaseStorage <<
"]"
264 <<
"\n" << std::flush;
270 this->model().globalStorage(storage);
273 if (this->gridView().comm().rank() == 0) {
274 std::cout <<
"Storage total: [" << storage <<
"]"
275 <<
"\n" << std::flush;
290 std::ostringstream oss;
292 <<
"_" << Model::name();
301 template <
class Context>
303 unsigned spaceIdx OPM_UNUSED,
304 unsigned timeIdx OPM_UNUSED)
const
305 {
return temperature_; }
310 template <
class Context>
314 unsigned timeIdx)
const
316 if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
324 template <
class Context>
327 unsigned timeIdx)
const
329 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
330 if (isFineMaterial_(pos))
331 return finePorosity_;
333 return coarsePorosity_;
339 template <
class Context>
340 const MaterialLawParams&
343 unsigned timeIdx)
const
345 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
346 if (isFineMaterial_(pos))
347 return fineMaterialParams_;
349 return coarseMaterialParams_;
358 template <
class Context>
360 unsigned spaceIdx OPM_UNUSED,
361 unsigned timeIdx OPM_UNUSED)
const
370 template <
class Context>
371 const HeatConductionLawParams &
374 unsigned timeIdx)
const
376 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
377 if (isFineMaterial_(pos))
378 return fineHeatCondParams_;
379 return coarseHeatCondParams_;
392 template <
class Context>
394 const Context& context,
396 unsigned timeIdx)
const
398 const auto& pos = context.pos(spaceIdx, timeIdx);
401 values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
402 else if (onOutlet_(pos))
403 values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
418 template <
class Context>
420 const Context& context,
422 unsigned timeIdx)
const
424 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
425 values.assignMassConservative(outletFluidState_, matParams);
434 template <
class Context>
436 const Context& context OPM_UNUSED,
437 unsigned spaceIdx OPM_UNUSED,
438 unsigned timeIdx OPM_UNUSED)
const
448 bool isFineMaterial_(
const GlobalPosition& pos)
const
449 {
return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
451 bool onInlet_(
const GlobalPosition& globalPos)
const
453 Scalar x = globalPos[0];
454 Scalar y = globalPos[1];
455 return x >= 60 - eps_ && y <= 10;
458 bool onOutlet_(
const GlobalPosition& globalPos)
const
460 Scalar x = globalPos[0];
461 Scalar y = globalPos[1];
462 return x < eps_ && y <= 10;
465 void initFluidStates_()
467 initFluidState_(inletFluidState_, coarseMaterialParams_,
469 initFluidState_(outletFluidState_, coarseMaterialParams_,
473 template <
class Flu
idState>
474 void initFluidState_(FluidState& fs,
const MaterialLawParams& matParams,
bool isInlet)
476 unsigned refPhaseIdx;
477 unsigned otherPhaseIdx;
480 fs.setTemperature(temperature_);
484 refPhaseIdx = liquidPhaseIdx;
485 otherPhaseIdx = gasPhaseIdx;
488 fs.setSaturation(liquidPhaseIdx, 1.0);
491 fs.setPressure(liquidPhaseIdx, 2e5);
494 fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
495 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
499 refPhaseIdx = gasPhaseIdx;
500 otherPhaseIdx = liquidPhaseIdx;
503 fs.setSaturation(gasPhaseIdx, 1.0);
506 fs.setPressure(gasPhaseIdx, 1e5);
509 fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
510 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
514 fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
518 MaterialLaw::capillaryPressures(pC, matParams, fs);
519 fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
520 + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
524 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem>
525 ComputeFromReferencePhase;
527 typename FluidSystem::template ParameterCache<Scalar> paramCache;
528 ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
533 void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
535 Scalar lambdaWater = 0.6;
536 Scalar lambdaGranite = 2.8;
538 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
539 * std::pow(lambdaWater, poro);
540 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
542 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
543 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
544 params.setVacuumLambda(lambdaDry);
550 Scalar coarsePorosity_;
551 Scalar finePorosity_;
553 MaterialLawParams fineMaterialParams_;
554 MaterialLawParams coarseMaterialParams_;
556 HeatConductionLawParams fineHeatCondParams_;
557 HeatConductionLawParams coarseHeatCondParams_;
559 Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
560 Opm::CompositionalFluidState<Scalar, FluidSystem> outletFluidState_;
std::string name() const
The problem name.
Definition: obstacleproblem.hh:288
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:359
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:312
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: obstacleproblem.hh:393
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: obstacleproblem.hh:189
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Declares the properties required for the NCP compositional multi-phase model.
void endTimeStep()
Called by the simulator after each time integration.
Definition: obstacleproblem.hh:251
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:325
const HeatConductionLawParams & heatConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:372
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: obstacleproblem.hh:419
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:302
Problem where liquid water is first stopped by a low-permeability lens and then seeps though it...
Definition: obstacleproblem.hh:56
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:341
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: obstacleproblem.hh:435
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394