{"action":"created","issue":{"url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/18","labels_url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/18/labels{/name}","comments_url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/18/comments","events_url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/18/events","html_url":"https://github.com/Fortran-FOSS-Programmers/FOODiE/issues/18","id":105049062,"number":18,"title":"Error in Adams-Bashforth integrator implementation","user":{"login":"milancurcic","id":4133310,"avatar_url":"https://avatars.githubusercontent.com/u/4133310?v=3","gravatar_id":"","url":"https://api.github.com/users/milancurcic","html_url":"https://github.com/milancurcic","followers_url":"https://api.github.com/users/milancurcic/followers","following_url":"https://api.github.com/users/milancurcic/following{/other_user}","gists_url":"https://api.github.com/users/milancurcic/gists{/gist_id}","starred_url":"https://api.github.com/users/milancurcic/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/milancurcic/subscriptions","organizations_url":"https://api.github.com/users/milancurcic/orgs","repos_url":"https://api.github.com/users/milancurcic/repos","events_url":"https://api.github.com/users/milancurcic/events{/privacy}","received_events_url":"https://api.github.com/users/milancurcic/received_events","type":"User","site_admin":false},"labels":[{"url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/labels/bug","name":"bug","color":"fc2929"}],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":1,"created_at":"2015-09-05T21:07:15Z","updated_at":"2015-09-06T11:39:29Z","closed_at":null,"body":"The `integrate` method of the `adams_bashforth_integrator` class is implemented as:\r\n\r\n\r\n```fortran\r\n  subroutine integrate(self, field, dt)\r\n  !---------------------------------------------------------------------------------------------------------------------------------\r\n  !< Integrate field with Adams-Bashforth class scheme.\r\n  !---------------------------------------------------------------------------------------------------------------------------------\r\n  class(adams_bashforth_integrator), intent(IN)    :: self  !< Actual AB integrator.\r\n  class(integrand),                  intent(INOUT) :: field !< Field to be integrated.\r\n  real(R_P),                         intent(in)    :: dt    !< Time step.\r\n  integer(I_P)                                     :: s     !< Steps counter.\r\n  !---------------------------------------------------------------------------------------------------------------------------------\r\n\r\n  !---------------------------------------------------------------------------------------------------------------------------------\r\n  do s=1, self%steps\r\n    field = field + field%t(n=s) * (dt * self%b(s))\r\n  enddo\r\n  call field%update_previous_steps\r\n  return\r\n  !---------------------------------------------------------------------------------------------------------------------------------\r\n  endsubroutine integrate\r\n```\r\n\r\nFor 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:\r\n```\r\nu(n+1) = u(n) + dt*[-0.5*F(n-1)+1.5*F(n)]\r\n```\r\nIn the code, this is calculated with the loop:\r\n```fortran\r\n  do s=1, self%steps\r\n    field = field + field%t(n=s) * (dt * self%b(s))\r\n  enddo\r\n```\r\nwhich is equivalent to:\r\n```fortran\r\nfield = field + field%t(n=1) * (dt * self%b(1))\r\nfield = field + field%t(n=2) * (dt * self%b(2))\r\n```\r\nThe 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.\r\n\r\nPlease 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:\r\n```fortran\r\n  do s=1, self%steps\r\n    field(self%steps) = field(self%steps) + field%t(n=s) * (dt * self%b(s))\r\n  enddo\r\n```\r\nbut 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?"},"comment":{"url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/comments/138074351","html_url":"https://github.com/Fortran-FOSS-Programmers/FOODiE/issues/18#issuecomment-138074351","issue_url":"https://api.github.com/repos/Fortran-FOSS-Programmers/FOODiE/issues/18","id":138074351,"user":{"login":"szaghi","id":1661942,"avatar_url":"https://avatars.githubusercontent.com/u/1661942?v=3","gravatar_id":"","url":"https://api.github.com/users/szaghi","html_url":"https://github.com/szaghi","followers_url":"https://api.github.com/users/szaghi/followers","following_url":"https://api.github.com/users/szaghi/following{/other_user}","gists_url":"https://api.github.com/users/szaghi/gists{/gist_id}","starred_url":"https://api.github.com/users/szaghi/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/szaghi/subscriptions","organizations_url":"https://api.github.com/users/szaghi/orgs","repos_url":"https://api.github.com/users/szaghi/repos","events_url":"https://api.github.com/users/szaghi/events{/privacy}","received_events_url":"https://api.github.com/users/szaghi/received_events","type":"User","site_admin":false},"created_at":"2015-09-06T11:39:29Z","updated_at":"2015-09-06T11:39:29Z","body":"@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."}}