Fix Flows/Flores for NNCs

This commit is contained in:
Tor Harald Sandve 2024-02-09 09:49:55 +01:00
parent edbe1a113a
commit e6c8eedeb8

View File

@ -496,7 +496,7 @@ private:
bf.area(), bf.area(),
bf.integrationPos()[dimWorld - 1], bf.integrationPos()[dimWorld - 1],
exFluidState}; exFluidState};
boundaryInfo_.push_back({myIdx, dir_id, bcdata}); boundaryInfo_.push_back({myIdx, dir_id, bfIndex, bcdata});
} }
} }
} }
@ -663,11 +663,9 @@ public:
// Flux term. // Flux term.
{ {
OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell); OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell);
short loc = 0;
for (const auto& nbInfo : nbInfos) { for (const auto& nbInfo : nbInfos) {
OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace); OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace);
// not for NNCs
if (nbInfo.res_nbinfo.dirId < 0)
continue;
unsigned globJ = nbInfo.neighbor; unsigned globJ = nbInfo.neighbor;
assert(globJ != globI); assert(globJ != globI);
adres = 0.0; adres = 0.0;
@ -676,30 +674,35 @@ public:
LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo); LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
adres *= nbInfo.res_nbinfo.faceArea; adres *= nbInfo.res_nbinfo.faceArea;
if (enableFlows) { if (enableFlows) {
for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
flowsInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = adres[phaseIdx].value(); flowsInfo_[globI][loc].flow[eqIdx] = adres[eqIdx].value();
} }
} }
if (enableFlores) { if (enableFlores) {
for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
floresInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = darcyFlux[phaseIdx].value(); floresInfo_[globI][loc].flow[eqIdx] = darcyFlux[eqIdx].value();
} }
} }
++loc;
} }
} }
} }
// Boundary terms. Only looping over cells with nontrivial bcs. // Boundary terms. Only looping over cells with nontrivial bcs.
for (const auto& bdyInfo : boundaryInfo_) { for (const auto& bdyInfo : boundaryInfo_) {
if (bdyInfo.bcdata.type == BCType::NONE)
continue;
ADVectorBlock adres(0.0); ADVectorBlock adres(0.0);
const unsigned globI = bdyInfo.cell; const unsigned globI = bdyInfo.cell;
const auto& nbInfos = neighborInfo_[globI];
const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0); const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI); LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
adres *= bdyInfo.bcdata.faceArea; adres *= bdyInfo.bcdata.faceArea;
const unsigned bfIndex = bdyInfo.bfIndex;
if (enableFlows) { if (enableFlows) {
for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
flowsInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value(); flowsInfo_[globI][nbInfos.size() + bfIndex].flow[eqIdx] = adres[eqIdx].value();
} }
} }
// TODO also store Flores? // TODO also store Flores?
@ -848,6 +851,9 @@ private:
// Boundary terms. Only looping over cells with nontrivial bcs. // Boundary terms. Only looping over cells with nontrivial bcs.
for (const auto& bdyInfo : boundaryInfo_) { for (const auto& bdyInfo : boundaryInfo_) {
if (bdyInfo.bcdata.type == BCType::NONE)
continue;
VectorBlock res(0.0); VectorBlock res(0.0);
MatrixBlock bMat(0.0); MatrixBlock bMat(0.0);
ADVectorBlock adres(0.0); ADVectorBlock adres(0.0);
@ -934,6 +940,7 @@ private:
{ {
unsigned int cell; unsigned int cell;
int dir; int dir;
unsigned int bfIndex;
BoundaryConditionData bcdata; BoundaryConditionData bcdata;
}; };
std::vector<BoundaryInfo> boundaryInfo_; std::vector<BoundaryInfo> boundaryInfo_;