My Project
Loading...
Searching...
No Matches
FlowProblemComp.hpp
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 Copyright 2024 SINTEF Digital
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 2 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20
21 Consult the COPYING file in the top-level source directory of this
22 module for the precise wording of the license and the list of
23 copyright holders.
24*/
30#ifndef OPM_FLOW_PROBLEM_COMP_HPP
31#define OPM_FLOW_PROBLEM_COMP_HPP
32
33
35#include <opm/material/fluidstates/CompositionalFluidState.hpp>
36
37#include <opm/material/thermal/EclThermalLawManager.hpp>
38
39
40#include <algorithm>
41#include <functional>
42#include <set>
43#include <string>
44#include <vector>
45
46namespace Opm {
47
54template <class TypeTag>
55class FlowProblemComp : public FlowProblem<TypeTag>
56{
57 // TODO: the naming of the Types will be adjusted
59
60 using typename FlowProblemType::Scalar;
61 using typename FlowProblemType::Simulator;
62 using typename FlowProblemType::GridView;
63 using typename FlowProblemType::FluidSystem;
64 using typename FlowProblemType::Vanguard;
65
66 // might not be needed
67 using FlowProblemType::dim;
68 using FlowProblemType::dimWorld;
69
70 using FlowProblemType::numPhases;
71 using FlowProblemType::numComponents;
72
73 using FlowProblemType::gasPhaseIdx;
74 using FlowProblemType::oilPhaseIdx;
75 using FlowProblemType::waterPhaseIdx;
76
77 using typename FlowProblemType::Indices;
78 using typename FlowProblemType::PrimaryVariables;
80 using typename FlowProblemType::Evaluation;
81 using typename FlowProblemType::MaterialLaw;
82 using typename FlowProblemType::RateVector;
83
84 using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
85
86public:
89
93 static void registerParameters()
94 {
96
97 // tighter tolerance is needed for compositional modeling here
98 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-7);
99 }
100
101
105 explicit FlowProblemComp(Simulator& simulator)
106 : FlowProblemType(simulator)
107 {
108 }
109
114 {
115 // TODO: there should be room to remove duplication for this function,
116 // but there is relatively complicated logic in the function calls in this function
117 // some refactoring is needed for this function
118 FlowProblemType::finishInit();
119
120 auto& simulator = this->simulator();
121
122 auto finishTransmissibilities = [updated = false, this]() mutable {
123 if (updated) {
124 return;
125 }
126 this->transmissibilities_.finishInit(
127 [&vg = this->simulator().vanguard()](const unsigned int it) { return vg.gridIdxToEquilGridIdx(it); });
128 updated = true;
129 };
130
132
133 const auto& eclState = simulator.vanguard().eclState();
134 const auto& schedule = simulator.vanguard().schedule();
135
136 // Set the start time of the simulation
137 simulator.setStartTime(schedule.getStartTime());
138 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
139
140 // We want the episode index to be the same as the report step index to make
141 // things simpler, so we have to set the episode index to -1 because it is
142 // incremented by endEpisode(). The size of the initial time step and
143 // length of the initial episode is set to zero for the same reason.
144 simulator.setEpisodeIndex(-1);
145 simulator.setEpisodeLength(0.0);
146
147 // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
148 // disables gravity, else the standard value of the gravity constant at sea level
149 // on earth is used
150 this->gravity_ = 0.0;
151 if (Parameters::Get<Parameters::EnableGravity>())
152 this->gravity_[dim - 1] = 9.80665;
153 if (!eclState.getInitConfig().hasGravity())
154 this->gravity_[dim - 1] = 0.0;
155
156 if (this->enableTuning_) {
157 // if support for the TUNING keyword is enabled, we get the initial time
158 // steping parameters from it instead of from command line parameters
159 const auto& tuning = schedule[0].tuning();
160 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
161 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
162 }
163
164 this->initFluidSystem_();
165
166 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
167 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
168 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
169 }
170
171 this->readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](const unsigned idx) {
172 std::array<int, dim> coords;
173 simulator.vanguard().cartesianCoordinate(idx, coords);
174 for (auto& c : coords) {
175 ++c;
176 }
177 return coords;
178 });
179 FlowProblemType::readMaterialParameters_();
180 FlowProblemType::readThermalParameters_();
181
182 const auto& initconfig = eclState.getInitConfig();
183 if (initconfig.restartRequested())
184 readEclRestartSolution_();
185 else
186 this->readInitialCondition_();
187
188 FlowProblemType::updatePffDofData_();
189
191 const auto& vanguard = this->simulator().vanguard();
192 const auto& gridView = vanguard.gridView();
193 int numElements = gridView.size(/*codim=*/0);
194 this->polymer_.maxAdsorption.resize(numElements, 0.0);
195 }
196
197 /* readBoundaryConditions_();
198
199 // compute and set eq weights based on initial b values
200 computeAndSetEqWeights_();
201
202 if (enableDriftCompensation_) {
203 drift_.resize(this->model().numGridDof());
204 drift_ = 0.0;
205 } */
206
207 // TODO: check wether the following can work with compostional
208 if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
209 simulator.setTimeStepSize(0.0);
211 }
212
213 // after finishing the initialization and writing the initial solution, we move
214 // to the first "real" episode/report step
215 // for restart the episode index and start is already set
216 if (!initconfig.restartRequested()) {
217 simulator.startNextEpisode(schedule.seconds(1));
218 simulator.setEpisodeIndex(0);
219 simulator.setTimeStepIndex(0);
220 }
221 }
222
228 template <class Context>
229 void boundary(BoundaryRateVector& values,
230 const Context& context,
231 unsigned spaceIdx,
232 unsigned /* timeIdx */) const
233 {
235 if (!context.intersection(spaceIdx).boundary())
236 return;
237
238 values.setNoFlow();
239
240 if (this->nonTrivialBoundaryConditions()) {
241 throw std::logic_error("boundary condition is not supported by compostional modeling yet");
242 }
243 }
244
251 template <class Context>
252 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
253 {
254 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
255 const auto& initial_fs = initialFluidStates_[globalDofIdx];
257 using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
258 for (unsigned p = 0; p < numPhases; ++p) { // TODO: assuming the phaseidx continuous
259 ComponentVector evals;
260 auto& last_eval = evals[numComponents - 1];
261 last_eval = 1.;
262 for (unsigned c = 0; c < numComponents - 1; ++c) {
263 const auto val = initial_fs.moleFraction(p, c);
264 const Evaluation eval = Evaluation::createVariable(val, c + 1);
265 evals[c] = eval;
266 last_eval -= eval;
267 }
268 for (unsigned c = 0; c < numComponents; ++c) {
269 fs.setMoleFraction(p, c, evals[c]);
270 }
271
272 // pressure
273 const auto p_val = initial_fs.pressure(p);
274 fs.setPressure(p, Evaluation::createVariable(p_val, 0));
275
276 const auto sat_val = initial_fs.saturation(p);
277 fs.setSaturation(p, sat_val);
278
279 const auto temp_val = initial_fs.temperature(p);
280 fs.setTemperature(temp_val);
281 }
282
283 {
284 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
285 paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
286 paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
287 fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
288 fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
289 fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
290 fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
291 }
292
293 // Set initial K and L
294 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
295 const Evaluation Ktmp = fs.wilsonK_(compIdx);
296 fs.setKvalue(compIdx, Ktmp);
297 }
298
299 const Evaluation& Ltmp = -1.0;
300 fs.setLvalue(Ltmp);
301
302 values.assignNaive(fs);
303 }
304
305 void addToSourceDense(RateVector&, unsigned, unsigned) const override
306 {
307 // we do nothing for now
308 }
309
310 const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
311 { return initialFluidStates_[globalDofIdx]; }
312
313 std::vector<InitialFluidState>& initialFluidStates()
314 { return initialFluidStates_; }
315
316 const std::vector<InitialFluidState>& initialFluidStates() const
317 { return initialFluidStates_; }
318
319 // TODO: do we need this one?
320 template<class Serializer>
321 void serializeOp(Serializer& serializer)
322 {
323 serializer(static_cast<FlowProblemType&>(*this));
324 }
325protected:
326
327 void updateExplicitQuantities_(int /* episodeIdx*/, int /* timeStepSize */, bool /* first_step_after_restart */) override
328 {
329 // we do nothing here for now
330 }
331
332 void readEquilInitialCondition_() override
333 {
334 throw std::logic_error("Equilibration is not supported by compositional modeling yet");
335 }
336
337 void readEclRestartSolution_()
338 {
339 throw std::logic_error("Restarting is not supported by compositional modeling yet");
340 }
341
342 void readExplicitInitialCondition_() override
343 {
344 readExplicitInitialConditionCompositional_();
345 }
346
347 void readExplicitInitialConditionCompositional_()
348 {
349 const auto& simulator = this->simulator();
350 const auto& vanguard = simulator.vanguard();
351 const auto& eclState = vanguard.eclState();
352 const auto& fp = eclState.fieldProps();
353 const bool has_pressure = fp.has_double("PRESSURE");
354 if (!has_pressure)
355 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
356 "keyword if the model is initialized explicitly");
357
358 const bool has_xmf = fp.has_double("XMF");
359 const bool has_ymf = fp.has_double("YMF");
360 const bool has_zmf = fp.has_double("ZMF");
361 if ( !has_zmf && !(has_xmf && has_ymf) ) {
362 throw std::runtime_error("The ECL input file requires the presence of ZMF or XMF and YMF "
363 "keyword if the model is initialized explicitly");
364 }
365
366 if (has_zmf && (has_xmf || has_ymf)) {
367 throw std::runtime_error("The ECL input file can not handle explicit initialization "
368 "with both ZMF and XMF or YMF");
369 }
370
371 if (has_xmf != has_ymf) {
372 throw std::runtime_error("The ECL input file needs XMF and YMF combined to do the explicit "
373 "initializtion when using XMF or YMF");
374 }
375
376 const bool has_temp = fp.has_double("TEMPI");
377
378 // const bool has_gas = fp.has_double("SGAS");
379 assert(fp.has_double("SGAS"));
380
381 std::size_t numDof = this->model().numGridDof();
382
383 initialFluidStates_.resize(numDof);
384
385 std::vector<double> waterSaturationData;
386 std::vector<double> gasSaturationData;
387 std::vector<double> soilData;
388 std::vector<double> pressureData;
389 std::vector<double> tempiData;
390
391 const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
392 const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
393 const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
394
395 if (water_active && Indices::numPhases > 1)
396 waterSaturationData = fp.get_double("SWAT");
397 else
398 waterSaturationData.resize(numDof);
399
400 pressureData = fp.get_double("PRESSURE");
401
402 if (has_temp) {
403 tempiData = fp.get_double("TEMPI");
404 } else {
405 ; // TODO: throw?
406 }
407
408 if (gas_active) // && FluidSystem::phaseIsActive(oilPhaseIdx))
409 gasSaturationData = fp.get_double("SGAS");
410 else
411 gasSaturationData.resize(numDof);
412
413 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
414 auto& dofFluidState = initialFluidStates_[dofIdx];
415 // dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
416
418 assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
419 dofFluidState.setTemperature(temperatureLoc);
420
421 if (gas_active) {
422 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
424 }
425 if (oil_active) {
426 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
427 1.0
430 }
431
433 // set phase pressures
435 const Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
436
437 // TODO: zero capillary pressure for now
438 const std::array<Scalar, numPhases> pc = {0};
439 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
440 if (!FluidSystem::phaseIsActive(phaseIdx))
441 continue;
442
443 if (Indices::oilEnabled)
444 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
445 else if (Indices::gasEnabled)
446 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
447 else if (Indices::waterEnabled)
448 // single (water) phase
449 dofFluidState.setPressure(phaseIdx, pressure);
450 }
451
452 if (has_xmf && has_ymf) {
453 const auto& xmfData = fp.get_double("XMF");
454 const auto& ymfData = fp.get_double("YMF");
455 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
456 const std::size_t data_idx = compIdx * numDof + dofIdx;
457 const Scalar xmf = xmfData[data_idx];
458 const Scalar ymf = ymfData[data_idx];
459
460 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
461 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
462 }
463 }
464
465 if (has_zmf) {
466 const auto& zmfData = fp.get_double("ZMF");
467 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
468 const std::size_t data_idx = compIdx * numDof + dofIdx;
469 const Scalar zmf = zmfData[data_idx];
470 if (gas_active) {
471 const auto ymf = (dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ? zmf : Scalar{0};
472 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
473 }
474 if (oil_active) {
475 const auto xmf = (dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ? zmf : Scalar{0};
476 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
477 }
478 }
479 }
480 }
481 }
482
483private:
484
485 void handleSolventBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
486 {
487 throw std::logic_error("solvent is disabled for compositional modeling and you're trying to add solvent to BC");
488 }
489
490 void handlePolymerBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
491 {
492 throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC");
493 }
494
495 std::vector<InitialFluidState> initialFluidStates_;
496
497 bool enableVtkOutput_;
498};
499
500} // namespace Opm
501
502#endif // OPM_FLOW_PROBLEM_COMP_HPP
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblemComp.hpp:56
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition FlowProblemComp.hpp:113
FlowProblemComp(Simulator &simulator)
Definition FlowProblemComp.hpp:105
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition FlowProblemComp.hpp:229
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition FlowProblemComp.hpp:252
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblemComp.hpp:93
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:92
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:815
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:670
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:179
void writeOutput(bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:491
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235