Error in Adams-Bashforth integrator implementation

[milancurcic]

The integrate method of the adams_bashforth_integrator class is implemented as:

  subroutine integrate(self, field, dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with Adams-Bashforth class scheme.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(adams_bashforth_integrator), intent(IN)    :: self  !< Actual AB integrator.
  class(integrand),                  intent(INOUT) :: field !< Field to be integrated.
  real(R_P),                         intent(in)    :: dt    !< Time step.
  integer(I_P)                                     :: s     !< Steps counter.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  do s=1, self%steps
    field = field + field%t(n=s) * (dt * self%b(s))
  enddo
  call field%update_previous_steps
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate

For simplicity, let's consider the 2nd order Adams-Bashforth (AB2). field%state then has the second dimension of size 2, which corresponds to the number of stored time levels. The solution to the ODE du/dt = F(u,t) using AB2 is:

u(n+1) = u(n) + dt*[-0.5*F(n-1)+1.5*F(n)]

In the code, this is calculated with the loop:

  do s=1, self%steps
    field = field + field%t(n=s) * (dt * self%b(s))
  enddo

which is equivalent to:

field = field + field%t(n=1) * (dt * self%b(1))
field = field + field%t(n=2) * (dt * self%b(2))

The problem that arises here is that field%t(n=1) * (dt * self%b(1)) is added to field using the overloaded operator + which performs element-wise sum for each respective array element. So, tendencies from each time level (n-2, n-1) are always added to the same time level, instead of time level n. The implementation is correct only for number of steps = 1, where the scheme degenerates to forward Euler, and there is only time level n in both the state and the tendency.

Please discuss what would be the best way to solve this. Somehow we need to add the tendency from previous time levels to the field at the present time level, something like:

  do s=1, self%steps
    field(self%steps) = field(self%steps) + field%t(n=s) * (dt * self%b(s))
  enddo

but this is not allowed because integrand is an abstract derived type, and we cannot access its state at this level in the code. Perhaps a solution would be to implement another method, perhaps field%update(), that would act to add tendencies from different time levels to state at the desired time level?

[szaghi]

@milancurcic thank you for the check. Yes this is a bug/problem. I cannot access to a computer before two days, but in the meanwhile if you find a solution, feel free to push it to the upstream.

[szaghi]

@milancurcic As I said, today and maybe tomorrow I gave not a computer under my hands, but I think a simple solution is to modify the overloaded operators like this

local_product%state = lhs%state * rhs%state

line https://github.com/Fortran-FOSS-Programmers/FOODiE/blob/master/src/tests/oscillation/type_oscillation.f90#L140

to operate only on the current time level

local_product%state = lhs%state
local_product%state(:, lhs%steps) = lhs%state(:, lhs%steps) * rhs%state(:, lhs%steps)

The reason is: operators should operate only on the current time level while the other time steps stored should be used only by the time derivative that, in turn, is used be the multi-step integrator. Does this sound reasonable?
I cannot check this until 2 days, sorry.

[milancurcic]

@szaghi Do you mean the sum instead of the product?

local_sum%steps = lhs%steps
local_sum%state(:, lhs%steps) = lhs%state(:, lhs%steps) + rhs%state(:, lhs%steps)

I think this may be a step in the right direction. However, it would not work exactly as it is, because the tendency fields can have non-zero state elements at indices different from lhs%steps. For example, if adding tendency from time level n-2 to n, the required sum would be:

local_sum%state(:, lhs%steps) = lhs%state(:, lhs%steps) + rhs%state(:, 1)

and if adding tendency from time level n-1 to n, the code would be:

local_sum%state(:, lhs%steps) = lhs%state(:, lhs%steps) + rhs%state(:, 2)

In some way, this function needs to be aware from which time level is the tendency being added to the field. One thing that I can think of right now is that the multi-step tendency field (i.e. result of the time_derivative function) always has only 1 non-zero element, because the function is evaluated for a specific time-level, while other time levels are initialized to 0 and remain such. Perhaps then we could do something like:

local_sum%state(:, lhs%steps) = lhs%state(:, lhs%steps) + sum(rhs%state,dim=2)

because the sum along the second dimension of rhs%state would yield exactly the value of the element that needs to be added to the field.

What do you think?

[szaghi]

@milancurcic

Nope, what I was thinking is something different. Let us consider a generic operator (sum, multiply, etc...). The overloaded operator should be aware of only 1 time step, the current one, namely self%state(:, self%steps). In this respect, the overloading operator procedure should be something like:

result%state(:, lhs%steps) = lhs%state(:, lhs%steps).op.rhs%state(:, lhs%steps)

Then, let us consider the problem of tendencies summation. It is done via the time derivative function that must return a type_integrand instance. This function, the time derivative one, should be the only (with the cyclic steps update procedure) public method allowed to operate on the previous steps. In particular, such a method should return its result in the current time step index of state storage array in order to sum correctly the tendency. For example:

field = field + field%t(n=1) * (dt * self%b(1))

should be computed with:

• field + ... sum operating on the field%state(:, field%steps) index;
• field%t(n=1) return a field object where the residuals results are stored into the result%state(:, field%steps) thus the previous summation (and multiplication) should work as expected.

I just thinking to an imaginary separation of the current time level from the previous one... maybe I can prove the concept implementing an actual separation storing the previous time level into a different array as a member of the type_osciallation. Maybe today I will access to a computer...

Thank you Milan for you great help!

[szaghi]

@milancurcic

I found a workstation and few minutes to test my idea... it seems to work. Briefly, I think that the problem (as you pointed out) was on the summation of tendencies due to the concurrent state(:, steps) = ... that were coming from both the time derivative procedure call and the overloaded-operators call. I have simply modified the tests implementation putting the previous times steps into a dedicated buffer and making the overloaded-operators acting on only the state buffer. This avoids the concurrent update of the state buffer, but at the cost of one more state copy (the effect of which shall be investigated into the performance analysis that will come soon).

The oscillation test now produces the expected results.

1st order

2nd order

3rd order

[milancurcic]

@szaghi Wonderful, thanks! It looks like AB2 and AB3 produce the expected solution.