Solving ODE in each node

[rogeraubert]

Hello,
I am trying to figure out, how to solve a diffusion-reaction PDE, that is coupling a diffusion kernel with an ODE for each point in space. Do you have a short example for something similar? I think, the Crane application is doing this, but they don't use advanced ODE solvers, they rather integrate the reaction into the PDE, which is not optimal for stiff reactions.

At first glance I thought that this example would demonstrate such a coupling:
https://mooseframework.inl.gov/getting_started/examples_and_tutorials/examples/ex18_scalar_kernel.html
But apparently the example solves an ODE for two scalar variables x(t) and y(t) independent from solving the PDE for the field "diffused". Am I wrong?

This brings me to another question. I am confused what is meant by a SCALAR variable in MOOSE. To my understanding, there is the general notion of either scalar or vector fields (having a one-dimensional or three-dimensional value in each point in space). MOOSE seems to use SCALAR for a single value (not a field). This is absolutely correct, of course, but I wonder if scalar fields can also be defined as SCALAR variables in MOOSE? The difference was not quite obvious to me when I went through various tutorials and workshop slides.

Thanks
Roger

[dschwen]

SCALAR is a variable type that does not live on the mesh. Think of a "global" quantity. For ODEs at each node look up NodalKernel.

[dschwen]

Btw. asking two questions in one post is not really suitable for GitHub discussions. If multiple people address different questions in your post. We rely on the asker to mark an answer so that the discussion isn't listed as "unanswered" in the overview.

[rogeraubert]

Sorry for that. I wanted to mark my question as answered when somebody replies that CVODE can be used to implement a TimeIntegrator for MOOSE, or why it cannot or should not be used.

[rogeraubert]

Thank you. As far as I see, the ReactionNodalKernel is just adding a term -k*u, and the time integration is still performed by low-order methods like CrankNicolson or ImplicitEuler. Shouldn't one use variable-order variable-timestep methods like CVODE for integrating systems of stiff reactions?

(The other question about SCALAR is now very clear)

Thanks
Roger

[keniley1]

I'm the developer of Crane, which you mentioned in your original post. You're right about what Crane is doing - when run as part of a set of diffusion-reaction equations in 1D-3D, the reaction terms are integrated into the full PDE and thus use the same time integration schemes as the diffusion terms. I've used Crane (coupled to Zapdos) pretty extensively for systems of several hundred reactions in plasma physics problems without issue, though I haven't actually measured the stiffness of my system of equations.

I'm willing to look into alternative formulations if they exist though! Do you have any examples of what you're looking for? Correct me if I'm wrong, but I think what you're trying to do is first solving the diffusion equation with finite elements for a given timestep, and then solving a set of ODEs at each node using more complex time integration schemes over the same timestep, and then iterating back and forth until the simulation completes -- is this correct? I think there might be a way to do something like this with the external_petsc_problem example in MOOSE. Not sure if it's possible, but maybe I can run PETSc as a MultiApp like in that example to handle the ODEs at each node? Maybe someone on the MOOSE/PETSc team can chime in with advice here -- I've never done something like this myself. I know that CVODE specifically is lumped into TS in PETSc, and last I checked the TS library is not accessible to MOOSE applications.

[rogeraubert]

Dear Shane,
your reply helped me very much to think about how to proceed. I don't have a solution, but I am only exploring the options.

What you describe as solving the FE diffusion time step separately from the (stiff) reaction ODEs and iterating, where only the reaction ODEs make use of advanced time integrators for stiff systems, is very close to the OpenFoam approach:
https://develop.openfoam.com/Development/openfoam/-/tree/master/applications/solvers/combustion/reactingFoam and https://develop.openfoam.com/Development/openfoam/-/blob/master/applications/solvers/combustion/reactingFoam/YEqn.H
Line #13 of YEqn.H updates the reaction rate (using Seulex or other integrators) and line #29 adds the rate to the transport equation.

OpenFoam seems to be successful with this approach although there may be situations where decoupling the diffusion from the reaction is not optimal.

I think the MOOSE framework can be even more efficient by solving the system fully implicitly with AD Jacobians. https://github.com/LLNL/sundials/blob/master/doc/cvode/cv_examples.pdf includes some examples how CVODE can integrate a discretized PDE, although I didn't find a stiff PDE in those examples. But I think it should be possible to add CVODE as TimeIntegrator to Moose and explore the advantages or disadvantages for solving reaction-diffusion PDEs fully implicitly, or are there clear reasons not to try this?

Of course, also BDF2 in MOOSE can integrate stiff systems without stability issues, but it is missing the linked order/timestep adaptivity from CVODE. The question is, which type of reactions in relation to which mesh sizes benefit most from the variable order schemes of CVODE. I'd like to have the choice of more time integrators than those currently available in MOOSE.

Thanks
Roger

[rogeraubert]

I understand the previous answers such that MOOSE doesn't allow the implementation of variable-order variable-timestep TimeIntegrators that are suitable for efficient integration of stiff reaction-diffusion systems. Nothing more efficient that BDF2 seems to be possible.