My Project
Loading...
Searching...
No Matches
richardslensproblem.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_RICHARDS_LENS_PROBLEM_HH
29#define EWOMS_RICHARDS_LENS_PROBLEM_HH
30
31#include <opm/models/richards/richardsmodel.hh>
32
33#include <opm/material/components/SimpleH2O.hpp>
34#include <opm/material/fluidsystems/LiquidPhase.hpp>
35#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
36#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
37#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
39
40#include <dune/grid/yaspgrid.hh>
41#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
42
43#include <dune/common/version.hh>
44#include <dune/common/fvector.hh>
45#include <dune/common/fmatrix.hh>
46
47namespace Opm {
48template <class TypeTag>
49class RichardsLensProblem;
50
51} // namespace Opm
52
53namespace Opm::Properties {
54
55// Create new type tags
56namespace TTag {
57struct RichardsLensProblem { using InheritsFrom = std::tuple<Richards>; };
58} // end namespace TTag
59
60// Use 2d YaspGrid
61template<class TypeTag>
62struct Grid<TypeTag, TTag::RichardsLensProblem> { using type = Dune::YaspGrid<2>; };
63
64// Set the physical problem to be solved
65template<class TypeTag>
66struct Problem<TypeTag, TTag::RichardsLensProblem> { using type = Opm::RichardsLensProblem<TypeTag>; };
67
68// Set the wetting phase
69template<class TypeTag>
70struct WettingFluid<TypeTag, TTag::RichardsLensProblem>
71{
72private:
73 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
74
75public:
76 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
77};
78
79// Set the material Law
80template<class TypeTag>
81struct MaterialLaw<TypeTag, TTag::RichardsLensProblem>
82{
83private:
84 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
85 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
86 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
87
88 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
89 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
90 /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
91 /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
92
93 // define the material law which is parameterized by effective
94 // saturations
95 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
96
97public:
98 // define the material law parameterized by absolute saturations
99 using type = Opm::EffToAbsLaw<EffectiveLaw>;
100};
101
102// Enable gravitational acceleration
103template<class TypeTag>
104struct EnableGravity<TypeTag, TTag::RichardsLensProblem> { static constexpr bool value = true; };
105
106// Use central differences to approximate the Jacobian matrix
107template<class TypeTag>
108struct NumericDifferenceMethod<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 0; };
109
110// Set the maximum number of newton iterations of a time step
111template<class TypeTag>
112struct NewtonMaxIterations<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 28; };
113
114// Set the "desireable" number of newton iterations of a time step
115template<class TypeTag>
116struct NewtonTargetIterations<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 18; };
117
118// Do not write the intermediate results of the newton method
119template<class TypeTag>
120struct NewtonWriteConvergence<TypeTag, TTag::RichardsLensProblem> { static constexpr bool value = false; };
121
122// The default for the end time of the simulation
123template<class TypeTag>
124struct EndTime<TypeTag, TTag::RichardsLensProblem>
125{
126 using type = GetPropType<TypeTag, Scalar>;
127 static constexpr type value = 3000;
128};
129
130// The default for the initial time step size of the simulation
131template<class TypeTag>
132struct InitialTimeStepSize<TypeTag, TTag::RichardsLensProblem>
133{
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 100;
136};
137
138// The default DGF file to load
139template<class TypeTag>
140struct GridFile<TypeTag, TTag::RichardsLensProblem> { static constexpr auto value = "./data/richardslens_24x16.dgf"; };
141
142} // namespace Opm::Properties
143
144namespace Opm {
145
162template <class TypeTag>
163class RichardsLensProblem : public GetPropType<TypeTag, Properties::BaseProblem>
164{
165 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
166
167 using GridView = GetPropType<TypeTag, Properties::GridView>;
168 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
169 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
170 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
171 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
172 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
173 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
174 using Model = GetPropType<TypeTag, Properties::Model>;
175 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
176 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
177
178 using Indices = GetPropType<TypeTag, Properties::Indices>;
179 enum {
180 // copy some indices for convenience
181 pressureWIdx = Indices::pressureWIdx,
182 contiEqIdx = Indices::contiEqIdx,
183 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
184 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
185 numPhases = FluidSystem::numPhases,
186
187 // Grid and world dimension
188 dimWorld = GridView::dimensionworld
189 };
190
191 // get the material law from the property system
192 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
194 using MaterialLawParams = typename MaterialLaw::Params;
195
196 using CoordScalar = typename GridView::ctype;
197 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
198 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
199 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
200
201public:
205 RichardsLensProblem(Simulator& simulator)
206 : ParentType(simulator)
207 , pnRef_(1e5)
208 {
209 dofIsInLens_.resize(simulator.model().numGridDof());
210 }
211
216 {
217 ParentType::finishInit();
218
219 eps_ = 3e-6;
220 pnRef_ = 1e5;
221
222 lensLowerLeft_[0] = 1.0;
223 lensLowerLeft_[1] = 2.0;
224
225 lensUpperRight_[0] = 4.0;
226 lensUpperRight_[1] = 3.0;
227
228 // parameters for the Van Genuchten law
229 // alpha and n
230 lensMaterialParams_.setVgAlpha(0.00045);
231 lensMaterialParams_.setVgN(7.3);
232 lensMaterialParams_.finalize();
233
234 outerMaterialParams_.setVgAlpha(0.0037);
235 outerMaterialParams_.setVgN(4.7);
236 outerMaterialParams_.finalize();
237
238 // parameters for the linear law
239 // minimum and maximum pressures
240 // lensMaterialParams_.setEntryPC(0);
241 // outerMaterialParams_.setEntryPC(0);
242 // lensMaterialParams_.setMaxPC(0);
243 // outerMaterialParams_.setMaxPC(0);
244
245 lensK_ = this->toDimMatrix_(1e-12);
246 outerK_ = this->toDimMatrix_(5e-12);
247
248 // determine which degrees of freedom are in the lens
249 Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
250 for (const auto& elem : elements(this->gridView())) {
251 stencil.update(elem);
252 for (unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
253 unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
254 const auto& dofPos = stencil.subControlVolume(dofIdx).center();
255 dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
256 }
257 }
258 }
259
264
268 std::string name() const
269 {
270 std::ostringstream oss;
271 oss << "lens_richards_"
272 << Model::discretizationName();
273 return oss.str();
274 }
275
280 {
281#ifndef NDEBUG
282 this->model().checkConservativeness();
283
284 // Calculate storage terms
285 EqVector storage;
286 this->model().globalStorage(storage);
287
288 // Write mass balance information for rank 0
289 if (this->gridView().comm().rank() == 0) {
290 std::cout << "Storage: " << storage << std::endl << std::flush;
291 }
292#endif // NDEBUG
293 }
294
298 template <class Context>
299 Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
300 { return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
301
302 Scalar temperature(unsigned /*globalSpaceIdx*/, unsigned /*timeIdx*/) const
303 { return 273.15 + 10; } // -> 10°C
304
308 template <class Context>
309 const DimMatrix& intrinsicPermeability(const Context& context,
310 unsigned spaceIdx,
311 unsigned timeIdx) const
312 {
313 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
314 if (isInLens_(pos))
315 return lensK_;
316 return outerK_;
317 }
318
322 template <class Context>
323 Scalar porosity(const Context& /*context*/,
324 unsigned /*spaceIdx*/,
325 unsigned /*timeIdx*/) const
326 { return 0.4; }
327
331 template <class Context>
332 const MaterialLawParams& materialLawParams(const Context& context,
333 unsigned spaceIdx,
334 unsigned timeIdx) const
335 {
336 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
337 return materialLawParams(globalSpaceIdx, timeIdx);
338 }
339
340 const MaterialLawParams& materialLawParams(unsigned globalSpaceIdx,
341 unsigned /*timeIdx*/) const
342 {
343 if (dofIsInLens_[globalSpaceIdx])
344 return lensMaterialParams_;
345 return outerMaterialParams_;
346 }
347
353 template <class Context>
354 Scalar referencePressure(const Context& context,
355 unsigned spaceIdx,
356 unsigned timeIdx) const
357 { return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
358
359 // the Richards model does not have an element context available at all places
360 // where the reference pressure is required...
361 Scalar referencePressure(unsigned /*globalSpaceIdx*/,
362 unsigned /*timeIdx*/) const
363 { return pnRef_; }
364
366
371
375 template <class Context>
376 void boundary(BoundaryRateVector& values,
377 const Context& context,
378 unsigned spaceIdx,
379 unsigned timeIdx) const
380 {
381 const auto& pos = context.pos(spaceIdx, timeIdx);
382
383 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
384 const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
385
386 Scalar Sw = 0.0;
387 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
388 fs.setSaturation(wettingPhaseIdx, Sw);
389 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
390
391 PhaseVector pC;
392 MaterialLaw::capillaryPressures(pC, materialParams, fs);
393 fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
394 fs.setPressure(nonWettingPhaseIdx, pnRef_);
395
396 typename FluidSystem::template ParameterCache<Scalar> paramCache;
397 paramCache.updateAll(fs);
398 fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx));
399 //fs.setDensity(nonWettingPhaseIdx, FluidSystem::density(fs, paramCache, nonWettingPhaseIdx));
400
401 fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx));
402 //fs.setViscosity(nonWettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, nonWettingPhaseIdx));
403
404 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
405 }
406 else if (onInlet_(pos)) {
407 RateVector massRate(0.0);
408
409 // inflow of water
410 massRate[contiEqIdx] = -0.04; // kg / (m * s)
411
412 values.setMassRate(massRate);
413 }
414 else
415 values.setNoFlow();
416 }
417
419
424
428 template <class Context>
429 void initial(PrimaryVariables& values,
430 const Context& context,
431 unsigned spaceIdx,
432 unsigned timeIdx) const
433 {
434 const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
435
436 Scalar Sw = 0.0;
437 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
438 fs.setSaturation(wettingPhaseIdx, Sw);
439 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
440
441 PhaseVector pC;
442 MaterialLaw::capillaryPressures(pC, materialParams, fs);
443 values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
444 }
445
452 template <class Context>
453 void source(RateVector& rate,
454 const Context& /*context*/,
455 unsigned /*spaceIdx*/,
456 unsigned /*timeIdx*/) const
457 { rate = Scalar(0.0); }
458
460
461private:
462 bool onLeftBoundary_(const GlobalPosition& pos) const
463 { return pos[0] < this->boundingBoxMin()[0] + eps_; }
464
465 bool onRightBoundary_(const GlobalPosition& pos) const
466 { return pos[0] > this->boundingBoxMax()[0] - eps_; }
467
468 bool onLowerBoundary_(const GlobalPosition& pos) const
469 { return pos[1] < this->boundingBoxMin()[1] + eps_; }
470
471 bool onUpperBoundary_(const GlobalPosition& pos) const
472 { return pos[1] > this->boundingBoxMax()[1] - eps_; }
473
474 bool onInlet_(const GlobalPosition& pos) const
475 {
476 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
477 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
478 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
479 }
480
481 bool isInLens_(const GlobalPosition& pos) const
482 {
483 for (unsigned i = 0; i < dimWorld; ++i) {
484 if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
485 return false;
486 }
487 return true;
488 }
489
490 GlobalPosition lensLowerLeft_;
491 GlobalPosition lensUpperRight_;
492
493 DimMatrix lensK_;
494 DimMatrix outerK_;
495 MaterialLawParams lensMaterialParams_;
496 MaterialLawParams outerMaterialParams_;
497
498 std::vector<bool> dofIsInLens_;
499
500 Scalar eps_;
501 Scalar pnRef_;
502};
503} // namespace Opm
504
505#endif
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain.
Definition richardslensproblem.hh:164
void finishInit()
Definition richardslensproblem.hh:215
std::string name() const
Definition richardslensproblem.hh:268
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition richardslensproblem.hh:332
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition richardslensproblem.hh:354
RichardsLensProblem(Simulator &simulator)
Definition richardslensproblem.hh:205
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition richardslensproblem.hh:309
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition richardslensproblem.hh:453
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition richardslensproblem.hh:376
void endTimeStep()
Definition richardslensproblem.hh:279
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition richardslensproblem.hh:299
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition richardslensproblem.hh:429
Scalar porosity(const Context &, unsigned, unsigned) const
Definition richardslensproblem.hh:323
Definition richardslensproblem.hh:57