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>
84 using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
115 FlowProblemType::finishInit();
117 auto& simulator = this->simulator();
119 auto finishTransmissibilities = [updated =
false,
this]()
mutable {
124 [&vg = this->simulator().vanguard()](
const unsigned int it) {
return vg.gridIdxToEquilGridIdx(it); });
128 finishTransmissibilities();
130 const auto& eclState = simulator.vanguard().eclState();
131 const auto& schedule = simulator.vanguard().schedule();
134 simulator.setStartTime(schedule.getStartTime());
135 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
141 simulator.setEpisodeIndex(-1);
142 simulator.setEpisodeLength(0.0);
147 this->gravity_ = 0.0;
148 if (Parameters::Get<Parameters::EnableGravity>())
149 this->gravity_[dim - 1] = 9.80665;
150 if (!eclState.getInitConfig().hasGravity())
151 this->gravity_[dim - 1] = 0.0;
156 const auto& tuning = schedule[0].tuning();
163 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
164 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
168 this->
readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](
const unsigned idx) {
169 std::array<int, dim> coords;
170 simulator.vanguard().cartesianCoordinate(idx, coords);
171 for (auto& c : coords) {
179 const auto& initconfig = eclState.getInitConfig();
180 if (initconfig.restartRequested())
187 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
188 const auto& vanguard = this->simulator().vanguard();
189 const auto& gridView = vanguard.gridView();
190 int numElements = gridView.size(0);
191 this->
polymer_.maxAdsorption.resize(numElements, 0.0);
205 if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
206 simulator.setTimeStepSize(0.0);
213 if (!initconfig.restartRequested()) {
214 simulator.startNextEpisode(schedule.seconds(1));
215 simulator.setEpisodeIndex(0);
216 simulator.setTimeStepIndex(0);
225 template <
class Context>
227 const Context& context,
231 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
232 if (!context.intersection(spaceIdx).boundary())
237 if (this->nonTrivialBoundaryConditions()) {
238 throw std::logic_error(
"boundary condition is not supported by compostional modeling yet");
248 template <
class Context>
251 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
252 const auto& initial_fs = initialFluidStates_[globalDofIdx];
253 Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
254 using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
255 for (
unsigned p = 0; p < numPhases; ++p) {
256 ComponentVector evals;
257 auto& last_eval = evals[numComponents - 1];
259 for (
unsigned c = 0; c < numComponents - 1; ++c) {
260 const auto val = initial_fs.moleFraction(p, c);
261 const Evaluation eval = Evaluation::createVariable(val, c + 1);
265 for (
unsigned c = 0; c < numComponents; ++c) {
266 fs.setMoleFraction(p, c, evals[c]);
270 const auto p_val = initial_fs.pressure(p);
271 fs.setPressure(p, Evaluation::createVariable(p_val, 0));
273 const auto sat_val = initial_fs.saturation(p);
274 fs.setSaturation(p, sat_val);
276 const auto temp_val = initial_fs.temperature(p);
277 fs.setTemperature(temp_val);
281 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
282 paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
283 paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
284 fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
285 fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
286 fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
287 fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
291 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
293 fs.setKvalue(compIdx, Ktmp);
299 values.assignNaive(fs);
308 {
return initialFluidStates_[globalDofIdx]; }
311 {
return initialFluidStates_; }
314 {
return initialFluidStates_; }
317 template<
class Serializer>
331 throw std::logic_error(
"Equilibration is not supported by compositional modeling yet");
336 throw std::logic_error(
"Restarting is not supported by compositional modeling yet");
341 readExplicitInitialConditionCompositional_();
346 const auto& simulator = this->simulator();
347 const auto& vanguard = simulator.vanguard();
348 const auto& eclState = vanguard.eclState();
349 const auto& fp = eclState.fieldProps();
351 assert(fp.has_double(
"XMF"));
353 assert(fp.has_double(
"YMF"));
354 const bool has_temp = fp.has_double(
"TEMPI");
355 const bool has_pressure = fp.has_double(
"PRESSURE");
358 assert(fp.has_double(
"SGAS"));
360 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
361 "keyword if the model is initialized explicitly");
363 std::size_t numDof = this->model().numGridDof();
365 initialFluidStates_.resize(numDof);
367 std::vector<double> xmfData;
368 std::vector<double> ymfData;
369 std::vector<double> waterSaturationData;
370 std::vector<double> gasSaturationData;
371 std::vector<double> soilData;
372 std::vector<double> pressureData;
373 std::vector<double> tempiData;
375 if (waterPhaseIdx > 0 && Indices::numPhases > 1)
376 waterSaturationData = fp.get_double(
"SWAT");
378 waterSaturationData.resize(numDof);
380 pressureData = fp.get_double(
"PRESSURE");
382 assert(waterPhaseIdx < 0);
384 xmfData = fp.get_double(
"XMF");
385 ymfData = fp.get_double(
"YMF");
387 tempiData = fp.get_double(
"TEMPI");
393 gasSaturationData = fp.get_double(
"SGAS");
395 gasSaturationData.resize(numDof);
399 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
400 auto& dofFluidState = initialFluidStates_[dofIdx];
403 Scalar temperatureLoc = tempiData[dofIdx];
404 assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
405 dofFluidState.setTemperature(temperatureLoc);
407 if (gasPhaseIdx > 0) {
408 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
409 gasSaturationData[dofIdx]);
411 if (oilPhaseIdx > 0) {
412 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
414 - waterSaturationData[dofIdx]
415 - gasSaturationData[dofIdx]);
421 const Scalar pressure = pressureData[dofIdx];
424 const std::array<Scalar, numPhases> pc = {0};
425 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
426 if (!FluidSystem::phaseIsActive(phaseIdx))
429 if (Indices::oilEnabled)
430 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
431 else if (Indices::gasEnabled)
432 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
433 else if (Indices::waterEnabled)
435 dofFluidState.setPressure(phaseIdx, pressure);
438 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
439 const std::size_t data_idx = compIdx * numDof + dofIdx;
440 const Scalar xmf = xmfData[data_idx];
441 const Scalar ymf = ymfData[data_idx];
443 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
444 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
451 void handleSolventBC(
const BCProp::BCFace& , RateVector& )
const override
453 throw std::logic_error(
"solvent is disabled for compositional modeling and you're trying to add solvent to BC");
456 void handlePolymerBC(
const BCProp::BCFace& , RateVector& )
const override
458 throw std::logic_error(
"polymer is disabled for compositional modeling and you're trying to add polymer to BC");
461 std::vector<InitialFluidState> initialFluidStates_;
463 bool enableVtkOutput_;
PolymerSolutionContainer< Scalar > polymer_
Definition: FlowGenericProblem.hpp:357
Scalar initialTimeStepSize_
Definition: FlowGenericProblem.hpp:368
bool enableTuning_
Definition: FlowGenericProblem.hpp:367
void readRockParameters_(const std::vector< Scalar > &cellCenterDepths, std::function< std::array< int, 3 >(const unsigned)> ijkIndex)
Definition: FlowGenericProblem_impl.hpp:138
std::vector< Scalar > maxOilSaturation_
Definition: FlowGenericProblem.hpp:358
void initFluidSystem_()
Definition: FlowGenericProblem_impl.hpp:482
Scalar maxTimeStepAfterWellEvent_
Definition: FlowGenericProblem.hpp:369
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblemComp.hpp:56
const std::vector< InitialFluidState > & initialFluidStates() const
Definition: FlowProblemComp.hpp:313
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: FlowProblemComp.hpp:110
FlowProblemComp(Simulator &simulator)
Definition: FlowProblemComp.hpp:102
void readExplicitInitialCondition_() override
Definition: FlowProblemComp.hpp:339
void readExplicitInitialConditionCompositional_()
Definition: FlowProblemComp.hpp:344
std::vector< InitialFluidState > & initialFluidStates()
Definition: FlowProblemComp.hpp:310
void readEclRestartSolution_()
Definition: FlowProblemComp.hpp:334
const InitialFluidState & initialFluidState(unsigned globalDofIdx) const
Definition: FlowProblemComp.hpp:307
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition: FlowProblemComp.hpp:226
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: FlowProblemComp.hpp:249
void readEquilInitialCondition_() override
Definition: FlowProblemComp.hpp:329
void serializeOp(Serializer &serializer)
Definition: FlowProblemComp.hpp:318
void updateExplicitQuantities_(int, int, bool) override
Definition: FlowProblemComp.hpp:324
void addToSourceDense(RateVector &, unsigned, unsigned) const override
Definition: FlowProblemComp.hpp:302
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:91
GetPropType< TypeTag, Properties::Evaluation > Evaluation
Definition: FlowProblem.hpp:153
@ numComponents
Definition: FlowProblem.hpp:113
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:813
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:668
GetPropType< TypeTag, Properties::Vanguard > Vanguard
Definition: FlowProblem.hpp:104
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: FlowProblem.hpp:98
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: FlowProblem.hpp:142
@ dim
Definition: FlowProblem.hpp:107
GetPropType< TypeTag, Properties::Indices > Indices
Definition: FlowProblem.hpp:154
@ gasPhaseIdx
Definition: FlowProblem.hpp:131
@ dimWorld
Definition: FlowProblem.hpp:108
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: FlowProblem.hpp:143
@ numPhases
Definition: FlowProblem.hpp:112
void readThermalParameters_()
Definition: FlowProblem.hpp:1349
@ waterPhaseIdx
Definition: FlowProblem.hpp:133
GetPropType< TypeTag, Properties::GridView > GridView
Definition: FlowProblem.hpp:99
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: FlowProblem.hpp:178
void updatePffDofData_()
Definition: FlowProblem.hpp:1461
GetPropType< TypeTag, Properties::PrimaryVariables > PrimaryVariables
Definition: FlowProblem.hpp:141
Vanguard::TransmissibilityType transmissibilities_
Definition: FlowProblem.hpp:1644
virtual void readInitialCondition_()
Definition: FlowProblem.hpp:1389
void writeOutput(bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition: FlowProblem.hpp:489
@ oilPhaseIdx
Definition: FlowProblem.hpp:132
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: FlowProblem.hpp:101
GetPropType< TypeTag, Properties::MaterialLaw > MaterialLaw
Definition: FlowProblem.hpp:151
void readMaterialParameters_()
Definition: FlowProblem.hpp:1315
Definition: blackoilboundaryratevector.hh:37
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