# Swapping arrays during time-stepping **URL:** \ **Category:** Poll\ **Created:** [February 6, 2024, 10:29pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354 "2024-02-06T22:29:36Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [February 6, 2024, 10:29pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/1 "2024-02-06T22:29:36Z")
For time-stepping large PDE problems, it’s common to read from one array (the previous time step), and write to a second array (the result at the new time step). Afterward the arrays are swapped. This is also known as the “flip-flop” approach, or ping-ponging. The separate arrays also make the update step easier to parallelise (no race conditions). What is your favourite implementation approach? (Scroll down for examples) _Poll ([view on site](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/1))_ ### Separate arrays with parity flag ```fortran real, allocatable :: f1(:), f2(:) integer :: flipflop ! ... initialize f flipflop = 1 do step = 1, nsteps select case(flipflop) case(1) call sweep(fin=f1,fout=f2) case(2) call sweep(fin=f2,fout=f1) end select flipflop = 3 - flipflop end do ``` If zero-based indexing is used, the parity flag can be inverted using ```fortran ! 0-based indexing flipflop = 1 - flipflop ``` ### Array with extra dimension This one relies on array slicing or dummy argument association. Integer flags are used to denote the old and new arrays. ```auto real, allocatable :: f(:,:) integer :: tmp, old, new allocate(f(n,2)) ! ... initialize f old = 1 new = 2 do step = 1, nsteps call sweep(fin=f(:,kold),fout=f(:,knew)) tmp = old old = new new = tmp end do ``` Alternatively, the integer flags can be calculated before the update like this: ```fortran old = mod(step+1,2) + 1 new = mod(step+2,2) + 1 ``` or after the update, like this: ```fortran old = 3 - old new = 3 - new ``` ### Pointer swap ```auto real, allocatable, target :: f1(:), f2(:) real, pointer :: fold(:), fnew(:), ftmp(:) fold => f1 fnew => f2 do step = 1, nsteps call sweep(fin=fold,fout=fnew) ftmp => fold fold => fnew fnew => ftmp end do ``` ### Overwrite old array ```fortran real, allocatable :: fold(:), fnew(:) do step = 1, nsteps call sweep(fin=fold,fout=fnew) fold = fnew end do ``` ---
**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\ **Post date:** [February 6, 2024, 11:23pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/2 "2024-02-06T23:23:56Z")
My first pass is somewhere between overwrite and swap, as I like to be able to write something like this at a high level. ```auto state = state + d_dt(state) * dt ``` For very large problems maybe this isn’t efficient, but I bet the compilers these days are pretty good at playing some tricks to optimize a statement like that. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [February 6, 2024, 11:28pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/3 "2024-02-06T23:28:12Z")
Just for completeness, I’ve also seen this done in the following way: ```auto do step = 1, nsteps/2 call sweep(fin=f(:,kold),fout=f(:,knew)) call sweep(fin=f(:,knew),fout=f(:,kold)) end do ``` ---
**Author:** ![amasaki203](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/amasaki203/32/2730_2.png) [@amasaki203](https://fortran-lang.discourse.group/u/amasaki203)\ **Post date:** [February 7, 2024, 12:02am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/4 "2024-02-07T00:02:38Z")
For most problems, I prefer the pointer swap approach rather than managing indexes or increasing if and case branches. When using this method, I always add a `contiguous` attribute to the declaration of the pointer variables to expect this will benefit from compiler optimizations. Thus, like that: ```fortran real, pointer, contiguous :: fold(:), fnew(:), ftmp(:) ``` For small-scale problems such as one-dimensional, I use the overwrite old array approach. ---
**Author:** ![oscardssmith](https://avatars.discourse-cdn.com/v4/letter/o/b9e5f3/32.png) [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)\ **Post date:** [February 7, 2024, 1:07am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/5 "2024-02-07T01:07:28Z")
The main other approach that is often better (assuming the arrays are small enough) is to use a higher order integrator that uses more than 2 state arrays. Once you are using a higher order integrator, you need to store a bunch of stages anyway, and the output is a linear combination of these stages, so it’s a somewhat moot point since you no longer have 2 arrays to swap. ---
**Author:** ![IanMartinAjzenszmidt](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ianmartinajzenszmidt/32/1992_2.png) [@IanMartinAjzenszmidt](https://fortran-lang.discourse.group/u/IanMartinAjzenszmidt)\ **Post date:** [February 7, 2024, 4:45am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/6 "2024-02-07T04:45:49Z")
The “flip-flop” or “ping-ponging” approach to time-stepping in large PDE (Partial Differential Equation) problems is indeed a common and effective strategy, especially in parallel computing environments. My preferred implementation approach would depend on a few factors, such as the specific PDE being solved, the computing environment (e.g., CPU vs. GPU), and the programming language being used. However, a general approach I favor involves a combination of clarity, efficiency, and adaptability to parallel computing. Here’s a broad outline: 1. **Use of High-Level Abstractions (when applicable):** In environments where high-level abstractions are available and performance is not severely impacted (like Python with NumPy or Julia), using these can greatly simplify code and make it more readable. This is especially true for operations that are inherently parallel, like array operations. 2. **Explicit Buffer Swapping:** Maintain two distinct buffers or arrays, one for the current time step and one for the next. After each time step, swap these buffers instead of copying data between them. This can be as simple as swapping pointers or references to the arrays, which is much more efficient than copying array contents. 3. **Parallelism-Friendly Code:** Write the update step in a way that is amenable to parallel execution. This usually means avoiding dependencies between different parts of the array during the update. For instance, if using a language like C++ with OpenMP or Python with Numba, you would ensure that the loop iterations are independent so they can be safely parallelized. 4. **Memory Considerations:** Especially in GPU computing, but also relevant for CPU, consider how memory access patterns affect performance. Strive for coalesced memory accesses in GPUs and cache-friendly patterns in CPUs. 5. **Error Handling and Validation:** Ensure that your implementation includes robust error handling and validation checks. This is crucial for debugging and for ensuring that parallelization does not introduce subtle bugs. 6. **Profiling and Optimization:** Use profiling tools to understand where bottlenecks lie and optimize those parts of the code. This might involve different strategies for different architectures (e.g., SIMD optimizations for CPUs, warp-shuffle operations for GPUs). 7. **Scalability and Adaptability:** Design the implementation to be scalable to different problem sizes and adaptable to different PDEs or boundary conditions. This often means abstracting the core computational routines from the specifics of the PDE being solved. 8. **Documentation and Comments:** Maintain clear documentation and comments, especially for parts of the code that handle the complexities of parallelization and memory management. In summary, my favorite approach balances performance with readability and maintainability, leveraging high-level abstractions when possible but also diving into lower-level optimizations where they are most beneficial. This approach should be adaptable to various types of PDEs and computing environments. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [February 7, 2024, 7:26am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/7 "2024-02-07T07:26:58Z")
I use array overwriting which for time stepping in large non-linear FEM problems is just transparent compared to the time required for solving the problem. With a two time steps approach it would globally speaking look like: 1. M \* dX = Residual( loads\_{t+dt} , X\_t, X\_{t-dt} ) 2. X\_{t+dt} = coef\_1 \* X\_t + coef\_2 \* X\_{t-dt} + coef\_3 \* dX 3. X\_{t-dt} = X\_t 4. X\_t = X\_{t+dt} Where in (1) the matrix M can be linear or non-linear so Newton-like iterations might be needed within. When parallelizing in a SPMD framework, with certain d.o.f. being shared across processes, one “just” needs to check that dX is synchronized at the end of (1), (2)-to-(4) can be done without any regards to communication. ---
**Author:** ![irukoa](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/irukoa/32/3306_2.png) [@irukoa](https://fortran-lang.discourse.group/u/irukoa)\ **Post date:** [February 7, 2024, 12:41pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/8 "2024-02-07T12:41:59Z")
To add: in problems involving unitary dynamics it is also common to employ exponential integrators. In these cases I go with something like ```fortran dt = (tend - tstart)/real(nt - 1) do it = 1, nt t = tstart + dt*(it - 1) tev = matmul(tev, matrix_exp(-cmplx_i*H(t)*dt)) enddo ``` This approach also requires no swapping. ---
**Author:** ![VladimirF](https://avatars.discourse-cdn.com/v4/letter/v/e9a140/32.png) [@VladimirF](https://fortran-lang.discourse.group/u/VladimirF)\ **Post date:** [February 8, 2024, 3:45pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/9 "2024-02-08T15:45:16Z")
In my main code I swap the pointers, but with `move_alloc()`, I do not use pointers unless I _really_ have to. ---
**Author:** ![VladimirF](https://avatars.discourse-cdn.com/v4/letter/v/e9a140/32.png) [@VladimirF](https://fortran-lang.discourse.group/u/VladimirF)\ **Post date:** [February 8, 2024, 3:50pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/10 "2024-02-08T15:50:33Z")
There are low-storage higher order Runge-Kutta solvers available. One has one or maybe two extra array per variable, but one has an array with a result of the stage and the array with the original value. And you need to somehow swap these. It doesn’t matter if it i 1st order Euler, or 5th order RK, the problem still remains the same. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [February 8, 2024, 4:21pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/11 "2024-02-08T16:21:31Z")
So like this? ```auto real, allocatable, dimension(:) :: ftmp, fold, fnew call move_alloc(from=fold, to=ftmp) call move_alloc(from=fnew, to=fold) call move_alloc(from=ftmp, to=fnew) ``` Seems like the language could introduce a generic `swap_alloc` intrinsic function. ---
**Author:** ![VladimirF](https://avatars.discourse-cdn.com/v4/letter/v/e9a140/32.png) [@VladimirF](https://fortran-lang.discourse.group/u/VladimirF)\ **Post date:** [February 8, 2024, 4:40pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/12 "2024-02-08T16:40:30Z")
Like that, yes, I have a subroutine for that. ---
**Author:** ![nncarlson](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/nncarlson/32/1175_2.png) [@nncarlson](https://fortran-lang.discourse.group/u/nncarlson)\ **Post date:** [February 8, 2024, 10:14pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/13 "2024-02-08T22:14:42Z")
Yes, this. Fast like swapping pointers, but can result in faster code when operating on the arrays (depending on context) because they won’t have the target attribute and the optimizer shouldn’t have to worry about possible aliasing. That’s been my experience, fwiw. ---
**Author:** ![manuel.arenaz](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/manuel.arenaz/32/3217_2.png) [@manuel.arenaz](https://fortran-lang.discourse.group/u/manuel.arenaz)\ **Post date:** [March 13, 2024, 6:52pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/14 "2024-03-13T18:52:57Z")
Interesting point… I have seen several of them in real codes ands others are new to me. With simplicity and maintainability in mind, assuming all of them provide a reasonable performance, I would vote for Overwriting arrays. Nice thread @ivanpribec ! ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [March 14, 2024, 7:20am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/15 "2024-03-14T07:20:53Z")
Thanks @manuel.arenaz. I also find it useful to have these kind of application-specific patterns collected in one place. * * * A few more remarks on this topic: - The method suggested by @RonShepard is essentially loop unrolling (see [PWR049](https://github.com/codee-com/open-catalog/tree/main/Checks/PWR049#example-3-loop-unrolling) in the Codee Open Catalog). Personally I’d prefer to modify the loop counter like this: - The various approaches can be combined. For example use of parity and pointers: - The approach with `move_alloc` works perfectly fine also with arrays mapped to OpenMP target devices and also with OpenACC and CUDA managed memory. - Swapping pointer arrays or swapping allocatable arrays using `move_alloc` is very similar in terms of the [amount of instructions required](https://godbolt.org/z/qG6TM653T). (Not that it really matters, because the bulk of work will be in the other subroutines within the loop body.) In other programming languages swapping looks like this: ```auto // C void swap(double **a, double** b) // swaps the pointers! { double *tmp = *a; *a = *b; *b = tmp; } // C++ #include std::swap(a, b); ``` ```python # Python (and also Julia) a, b = b, a ``` ```txt % MATLAB % https://www.mathworks.com/matlabcentral/answers/362680-how-to-swap-values-of-two-variables % Option 1 [b, a] = deal(a,b); % Option 2 [b, a] = swap(a,b) function [b, a] = swap(a, b) % This function has no body! ``` * * * Edit, 05/04/2024: I added the link which helped me realize that doing two time steps is equivalent to loop unrolling. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [March 14, 2024, 3:39pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/16 "2024-03-14T15:39:20Z")
> [@ivanpribec](#): > > ```auto > if (modulo(step,2) == 0) then > ,,, > > call sweep(pold, pnew) > > ``` This approach also works when you have multiple array swaps, e.g. when you are using a circular buffer workspace. It can be done with the last index of a multidimensional array, with array offsets, with pointers, or with allocatable arrays and `move_alloc()`. Some care must be taken when using array offsets or pointers to hide the possible memory aliases. The compiler can produce the optimal code when it knows that aliases cannot occur. One way to do this is with a subroutine call, as above. Once inside the subroutine, the compiler knows that its dummy arguments cannot be aliased. I’m curious if the same semantics occurs when the subroutine call is replaced with an ASSOCIATE block? Can the compiler know that the associate block renames are not aliased? ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [September 29, 2024, 4:03pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/17 "2024-09-29T16:03:06Z")
Here is a more object-oriented approach adapted from a pattern encountered with swapping kernel buffers in OpenCL ([How to effectively swap OpenCL memory buffers? - Stack Overflow](https://stackoverflow.com/a/21984042/4283055)): ```auto program timesweep implicit none interface subroutine sweep(n,fin,fout) integer, value :: n real, intent(in) :: fin(n) real, intent(out) :: fout(n) end subroutine end interface type :: kernel_t integer :: n real, pointer, contiguous :: fin(:), fout(:) procedure(sweep), pointer, nopass :: kernel => sweep end type integer :: n, nsteps, step real, allocatable, target :: f1(:), f2(:) type(kernel_t) :: odd, even nsteps = 1000 n = 100 allocate(f1(n),f2(n)) ! ... Initialize f1, f2 ... even = kernel_t(n,f1,f2) odd = kernel_t(n,f2,f1) do step = 0, nsteps-1 if (modulo(step,2) == 0) then call exec(even) else call exec(odd) end if ! in F2023, also ! call exec(mod(step,2) == 0 ? even : odd) end do contains subroutine exec(kern) type(kernel_t), intent(in) :: kern call kern%kernel(n=kern%n,fin=kern%fin,fout=kern%fout) end subroutine end program ``` The nice thing about this approach is it is easy to change the sweep method - you just need to supply a different procedure pointer. This decision could also happen at runtime. Unless the kernel would provide some other methods/behaviors, I don’t think this pattern is really worth the extra effort. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [December 18, 2025, 7:30am UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/18 "2025-12-18T07:30:42Z")
In the context of coarray programming `move_alloc` implies a synchronization of images. I haven’t considered other options yet. ---
**Author:** ![Machalot](https://avatars.discourse-cdn.com/v4/letter/m/ebca7d/32.png) [@Machalot](https://fortran-lang.discourse.group/u/Machalot)\ **Post date:** [December 19, 2025, 6:01pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/19 "2025-12-19T18:01:00Z")
> [@ivanpribec](#): > > ### Separate arrays with parity flag > > ```auto > real, allocatable :: f1(:), f2(:) > integer :: flipflop > > ! ... initialize f > > flipflop = 1 > > do step = 1, nsteps > > select case(flipflop) > case(1) > call sweep(fin=f1,fout=f2) > case(2) > call sweep(fin=f2,fout=f1) > end select > flipflop = 3 - flipflop > > end do > > ``` > > If zero-based indexing is used, the parity flag can be inverted using > > ```auto > ! 0-based indexing > flipflop = 1 - flipflop > > ``` Is there a reason to update `flipflop` using arithmetic rather than just setting its value explicitly in each `case`, like this? ```Fortran select case(flipflop) case(1) call sweep(fin=f1,fout=f2) flipflop = 2 case(2) call sweep(fin=f2,fout=f1) flipflop = 1 end select ``` ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [December 19, 2025, 7:37pm UTC](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354/20 "2025-12-19T19:37:36Z")
> [@ivanpribec](#): > > In the context of coarray programming `move_alloc` implies a synchronization of images. Can you elaborate on this a little? It seems like `move_alloc()` would be an operation local to an image. At what point does the synchronization between images occur within `move_alloc()`? [Next page](https://fortran-lang.discourse.group/t/swapping-arrays-during-time-stepping/7354.md?page=2)