30#ifndef OPM_FLOW_PROBLEM_COMP_HPP
31#define OPM_FLOW_PROBLEM_COMP_HPP
35#include <opm/material/fluidstates/CompositionalFluidState.hpp>
37#include <opm/material/thermal/EclThermalLawManager.hpp>
54template <
class TypeTag>
55class FlowProblemComp :
public FlowProblem<TypeTag>
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;
67 using FlowProblemType::dim;
68 using FlowProblemType::dimWorld;
70 using FlowProblemType::numPhases;
71 using FlowProblemType::numComponents;
73 using FlowProblemType::gasPhaseIdx;
74 using FlowProblemType::oilPhaseIdx;
75 using FlowProblemType::waterPhaseIdx;
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;
84 using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
98 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-7);
118 FlowProblemType::finishInit();
120 auto& simulator = this->simulator();
122 auto finishTransmissibilities = [updated =
false,
this]()
mutable {
126 this->transmissibilities_.finishInit(
127 [&vg = this->simulator().vanguard()](
const unsigned int it) {
return vg.gridIdxToEquilGridIdx(it); });
131 finishTransmissibilities();
133 const auto& eclState = simulator.vanguard().eclState();
134 const auto& schedule = simulator.vanguard().schedule();
137 simulator.setStartTime(schedule.getStartTime());
138 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
144 simulator.setEpisodeIndex(-1);
145 simulator.setEpisodeLength(0.0);
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;
156 if (this->enableTuning_) {
159 const auto& tuning = schedule[0].tuning();
160 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
161 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
164 this->initFluidSystem_();
166 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
167 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
168 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
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) {
179 FlowProblemType::readMaterialParameters_();
180 FlowProblemType::readThermalParameters_();
182 const auto& initconfig = eclState.getInitConfig();
183 if (initconfig.restartRequested())
184 readEclRestartSolution_();
186 this->readInitialCondition_();
188 FlowProblemType::updatePffDofData_();
191 const auto& vanguard = this->simulator().vanguard();
192 const auto& gridView = vanguard.gridView();
193 int numElements = gridView.size(0);
194 this->polymer_.maxAdsorption.resize(numElements, 0.0);
208 if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
209 simulator.setTimeStepSize(0.0);
216 if (!initconfig.restartRequested()) {
217 simulator.startNextEpisode(schedule.seconds(1));
218 simulator.setEpisodeIndex(0);
219 simulator.setTimeStepIndex(0);
228 template <
class Context>
230 const Context& context,
234 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
235 if (!context.intersection(spaceIdx).boundary())
240 if (this->nonTrivialBoundaryConditions()) {
241 throw std::logic_error(
"boundary condition is not supported by compostional modeling yet");
251 template <
class Context>
252 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
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) {
259 ComponentVector evals;
260 auto& last_eval = evals[numComponents - 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);
268 for (
unsigned c = 0; c < numComponents; ++c) {
269 fs.setMoleFraction(p, c, evals[c]);
273 const auto p_val = initial_fs.pressure(p);
274 fs.setPressure(p, Evaluation::createVariable(p_val, 0));
276 const auto sat_val = initial_fs.saturation(p);
277 fs.setSaturation(p, sat_val);
279 const auto temp_val = initial_fs.temperature(p);
280 fs.setTemperature(temp_val);
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));
294 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
295 const Evaluation Ktmp = fs.wilsonK_(compIdx);
296 fs.setKvalue(compIdx, Ktmp);
299 const Evaluation& Ltmp = -1.0;
302 values.assignNaive(fs);
305 void addToSourceDense(RateVector&,
unsigned,
unsigned)
const override
310 const InitialFluidState& initialFluidState(
unsigned globalDofIdx)
const
311 {
return initialFluidStates_[globalDofIdx]; }
313 std::vector<InitialFluidState>& initialFluidStates()
314 {
return initialFluidStates_; }
316 const std::vector<InitialFluidState>& initialFluidStates()
const
317 {
return initialFluidStates_; }
320 template<
class Serializer>
321 void serializeOp(Serializer& serializer)
323 serializer(
static_cast<FlowProblemType&
>(*
this));
327 void updateExplicitQuantities_(
int ,
int ,
bool )
override
332 void readEquilInitialCondition_()
override
334 throw std::logic_error(
"Equilibration is not supported by compositional modeling yet");
337 void readEclRestartSolution_()
339 throw std::logic_error(
"Restarting is not supported by compositional modeling yet");
342 void readExplicitInitialCondition_()
override
344 readExplicitInitialConditionCompositional_();
347 void readExplicitInitialConditionCompositional_()
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");
355 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
356 "keyword if the model is initialized explicitly");
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");
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");
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");
376 const bool has_temp = fp.has_double(
"TEMPI");
379 assert(fp.has_double(
"SGAS"));
381 std::size_t numDof = this->model().numGridDof();
383 initialFluidStates_.resize(numDof);
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;
391 const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
392 const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
393 const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
395 if (water_active && Indices::numPhases > 1)
396 waterSaturationData = fp.get_double(
"SWAT");
398 waterSaturationData.resize(numDof);
400 pressureData = fp.get_double(
"PRESSURE");
403 tempiData = fp.get_double(
"TEMPI");
409 gasSaturationData = fp.get_double(
"SGAS");
411 gasSaturationData.resize(numDof);
413 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
414 auto& dofFluidState = initialFluidStates_[dofIdx];
417 Scalar temperatureLoc = tempiData[dofIdx];
418 assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
419 dofFluidState.setTemperature(temperatureLoc);
422 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
423 gasSaturationData[dofIdx]);
426 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
428 - waterSaturationData[dofIdx]
429 - gasSaturationData[dofIdx]);
435 const Scalar pressure = pressureData[dofIdx];
438 const std::array<Scalar, numPhases> pc = {0};
439 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
440 if (!FluidSystem::phaseIsActive(phaseIdx))
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)
449 dofFluidState.setPressure(phaseIdx, pressure);
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];
460 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
461 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
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];
471 const auto ymf = (dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ? zmf : Scalar{0};
472 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
475 const auto xmf = (dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ? zmf : Scalar{0};
476 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
485 void handleSolventBC(
const BCProp::BCFace& , RateVector& )
const override
487 throw std::logic_error(
"solvent is disabled for compositional modeling and you're trying to add solvent to BC");
490 void handlePolymerBC(
const BCProp::BCFace& , RateVector& )
const override
492 throw std::logic_error(
"polymer is disabled for compositional modeling and you're trying to add polymer to BC");
495 std::vector<InitialFluidState> initialFluidStates_;
497 bool enableVtkOutput_;
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