[Michael@UW]

Hi deal Foamers,

I have a couple of questions about the ODE solver in OpenFOAM.

Question 1:

In the file test-ODE.C


I don't see the purpose of the following two lines:

    
// Print the evolution of the solution and the time-step
    scalarField dyStart(ode.nEqns());
    ode.derivatives(xStart, yStart, dyStart);

I understand the derivatives are calculated through the function ode.derivatives and stored in dyStart, but dyStart
is never used again later.
The ODE is solved by calling the function

odeSolver->solve(x, xEnd, y, dxEst);

If we look at the solve function for an instantiated solve, e.g. RKCK45, odes_.derivatives is called there. So it makes no sense to me to call ode.derivatives above.

Question 2:
In the ODESolver.C, the following two solve functions call each other recursively, so there seems to be a dead loop. Can anyone explain how they work?

void Foam::ODESolver::solve
(
    scalar& x,
    scalarField& y,
    scalar& dxTry
) const
{
    stepState step(dxTry);
    solve(x, y, step);
    dxTry = step.dxTry;
}
 
 
void Foam::ODESolver::solve
(
    scalar& x,
    scalarField& y,
    stepState& step
) const
{
    scalar x0 = x;
    solve(x, y, step.dxTry);
    step.dxDid = x - x0;
}

Feel free to answer or discuss any of these questions. Thank you!

Michael

[dlahaye]

Many thanks for raising these issues!

Q1: not sure. Possibly it helps to write down the problem more formally, as \dot{y} = f(y(x), x) supplied with initial guess y(0) = y0 and try to reason on an algorithm from here. You raise the issue of derivatives. I assume that derivatives are used by implicit algorithms only. RK45 is an explicit algorithm. Not sure whether this helps.

Q2: not sure either. It is possible that the size of the step is adapted inbetween calls of the solve function?

I am happy to discuss further.

[Michael@UW]

[QUOTE=dlahaye;860062]Many thanks for raising these issues!

Thank you for the discussion.

Originally Posted by dlahaye

Q1: not sure. Possibly it helps to write down the problem more formally, as \dot{y} = f(y(x), x) supplied with initial guess y(0) = y0 and try to reason on an algorithm from here. You raise the issue of derivatives. I assume that derivatives are used by implicit algorithms only. RK45 is an explicit algorithm. Not sure whether this helps.

It seems to make sense, for example, if we want to call

Foam::scalar Foam::Euler::solve
(
    const scalar x0,
    const scalarField& y0,
    const scalarField& dydx0,
    const scalar dx,
    scalarField& y
) const
{
    // Calculate error estimate from the change in state:
    forAll(err_, i)
    {
        err_[i] = dx*dydx0[i];
    }
 
    // Update the state
    forAll(y, i)
    {
        y[i] = y0[i] + err_[i];
    }
 
    return normalizeError(y0, y, err_);
}

directly which requires the derivatives dydx0, of course we also need to specify the solver type (Euler) explicitly.

By contrast, if we call

void Foam::Euler::solve
(
    scalar& x,
    scalarField& y,
    scalar& dxTry
) const
{
    adaptiveSolver::solve(odes_, x, y, dxTry);
}

which will eventually call the first solve function, but the difference is that the derivatives are calculated during this call process somewhere, so there is no need to call ode.derivatives outside these functions. But for the ODE-test.C, it is confusing, since the derivatives are not used explicitly.

Originally Posted by dlahaye

Q2: not sure either. It is possible that the size of the step is adapted inbetween calls of the solve function?

Further looking at the code of ODESolver class, it is an abstract class, so the inter-call between

void Foam::ODESolver::solve
(
    scalar& x,
    scalarField& y,
    scalar& dxTry
) const
{
    stepState step(dxTry);
    solve(x, y, step);
    dxTry = step.dxTry;
}

and

void Foam::ODESolver::solve
(
    scalar& x,
    scalarField& y,
    stepState& step
) const
{
    scalar x0 = x;
    solve(x, y, step.dxTry);
    step.dxDid = x - x0;
}

should never happen, since all instantiations have implemented these two virtual functions, so if the first solve function is called, the second function of the instance instead of that of the base (ODESolver) will be called.

Originally Posted by dlahaye

I am happy to discuss further.

Your response helped me to think more deeply. Hope more ideas can come in and clarify these issues.

[dlahaye]

Dear Michael,

I am confused about your most recent reply. I miss the head, body and tail of the points you raise.

Is it valuable for you to re-iterative on the points you raise? Possibly it allows me to provide input.

Thanks, Domenico.

[Michael@UW]

Sorry for the confusion. I came out of these thoughts when I was trying to comment on your comments, but I did not address your comments actually. Let me rephrase my thoughts.

(Q1) First of all, I don't think the derivatives are needed in the code of ODE-Test.C, because odeSolver->solve(x, xEnd, y, dxEst) takes care of the derivatives inside the function.

However, if you call another version of solve function, such as odeSolver->solve(x0, y0, dydx0, dx, y), then we need to call ode.derivatives(x0, y0, dxdy0) explicitly. This responds to your comment: "Possibly it helps to write down the problem more formally, as \dot{y} = f(y(x), x) supplied with initial guess y(0) = y0 and try to reason on an algorithm from here."

By the way, both the implicit and explicit ODE solver use the derivatives for the integration.

Again, as a test code, the appearance of the unnecessary ode.derivatives is confusing.

(Q2) "It is possible that the size of the step is adapted in-between calls of the solve function" Yes, the size of the step can be adapted by these in-between calls, but it does not explain how to exit the dead-loop call.

My explanation is that these two functions are virtual functions in the base class (ODESolver), so if one of the two solve functions is (must be) overridden in the derived class, this call-loop will be broken. The design seems subtle and I cannot fully understand.

Hope this makes more sense to you. Please feel free to share your inputs and post more discussion.