Skip to content

Store linear invariants matrix in (Conservative)PDSProblem#163

Merged
ranocha merged 13 commits intoNumericalMathematics:mainfrom
jmbender:jb/linearInvariants
May 30, 2025
Merged

Store linear invariants matrix in (Conservative)PDSProblem#163
ranocha merged 13 commits intoNumericalMathematics:mainfrom
jmbender:jb/linearInvariants

Conversation

@jmbender
Copy link
Collaborator

@jmbender jmbender commented May 26, 2025

The use of the Sandu projection in #162 requires a matrix containing the linear invariants. A PDSFuntion or ConservativePDSFunction should be able to store this matrix.

  • Added the field lin_invariants to the structs PDSFunction and ConservativePDSFunction
  • Store the corresponding matrix of linear invariants in each problem contained in PDSProblemLibrary.jl
  • Add tests

@coveralls
Copy link

coveralls commented May 26, 2025

Pull Request Test Coverage Report for Build 15341933359

Details

  • 11 of 11 (100.0%) changed or added relevant lines in 1 file are covered.
  • No unchanged relevant lines lost coverage.
  • Overall coverage remained the same at 97.768%

Totals Coverage Status
Change from base Build 15070539288: 0.0%
Covered Lines: 1971
Relevant Lines: 2016

💛 - Coveralls

Copy link
Member

@ranocha ranocha left a comment

Choose a reason for hiding this comment

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

I guess you are working with @SKopecz on this, correct? What are the plans for this feature - shall modified Patankar methods be able to use a single given linear invariant to rescale the setup to conserve it instead of the default linear invariant with weight vector $(1, \dots, 1)^T$?

Depending on the goal, we may also consider using https://github.com/JuliaArrays/FillArrays.jl (Trues or Ones) as default for conservative problems.

@jmbender
Copy link
Collaborator Author

@ranocha Yes, I'm working with @SKopecz on it. Our motivation was to implement the Sandu optimizer. Therefore, linear invariants are needed. At first, we hadn't considered MPRK methods in this context.

However, your remark about the default values could make sense. Do you mean that for conservative problems, we set the linear invariants to one, but this can be changed if you have a problem that preserves additional linear invariants?

@SKopecz SKopecz changed the title Jb/linear invariants WIP: Store linear invariants matrix in (Conservative)PDSProblem May 27, 2025
@SKopecz SKopecz marked this pull request as draft May 27, 2025 09:52
@SKopecz
Copy link
Collaborator

SKopecz commented May 27, 2025

What are the plans for this feature - shall modified Patankar methods be able to use a single given linear invariant to rescale the setup to conserve it instead of the default linear invariant with weight vector ( 1 , … , 1 ) T ?

The plan is to provide the matrix of linear invariants directly with the problems contained in PDSProblemLibrary.jl. These matrices than can be used to test the Sandu projection introduced in #162. So lin_invariant should be an SMatrix instead of an SVector.

To implement what you suggested, we would only need a vector. I think this should be implemented separately.

@ranocha
Copy link
Member

ranocha commented May 27, 2025

I see, thanks! This makes sense. Would it make sense to rename the argument to linear_invariants (plural) in this case to indicate that it should be a matrix (SMatrix for the small problems) collecting the coefficients of all linear invariants?

@SKopecz
Copy link
Collaborator

SKopecz commented May 27, 2025

I think so too.

@codecov
Copy link

codecov bot commented May 27, 2025

Codecov Report

All modified and coverable lines are covered by tests ✅

📢 Thoughts on this report? Let us know!

@jmbender jmbender requested a review from ranocha May 28, 2025 13:29
@SKopecz SKopecz marked this pull request as ready for review May 28, 2025 13:31
@SKopecz SKopecz changed the title WIP: Store linear invariants matrix in (Conservative)PDSProblem Store linear invariants matrix in (Conservative)PDSProblem May 28, 2025
Copy link
Member

@ranocha ranocha left a comment

Choose a reason for hiding this comment

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

Thanks! This looks already quite good to me.

@@ -43,14 +43,15 @@ The functions `P` and `D` can be used either in the out-of-place form with signa
struct PDSProblem{iip} <: AbstractPDSProblem end
Copy link
Member

Choose a reason for hiding this comment

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

Since linear_invariants is a new feature, it should be documented at least in the docstring above. You can mention, for now, that it is an experimental feature so that we may change its behavior without having to release a breaking version of PositiveIntegrators.jl. If you do so, could you please also open an issue allowing us to track the status of this experimental feature?

Copy link
Collaborator

@SKopecz SKopecz May 29, 2025

Choose a reason for hiding this comment

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

linear_invariants needs to be added at the beginning of the docstring as well.

PDSProblem(P, D, u0, tspan, p = NullParameters();
               p_prototype = nothing,
               analytic = nothing,
               std_rhs = nothing,
               linear_invariants = nothing)

In addition, linear_invariants needs to be mentioned in the docstring of ConservativePDSProblem.

Copy link
Collaborator

@SKopecz SKopecz May 29, 2025

Choose a reason for hiding this comment

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

Since linear_invariants is a new feature, it should be documented at least in the docstring above. You can mention, for now, that it is an experimental feature so that we may change its behavior without having to release a breaking version of PositiveIntegrators.jl.

It's not clear to me why the introduction of linear_invariants could be breaking. It's an optional feature and every code that ran in the past will also run in the future. Is it breaking because the PDSFunction and ConservativePDSFunction structs changed? But these structs are only used internally and are not documented.

Copy link
Member

Choose a reason for hiding this comment

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

What Hendrik means is not that adding this feature is breaking, but rather that in the future it might be necessary or desirable to change this feature again (e.g. changing the API), which could break code then. Adding a note that this is experimental, allows us to do that without the need to care too much about creating a breaking release.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

linear_invariants needs to be added at the beginning of the docstring as well.

PDSProblem(P, D, u0, tspan, p = NullParameters();
               p_prototype = nothing,
               analytic = nothing,
               std_rhs = nothing,
               linear_invariants = nothing)

In addition, linear_invariants needs to be mentioned in the docstring of ConservativePDSProblem.

done

src/proddest.jl Outdated
the production-destruction representation of the ODE, will use this function
instead to compute the solution. If not specified,
a default implementation calling `P` and `D` is used.
-`linear_invariants`: Specifies the linear invariants of the problem as a `StaticArrays.SMatrix`.
Copy link
Collaborator

Choose a reason for hiding this comment

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

linear_invariants must not be an SMatrix. I'd suggest the following:

`linear_invariants`: The rows of this matrix contain the linear invariants of the ODE. Certain solvers or callbacks require this matrix.

The only solver/solution technique I can think of right now that requires the matrix of linear invariants is the Sandu projection. But I'd suggest to keep the docstring fairly general.

Furthermore, currently linear_invariants must not even be a matrix. Should we throw an error if it's not?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I changed the text to the line above + note to experimental feature

Copy link
Member

@ranocha ranocha left a comment

Choose a reason for hiding this comment

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

Thanks a lot!

@ranocha ranocha merged commit 1e80944 into NumericalMathematics:main May 30, 2025
10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants