OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4182 from totto82/bc_dirichlet

Browse Source 
implement dirichlet boundary conditions
This commit is contained in:
Markus Blatt 2022-11-11 22:23:14 +01:00 committed by GitHub
commit 81b296a682
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 99 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
ebos/eclproblem.hh
View File

@ -1672,8 +1672,11 @@ public:
unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(indexInInside);
const auto& dirichlet = dirichlet_(dir)[globalDofIdx];
if (freebc_(dir)[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
else if (std::get<0>(dirichlet) != BCComponent::NONE)
values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
else
values.setMassRate(massratebc_(dir)[globalDofIdx], pvtRegionIdx);
}
@ -1935,6 +1938,89 @@ public:
bool nonTrivialBoundaryConditions() const
{ return nonTrivialBoundaryConditions_; }
const InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
{
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
const auto& dirichlet = dirichlet_(dir)[globalDofIdx];
if(std::get<0>(dirichlet) == BCComponent::NONE)
return initialFluidStates_[globalDofIdx];
InitialFluidState fluidState;
const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
fluidState.setPvtRegionIndex(pvtRegionIdx);
double pressure = initialFluidStates_[globalDofIdx].pressure(oilPhaseIdx);
const auto pressure_input = std::get<1>(dirichlet);
if(pressure_input)
pressure = *pressure_input;
std::array<Scalar, numPhases> pc = {0};
const auto& matParams = materialLawParams(globalDofIdx);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(pc);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
if (Indices::oilEnabled)
fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
else if (Indices::gasEnabled)
fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
else if (Indices::waterEnabled)
//single (water) phase
fluidState.setPressure(phaseIdx, pressure);
}
switch (std::get<0>(dirichlet)) {
case BCComponent::OIL:
if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
throw std::logic_error("oil is not active and you're trying to add oil BC");
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
break;
case BCComponent::GAS:
if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
throw std::logic_error("gas is not active and you're trying to add gas BC");
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
break;
case BCComponent::WATER:
if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
throw std::logic_error("water is not active and you're trying to add water BC");
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
break;
case BCComponent::SOLVENT:
case BCComponent::POLYMER:
case BCComponent::NONE:
throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
break;
}
double temperature = initialFluidStates_[globalDofIdx].temperature(oilPhaseIdx);
const auto temperature_input = std::get<2>(dirichlet);
if(temperature_input)
temperature = *temperature_input;
fluidState.setTemperature(temperature);
fluidState.setRs(0.0);
fluidState.setRv(0.0);
if (FluidSystem::enableVaporizedWater())
fluidState.setRvw(0.0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
fluidState.setInvB(phaseIdx, b);
const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
fluidState.setDensity(phaseIdx, rho);
}
return fluidState;
}
/*!
* \brief Propose the size of the next time step to the simulator.
*
@ -2045,7 +2131,9 @@ public:
return { false, RateVector(0.0) };
}
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
return { freebc_(dir)[globalSpaceIdx], massratebc_(dir)[globalSpaceIdx] };
const auto& dirichlet = dirichlet_(dir)[globalSpaceIdx];
bool free = freebc_(dir)[globalSpaceIdx] || std::get<0>(dirichlet) != BCComponent::NONE;
return { free, massratebc_(dir)[globalSpaceIdx] };
}
private:
@ -2742,6 +2830,7 @@ private:
massratebc_.resize(numElems, 0.0);
freebc_.resize(numElems, false);
dirichlet_.resize(numElems, {BCComponent::NONE, 0.0,0.0});
auto loopAndApply = [&cartesianToCompressedElemIdx,
&vanguard](const auto& bcface,
@ -2808,6 +2897,13 @@ private:
if (initconfig.restartRequested()) {
throw std::logic_error("restart is not compatible with using free boundary conditions");
}
} else if (type == BCType::DIRICHLET) {
const auto component = bcface.component;
const auto pressure = bcface.pressure;
const auto temperature = bcface.temperature;
std::vector<std::tuple<BCComponent, std::optional<double>, std::optional<double>>>& data = dirichlet_(bcface.dir);
loopAndApply(bcface,
[&data,component,pressure,temperature](int elemIdx) { data[elemIdx] = {component, pressure, temperature}; });
} else {
throw std::logic_error("invalid type for BC. Use FREE or RATE");
}
@ -2946,6 +3042,7 @@ private:
BCData<bool> freebc_;
BCData<RateVector> massratebc_;
BCData<std::tuple<BCComponent, std::optional<double>, std::optional<double>>> dirichlet_;
bool nonTrivialBoundaryConditions_ = false;
};