GCC Code Coverage Report

Line	Branch	Exec	Source
1			//
2			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
3			// SPDX-License-Identifier: GPL-3.0-or-later
4			//
5			/*!
6			* \file
7			* \ingroup MPNCTests
8			* \brief Problem where hot, pure liquid water is injected from the left hand side into an initially
9			* isothermal domain.
10			*
11			* The water is fully evaporated by a strong heat source.
12			* A local thermal non-equilibrium model is used: i.e. two different (fluid, solid)
13			* temperatures are primary variables.
14			*
15			* \author Philipp Nuske
16			*/
17
18			#ifndef DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH
19			#define DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH
20
21			#include <dumux/common/properties.hh>
22			#include <dumux/common/parameters.hh>
23			#include <dumux/common/boundarytypes.hh>
24			#include <dumux/common/numeqvector.hh>
25
26			#include <dumux/porousmediumflow/problem.hh>
27			#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh>
28
29			namespace Dumux {
30
31			/*!
32			* \ingroup MPNCTests
33			* \brief Problem where water is injected from the left hand side into a porous media filled domain,
34			* which is supplied with energy from the right hand side to evaporate the water.
35			*/
36			template<class TypeTag>
37			class CombustionProblemOneComponent: public PorousMediumFlowProblem<TypeTag>
38			{
39			using ParentType = PorousMediumFlowProblem<TypeTag>;
40			using Scalar = GetPropType<TypeTag, Properties::Scalar>;
41			using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
42			using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
43			using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
44			using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
45
46			using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
47			using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
48			using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
49
50			using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
51			using SubControlVolume = typename FVElementGeometry::SubControlVolume;
52			using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
53			using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
54			using Element = typename GridView::template Codim<0>::Entity;
55			using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
56			using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
57			using FluidState = GetPropType<TypeTag, Properties::FluidState>;
58			using ParameterCache = typename FluidSystem::ParameterCache;
59
60			using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
61			using Indices = typename ModelTraits::Indices;
62
63			static constexpr int dimWorld = GridView::dimensionworld;
64			static constexpr int numComponents = ModelTraits::numFluidComponents();
65			static constexpr int s0Idx = Indices::s0Idx;
66			static constexpr int p0Idx = Indices::p0Idx;
67			static constexpr int conti00EqIdx = Indices::conti0EqIdx;
68			static constexpr int energyEq0Idx = Indices::energyEqIdx;
69			static constexpr int numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
70			static constexpr int numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
71			static constexpr int energyEqSolidIdx = energyEq0Idx + numEnergyEqFluid + numEnergyEqSolid - 1;
72			static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx;
73			static constexpr int nPhaseIdx = FluidSystem::nPhaseIdx;
74			static constexpr int wCompIdx = FluidSystem::H2OIdx;
75			static constexpr int nCompIdx = FluidSystem::N2Idx;
76
77			static constexpr auto numPhases = ModelTraits::numFluidPhases();
78
79			public:
80		✗	CombustionProblemOneComponent(std::shared_ptr<const GridGeometry> gridGeometry)
81		✗	: ParentType(gridGeometry)
82			{
83		✗	outputName_ = getParam<std::string>("Problem.Name");
84		✗	nRestart_ = getParam<Scalar>("Constants.nRestart");
85		✗	TInitial_ = getParam<Scalar>("InitialConditions.TInitial");
86		✗	TRight_ = getParam<Scalar>("InitialConditions.TRight");
87		✗	pnInitial_ = getParam<Scalar>("InitialConditions.pnInitial");
88		✗	TBoundary_ = getParam<Scalar>("BoundaryConditions.TBoundary");
89		✗	SwBoundary_ = getParam<Scalar>("BoundaryConditions.SwBoundary");
90		✗	SwOneComponentSys_= getParam<Scalar>("BoundaryConditions.SwOneComponentSys");
91		✗	massFluxInjectedPhase_ = getParam<Scalar>("BoundaryConditions.massFluxInjectedPhase");
92		✗	heatFluxFromRight_ = getParam<Scalar>("BoundaryConditions.heatFluxFromRight");
93		✗	coldTime_ =getParam<Scalar>("BoundaryConditions.coldTime");
94		✗	time_ = 0.0;
95		✗	}
96
97		✗	void setTime(Scalar time)
98		✗	{ time_ = time; }
99
100			void setGridVariables(std::shared_ptr<GridVariables> gridVariables)
101		✗	{ gridVariables_ = gridVariables; }
102
103			const GridVariables& gridVariables() const
104		✗	{ return *gridVariables_; }
105
106			/*!
107			* \name Problem Params
108			*/
109			// \{
110			/*!
111			* \brief The problem name.
112			*
113			* This is used as a prefix for files generated by the simulation.
114			*/
115			const std::string name() const
116		✗	{ return outputName_;}
117
118			// \}
119
120			/*!
121			* \brief Evaluates the source term for all balance equations within a given
122			* sub-control volume.
123			*
124			* \param element The finite element
125			* \param fvGeometry The finite volume geometry of the element
126			* \param elemVolVars The volume variables of the element
127			* \param scv The local index of the sub-control volume
128			*
129			* Positive values mean that mass is created, negative ones mean that it vanishes.
130			*/
131		✗	NumEqVector source(const Element &element,
132			const FVElementGeometry& fvGeometry,
133			const ElementVolumeVariables& elemVolVars,
134			const SubControlVolume &scv) const
135			{
136		✗	NumEqVector values(0.0);
137
138		✗	const auto& globalPos = scv.dofPosition();
139
140		✗	const Scalar volume = scv.volume();
141		✗	const Scalar numScv = fvGeometry.numScv(); // box: numSCV, cc:1
142
143		✗	if (time_ > coldTime_ )
144			{
145		✗	if (onRightBoundaryPorousMedium_(globalPos))
146			{
147			// Testing the location of a vertex, but function is called for each associated scv. Compensate for that
148		✗	values[energyEqSolidIdx] = heatFluxFromRight_ / volume / numScv;
149			}
150			}
151		✗	return values;
152			}
153
154			/*!
155			* \brief Specifies which kind of boundary condition should be
156			* used for which equation on a given boundary segment.
157			*
158			* \param globalPos The position for which the bc type should be evaluated
159			*/
160		✗	BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
161			{
162		✗	BoundaryTypes bcTypes;
163			// Default: Neumann
164		✗	bcTypes.setAllNeumann();
165
166		✗	if(onRightBoundary_(globalPos) ) {
167			bcTypes.setAllDirichlet();
168			}
169		✗	return bcTypes;
170			}
171
172			/*!
173			* \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
174			*
175			* \param globalPos The global position
176			*
177			* For this method, the \a values parameter stores primary variables.
178			*/
179			PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
180			{
181		✗	return initial_(globalPos);
182			}
183
184			/*!
185			* \brief Evaluates the boundary conditions for a Neumann boundary segment.
186			*
187			* \param element The finite element
188			* \param fvGeometry The finite-volume geometry in the box scheme
189			* \param elemVolVars The volume variables of the element
190			* \param elemFluxVarsCache Flux variables caches for all faces in stencil
191			* \param scvf The sub-control volume face
192			*
193			* For this method, the \a values parameter stores the mass flux
194			* in normal direction of each phase. Negative values mean influx.
195			*/
196		✗	NumEqVector neumann(const Element &element,
197			const FVElementGeometry& fvGeometry,
198			const ElementVolumeVariables& elemVolVars,
199			const ElementFluxVariablesCache& elemFluxVarsCache,
200			const SubControlVolumeFace& scvf) const
201			{
202		✗	NumEqVector values(0.0);
203
204		✗	const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
205		✗	const auto& scvIdx = scvf.insideScvIdx();
206		✗	const Scalar massFluxInjectedPhase = massFluxInjectedPhase_;
207
208		✗	FluidState fluidState;
209
210		✗	const Scalar pn = elemVolVars[scvIdx].pressure(nPhaseIdx);
211		✗	const Scalar pw = elemVolVars[scvIdx].pressure(wPhaseIdx);
212
213		✗	fluidState.setPressure(nPhaseIdx, pn);
214		✗	fluidState.setPressure(wPhaseIdx, pw);
215
216		✗	fluidState.setTemperature(TBoundary_);
217			ParameterCache dummyCache;
218		✗	fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
219		✗	fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
220			// compute density of injection phase
221		✗	const Scalar density = FluidSystem::density(fluidState,
222			dummyCache,
223			wPhaseIdx);
224		✗	fluidState.setDensity(wPhaseIdx, density);
225		✗	const Scalar molarDensity = FluidSystem::molarDensity(fluidState,
226			dummyCache,
227			wPhaseIdx);
228		✗	fluidState.setMolarDensity(wPhaseIdx, molarDensity);
229
230		✗	for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) {
231		✗	const Scalar h = FluidSystem::enthalpy(fluidState,
232			dummyCache,
233			phaseIdx);
234		✗	fluidState.setEnthalpy(phaseIdx, h);
235			}
236
237		✗	const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(wPhaseIdx);
238
239		✗	if (onLeftBoundary_(globalPos))
240			{
241		✗	values[conti00EqIdx + wCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, wCompIdx);
242		✗	values[conti00EqIdx + nCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, nCompIdx);
243		✗	values[energyEq0Idx] = - massFluxInjectedPhase * fluidState.enthalpy(wPhaseIdx);
244			}
245		✗	return values;
246			}
247
248			/*!
249			* \brief Evaluates the initial value for a control volume.
250			*
251			* For this method, the \a values parameter stores primary
252			* variables.
253			*
254			* \param globalPos The global position
255			*/
256			PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
257			{
258		✗	return initial_(globalPos);
259			}
260
261			private:
262			// the internal method for the initial condition
263		✗	PrimaryVariables initial_(const GlobalPosition &globalPos) const
264			{
265		✗	const Scalar curPos = globalPos[0];
266		✗	const Scalar slope = (SwBoundary_ - SwOneComponentSys_)/this->spatialParams().lengthPM();
267			Scalar S[numPhases];
268		✗	const Scalar thisSaturation = SwOneComponentSys_ + curPos * slope;
269
270		✗	S[wPhaseIdx] = SwBoundary_;
271		✗	if (inPM_(globalPos) )
272		✗	S[wPhaseIdx] = thisSaturation;
273		✗	S[nPhaseIdx] = 1. - S[wPhaseIdx];
274
275		✗	Scalar thisTemperature = TInitial_;
276		✗	if(onRightBoundary_(globalPos))
277		✗	thisTemperature = TRight_;
278
279		✗	FluidState fs;
280		✗	for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
281			{
282		✗	fs.setSaturation(phaseIdx, S[phaseIdx]);
283		✗	fs.setTemperature(thisTemperature);
284			}
285
286			//obtain pc according to saturation
287		✗	const int wettingPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
288		✗	const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
289		✗	const auto capPress = fm.capillaryPressures(fs, wettingPhaseIdx);
290			Scalar p[numPhases];
291
292			using std::abs;
293		✗	p[wPhaseIdx] = pnInitial_ - abs(capPress[wPhaseIdx]);
294		✗	p[nPhaseIdx] = p[wPhaseIdx] + abs(capPress[wPhaseIdx]);
295		✗	for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
296		✗	fs.setPressure(phaseIdx, p[phaseIdx]);
297
298		✗	fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
299		✗	fs.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
300		✗	fs.setMoleFraction(nPhaseIdx, wCompIdx, 1.0);
301		✗	fs.setMoleFraction(nPhaseIdx, nCompIdx, 0.0);
302
303		✗	const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx;
304
305			MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper;
306		✗	helper.solve(fs, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx});
307		✗	auto priVars = helper.getPrimaryVariables(fs);
308
309			// additionally set the temperature for thermal non-equilibrium for the phases
310		✗	priVars[energyEq0Idx] = thisTemperature;
311		✗	priVars[energyEqSolidIdx] = thisTemperature;
312		✗	return priVars;
313			}
314
315			/*!
316			* \brief Returns whether the tested position is on the left boundary of the domain.
317			*/
318			bool onLeftBoundary_(const GlobalPosition & globalPos) const
319		✗	{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_;}
320
321			/*!
322			* \brief Returns whether the tested position is on the right boundary of the domain.
323			*/
324			bool onRightBoundary_(const GlobalPosition & globalPos) const
325		✗	{ return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;}
326
327			/*!
328			* \brief Returns whether the tested position is in a specific region (right) in the domain
329			* \todo this needs to be more sophisticated in order to allow for meshes with nodes not directly on the boundary
330			*/
331			bool onRightBoundaryPorousMedium_(const GlobalPosition & globalPos) const
332			{
333			using std::abs;
334		✗	return (abs(globalPos[0] - (this->spatialParams().lengthPM())) < eps_ );
335			}
336
337			/*!
338			* \brief Returns whether the tested position is in the porous medium.
339			*/
340			bool inPM_(const GlobalPosition & globalPos) const
341		✗	{ return !this->spatialParams().inOutFlow(globalPos); }
342
343			private:
344			static constexpr Scalar eps_ = 1e-6;
345			int nTemperature_;
346			int nPressure_;
347			std::string outputName_;
348			int nRestart_;
349			Scalar TInitial_;
350			Scalar TRight_;
351
352			Scalar pnInitial_;
353
354			Dune::ParameterTree inputParameters_;
355			Scalar x_[numPhases][numComponents];
356
357			Scalar TBoundary_;
358			Scalar SwBoundary_;
359			Scalar SwOneComponentSys_;
360
361			Scalar massFluxInjectedPhase_;
362			Scalar heatFluxFromRight_;
363			Scalar lengthPM_;
364			Scalar coldTime_;
365
366			Scalar time_;
367			std::shared_ptr<GridVariables> gridVariables_;
368			};
369
370			} // end namespace Dumux
371
372			#endif
373