Merge pull request #4171 from akva2/bcdata_members

Make related BC data a single class member
This commit is contained in:
Bård Skaflestad
2022-10-18 13:56:45 +02:00
committed by GitHub

View File

@@ -1649,7 +1649,7 @@ public:
unsigned spaceIdx,
unsigned timeIdx) const
{
if(!context.intersection(spaceIdx).boundary())
if (!context.intersection(spaceIdx).boundary())
return;
if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
@@ -1670,47 +1670,11 @@ public:
unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
switch (indexInInside) {
case 0:
if (freebcXMinus_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcXMinus_[globalDofIdx], pvtRegionIdx);
break;
case 1:
if (freebcX_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcX_[globalDofIdx], pvtRegionIdx);
break;
case 2:
if (freebcYMinus_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcYMinus_[globalDofIdx], pvtRegionIdx);
break;
case 3:
if (freebcY_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcY_[globalDofIdx], pvtRegionIdx);
break;
case 4:
if (freebcZMinus_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcZMinus_[globalDofIdx], pvtRegionIdx);
break;
case 5:
if (freebcZ_[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebcZ_[globalDofIdx], pvtRegionIdx);
break;
default:
throw std::logic_error("invalid face index for boundary condition");
}
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(indexInInside);
if (freebc_(dir)[globalDofIdx])
values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
else
values.setMassRate(massratebc_(dir)[globalDofIdx], pvtRegionIdx);
}
}
@@ -2079,22 +2043,8 @@ public:
if (!nonTrivialBoundaryConditions_) {
return { false, RateVector(0.0) };
}
switch (directionId) {
case 0:
return { freebcXMinus_[globalSpaceIdx], massratebcXMinus_[globalSpaceIdx] };
case 1:
return { freebcX_[globalSpaceIdx], massratebcX_[globalSpaceIdx] };
case 2:
return { freebcYMinus_[globalSpaceIdx], massratebcYMinus_[globalSpaceIdx] };
case 3:
return { freebcY_[globalSpaceIdx], massratebcY_[globalSpaceIdx] };
case 4:
return { freebcZMinus_[globalSpaceIdx], massratebcZMinus_[globalSpaceIdx] };
case 5:
return { freebcZ_[globalSpaceIdx], massratebcZ_[globalSpaceIdx] };
default:
return { false, RateVector(0.0) };
}
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
return { freebc_(dir)[globalSpaceIdx], massratebc_(dir)[globalSpaceIdx] };
}
private:
@@ -2776,7 +2726,6 @@ private:
void readBoundaryConditions_()
{
nonTrivialBoundaryConditions_ = false;
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
@@ -2790,18 +2739,24 @@ private:
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
massratebcXMinus_.resize(numElems, 0.0);
massratebcX_.resize(numElems, 0.0);
massratebcYMinus_.resize(numElems, 0.0);
massratebcY_.resize(numElems, 0.0);
massratebcZMinus_.resize(numElems, 0.0);
massratebcZ_.resize(numElems, 0.0);
freebcX_.resize(numElems, false);
freebcXMinus_.resize(numElems, false);
freebcY_.resize(numElems, false);
freebcYMinus_.resize(numElems, false);
freebcZ_.resize(numElems, false);
freebcZMinus_.resize(numElems, false);
massratebc_.resize(numElems, 0.0);
freebc_.resize(numElems, false);
auto loopAndApply = [&cartesianToCompressedElemIdx,
&vanguard](const auto& bcface,
auto apply)
{
for (int i = bcface.i1; i <= bcface.i2; ++i) {
for (int j = bcface.j1; j <= bcface.j2; ++j) {
for (int k = bcface.k1; k <= bcface.k2; ++k) {
std::array<int, 3> tmp = {i,j,k};
auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
if (elemIdx >= 0)
apply(elemIdx);
}
}
}
};
for (const auto& bcface : bcconfig) {
const auto& type = bcface.bctype;
@@ -2835,76 +2790,16 @@ private:
break;
}
std::vector<RateVector>* data = nullptr;
switch (bcface.dir) {
case FaceDir::XMinus:
data = &massratebcXMinus_;
break;
case FaceDir::XPlus:
data = &massratebcX_;
break;
case FaceDir::YMinus:
data = &massratebcYMinus_;
break;
case FaceDir::YPlus:
data = &massratebcY_;
break;
case FaceDir::ZMinus:
data = &massratebcZMinus_;
break;
case FaceDir::ZPlus:
data = &massratebcZ_;
break;
case FaceDir::Unknown:
throw std::runtime_error("Unexpected unknown face direction");
}
std::vector<RateVector>& data = massratebc_(bcface.dir);
const Evaluation rate = bcface.rate;
for (int i = bcface.i1; i <= bcface.i2; ++i) {
for (int j = bcface.j1; j <= bcface.j2; ++j) {
for (int k = bcface.k1; k <= bcface.k2; ++k) {
std::array<int, 3> tmp = {i,j,k};
auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
if (elemIdx >= 0)
(*data)[elemIdx][compIdx] = rate;
}
}
}
loopAndApply(bcface,
[&data,compIdx,rate](int elemIdx)
{ data[elemIdx][compIdx] = rate; });
} else if (type == BCType::FREE) {
std::vector<bool>* data = nullptr;
switch (bcface.dir) {
case FaceDir::XMinus:
data = &freebcXMinus_;
break;
case FaceDir::XPlus:
data = &freebcX_;
break;
case FaceDir::YMinus:
data = &freebcYMinus_;
break;
case FaceDir::YPlus:
data = &freebcY_;
break;
case FaceDir::ZMinus:
data = &freebcZMinus_;
break;
case FaceDir::ZPlus:
data = &freebcZ_;
break;
case FaceDir::Unknown:
throw std::runtime_error("Unexpected unknown face direction");
}
for (int i = bcface.i1; i <= bcface.i2; ++i) {
for (int j = bcface.j1; j <= bcface.j2; ++j) {
for (int k = bcface.k1; k <= bcface.k2; ++k) {
std::array<int, 3> tmp = {i,j,k};
auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
if (elemIdx >= 0)
(*data)[elemIdx] = true;
}
}
}
std::vector<bool>& data = freebc_(bcface.dir);
loopAndApply(bcface,
[&data](int elemIdx) { data[elemIdx] = true; });
// TODO: either the real initial solution needs to be computed or read from the restart file
const auto& eclState = simulator.vanguard().eclState();
@@ -2960,7 +2855,8 @@ private:
return dtNext;
}
void computeAndSetEqWeights_() {
void computeAndSetEqWeights_()
{
std::vector<Scalar> sumInvB(numPhases, 0.0);
const auto& gridView = this->gridView();
ElementContext elemCtx(this->simulator());
@@ -3018,20 +2914,38 @@ private:
EclActionHandler actionHandler_;
std::vector<bool> freebcX_;
std::vector<bool> freebcXMinus_;
std::vector<bool> freebcY_;
std::vector<bool> freebcYMinus_;
std::vector<bool> freebcZ_;
std::vector<bool> freebcZMinus_;
template<class T>
struct BCData
{
std::array<std::vector<T>,6> data;
bool nonTrivialBoundaryConditions_;
std::vector<RateVector> massratebcX_;
std::vector<RateVector> massratebcXMinus_;
std::vector<RateVector> massratebcY_;
std::vector<RateVector> massratebcYMinus_;
std::vector<RateVector> massratebcZ_;
std::vector<RateVector> massratebcZMinus_;
void resize(size_t size, T defVal)
{
for (auto& d : data)
d.resize(size, defVal);
}
const std::vector<T>& operator()(FaceDir::DirEnum dir) const
{
if (dir == FaceDir::DirEnum::Unknown)
throw std::runtime_error("Tried to acess BC data for the 'Unknown' direction");
int idx = 0;
int div = static_cast<int>(dir);
while ((div /= 2) >= 1)
++idx;
assert(idx >= 0 && idx <= 5);
return data[idx];
}
std::vector<T>& operator()(FaceDir::DirEnum dir)
{
return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
}
};
BCData<bool> freebc_;
BCData<RateVector> massratebc_;
bool nonTrivialBoundaryConditions_ = false;
};
} // namespace Opm