search

gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
695 stars 95 forks source link

ode_start() in Generalized_Alpha results into singular system #986

Closed oriolcg closed 5 months ago

oriolcg commented 6 months ago

Hi @AlexandreMagueresse, I'm using the new version of the ODEs and I'm getting a singular exception when solving a problem where the mass matrix is only defined in part of the domain. The overall system is well posed as the stiffness and/or damping have values on the full domain. However, when initializing the solution, the GeneralizedAlpha solves a "mass-only" problem in ode_start https://github.com/gridap/Gridap.jl/blob/92f38a29c2cb47b4ecb6d29955a99e88cad892a3/src/ODEs/ODESolvers/GeneralizedAlpha2.jl#L90

Before it was taking an initial acceleration as input. I don't see an easy way to address this issue that would not involve sending different operator for the initialization (or just pass an initial acceleration as before). Do you have any ideas on this?

Thanks!

Tagging @janmodderman and @shagun751 to follow this discussion.

JordiManyer commented 6 months ago

Hi @oriolcg . Thanks for reporting this. We are currently preparing the new release in PR #985 . We also found a couple things that need to be changed. We'll take this into account.

AlexandreMagueresse commented 6 months ago

Hi @oriolcg . Thank you for pointing this out and sorry for the delay in replying. I agree that in this scenario, the n-th order derivative (velocity for GeneralizedAlpha1, acceleration for GeneralizedAlpha2) has to be provided in another way, i.e. as you suggest, either manually, or as the solution of another FEOperator.

I have restored the possibility to directly provide these initial derivatives manually when calling solve as before (see this commit 834289d96fed6f88d57b5438384f976d24e98a26). You can take a look at the following tests (these links directly point to the lines of interest): Order1ODETests.jl, Order2ODETests.jl, HeatEquationScalarTests.jl or SecondOrderEquationTests.jl. I will also add a note in docs.