Closed milancurcic closed 9 years ago
@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.
@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
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.
@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?
@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!
@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.
@szaghi Wonderful, thanks! It looks like AB2 and AB3 produce the expected solution.
The
integrate
method of theadams_bashforth_integrator
class is implemented as: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:In the code, this is calculated with the loop:
which is equivalent to:
The problem that arises here is that
field%t(n=1) * (dt * self%b(1))
is added tofield
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:
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, perhapsfield%update()
, that would act to add tendencies from different time levels to state at the desired time level?