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

add an abstraction layer for the matrix between the linearizer and linear solver

Browse Source 
this allows to assemble the Jacobian matrices directly into the native
format expected by linear solver. So far, only backends using
Dune::BCRSMatrix are provided, but there are work-in-progress patches
for dune-fem, vienna-CL and PETSc backends.
This commit is contained in:
Robert Kloefkorn 2018-11-12 12:40:04 +01:00 committed by Andreas Lauser
parent a5493919da
commit cfcf04ed00
1 changed files with 21 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
ebos/eclpeacemanwell.hh
View File

@ -343,14 +343,17 @@ public:
{ {
const SolutionVector& curSol = simulator_.model().solution(/*timeIdx=*/0); const SolutionVector& curSol = simulator_.model().solution(/*timeIdx=*/0);
typedef typename JacobianMatrix::block_type MatrixBlock;
unsigned wellGlobalDofIdx = AuxModule::localToGlobalDof(/*localDofIdx=*/0); unsigned wellGlobalDofIdx = AuxModule::localToGlobalDof(/*localDofIdx=*/0);
residual[wellGlobalDofIdx] = 0.0; residual[wellGlobalDofIdx] = 0.0;
auto& diagBlock = matrix[wellGlobalDofIdx][wellGlobalDofIdx]; MatrixBlock diagBlock(0.0);
diagBlock = 0.0;
for (unsigned i = 0; i < numModelEq; ++ i) for (unsigned i = 0; i < numModelEq; ++ i)
diagBlock[i][i] = 1.0; diagBlock[i][i] = 1.0;
MatrixBlock block(0.0);
if (wellStatus() == Shut) { if (wellStatus() == Shut) {
// if the well is shut, make the auxiliary DOFs a trivial equation in the // if the well is shut, make the auxiliary DOFs a trivial equation in the
// matrix: the main diagonal is already set to the identity matrix, the // matrix: the main diagonal is already set to the identity matrix, the
@ -358,16 +361,16 @@ public:
auto wellDofIt = dofVariables_.begin(); auto wellDofIt = dofVariables_.begin();
const auto& wellDofEndIt = dofVariables_.end(); const auto& wellDofEndIt = dofVariables_.end();
for (; wellDofIt != wellDofEndIt; ++ wellDofIt) { for (; wellDofIt != wellDofEndIt; ++ wellDofIt) {
matrix[wellGlobalDofIdx][wellDofIt->first] = 0.0; matrix.setBlock(wellGlobalDofIdx, wellDofIt->first, block);
matrix[wellDofIt->first][wellGlobalDofIdx] = 0.0; matrix.setBlock(wellDofIt->first, wellGlobalDofIdx, block);
} }
matrix[wellGlobalDofIdx][wellGlobalDofIdx] = diagBlock; matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock);
residual[wellGlobalDofIdx] = 0.0; residual[wellGlobalDofIdx] = 0.0;
return; return;
} }
else if (dofVariables_.empty()) { else if (dofVariables_.empty()) {
// the well does not feature any perforations on the local process // the well does not feature any perforations on the local process
matrix[wellGlobalDofIdx][wellGlobalDofIdx] = diagBlock; matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock);
residual[wellGlobalDofIdx] = 0.0; residual[wellGlobalDofIdx] = 0.0;
return; return;
} }
@ -389,10 +392,9 @@ public:
///////////// /////////////
// influence of grid on well // influence of grid on well
auto& curBlock = matrix[wellGlobalDofIdx][gridDofIdx];
elemCtx.updateStencil( dofVars.element ); elemCtx.updateStencil( dofVars.element );
curBlock = 0.0; // reset block from previous values
block = 0.0;
for (unsigned priVarIdx = 0; priVarIdx < numModelEq; ++priVarIdx) { for (unsigned priVarIdx = 0; priVarIdx < numModelEq; ++priVarIdx) {
// calculate the derivative of the well equation w.r.t. the current // calculate the derivative of the well equation w.r.t. the current
// primary variable using forward differences // primary variable using forward differences
@ -408,11 +410,13 @@ public:
Scalar dWellEq_dPV = Scalar dWellEq_dPV =
(wellResidual_(actualBottomHolePressure_, &tmpDofVars, gridDofIdx) - wellResid) (wellResidual_(actualBottomHolePressure_, &tmpDofVars, gridDofIdx) - wellResid)
/ eps; / eps;
curBlock[0][priVarIdx] = dWellEq_dPV; block[0][priVarIdx] = dWellEq_dPV;
// go back to the original primary variables // go back to the original primary variables
priVars[priVarIdx] -= eps; priVars[priVarIdx] -= eps;
} }
matrix.setBlock(wellGlobalDofIdx, gridDofIdx, block);
// //
///////////// /////////////
@ -462,10 +466,12 @@ public:
// the black-oil model for now anyway, so this should not be too much of a // the black-oil model for now anyway, so this should not be too much of a
// problem... // problem...
Opm::Valgrind::CheckDefined(q); Opm::Valgrind::CheckDefined(q);
auto& matrixEntry = matrix[gridDofIdx][wellGlobalDofIdx]; block = 0.0;
matrixEntry = 0.0;
for (unsigned eqIdx = 0; eqIdx < numModelEq; ++ eqIdx) for (unsigned eqIdx = 0; eqIdx < numModelEq; ++ eqIdx)
matrixEntry[eqIdx][0] = - Toolbox::value(q[eqIdx])/dofVars.totalVolume; block[eqIdx][0] = - Opm::getValue(q[eqIdx])/dofVars.totalVolume;
matrix.setBlock(gridDofIdx, wellGlobalDofIdx, block);
// //
///////////// /////////////
} }
@ -477,6 +483,8 @@ public:
*std::max<Scalar>(1e7, targetBottomHolePressure_); *std::max<Scalar>(1e7, targetBottomHolePressure_);
Scalar wellResidStar = wellResidual_(actualBottomHolePressure_ + eps); Scalar wellResidStar = wellResidual_(actualBottomHolePressure_ + eps);
diagBlock[0][0] = (wellResidStar - wellResid)/eps; diagBlock[0][0] = (wellResidStar - wellResid)/eps;
matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock);
} }