Refactoring of GridapODEs #965
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
…elated to transient FE spaces
…d variable names more consistent, added tests
…rs with optimised computation of residuals/jacobians
…n an AffineOperator whenever possible, removed BackwardEuler
…a generic Runge Kutta solver. Explicit Runge-Kutta working and tested, Diagonally-Implicit coming
…screteODEOperator whenever possible
…ests for 2nd order ODEs
…, wrote IMEXRK solver in the linear case, introduced more consistent notations across all ODE solvers, restructured IMEX operator to sum with operator of one degree less
Quasilinear and added Semilinear (for IMEX RK), fixed filters
…ts for second-order and IMEX TransientFEOperators
…or issues. Tests still in progress
…itrary-order time derivative, unified order of arguments, need to check AD with second-order
…acobian matrix is not constant across steps). This implied some restructuring
AlexandreMagueresse
unassigned oriolcg, JordiManyer, amartinhuertas, santiagobadia and tamaratambyah
oriolcg
approved these changes
Dec 15, 2023
oriolcg
requested changes
Dec 15, 2023
…ethods (start/update/finish), reuse of factorisations whenever possible, additional constructors for TransientFEOperators to match the existing ones, imported more tableaus and added generic constructors for low-order EXRK and DIRK schemes, introduced a StageOperator type with linear and nonlinear versions to save redundancy
|
Hi @santiagobadia @amartinhuertas @JordiManyer @oriolcg @tamaratambyah, I am finally done with all the changes I had in mind. I think the pull request is ready for review. The only thing I would like to add is a markdown file |
… merging cell vectors/matrices, fixed LaTeX typos
When the jacobian of a residual is obtained through automatic differentiation, the return type is BlockArray{<:SubArray} and the behaviour of testvalue does not allow broadcasting operations between BlockArray{<:AbstractMatrix} and BlockArray{<:SubArray}. This function returns a matrix of size a P-dimensional array where each dimension has length 0, i.e., (0, ..., 0).
…matrices if they have contributions
|
Hi @AlexandreMagueresse, I'm encountering the following error when using the I'm solving the transient Navier-Stokes equations where, indeed, the mass and stiffness matrices don't have the same sparsity pattern. Have you thought about this issue? Is there a way to address this? Or should I just use another operator? |
…tern as the overall jacobian. Similarly, allocate the jacobian of an IMEX operator with the sparsity structure of the sum of the implicit and explicit jacobians
|
Hi @oriolcg, thank you for finding that issue. My last commit should fix it. |
|
Hi @AlexandreMagueresse, now I'm getting another error when using AD with multifield. I'm not sure if this is an issue only on parallel or if in serial also appears... It's failing in the allocation of the jacobians when using a SemilinearOperator, I think it is expecting a |
|
@oriolcg AD is not supported by GridapDistributed yet (and never was). We have a working prototype, but it's not part of the library. Although the error you point out is indeed caused by accessing properties of a So long story short, this is an issue in distributed only, and it has nothing to do with this PR. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
As the ODE module of Gridap is expanding, its structure is becoming a little unclear and redundant. I propose the following major changes to the API. The general framework of the module is described here. I made other minor changes internally.
ODE operators
NonlinearODE,QuasilinearODE,SemilinearODEandLinearODENonlinearODE:QuasilinearODE:SemilinearODE:LinearODE:Transient FE spaces, cell fields and operators
TransientTrialFESpaceandTransientMultiFieldFESpaceare now subtypes ofAbstractTransientTrialFESpace, which is itself a subtype of `FESpace.TransientFEOperator(res, jacs, trial, test)andTransientFEOperator(res, trial, test; order)for the version with AD. The residual is expected to have the signatureresidual(t, u, v).TransientQuasilinearFEOperator(mass, res, jacs, trial, test)andTransientQuasilinearFEOperator(mass, res, trial, test; order)for the version with AD. The mass and residual are expected to have the signaturesmass(t, u, dtu, v)andresidual(t, u, v), i.e.massis written as a linear form of the highest-order time derivativedtu. Since in this setting the mass matrix is supposed to depend on lower-order time derivatives,uis provided for the nonlinearity of the mass matrix.TransientSemilinearFEOperator(mass, res, jacs, trial, test; constant_mass)andTransientSemilinearFEOperator(mass, res, trial, test; order, constant_mass)for the version with AD (no AD is needed for the mass term). The mass and residual are expected to have the signaturesmass(t, dtu, v)andresidual(t, u, v), where here againdtuis the highest-order derivative, i.e.massis specified as a linear form ofdu.TransientLinearFEOperator(forms, res, jacs, trial, test; constant_forms)andTransientLinearFEOperator(forms, res, trial, test; constant_forms)for the version with AD (in fact no AD is needed at all). The forms and residual are expected to have the signaturesform_k(t, dtku, v)andresidual(t, v), i.e.form_kis a linear form of theu.TransientFEOperator:res(t, u, v) = ∫(∂t(u) ⋅ v + ∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩTransientQuasilinearFEOperator:mass(t, u, dtu, v) = ∫(dtu ⋅ v) * dΩandres(t, u, v) = ∫(∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩTransientQuasilinearFEOperator:mass(t, dtu, v) = ∫(dtu ⋅ v) * dΩandres(t, u, v) = ∫(∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩTransientLinearFEOperator:mass(t, dtu, v) = ∫(dtu ⋅ v) ) * dΩ,stiffness(t, u, v) = ∫(∇(u) ⋅ ∇(v)) * dΩandres(t, u, v) = ∫(- f(t) ⋅ v) * dΩ. Important: note the minus sign here. We could change that if it is more intuitive to haveresequal to the forcing term rather than its opposite.jacobian!.ODEOpFromFEOp, which is a subtype ofODEOperator. In that senseTransientFEOperatordoes not fully implement the interface of anFEOperator, it is rather at the intersection ofFEOperatorandODEOperator.ODESolutioninterface,FESolutionis now an abstract type and the equivalent ofGenericODESolutionisGenericFESolution.solve(solver, operator, u0, t0, tF), it is nowsolve(solver, operator, t0, tF, u0). Similarly, when iterating over an FE solution, the output is(t_n, u_n)rather than(u_n, t_n). These changes follow from using the order(space, time, functional space)in the whole module.ODE solvers
DiscreteODEOperatoris created (it is a subtype ofNonlinearOperator). It represents the discretisation of the time derivative performed by the ODE solver. However, whenever possible, an ODE solver should instantiate aLinearDiscreteODEOperatorinstead, which behaves like anAffineOperator. This is always possible when an explicit ODE solver is applied to a quasilinear ODE, or when a (diagonally-) implicit ODE solver is applied to a linear ODE.DiscreteODEOperatorneeds to be updated across consecutive time steps, even when theDiscreteODEOperatoris actually aLinearDiscreteODEOperatorand the solver for this system of equation is aLinearSolver. In particular, thesolve!function has been overwritten to recompute the jacobian matrix and its factorisation. However, under special circumstances, the jacobian matrix may be constant, and it is possible to indicate so at theLinearDiscreteODEOperatorlevel.Notes
Additionally to these changes to the ODE module, I made the following changes.
axpy_entries!function toAlgebrathat enables adding two matrices, with an efficient specialisation for sparse matrices.MultiFieldCellField. The issue was that the return type of the jacobian provided byForwardDiffwasBlockMatrix{<:SubArray}, and its correspondingtestvaluewas wrong (it wasArray{T, 0}instead ofArray{T, P}, wherePis the dimension of the quantity being differentiated, i.e.2for jacobians).MultiFieldCellFieldandMultiFieldFEFunctionnow implementBase.lengthso that they can be zipped and mapped over.To discuss
I would like to discuss the following changes that directly impact the user.
FESolution, the iterator returns the time first and then theFEFunction. It used to return theFEFunctionfirst and then the time.TransientLinearFEOperator, the order is stiffness, then mass, then residual, where residual is the opposite of the forcing term. This has the advantage of keeping the same order as the other constructors (nonlinear, semilinear and quasilinear) but is quite unnatural from a user perspective. Another question is the order in which the jacobians are given, and whether the linear forms are constant. I am not happy about this but I have not been able to come up with a satisfying answer for now.I have some final minor questions.
const GridapODEs = ODEsuseful? See here.(space::FESpace)(t) = evaluate(space, t)is only defined for julia versions above 1.3. Do we need to keep compatibility with julia versions lower than 1.3? Also, why is this not being defined for earlier versions? See here.