Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Avoid assembling ghost row off-diagonals. #623

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Avoid assembling ghost row off-diagonals. #623

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions opm/models/discretization/common/fvbaselinearizer.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -349,10 +349,18 @@ private:

for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);

for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
if (dofIdx > 0 && stencil.element(dofIdx).partitionType() != Dune::PartitionType::InteriorEntity) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will break all discretizations besides cell-centered finite volume in parallel.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feared something like that... So what is the best solution? Make an cell-centered linearizer class that overrides these methods? I first experimented with modifying the stencil (which is specific to the discretization), but that did not work because it is used also to assemble the fluxes (obvious in hindsight...) which should be included.

Copy link
Member

@blattms blattms Aug 14, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And I think it might even not work for cell-centered:
If we are on an interior element, then we want to have the connections to other non-interior dofs in the sparsity pattern. But with you changes we do not get these are we. Hence even for cell-centered finite volumes we need to have an additional check whether the current element is in the interior and only if not execute your new code.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On that, I think you are not correct. The proposed code eliminates the ghost rows (except diagonal) but not ghost columns. The sparsity pattern created is not symmetric any more, note that the myIdx is the column number, and we do not eliminate based on that, but on the neighborIdx (using the corresponding element).

// Skip non-diagonal elements on non-interior rows.
continue;
}
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
sparsityPattern[myIdx].insert(neighborIdx);
// When the stencil is used later (in
// linearizeElement_(), via an element context) it
// assumes that the stencil is column-wise,
// i.e. the first index myIdx refers to the column
// number.
sparsityPattern[neighborIdx].insert(myIdx);
}
}
}
Expand Down Expand Up @@ -526,6 +534,10 @@ private:

// update the global Jacobian matrix
for (unsigned dofIdx = 0; dofIdx < elementCtx->numDof(/*timeIdx=*/0); ++ dofIdx) {
if (dofIdx > 0 && elementCtx->stencil(/*timeIdx=*/0).element(dofIdx).partitionType() != Dune::PartitionType::InteriorEntity) {
// Skip non-diagonal elements on non-interior rows.
continue;
}
unsigned globJ = elementCtx->globalSpaceIndex(/*spaceIdx=*/dofIdx, /*timeIdx=*/0);

jacobian_->addToBlock(globJ, globI, localLinearizer.jacobian(dofIdx, primaryDofIdx));
Expand Down