combustionproblem1c.hh 21.6 KB
Newer Older
1
/*****************************************************************************
2
 *   See the file COPYING for full copying permissions.                      *
3
4
5
6
7
8
9
10
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
11
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
12
13
14
15
16
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
17
/*!
Katharina Heck's avatar
Katharina Heck committed
18
19
 * \file
 * \ingroup MPNCTests
20
21
22
 * \brief Problem where hot, pure liquid water is injected from the left hand side into a initially
 *        isotherm domain. The water is fully evaporated by a strong heat source.
 *        A local thermal non-equilibrium model is used: i.e. two different (fluid, solid)
23
24
25
26
27
28
29
 *        temperatures are primary variables.
 *
 * \author Philipp Nuske
 */
#ifndef DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH
#define DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH

30
31
#include <dune/grid/onedgrid.hh>

Katharina Heck's avatar
Katharina Heck committed
32
#include <dumux/discretization/box/properties.hh>
33

Katharina Heck's avatar
Katharina Heck committed
34
35
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/mpnc/model.hh>
36
#include <dumux/porousmediumflow/mpnc/pressureformulation.hh>
37

38
39
#include <dumux/material/solidstates/compositionalsolidstate.hh>
#include <dumux/material/solidsystems/compositionalsolidphase.hh>
40
#include <dumux/material/components/constant.hh>
41

42
43
44
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>

Katharina Heck's avatar
Katharina Heck committed
45
#include "combustionspatialparams.hh"
46
#include "combustionfluidsystem.hh"
47

48
namespace Dumux {
49

50
template<class TypeTag>
51
52
class CombustionProblemOneComponent;

53
//! Custom model traits to deactivate diffusion for this test
54
55
template<int numP, int numC, MpNcPressureFormulation formulation, bool useM>
struct CombustionModelTraits : public MPNCModelTraits<numP, numC, formulation, useM>
56
57
58
59
{
    static constexpr bool enableMolecularDiffusion() { return false; }
};

60
61
namespace Properties {
NEW_TYPE_TAG(CombustionOneComponentTypeTag, INHERITS_FROM(MPNCNonequil));
62
NEW_TYPE_TAG(CombustionOneComponentBoxTypeTag, INHERITS_FROM(BoxModel, CombustionOneComponentTypeTag));
63
64

// Set the grid type
65
SET_TYPE_PROP(CombustionOneComponentTypeTag, Grid, Dune::OneDGrid);
66
67

// Set the problem property
68
SET_TYPE_PROP(CombustionOneComponentTypeTag,
69
70
              Problem,
              CombustionProblemOneComponent<TypeTag>);
71

72

73

74
SET_TYPE_PROP(CombustionOneComponentTypeTag,
Katharina Heck's avatar
Katharina Heck committed
75
              FluidSystem,
76
              FluidSystems::CombustionFluidsystem<typename GET_PROP_TYPE(TypeTag, Scalar)>);
77
78

//! Set the default pressure formulation: either pw first or pn first
79
80
81
82
83
SET_PROP(CombustionOneComponentTypeTag, PressureFormulation)
{
public:
    static const MpNcPressureFormulation value = MpNcPressureFormulation::mostWettingFirst;
};
84
85

// Set the type used for scalar values
86
SET_TYPE_PROP(CombustionOneComponentTypeTag, Scalar, double );
87
// quad / double
88

89
90
91
92
93
94
// We use different model traits for the equilibrium part because we want to deactivate diffusion
SET_PROP(CombustionOneComponentTypeTag, EquilibriumModelTraits)
{
private:
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
public:
95
96
    using type = CombustionModelTraits< FluidSystem::numPhases,
                                        FluidSystem::numComponents,
97
98
                                        GET_PROP_VALUE(TypeTag, PressureFormulation),
                                        GET_PROP_VALUE(TypeTag, UseMoles) >;
Katharina Heck's avatar
Katharina Heck committed
99
};
100

101
SET_PROP(CombustionOneComponentTypeTag, FluidState)
102
{
103
private:
104
105
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
106
public:
107
    using type = CompositionalFluidState<Scalar, FluidSystem>;
108
};
Katharina Heck's avatar
Katharina Heck committed
109
110
111
//#################
//changes from the default settings which also assume chemical non-equilibrium
//set the number of energyequations we want to use
112
113
SET_INT_PROP(CombustionOneComponentTypeTag, NumEnergyEqFluid, 1);
SET_INT_PROP(CombustionOneComponentTypeTag, NumEnergyEqSolid, 1);
Katharina Heck's avatar
Katharina Heck committed
114
115

// by default chemical non equilibrium is enabled in the nonequil model, switch that off here
116
SET_BOOL_PROP(CombustionOneComponentTypeTag, EnableChemicalNonEquilibrium, false);
Katharina Heck's avatar
Katharina Heck committed
117
//#################
118

119
120
121
SET_PROP(CombustionOneComponentTypeTag, SolidSystem)
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
122
123
124
125
    using ComponentOne = Dumux::Components::Constant<1, Scalar>;
    using ComponentTwo = Dumux::Components::Constant<2, Scalar>;
    static constexpr int numInertComponents = 2;
    using type = SolidSystems::CompositionalSolidPhase<Scalar, ComponentOne, ComponentTwo, numInertComponents>;
126
};
127
128
129
130
131
132
133
134
135
136
137
138
139

SET_PROP(CombustionOneComponentTypeTag, SolidState)
{
private:
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using SolidSystem = typename GET_PROP_TYPE(TypeTag, SolidSystem);
public:
    using type = CompositionalSolidState<Scalar, SolidSystem>;
};

// Set the spatial parameters
SET_TYPE_PROP(CombustionOneComponentTypeTag, SpatialParams, CombustionSpatialParams<TypeTag>);

140
141
}
/*!
Katharina Heck's avatar
Katharina Heck committed
142
 * \ingroup MPNCTests
143
 * \brief Problem where water is injected from the left hand side into a porous media filled domain,
144
 *        which is supplied with energy from the right hand side to evaporate the water.
145
 */
146
template<class TypeTag>
Katharina Heck's avatar
Katharina Heck committed
147
148
149
class CombustionProblemOneComponent: public PorousMediumFlowProblem<TypeTag>
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
150
151
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
Katharina Heck's avatar
Katharina Heck committed
152
153
    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
154
    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
155
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
156
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
157
158
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
Katharina Heck's avatar
Katharina Heck committed
159
160
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
161
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
Katharina Heck's avatar
Katharina Heck committed
162
    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
163
164
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
    using ParameterCache = typename FluidSystem::ParameterCache;
Katharina Heck's avatar
Katharina Heck committed
165
166
    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);

167
168
169
    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
    using Indices = typename ModelTraits::Indices;

170
    enum {dimWorld = GridView::dimensionworld};
171
172
    enum {numPhases = ModelTraits::numPhases()};
    enum {numComponents = ModelTraits::numComponents()};
173
174
    enum {s0Idx = Indices::s0Idx};
    enum {p0Idx = Indices::p0Idx};
175
176
    enum {conti00EqIdx = Indices::conti0EqIdx};
    enum {energyEq0Idx = Indices::energyEqIdx};
177
178
    enum {numEnergyEqFluid = ModelTraits::numEnergyEqFluid()};
    enum {numEnergyEqSolid = ModelTraits::numEnergyEqSolid()};
Katharina Heck's avatar
Katharina Heck committed
179
    enum {energyEqSolidIdx = energyEq0Idx + numEnergyEqFluid + numEnergyEqSolid - 1};
180
181
182
183
    enum {wPhaseIdx = FluidSystem::wPhaseIdx};
    enum {nPhaseIdx = FluidSystem::nPhaseIdx};
    enum {wCompIdx = FluidSystem::H2OIdx};
    enum {nCompIdx = FluidSystem::N2Idx};
Katharina Heck's avatar
Katharina Heck committed
184

185
    // formulations
186
187
188
    static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
    static constexpr auto mostWettingFirst = MpNcPressureFormulation::mostWettingFirst;
    static constexpr auto leastWettingFirst = MpNcPressureFormulation::leastWettingFirst;
189
190

public:
Katharina Heck's avatar
Katharina Heck committed
191
192
    CombustionProblemOneComponent(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
        : ParentType(fvGridGeometry)
193
    {
Katharina Heck's avatar
Katharina Heck committed
194
195
196
197
198
199
200
201
202
203
204
205
            outputName_ = getParam<std::string>("Problem.Name");
            nRestart_ = getParam<Scalar>("Constants.nRestart");
            TInitial_ = getParam<Scalar>("InitialConditions.TInitial");
            TRight_ = getParam<Scalar>("InitialConditions.TRight");
            pnInitial_ = getParam<Scalar>("InitialConditions.pnInitial");
            TBoundary_ = getParam<Scalar>("BoundaryConditions.TBoundary");
            SwBoundary_ = getParam<Scalar>("BoundaryConditions.SwBoundary");
            SwOneComponentSys_= getParam<Scalar>("BoundaryConditions.SwOneComponentSys");
            massFluxInjectedPhase_ = getParam<Scalar>("BoundaryConditions.massFluxInjectedPhase");
            heatFluxFromRight_ = getParam<Scalar>("BoundaryConditions.heatFluxFromRight");
            coldTime_ =getParam<Scalar>("BoundaryConditions.coldTime");
            time_ = 0.0;
206
207
    }

Katharina Heck's avatar
Katharina Heck committed
208
209
210
211
212
213
214
215
    void setTime(Scalar time)
    { time_ = time; }

    void setGridVariables(std::shared_ptr<GridVariables> gridVariables)
    { gridVariables_ = gridVariables; }

     const GridVariables& gridVariables() const
    { return *gridVariables_; }
216
217
218
219
220
221
222
223
224

    /*!
     * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
     *
     * This is not specific to the discretization.
     */
    Scalar temperature() const
    {   return TInitial_;}

Katharina Heck's avatar
Katharina Heck committed
225

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
    /*!
     * \name Problem Params
     */
    // \{
    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
    const std::string name() const
    {   return outputName_;}

    // \}

    /*!
Katharina Heck's avatar
Katharina Heck committed
241
     * \brief Evaluate the source term for all balance equations within a given
242
243
     *        sub-control-volume.
     *
Katharina Heck's avatar
Katharina Heck committed
244
245
     * \param values Stores the solution for the conservation equations in
     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
246
247
248
     * \param element The finite element
     * \param fvGeometry The finite volume geometry of the element
     * \param scvIdx The local index of the sub-control volume
Katharina Heck's avatar
Katharina Heck committed
249
250
     *
     * Positive values mean that mass is created, negative ones mean that it vanishes.
251
     */
Katharina Heck's avatar
Katharina Heck committed
252
    //! \copydoc Dumux::ImplicitProblem::source()
253
    NumEqVector source(const Element &element,
Katharina Heck's avatar
Katharina Heck committed
254
255
256
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolume &scv) const
257
    {
258
        NumEqVector values(0.0);
Katharina Heck's avatar
Katharina Heck committed
259
260

        const auto& globalPos = scv.dofPosition();
261

Katharina Heck's avatar
Katharina Heck committed
262
263
        const Scalar volume = scv.volume();
        const Scalar numScv = fvGeometry.numScv(); // box: numSCV, cc:1
264

Katharina Heck's avatar
Katharina Heck committed
265
        if (time_ > coldTime_ )
266
267
268
        {
            if (onRightBoundaryPorousMedium_(globalPos))
            {
Katharina Heck's avatar
Katharina Heck committed
269
                // Testing the location of a vertex, but function is called for each associated scv. Compensate for that
270
                values[energyEqSolidIdx] = heatFluxFromRight_ / volume / numScv;
271
            }
Katharina Heck's avatar
Katharina Heck committed
272
         }
273
        return values;
274
275
276
277
278
279
    }

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
Katharina Heck's avatar
Katharina Heck committed
280
281
     * \param values The boundary types for the conservation equations
     * \param globalPos The position for which the bc type should be evaluated
282
     */
Katharina Heck's avatar
Katharina Heck committed
283
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
284
    {
Katharina Heck's avatar
Katharina Heck committed
285
        BoundaryTypes bcTypes;
286
        // Default: Neumann
Katharina Heck's avatar
Katharina Heck committed
287
         bcTypes.setAllNeumann();
288

Katharina Heck's avatar
Katharina Heck committed
289
290
291
292
         if(onRightBoundary_(globalPos) ) {
            bcTypes.setAllDirichlet();
         }
        return bcTypes;
293
294
295
296
297
298
299
300
301
302
303
304
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
     *        boundary segment.
     *
     * \param priVars Stores the Dirichlet values for the conservation equations in
     *               \f$ [ \textnormal{unit of primary variable} ] \f$
     * \param globalPos The global position
     *
     * For this method, the \a values parameter stores primary variables.
     */
Katharina Heck's avatar
Katharina Heck committed
305
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
306
    {
Katharina Heck's avatar
Katharina Heck committed
307
       return  initial_(globalPos);
308
309
310
311
312
313
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
Katharina Heck's avatar
Katharina Heck committed
314
     * \param values The neumann values for the conservation equations
315
     * \param element The finite element
Katharina Heck's avatar
Katharina Heck committed
316
     * \param fvGeometry The finite-volume geometry in the box scheme
317
     * \param intersection The intersection between element and boundary
Katharina Heck's avatar
Katharina Heck committed
318
     * \param scvIdx The local vertex index
319
     * \param boundaryFaceIdx The index of the boundary face
Katharina Heck's avatar
Katharina Heck committed
320
321
322
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
323
     */
324
    NumEqVector neumann(const Element &element,
Katharina Heck's avatar
Katharina Heck committed
325
326
327
328
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolumeFace& scvf) const
    {
329
        NumEqVector values(0.0);
330

Katharina Heck's avatar
Katharina Heck committed
331
332
        const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
        const auto& scvIdx = scvf.insideScvIdx();
333
334
335
336
        const Scalar massFluxInjectedPhase = massFluxInjectedPhase_;

        FluidState fluidState;

337
338
        const Scalar pn = elemVolVars[scvIdx].pressure(nPhaseIdx);
        const Scalar pw = elemVolVars[scvIdx].pressure(wPhaseIdx);
339
340
341
342

        fluidState.setPressure(nPhaseIdx, pn);
        fluidState.setPressure(wPhaseIdx, pw);

343
        fluidState.setTemperature(TBoundary_);
344
345
346
347
348
        ParameterCache dummyCache;
        fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
        fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
        // compute density of injection phase
        const Scalar density = FluidSystem::density(fluidState,
Simon Scholz's avatar
Simon Scholz committed
349
350
                                                    dummyCache,
                                                    wPhaseIdx);
351
        fluidState.setDensity(wPhaseIdx, density);
352
        const Scalar molarDensity = FluidSystem::molarDensity(fluidState,
Simon Scholz's avatar
Simon Scholz committed
353
354
                                                              dummyCache,
                                                              wPhaseIdx);
355
        fluidState.setMolarDensity(wPhaseIdx, molarDensity);
356
357
358

        for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) {
            const Scalar h = FluidSystem::enthalpy(fluidState,
Simon Scholz's avatar
Simon Scholz committed
359
360
                                                   dummyCache,
                                                   phaseIdx);
361
362
363
364
365
366
367
            fluidState.setEnthalpy(phaseIdx, h);
        }

        const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(wPhaseIdx);

        if (onLeftBoundary_(globalPos))
        {
368
369
            values[conti00EqIdx + wCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, wCompIdx);
            values[conti00EqIdx + nCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, nCompIdx);
370
            values[energyEq0Idx] = - massFluxInjectedPhase * fluidState.enthalpy(wPhaseIdx);
371
        }
372
        return values;
373
374
375
376
377
378
379
380
381
382
383
384
    }

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * For this method, the \a values parameter stores primary
     * variables.
     *
     * \param values Stores the Neumann values for the conservation equations in
     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
     * \param globalPos The global position
     */
Katharina Heck's avatar
Katharina Heck committed
385
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
386
    {
Katharina Heck's avatar
Katharina Heck committed
387
       return initial_(globalPos);
388
    }
389
390

private:
391
    // the internal method for the initial condition
Katharina Heck's avatar
Katharina Heck committed
392
    PrimaryVariables initial_(const GlobalPosition &globalPos) const
393
    {
Katharina Heck's avatar
Katharina Heck committed
394
        PrimaryVariables priVars(0.0);
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
        const Scalar curPos = globalPos[0];
        const Scalar slope = (SwBoundary_-SwOneComponentSys_) / (this->spatialParams().lengthPM());
        Scalar S[numPhases];
        const Scalar thisSaturation = SwOneComponentSys_ + curPos * slope;

        S[wPhaseIdx] = SwBoundary_;
        if (inPM_(globalPos) ) {
            S[wPhaseIdx] = thisSaturation;
        }

        S[nPhaseIdx] = 1. - S[wPhaseIdx];

        //////////////////////////////////////
        // Set saturation
        //////////////////////////////////////
        for (int i = 0; i < numPhases - 1; ++i) {
411
            priVars[s0Idx + i] = S[i];
412
413
414
415
416
417
418
419
420
421
422
423
424
        }

        FluidState fluidState;

        Scalar thisTemperature = TInitial_;
        if(onRightBoundary_(globalPos))
        thisTemperature = TRight_;

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            fluidState.setSaturation(phaseIdx, S[phaseIdx]);

            fluidState.setTemperature(thisTemperature );

Katharina Heck's avatar
Katharina Heck committed
425
        }
426
427
428
        //////////////////////////////////////
        // Set temperature
        //////////////////////////////////////
Katharina Heck's avatar
Katharina Heck committed
429
430
        priVars[energyEq0Idx] = thisTemperature;
        priVars[energyEqSolidIdx] = thisTemperature;
431
432
433
        Scalar capPress[numPhases];

        //obtain pc according to saturation
Katharina Heck's avatar
Katharina Heck committed
434
435
        const auto &materialParams =
        this->spatialParams().materialLawParamsAtPos(globalPos);
436
        using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
Katharina Heck's avatar
Katharina Heck committed
437
         MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
438
439
440

        Scalar p[numPhases];

441
442
443
        using std::abs;
        p[wPhaseIdx] = pnInitial_ - abs(capPress[wPhaseIdx]);
        p[nPhaseIdx] = p[wPhaseIdx] + abs(capPress[wPhaseIdx]);
444
445
446
447
448
449
450
451
452
453

        for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
        fluidState.setPressure(phaseIdx, p[phaseIdx]);

        //////////////////////////////////////
        // Set pressure
        //////////////////////////////////////
        if(pressureFormulation == mostWettingFirst) {
            // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pw
454
            priVars[p0Idx] = p[wPhaseIdx];
455
456
457
458
        }
        else if(pressureFormulation == leastWettingFirst) {
            // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pn
459
            priVars[p0Idx] = p[nPhaseIdx];
460
        }
461
462
        else
            DUNE_THROW(Dune::InvalidStateException, "CombustionProblemOneComponent does not support the chosen pressure formulation.");
463
464
465
466
467
468
469

        fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
        fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);

        fluidState.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
        fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);

Katharina Heck's avatar
Katharina Heck committed
470
       int refPhaseIdx;
471

Katharina Heck's avatar
Katharina Heck committed
472
473
        // on right boundary: reference is gas
        refPhaseIdx = nPhaseIdx;
474

Katharina Heck's avatar
Katharina Heck committed
475
476
477
        if(inPM_(globalPos)) {
            refPhaseIdx = wPhaseIdx;
        }
478

Katharina Heck's avatar
Katharina Heck committed
479
480
481
482
        // obtain fugacities
        using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
        ParameterCache paramCache;
        ComputeFromReferencePhase::solve(fluidState,
483
                                         paramCache,
484
                                         refPhaseIdx);
Katharina Heck's avatar
Katharina Heck committed
485
486
487
488
489
490
491
492

        //////////////////////////////////////
        // Set fugacities
        //////////////////////////////////////
        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
            priVars[conti00EqIdx + compIdx] = fluidState.fugacity(refPhaseIdx,compIdx);
        }
        return priVars;
493
494
495
496
497
498
    }

    /*!
     * \brief Give back whether the tested position (input) is a specific region (left) in the domain
     */
    bool onLeftBoundary_(const GlobalPosition & globalPos) const
Katharina Heck's avatar
Katharina Heck committed
499
    {   return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;}
500
501
502
503
504

    /*!
     * \brief Give back whether the tested position (input) is a specific region (right) in the domain
     */
    bool onRightBoundary_(const GlobalPosition & globalPos) const
Katharina Heck's avatar
Katharina Heck committed
505
    {   return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;}
506
507
508
509
510
511

    /*!
     * \brief Give back whether the tested position (input) is a specific region (right) in the domain
     *  \todo this needs to be more sophisticated in order to allow for meshes with nodes not directly on the boundary
     */
    bool onRightBoundaryPorousMedium_(const GlobalPosition & globalPos) const
512
513
    {
        using std::abs;
Katharina Heck's avatar
Katharina Heck committed
514
        return (abs(globalPos[0] - (this->spatialParams().lengthPM())) < eps_ );
515
    }
516
517
518
519
520

    /*!
     * \brief Give back whether the tested position (input) is a specific region (right) in the domain
     */
    bool inPM_(const GlobalPosition & globalPos) const
521
    { return !this->spatialParams().inOutFlow(globalPos); }
522

523
private:
524
    static constexpr Scalar eps_ = 1e-6;
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
    int nTemperature_;
    int nPressure_;
    std::string outputName_;
    int nRestart_;
    Scalar TInitial_;
    Scalar TRight_;

    Scalar pnInitial_;

    Dune::ParameterTree inputParameters_;
    Scalar x_[numPhases][numComponents];

    Scalar TBoundary_;
    Scalar SwBoundary_;
    Scalar SwOneComponentSys_;

    Scalar massFluxInjectedPhase_;
    Scalar heatFluxFromRight_;
    Scalar lengthPM_;
    Scalar coldTime_;
545

Katharina Heck's avatar
Katharina Heck committed
546
547
    Scalar time_;
    std::shared_ptr<GridVariables> gridVariables_;
548
549
};

Katharina Heck's avatar
Katharina Heck committed
550
} //end namespace
551
552

#endif