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
58 using FlowProblemType = FlowProblem<TypeTag>;
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;
79 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
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
131 finishTransmissibilities();
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 {
234 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
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];
256 Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
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
417 Scalar temperatureLoc = tempiData[dofIdx];
418 assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
419 dofFluidState.setTemperature(temperatureLoc);
420
421 if (gas_active) {
422 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
423 gasSaturationData[dofIdx]);
424 }
425 if (oil_active) {
426 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
427 1.0
428 - waterSaturationData[dofIdx]
429 - gasSaturationData[dofIdx]);
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...
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
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:209
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