Pointwise ODE solved on Q-vector

[lmolin3]

Hi,
In a diffusion-reaction problem, after applying operator splitting, the reaction step following diffusion is essentially a pointwise ODE (u is my diffused quantity, w are some internal variables):

du/dt = g(u,w)
dw/dt = f(u,w)

Currently, I’m solving this system at each entry of the T-vector, but I’d like to perform the solve at the quadrature points instead (on the Q-vector).
The quadrature point operator would be a call to another function (something like solve_ode(u_q,w_q)), because I have code for that part autogenerated from another library.
I know how to get from T-->E and back (like e.g. in the hooke miniapp), and then I could use a QuadratureInterpolator for the E-->Q step. I was wondering if there's an existing MFEM functionality already implemented to compute the Q-->E transformation directly; or would I need to implement the integration loop using the DofToQuad map B, to compute the E-vector?. (i saw some examples in TransferKernels::Prolongation and ::Restriction, and in the Apply3D kernels in the hooke miniapp).

Alternatively, if I define a quadrature function to perform my ode solve and provide it to a DifferentiableOperator (input/output would use VALUE as FieldOperator), would this handle all the above when calling Mult()? Or is there any issue with that or a better alternative?

[termi-official]

I’m solving this system at each entry of the T-vector, but I’d like to perform the solve at the quadrature points instead (on the Q-vector).

Any reason why you would prefer the Q vector here, or why you would like to solve something at the quadrature points?

In reaction diffusion splittings using e.g. with the Lie-Trotter method with Lagrangian ansatz functions, the reaction dofs are typically not on the quadrature points but really located the nodes of your ansatz functions. If you would use a quadrature rule now to locate the variables, then you have redundant solves, as the reactions at the element boundary are solved multiple times form. Does that make sense or do you apply something more special?

[lmolin3]

Thanks for the clarification, that makes sense. For now, we're using CG basis functions (though DG could be an option in the future) and a second-order splitting, with the mid step still being the decoupled reaction one.

I agree that solving the reaction ODEs at the DOFs is the most natural and efficient approach, since it avoids redundant evaluations and aligns with the nodal representation of the field. For now, we're fine with what’s implemented on T-vectors, but I was just curious about solving on the Q-vector for more special cases, like spatially heterogeneous parameters defined per quadrature point, or other scenarios that might benefit from using Q-vectors, and might come up in the future.

[termi-official]

I typically also move the parameter fields on the nodes in one way or another (e.g. by having the raw data on the nodes or by some interpolation/projection onto the nodal space). You typically just need the Q vectors for e.g. IMEX or explicit schemes - or at least I never found a case where I used Q vectors in problems which are RD split. At least IMEX schemes tend to perform worse (see e.g. https://www.sciencedirect.com/science/article/pii/S0377042715004677 for a formal analysis on Bidomain problems, which are degenerate reaction diffusion problems, but the analysis should carry over to your formulation). If you have DG ansatz, then you can still can use the provided nodal space.