Community
cancel
Showing results for
Did you mean:
Highlighted
Beginner
46 Views

## Time splitting issue

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

Tags (1)

Accepted Solutions
Highlighted
Black Belt
42 Views

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

22 Replies
Highlighted
Valued Contributor III
42 Views
```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!

:-)

Highlighted
Valued Contributor III
42 Views

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.

Highlighted
Beginner
42 Views

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 ```

Highlighted
Valued Contributor III
42 Views

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.....

Highlighted
Beginner
42 Views

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).

Highlighted
Black Belt
42 Views

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

Highlighted
Beginner
42 Views

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).

Highlighted
Beginner
42 Views

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

Highlighted
Beginner
42 Views

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?

Highlighted
Black Belt
43 Views

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

Highlighted
Black Belt
42 Views

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

Highlighted
Black Belt
42 Views

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:

Jim Dempsey

Highlighted
Beginner
42 Views

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

Highlighted
Black Belt
42 Views

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

Highlighted
Beginner
42 Views

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.

Highlighted
Valued Contributor III
42 Views

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..... `

Highlighted
Beginner
42 Views

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

Highlighted
Beginner
42 Views

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)

Highlighted
Beginner
42 Views

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