Time splitting issue

[charlotte_M_]

 Hello Everyone,

I'm trying to solve a problem use time splitting method.  For this i have to solve linear systems  of equation using tridag (fortran ) routine.

I have two problems

The first one is that the coefficients  use in tridag are 2 D array while this routine use 1 D array.

What i did to overcome this is the following (just a example):

  do  k=1,NT         ! Main loop over the time step

!Compute  the solution of the first equation 

do i=1, N1

do j= 1, N2

   A(j) = x(i)*y(j)

  B(j)= x(i)*y(j)

  C(j)=x(i)*y(j)

 End do 

 Call tridag( A,B,C,F0, F1, N)     !F0  is a two dimensional function , F1, a one dimensional function

 End do 

!Once done i will use F1 as the initial condition for the second equation

 Do j =1, N2

 Do i=1,N1

  E(i)= xx(i)*y(jj)

 F(i) =xx(i)*yy(j)

 G(i)= xx(i)*yy(i)

End do 

  Call tridag(E,F,G,F1,F2,n1)

end do 

!Next step Write result in the file, for different time

!Finally   i update my initial condition

   DO i=1,N1

 DO j=1,N2

   F0(i,j) =F2(i,j)

  end do

end do

END DO         !ENd of the time loop

****** The second  problem is that , i have the same result for different  time that i print the solution.

 I am doing something wrong?

 sorry for the long text

  Thank you

[andrew_4619]

do  k=1,NT         ! Main loop over the time step
   !Compute  the solution of the first equation 
   do i=1, N1
      do j= 1, N2
         A(j) = x(i)*y(j)
         B(j)= x(i)*y(j)
         C(j)=x(i)*y(j)
      End do 
      Call tridag( A,B,C,F0, F1, N) 
      !F0  is a two dimensional function , F1, a one dimensional function
   End do 

   !Once done i will use F1 as the initial condition for the second equation
   Do j =1, N2
      Do i=1,N1
         E(i)= xx(i)*y(jj)
         F(i) =xx(i)*yy(j)
         G(i)= xx(i)*yy(i)
      End do 
      Call tridag(E,F,G,F1,F2,n1)
   end do 

   !Next step Write result in the file, for different time
   !Finally   i update my initial condition
   DO i=1,N1
      DO j=1,N2
         F0(i,j) =F2(i,j)
      end do
   end do
END DO         !ENd of the time loop

! if you put the code in a Fortran window ( the {code} icon) from the tool bar  it is more readable!

:-)

[andrew_4619]

I am slightly confused there seem to be several sources (variants) of TRIDAG from a quick google search but all the ones I looked at have 1D array inputs, perhaps you could ID the source of TRIDAG you are using or post the first few lines with the definitions of the items passed in/out.

[charlotte_M_]

Hello Andrew_4619,

Thank you for the answer and the fortran code

This is the routine.

SUBROUTINE tridag(a,b,c,r,u,n) 
PARAMETER (nmax=100) 
REAL gam(nmax),a(n),b(n),c(n),u(n),r(n) 
INTEGER j,n 
if (b(1)==0.) pause 'b(1)=0 in tridag' 
bet=b(1) 
u(1)=r(1)/bet 
do j=2,n 
    gam(j)=c(j-1)/bet 
    bet=b(j)-a(j)*gam(j) 
    if (bet==0.) pause 'bet=0 in tridag' 
    u(j)=(r(j)-a(j)*u(j-1))/bet 
end do 
do j=n-1,1,-1 
    u(j)=u(j)-gam(j+1)*u(j+1) 
end do 
END SUBROUTINE tridag 

[andrew_4619]

OK I do not know what the program is trying to do but you are passing a 2D array to a 1D array so in the first call to TRIDAG you are using only the first N items of the array which is in column order ie F0(1,1), F0(2,1),....F(N,1) or if the size first index of F0  is <N you will wrap into column 2.

I am guessing you want to pass a slice of F0 e.g. F0(:,i) (passes a 1D array that is the column number i) or maybe F0(j,:) for row number J. Note that the second case is not effecient as the slice represents non-adjacent memory locations so the compiler will create a temporary array to hold this data.....

[charlotte_M_]

 For the first  equation What i want is for a given value of i, solve the equation for all the value of j.

The process has to be repeat for all the value of i ( loop over i).

The initial distribution is F0(i,j).

[jimdempseyatthecove]

Your comment in your code posting state that F0 and F1 are functions, not arrays. The subroutine tridag specifies that dummy arguments r and u are REAL arrays.

If your comment is wrong, or should I say partly wrong, in that F0 and F1 are arrays, with F0 being a 2D array, then Andrew's comment about what are you intending to pass into tridag as r? IOW, is it a row in F0 or a column in F0?

Your code having the outer DO k=1,NT would seem to indicate that F0 may contain different time steps as either by row or by column. We do not have information as to how F0 is arranged.

Jim Dempsey

[charlotte_M_]

Sorry my comment is wrong. F0(i,j) and F1 are both real arrays. F0  is  2 D array.

F0 is r in tridag.

Let's say F0(i,j) = exp(- i**2)*exp(-j).

[charlotte_M_]

 Searching online, i think my  problem is how do i pass this 2 D array to tridag routine? 

And once i have the result from tridag, how can i have back the array.

Thank you so much

[charlotte_M_]

In all the case i try to pass in tridag the index of the second loop .. It means a column in F0.

However  for the second equation,  the column becomes the row and the row the column.

 Little mistake, the first call is CALL tridag( A,B,C,F0,F1,N2) 

F0 is allocate array of dimension (N1, N2)..

 So yes i want to pass a slice of F0

For the first equation i want to pass F0(i,:) and the second F0(:,j)

Any suggestion?

[jimdempseyatthecove · Accepted solution]

Pass F0(i,:) and F0(:,j).

Note, in the case of F0(i,:), the compiler will generate a temporary (because the slice has a stride of N1, size of N2).
For the F0(:,J), the compiler should use the data in place (because the slice has a stride of 1, size of N1).

Jim Dempsey

[jimdempseyatthecove]

By the way. Consider using:

SUBROUTINE tridag(a,b,c,r,u,n) 
REAL gam(n),a(n),b(n),c(n),u(n),r(n) ! have gam assume the size of the other arrays (n)

This will prevent a buffer overflow should n (N1 or N2 from caller) exceed 100.

If n is very large (say over 1,000,000) consider making gam allocatable:

SUBROUTINE tridag(a,b,c,r,u,n) 
REAL a(n),b(n),c(n),u(n),r(n)
REAL, allocatable :: gam(:)
allocate(gam(n))
...
deallocate(gam) ! optional on newer specifications

Jim Dempsey

[jimdempseyatthecove]

Once you get the code working, and if the entire process takes longer to execute than you would like, readers here can show you how to get better performance.

If N1 and N2 get reasonably large, then parallelization and/or vectorizaton improvements may be available.
And, even if parallelization is not warranted, two threads can have one thread computing time steps and the other thread performing the file writes. IOW runtime goes from:

(computation time) + (File Write time)

to:

max((computation time), (File Write time)) + (a tad for threading overhead)

Jim Dempsey

[charlotte_M_]

 Thank you so much  Jim Dempsey for your answer.

For the first equation i passed F0(i,:).  Everything was ok.

The solution given by TRIDAG is F1.

Now i want to use F1 as F0 for the next call (second equation).  F1 becomes my initial array (IOW r)

But F1 is not a 2D array, how can  i pass F1(j,:)

I made the modification for the gam ..Thank you so much

[jimdempseyatthecove]

In looking at your original sketch code (#1), you are not using the F1 results after the call to TRIDAG and its end do. IOW you are overwriting the results. Therefor, I assume you intended to have made F1 2D as well, then process the row/col of F1 array in the second nested loop.

Jim Dempsey

[charlotte_M_]

jimdempseyatthecove wrote:

In looking at your original sketch code (#1), you are not using the F1 results after the call to TRIDAG and its end do. IOW you are overwriting the results. Therefor, I assume you intended to have made F1 2D as well, then process the row/col of F1 array in the second nested loop.

Jim Dempsey

Yes  exactly, after making F1 2D array , i pass F1(:,j) for the second call.

[andrew_4619]

well in the first case you have Call tridag( A,B,C,F0(i,:), F1, N2) 

but your second points confuses me, is it that at this first call you need the values of F1 for each i to use later?

If so make F1 2D to the same dims as F0 and then have  as Call tridag( A,B,C,F0(i,:), F1(i,:), N2)

then later you would have: Call tridag(E,F,G,F1(:,J),F2,n1)

I am guessing as I have made no attempt to understand teh nature of your calculation..... 

[charlotte_M_]

  F1 is the solution of TRIDAG, no need to know it.

When i call TRIDAG  for the first time , the unknown  is F1, so TRIDAG will give me F1

[charlotte_M_]

  This is what i did :

First call :  Call TRIDAG  ( A,B,C,F0(i,:), F1(i,:) , N2)

Then for the second 

       call TRIDAG (E,F,G, F1(:,j), F2(:,j) , N1)

[charlotte_M_]

Thank you so much Jim Dempsey and andrew_4619.

Everything is working now. And i have a better understanding.

The next step will be parallelization..

So see you soon

[andrew_4619]

! is this the real code or just example?
      do j= 1, N2
         A(j) = x(i)*y(j)
         B(j)= x(i)*y(j)
         C(j)=x(i)*y(j)
      End do 
!assuming A,B,C any Y are all declared the same size
!the loop can be eliminated and the code is:
A=x(i)*Y
B=A
C=A
!the compiler is likely to vectorise these whole array ops (1 stage parallelisation)
!
!if they are different length you would say:
A(1:N2) = x(i)*Y(1:N2)
A(1:N2) =B(1:N2)
C(1:N2) =C(1:N2)
! that might not be quite so efficient (depends on N2 which is possibly unknown
! at compile time) as the optimiser are to make some decisions but is still 
! likely to be faster.