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;
80 using typename FlowProblemType::Evaluation;
81 using typename FlowProblemType::MaterialLaw;
82 using typename FlowProblemType::RateVector;
98 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1
e-7);
118 FlowProblemType::finishInit();
120 auto& simulator = this->simulator();
126 this->transmissibilities_.finishInit(
127 [&
vg = this->simulator().vanguard()](
const unsigned int it) {
return vg.gridIdxToEquilGridIdx(it); });
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();
184 readEclRestartSolution_();
186 this->readInitialCondition_();
188 FlowProblemType::updatePffDofData_();
191 const auto& vanguard = this->simulator().vanguard();
192 const auto& gridView = vanguard.gridView();
194 this->polymer_.maxAdsorption.resize(
numElements, 0.0);
208 if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
209 simulator.setTimeStepSize(0.0);
217 simulator.startNextEpisode(schedule.seconds(1));
218 simulator.setEpisodeIndex(0);
219 simulator.setTimeStepIndex(0);
228 template <
class Context>
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>
257 using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
258 for (
unsigned p = 0;
p < numPhases; ++
p) {
259 ComponentVector
evals;
262 for (
unsigned c = 0;
c < numComponents - 1; ++
c) {
264 const Evaluation eval = Evaluation::createVariable(
val,
c + 1);
268 for (
unsigned c = 0;
c < numComponents; ++
c) {
274 fs.setPressure(
p, Evaluation::createVariable(
p_val, 0));
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));
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
313 std::vector<InitialFluidState>& initialFluidStates()
314 {
return initialFluidStates_; }
316 const std::vector<InitialFluidState>& initialFluidStates()
const
317 {
return initialFluidStates_; }
320 template<
class 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();
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");
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");
367 throw std::runtime_error(
"The ECL input file can not handle explicit initialization "
368 "with both ZMF and XMF or 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");
381 std::size_t numDof = this->model().numGridDof();
383 initialFluidStates_.resize(numDof);
391 const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
392 const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
393 const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
438 const std::array<Scalar, numPhases>
pc = {0};
440 if (!FluidSystem::phaseIsActive(
phaseIdx))
443 if (Indices::oilEnabled)
445 else if (Indices::gasEnabled)
447 else if (Indices::waterEnabled)
453 const auto&
xmfData =
fp.get_double(
"XMF");
454 const auto&
ymfData =
fp.get_double(
"YMF");
466 const auto&
zmfData =
fp.get_double(
"ZMF");
471 const auto ymf = (
dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ?
zmf : Scalar{0};
475 const auto xmf = (
dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ?
zmf : Scalar{0};
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_;
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
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:815
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