My Project
Loading...
Searching...
No Matches
outflowproblem.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*/
27#ifndef EWOMS_OUTFLOW_PROBLEM_HH
28#define EWOMS_OUTFLOW_PROBLEM_HH
29
30#include <opm/models/pvs/pvsproperties.hh>
31
32#include <opm/material/fluidstates/CompositionalFluidState.hpp>
33#include <opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp>
34
35#include <dune/grid/yaspgrid.hh>
36#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
37
38#include <dune/common/version.hh>
39#include <dune/common/fvector.hh>
40#include <dune/common/fmatrix.hh>
41
42namespace Opm {
43template <class TypeTag>
44class OutflowProblem;
45}
46
47namespace Opm::Properties {
48
49namespace TTag {
50
52
53} // namespace TTag
54
55// Set the grid type
56template<class TypeTag>
57struct Grid<TypeTag, TTag::OutflowBaseProblem> { using type = Dune::YaspGrid<2>; };
58
59// Set the problem property
60template<class TypeTag>
61struct Problem<TypeTag, TTag::OutflowBaseProblem> { using type = Opm::OutflowProblem<TypeTag>; };
62
63// Set fluid system
64template<class TypeTag>
65struct FluidSystem<TypeTag, TTag::OutflowBaseProblem>
66{
67private:
68 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
69
70public:
71 // Two-component single phase fluid system
72 using type = Opm::H2ON2LiquidPhaseFluidSystem<Scalar>;
73};
74
75// Disable gravity
76template<class TypeTag>
77struct EnableGravity<TypeTag, TTag::OutflowBaseProblem> { static constexpr bool value = false; };
78
79// Also write mass fractions to the output
80template<class TypeTag>
81struct VtkWriteMassFractions<TypeTag, TTag::OutflowBaseProblem> { static constexpr bool value = true; };
82
83// The default for the end time of the simulation
84template<class TypeTag>
85struct EndTime<TypeTag, TTag::OutflowBaseProblem>
86{
87 using type = GetPropType<TypeTag, Scalar>;
88 static constexpr type value = 100;
89};
90
91// The default for the initial time step size of the simulation
92template<class TypeTag>
93struct InitialTimeStepSize<TypeTag, TTag::OutflowBaseProblem>
94{
95 using type = GetPropType<TypeTag, Scalar>;
96 static constexpr type value = 1;
97};
98
99// The default DGF file to load
100template<class TypeTag>
101struct GridFile<TypeTag, TTag::OutflowBaseProblem> { static constexpr auto value = "./data/outflow.dgf"; };
102
103} // namespace Opm::Properties
104
105namespace Opm {
123template <class TypeTag>
124class OutflowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
125{
126 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
127
128 using GridView = GetPropType<TypeTag, Properties::GridView>;
129 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
130 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
131 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
132 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
133 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
134 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
135 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
136 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
137
138 // copy some indices for convenience
139 enum {
140 // Grid and world dimension
141 dim = GridView::dimension,
142 dimWorld = GridView::dimensionworld,
143
144 numPhases = FluidSystem::numPhases,
145
146 // component indices
147 H2OIdx = FluidSystem::H2OIdx,
148 N2Idx = FluidSystem::N2Idx
149 };
150
151 using CoordScalar = typename GridView::ctype;
152 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
153
154 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
155
156public:
160 OutflowProblem(Simulator& simulator)
161 : ParentType(simulator)
162 , eps_(1e-6)
163 { }
164
169 {
170 ParentType::finishInit();
171
172 temperature_ = 273.15 + 20;
173 FluidSystem::init(/*minT=*/temperature_ - 1, /*maxT=*/temperature_ + 2,
174 /*numT=*/3,
175 /*minp=*/0.8e5, /*maxp=*/2.5e5, /*nump=*/500);
176
177 // set parameters of porous medium
178 perm_ = this->toDimMatrix_(1e-10);
179 porosity_ = 0.4;
180 tortuosity_ = 0.28;
181 }
182
187
191 std::string name() const
192 { return "outflow"; }
193
198 {
199#ifndef NDEBUG
200 this->model().checkConservativeness();
201
202 // Calculate storage terms
203 EqVector storage;
204 this->model().globalStorage(storage);
205
206 // Write mass balance information for rank 0
207 if (this->gridView().comm().rank() == 0) {
208 std::cout << "Storage: " << storage << std::endl << std::flush;
209 }
210#endif // NDEBUG
211 }
212
218 template <class Context>
219 Scalar temperature(const Context& /*context*/,
220 unsigned /*spaceIdx*/,
221 unsigned /*timeIdx*/) const
222 { return temperature_; } // in [K]
223
229 template <class Context>
230 const DimMatrix& intrinsicPermeability(const Context& /*context*/,
231 unsigned /*spaceIdx*/,
232 unsigned /*timeIdx*/) const
233 { return perm_; }
234
240 template <class Context>
241 Scalar porosity(const Context& /*context*/,
242 unsigned /*spaceIdx*/,
243 unsigned /*timeIdx*/) const
244 { return porosity_; }
245
246#if 0
251 template <class Context>
252 Scalar tortuosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
253 { return tortuosity_; }
254
259 template <class Context>
260 Scalar dispersivity(const Context& context,
261 unsigned spaceIdx, unsigned timeIdx) const
262 { return 0; }
263#endif
264
266
271
275 template <class Context>
276 void boundary(BoundaryRateVector& values, const Context& context,
277 unsigned spaceIdx, unsigned timeIdx) const
278 {
279 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
280
281 if (onLeftBoundary_(globalPos)) {
282 Opm::CompositionalFluidState<Scalar, FluidSystem,
283 /*storeEnthalpy=*/false> fs;
284 initialFluidState_(fs, context, spaceIdx, timeIdx);
285 fs.setPressure(/*phaseIdx=*/0, fs.pressure(/*phaseIdx=*/0) + 1e5);
286
287 Scalar xlN2 = 2e-4;
288 fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, xlN2);
289 fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1 - xlN2);
290
291 typename FluidSystem::template ParameterCache<Scalar> paramCache;
292 paramCache.updateAll(fs);
293 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
294 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
295 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
296 }
297
298 // impose an freeflow boundary condition
299 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
300 }
301 else if (onRightBoundary_(globalPos)) {
302 Opm::CompositionalFluidState<Scalar, FluidSystem,
303 /*storeEnthalpy=*/false> fs;
304 initialFluidState_(fs, context, spaceIdx, timeIdx);
305
306 // impose an outflow boundary condition
307 values.setOutFlow(context, spaceIdx, timeIdx, fs);
308 }
309 else
310 // no flow on top and bottom
311 values.setNoFlow();
312 }
313
315
320
324 template <class Context>
325 void initial(PrimaryVariables& values,
326 const Context& context,
327 unsigned spaceIdx,
328 unsigned timeIdx) const
329 {
330 Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/false> fs;
331 initialFluidState_(fs, context, spaceIdx, timeIdx);
332
333 values.assignNaive(fs);
334 }
335
342 template <class Context>
343 void source(RateVector& rate,
344 const Context& /*context*/,
345 unsigned /*spaceIdx*/,
346 unsigned /*timeIdx*/) const
347 { rate = Scalar(0.0); }
348
350
351private:
352 bool onLeftBoundary_(const GlobalPosition& pos) const
353 { return pos[0] < eps_; }
354
355 bool onRightBoundary_(const GlobalPosition& pos) const
356 { return pos[0] > this->boundingBoxMax()[0] - eps_; }
357
358 template <class FluidState, class Context>
359 void initialFluidState_(FluidState& fs, const Context& context,
360 unsigned spaceIdx, unsigned timeIdx) const
361 {
362 Scalar T = temperature(context, spaceIdx, timeIdx);
363 // Scalar rho = FluidSystem::H2O::liquidDensity(T, /*pressure=*/1.5e5);
364 // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
365 // this->boundingBoxMax()[dim - 1];
366 // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
367 // this->boundingBoxMax()[dim - 1];
368
369 fs.setSaturation(/*phaseIdx=*/0, 1.0);
370 fs.setPressure(/*phaseIdx=*/0, 1e5 /* + rho*z */);
371 fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1.0);
372 fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, 0);
373 fs.setTemperature(T);
374
375 typename FluidSystem::template ParameterCache<Scalar> paramCache;
376 paramCache.updateAll(fs);
377 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
378 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
379 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
380 }
381 }
382
383 const Scalar eps_;
384
385 MaterialLawParams materialParams_;
386 DimMatrix perm_;
387 Scalar temperature_;
388 Scalar porosity_;
389 Scalar tortuosity_;
390};
391} // namespace Opm
392
393#endif
Problem where dissolved nitrogen is transported with the water phase from the left side to the right.
Definition outflowproblem.hh:125
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition outflowproblem.hh:325
void finishInit()
Definition outflowproblem.hh:168
void endTimeStep()
Definition outflowproblem.hh:197
OutflowProblem(Simulator &simulator)
Definition outflowproblem.hh:160
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition outflowproblem.hh:343
Scalar temperature(const Context &, unsigned, unsigned) const
Definition outflowproblem.hh:219
Scalar porosity(const Context &, unsigned, unsigned) const
Definition outflowproblem.hh:241
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition outflowproblem.hh:276
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition outflowproblem.hh:230
std::string name() const
Definition outflowproblem.hh:191
Definition outflowproblem.hh:51