Assertion error in solve! for transient operators

[oriolcg]

Hi @amartinhuertas , @fverdugo . I'm getting an assertion error here:

GridapPETSc.jl/src/PETScNonlinearSolvers.jl

Line 161 in 119a9d0

@assert cache.op === op

when using PETScNonlinearSolver in transient problem. I get the error at the second time step, which makes me think that indeed op and cache.op are not the same. Should that be the case? If yes, what would be the best solution to this issue?

[hyun825]

Hi! I am getting the same error at the second time step after convergence of nonlinear solver in the first time step.
ERROR: AssertionError: cache.op === op

[carlodev]

Using the theta method I found out that cache.op.tθ is not updated at the second time step, while op.tθ it is; and easy fix that seems to work is:
cache.op = op

[oriolcg]

I would say that cache.op shouldn't be the same as op . I believe this assert should be removed (or not applied) for time-dependent problems.

[carlodev]

I tested the PETScNonlinearSolver to solve the Poisson problem, precisely the one in the transient tutorial. I can't see any difference in comparing the solution with the Gridap internal NLSolver. I know it is a Linear Problem, but for testing purpose seems reasonable. Just commenting cache.op === op does not produce the same result, while adding also cache.op = op produce the same results as using the internal NLSolver. But honestly, the difference between op and cache.op is not really clear to me.

[oriolcg]

Where do you add this operation: cache.op = op ?

[carlodev]

I added cache.op = op https://github.com/carlodev/GridapPETSc.jl/blob/2b345e3c3c51483139d04ad6581d9ccce3e23241/src/PETScNonlinearSolvers.jl#L162
While this is the test I made https://github.com/carlodev/GridapPETSc.jl/blob/810d2e392bf8aa0b78af7bba09e6bedf33b48630/test/Nonlinear_timedependentTests.jl

[oriolcg]

I see. The problem is that the only way PETScNonlinearSolver knows about the operator (i.e. residual) is through the cache. If this cache is not updated at each time step, the residual is computed taking into account outdated data. I think what you do here is a solution, but not sure if the most elegant/efficient. Any thoughts @amartinhuertas @fverdugo ?

[oriolcg]

Hi @amartinhuertas and @fverdugo, what is your view on this issue? should we replace cache.op by op in the solve! function?

[amartinhuertas]

Hi @oriolcg (and the others) ! Thanks for taking care of this.

I think that cache.op = op is working in your particular case because there is "compatibility" among the nonlinear operators being generated at each time step, e.g., the domain (trial FE space) and range (test FE space) of the operators, and how they are distributed among processors, match. This is certainly going to be the case in all transient simulations (as far as one does not have adaptive mesh refinement + mesh redistribution in a parallel context, we dont have this yet anyway).

However, the concern I have of just doing cache.op = op is that one opens the door to passing an op with domain/ranges not compatible with cache.op and this is for sure going to cause trouble if not further action is taken.

Ideally, what I think we should do is:

1. If cache.op and op are the same, do nothing.
2. If cache.op and op are different:
   2.1) If they have compatible domain/ranges, cache.op = op.
   2.2) If not, destroy the current cache and generate a new one.

[oriolcg]

Hi @amartinhuertas ,

I agree with this scheme, but I see one implementation issue. It is not clear to me how to check the domain/range size of op since the abstract type NonlinearOperator only has residual and jacobian as required methods. To check the size it would require to allocate the jacobian, which I think it's not a good solution. Any suggestion on this?

[oriolcg]

Hi @amartinhuertas ,

I agree with this scheme, but I see one implementation issue. It is not clear to me how to check the domain/range size of op since the abstract type NonlinearOperator only has residual and jacobian as required methods. To check the size it would require to allocate the jacobian, which I think it's not a good solution. Any suggestion on this?

@amartinhuertas @fverdugo , any view on this?

[amartinhuertas]

@amartinhuertas @fverdugo , any view on this?

As far I can see, and you point out, there are expressivity limitations in the abstract concept of a NonLinearOperator as now conceived in Gridap. Thus, if we want to implement the algorithm in #81 (comment), I am afraid we should provide a richer abstract interface for NonLinearOperators, in the spirit, of e.g., https://www.sciencedirect.com/science/article/pii/S089812211630205X