opm-simulators icon indicating copy to clipboard operation
opm-simulators copied to clipboard

Make CPR reuse setup option only apply to CPR

Open jakobtorben opened this issue 1 year ago • 3 comments

The cpr-reuse-setup option is currently not limited to CPR and is being used for all preconditioners. The default option is to recreate the linear solver every 30th linear iteration. The motivation behind this option is that the AMG hierarchy depends on the values of the matrix, such that as the values of the matrix changes, the preconditioner should be updated. For preconditioners where the construction of the preconditioner only depends on the matrix structure, it is never needed to recreate the preconditioner, as the matrix structure should be constant throughout the simulation.

I therefore suggest that we make this option only apply to CPR. This will also make the effect of the option match the name.

jakobtorben avatar Jun 27 '24 10:06 jakobtorben

jenkins build this please

multitalentloes avatar Jun 27 '24 10:06 multitalentloes

image Here is a comparison on Norne where I used the ILU0 preconditioner, that only depends on the matrix structure. To the left is the simulation result when never recreating the linear solver and to the right is the result when recreating the preconditioner every linear solve. As expected, the number of timesteps, Newton iterations and linear iterations is the same.

jakobtorben avatar Jun 27 '24 10:06 jakobtorben

Note that this a quick fix to make the option accurate and avoid recreating the linear solver for other preconditioners. In a later stage we should introduce a flag that determines if a preconditioner should be recreated and see what can be resued between linear solves.

jakobtorben avatar Jun 27 '24 11:06 jakobtorben

@jakobtorben : Just a quick clarification on this particular point

For preconditioners where the construction of the preconditioner only depends on the matrix structure, it is never needed to recreate the preconditioner, as the matrix structure should be constant throughout the simulation.

Do you consider upwind directions to be a part of the matrix structure? Those directions could change during the simulation run.

bska avatar Aug 19 '24 09:08 bska

@bska I am not familiar with how the upwind directions affect the matrix structure in OPM. What I meant here is that certain preconditioners only depend on the matrix sparsity pattern, such that they only need to be recreated if this changes. If the upwind directions changes the sparsity pattern of a matrix, then such preconditioners should be recreated.

jakobtorben avatar Aug 20 '24 07:08 jakobtorben

@bska I am not familiar with how the upwind directions affect the matrix structure in OPM.

I think they don't. For any two neighbour cells i and j, the flux between them depends on both cells' pressures, so both (i, j) and (j, i) elements of the Jacobian are nonzero. Due to upwinding, there will be some zeroes inside one of these blocks and not in the other, but the blocks as a whole are always nonzero because of the pressure dependency.

So that should not be a hindrance. However, the PR must still be changed as I requested earlier before it can be merged.

atgeirr avatar Aug 20 '24 08:08 atgeirr

What I meant here is that certain preconditioners only depend on the matrix sparsity pattern, such that they only need to be recreated if this changes

Right, that's what I thought. Then this discussion probably comes down to whether the sparsity pattern comprises the potential non-zero matrix elements or the actual non-zero matrix elements. The potential non-zero elements are static and depend only on the grid and the well/reservoir connections. The actual non-zero elements also depend on the dynamic state and from which side of a connection we pick the property values (i.e., the upwinding).

bska avatar Aug 20 '24 08:08 bska

What I meant here is that certain preconditioners only depend on the matrix sparsity pattern, such that they only need to be recreated if this changes

But these (Jacobi SOR SSOR) are only very seldom used. All the rest actually depend on the values (ILU*, AMG, etc.)

blattms avatar Aug 20 '24 08:08 blattms

Ok, thanks for clarifying. I will make the necessary changes to get this in now. Going for the bool hasPerfectUpdate() solution, that defaults to False. @blattms The preconditioners where the update depends on the values (such as ILU, AMG) will still be updated but not recreated.

jakobtorben avatar Aug 20 '24 08:08 jakobtorben

Ok, thanks for clarifying. I will make the necessary changes to get this in now. Going for the bool hasPerfectUpdate() solution, that defaults to False.

Sounds good. There are some things to consider along that line though:

  • This should probably be a constexpr property for each preconditioner, since it only should depend on the class, not the concrete instance. However we may then need the virtual bool hasPerfectUpdate() in addition to actually use it in contexts where we only have a base class pointer or reference.
  • For the wrappers that take Dune preconditioners we must think carefully about how to define the property. The RebuildOnUpdatePreconditioner should be true, for the DummyUpdatePreconditioner I am not 100% sure: Probably true here as well, as we anyway leave it to the code that employs this wrapper to only use it when an empty update() is indeed correct.

atgeirr avatar Aug 20 '24 09:08 atgeirr

Just a note: There is a difference between ILU and AMG here. For the former the sparsity pattern is fixed. For the latter the pattern on coarse levels actually depends on the values of the fine level.

blattms avatar Aug 20 '24 12:08 blattms

For the former the sparsity pattern is fixed. For the latter the pattern on coarse levels actually depends on the values of the fine level.

Good point! After all, this is why we have the feature to rebuild it in the first place.

atgeirr avatar Aug 20 '24 12:08 atgeirr

Here is a proposed solution that adds a method to preconditionersWithUpdate, which defines if the preconditioner has a perfect update. I opted for forcing the authors of preconditioners to define this property. For all CPR/AMG and two-level preconditioners I set it to false. There is probably a possibility of reusing more of the preconditioner/smoother in AMG, but I did not address this complexity here. I am open to other name suggestions as hasPerfectUpdate is not so clear.

I tested all available preconditioners in OPM, here is the result of calling hasPerfectUpdate():

ilu0 true dilu true DuneILU true GS true ILUn true Jac true ParOverILU0 true SOR true SSOR true isai true

CUDILU true CUDILUFloat true CUILU0 true CUILU0Float true CUJac true OPMCUILU0 true

cprw false cpr false cprt false cpr_quasiimpes false cpr_trueimpes false cpr_trueimpesanalytic false

jakobtorben avatar Aug 20 '24 15:08 jakobtorben

jenkins build this please

jakobtorben avatar Aug 20 '24 15:08 jakobtorben

jenkins build this please

bska avatar Aug 20 '24 16:08 bska

I tested all available preconditioners in OPM, here is the result of calling hasPerfectUpdate():

AMG is missing from the list, but I think it should be false, and that this code would indeed return false.

atgeirr avatar Aug 21 '24 06:08 atgeirr

Thanks for the work and refactoring, merging!

I you find a better name for hasPerfectUpdate() I will listen. Also, better name for update() might make it clearer what is expected.

atgeirr avatar Aug 21 '24 06:08 atgeirr

I was not able to get AMG working as an independent preconditioner in OPM. It is used in CPR preconditioners (twolevelmethodcpr) where it is set to false.

I did not come up with any better name, so I just left it as is.

Thanks for merging @atgeirr!

jakobtorben avatar Aug 21 '24 07:08 jakobtorben