Memory management in time-domain simulation

Question

I develop a parallel time-domain DAE simulation software using OpenMP and Fortran 2003. The main algorithm is:

while t<time_horizon do
   X_1=X
   while not converged do
      !$OMP PARALLEL DO
         Solve sub-domain DAE using a portion on X and X_1 to get DeltaX
         Correct X=X+DeltaX
      !$OMP END PARALLEL DO
   end
   t=t+h
end

X is the state variable vector and X_1 is the the vector at the previous time-step used for the discretization algorithm (Trapezoidal method).

The sequential implementation profiling shows that the memory copy takes 15% of the computational time, a significant amount...

I tried:

!$omp parallel workshare
   X_1=X
!$omp end parallel worshare

And:

!$omp parallel
   i=omp_get_thread_num()
   X_1(i*length/num_threads+1:(i+1)*length/num_threads)= X(i*length/num_threads+1:(i+1)*length/num_threads)
!$omp end parallel

This helped a little bit (dropped to 13%) but not much... Any thoughts? Any other way to keep the history? I keep the X of all subdomains as one vector as it's easier to handle (apply infinite norm, copy, etc.). Should I break it into one for each subdomain?

Thanks in advance!

Answer 1

You do not need to copy the state vector. Assuming you are storing two copies X1 and X2, you can have the parallel region first work on X1, and then work on X2 on the next iteration. In Fortran, one easy way to do so is to move the parallel region into a subroutine, and to call it from within the main loop:

while t<time_horizon do
   call evolve(X1, X2)
   t=t+h
   if (t>=time_horizon) then
      "swap x1 and x2"
      exit
   end if
   call evolve(X2, X1)
   t=t+h
end while