My Project
Loading...
Searching...
No Matches
diffusionproblem.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_POWER_INJECTION_PROBLEM_HH
29#define EWOMS_POWER_INJECTION_PROBLEM_HH
30
31#include <opm/models/ncp/ncpproperties.hh>
32
33#include <opm/models/io/cubegridvanguard.hh>
34
35#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
36#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
37#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
38#include <opm/material/fluidstates/CompositionalFluidState.hpp>
39#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
40
41#include <dune/grid/yaspgrid.hh>
42#include <dune/common/version.hh>
43#include <dune/common/fvector.hh>
44#include <dune/common/fmatrix.hh>
45
46#include <sstream>
47#include <string>
48
49namespace Opm {
50template <class TypeTag>
51class DiffusionProblem;
52}
53
54namespace Opm::Properties {
55
56namespace TTag {
57
59
60} // namespace TTag
61
62// Set the grid implementation to be used
63template<class TypeTag>
64struct Grid<TypeTag, TTag::DiffusionBaseProblem> { using type = Dune::YaspGrid</*dim=*/1>; };
65
66// set the Vanguard property
67template<class TypeTag>
68struct Vanguard<TypeTag, TTag::DiffusionBaseProblem> { using type = Opm::CubeGridVanguard<TypeTag>; };
69
70// Set the problem property
71template<class TypeTag>
72struct Problem<TypeTag, TTag::DiffusionBaseProblem> { using type = Opm::DiffusionProblem<TypeTag>; };
73
74// Set the fluid system
75template<class TypeTag>
76struct FluidSystem<TypeTag, TTag::DiffusionBaseProblem>
77{
78private:
79 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
80
81public:
82 using type = Opm::H2ON2FluidSystem<Scalar>;
83};
84
85// Set the material Law
86template<class TypeTag>
87struct MaterialLaw<TypeTag, TTag::DiffusionBaseProblem>
88{
89private:
90 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
92
93 static_assert(FluidSystem::numPhases == 2,
94 "A fluid system with two phases is required "
95 "for this problem!");
96
97 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
98 /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
99 /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>;
100
101public:
102 using type = Opm::LinearMaterial<Traits>;
103};
104
105// Enable molecular diffusion for this problem
106template<class TypeTag>
107struct EnableDiffusion<TypeTag, TTag::DiffusionBaseProblem> { static constexpr bool value = true; };
108
109// Disable gravity
110template<class TypeTag>
111struct EnableGravity<TypeTag, TTag::DiffusionBaseProblem> { static constexpr bool value = false; };
112
113// define the properties specific for the diffusion problem
114template<class TypeTag>
115struct DomainSizeX<TypeTag, TTag::DiffusionBaseProblem>
116{
117 using type = GetPropType<TypeTag, Scalar>;
118 static constexpr type value = 1.0;
119};
120template<class TypeTag>
121struct DomainSizeY<TypeTag, TTag::DiffusionBaseProblem>
122{
123 using type = GetPropType<TypeTag, Scalar>;
124 static constexpr type value = 1.0;
125};
126template<class TypeTag>
127struct DomainSizeZ<TypeTag, TTag::DiffusionBaseProblem>
128{
129 using type = GetPropType<TypeTag, Scalar>;
130 static constexpr type value = 1.0;
131};
132
133template<class TypeTag>
134struct CellsX<TypeTag, TTag::DiffusionBaseProblem> { static constexpr unsigned value = 250; };
135template<class TypeTag>
136struct CellsY<TypeTag, TTag::DiffusionBaseProblem> { static constexpr unsigned value = 1; };
137template<class TypeTag>
138struct CellsZ<TypeTag, TTag::DiffusionBaseProblem> { static constexpr unsigned value = 1; };
139
140// The default for the end time of the simulation
141template<class TypeTag>
142struct EndTime<TypeTag, TTag::DiffusionBaseProblem>
143{
144 using type = GetPropType<TypeTag, Scalar>;
145 static constexpr type value = 1e6;
146};
147
148// The default for the initial time step size of the simulation
149template<class TypeTag>
150struct InitialTimeStepSize<TypeTag, TTag::DiffusionBaseProblem>
151{
152 using type = GetPropType<TypeTag, Scalar>;
153 static constexpr type value = 1000;
154};
155
156} // namespace Opm::Properties
157
158namespace Opm {
169template <class TypeTag>
170class DiffusionProblem : public GetPropType<TypeTag, Properties::BaseProblem>
171{
172 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
173
174 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
175 using GridView = GetPropType<TypeTag, Properties::GridView>;
176 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
177 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
178 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
179 using Model = GetPropType<TypeTag, Properties::Model>;
180
181 enum {
182 // number of phases
183 numPhases = FluidSystem::numPhases,
184
185 // phase indices
186 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
187 gasPhaseIdx = FluidSystem::gasPhaseIdx,
188
189 // component indices
190 H2OIdx = FluidSystem::H2OIdx,
191 N2Idx = FluidSystem::N2Idx,
192
193 // Grid and world dimension
194 dim = GridView::dimension,
195 dimWorld = GridView::dimensionworld
196 };
197
198 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
199 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
200 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
201
202 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
203 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
204
205 using CoordScalar = typename GridView::ctype;
206 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
207
208 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
209
210public:
214 DiffusionProblem(Simulator& simulator)
215 : ParentType(simulator)
216 { }
217
222 {
223 ParentType::finishInit();
224
225 FluidSystem::init();
226
227 temperature_ = 273.15 + 20.0;
228
229 materialParams_.finalize();
230
231 K_ = this->toDimMatrix_(1e-12); // [m^2]
232
233 setupInitialFluidStates_();
234 }
235
240
244 std::string name() const
245 { return std::string("diffusion_") + Model::name(); }
246
251 {
252#ifndef NDEBUG
253 this->model().checkConservativeness();
254
255 // Calculate storage terms
256 EqVector storage;
257 this->model().globalStorage(storage);
258
259 // Write mass balance information for rank 0
260 if (this->gridView().comm().rank() == 0) {
261 std::cout << "Storage: " << storage << std::endl << std::flush;
262 }
263#endif // NDEBUG
264 }
265
267
272
276 template <class Context>
277 const DimMatrix& intrinsicPermeability(const Context& /*context*/,
278 unsigned /*spaceIdx*/,
279 unsigned /*timeIdx*/) const
280 { return K_; }
281
285 template <class Context>
286 Scalar porosity(const Context& /*context*/,
287 unsigned /*spaceIdx*/,
288 unsigned /*timeIdx*/) const
289 { return 0.35; }
290
294 template <class Context>
295 const MaterialLawParams&
296 materialLawParams(const Context& /*context*/,
297 unsigned /*spaceIdx*/,
298 unsigned /*timeIdx*/) const
299 { return materialParams_; }
300
304 template <class Context>
305 Scalar temperature(const Context& /*context*/,
306 unsigned /*spaceIdx*/,
307 unsigned /*timeIdx*/) const
308 { return temperature_; }
309
311
316
322 template <class Context>
323 void boundary(BoundaryRateVector& values,
324 const Context& /*context*/,
325 unsigned /*spaceIdx*/,
326 unsigned /*timeIdx*/) const
327 { values.setNoFlow(); }
328
330
335
339 template <class Context>
340 void initial(PrimaryVariables& values,
341 const Context& context,
342 unsigned spaceIdx,
343 unsigned timeIdx) const
344 {
345 const auto& pos = context.pos(spaceIdx, timeIdx);
346 if (onLeftSide_(pos))
347 values.assignNaive(leftInitialFluidState_);
348 else
349 values.assignNaive(rightInitialFluidState_);
350 }
351
358 template <class Context>
359 void source(RateVector& rate,
360 const Context& /*context*/,
361 unsigned /*spaceIdx*/,
362 unsigned /*timeIdx*/) const
363 { rate = Scalar(0.0); }
364
366
367private:
368 bool onLeftSide_(const GlobalPosition& pos) const
369 { return pos[0] < (this->boundingBoxMin()[0] + this->boundingBoxMax()[0]) / 2; }
370
371 void setupInitialFluidStates_()
372 {
373 // create the initial fluid state for the left half of the domain
374 leftInitialFluidState_.setTemperature(temperature_);
375
376 Scalar Sl = 0.0;
377 leftInitialFluidState_.setSaturation(liquidPhaseIdx, Sl);
378 leftInitialFluidState_.setSaturation(gasPhaseIdx, 1 - Sl);
379
380 Scalar p = 1e5;
381 leftInitialFluidState_.setPressure(liquidPhaseIdx, p);
382 leftInitialFluidState_.setPressure(gasPhaseIdx, p);
383
384 Scalar xH2O = 0.01;
385 leftInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
386 leftInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
387
388 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
389 typename FluidSystem::template ParameterCache<Scalar> paramCache;
390 CFRP::solve(leftInitialFluidState_, paramCache, gasPhaseIdx,
391 /*setViscosity=*/false, /*setEnthalpy=*/false);
392
393 // create the initial fluid state for the right half of the domain
394 rightInitialFluidState_.assign(leftInitialFluidState_);
395 xH2O = 0.0;
396 rightInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
397 rightInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
398 CFRP::solve(rightInitialFluidState_, paramCache, gasPhaseIdx,
399 /*setViscosity=*/false, /*setEnthalpy=*/false);
400 }
401
402 DimMatrix K_;
403 MaterialLawParams materialParams_;
404
405 Opm::CompositionalFluidState<Scalar, FluidSystem> leftInitialFluidState_;
406 Opm::CompositionalFluidState<Scalar, FluidSystem> rightInitialFluidState_;
407 Scalar temperature_;
408};
409
410} // namespace Opm
411
412#endif
1D problem which is driven by molecular diffusion.
Definition diffusionproblem.hh:171
Scalar temperature(const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:305
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition diffusionproblem.hh:340
void finishInit()
Definition diffusionproblem.hh:221
Scalar porosity(const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:286
void boundary(BoundaryRateVector &values, const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:323
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:359
DiffusionProblem(Simulator &simulator)
Definition diffusionproblem.hh:214
void endTimeStep()
Definition diffusionproblem.hh:250
std::string name() const
Definition diffusionproblem.hh:244
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:296
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition diffusionproblem.hh:277
Definition diffusionproblem.hh:58