Solving Diffusion Equation with Implicit Method: Problem with DGESV

[debsankar]

I am trying to solve Diffusion eq with implicit method using DGESV subroutine from lapack...
But my periodic boundary condition is not working..
Here is my code: attached...

please help.

 

[gsal]

General suggestions:

When typing program source code, never let tabs persist; it is generally not well liked among programmers.

Tabs are nice for quick spacing, but I suggest you use one of those editors that offer the feature of emulating tabs with a given number of blank spaces, instead of the actual tab character. It is just that tabs are rendered differently by different editors or people set the lengths different according to their own taste, etc...there is just no consistency, so stay away from persistent tabs.

Please use indentation.

Your program uses two external files, you did not provide them.

Be more specific about your request...you just said "please help" but mentioned nothing about any specific issues you are having: does it compile? does it run? does it give any results? how are they different from what you expect? etc, etc.

 

[Mark44]

And to add to what gsal said, attaching you code as a zip file is not the best way to do things. For such a small file, just paste the text right in the text window. That way we can see your code without have to go to the effort of opening a zip file, and opening the source code file in an editor.

 

[Solkar]

Your program uses two external files, you did not provide them.

Those files are accessed write-only,

But I agree with the rest of gsal's and Mark44's critics.
There are even many more coding issues, but let's skip those and only look at the core

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!	CODE TO SOLVE DIFFUSION EQUATION USING SHERMAN-MORRISON METHOD
!	Df= diffusion co-eff, delt and delx are grid spacing of time 
!	and space in discretized equations.

Sherman-Morrison is good for calculating inverses cheaply while iterating, but with

do k=1,T
    do i=1,L
        p(i)=x(i)
    ! write(*,*) p(i)
    enddo
	
    call DGESV( L, nhrs, A, L, pivot, p, L, OK )

    do i=1,L
        x(i)=p(i)
    ! write(*,*) p(i)
    enddo
enddo

write(*,*) OK

(indentation fixed)

you T times solve straight ahead
A p_{k+1} = p_{k}
for p_{k+1},

which is, btw, expensive compared to calling DGETRF once and subsequently only calling DGETRS for solving; lapack will, without optimizations, try to LUP-decompose an already LU-decomposed matrix your way.

---

What is also amazing is, that you check your "OK" (conventionally called "info", btw) only once at the end, but never in the loop.

---

Anyway - when I compile and run that I get two files with a bunch of fp numbers, and a single DGESV info of "0" displayed, which, as usual, means "all is well".

How are we supposed to see without further explanation that those fp numbers are flawed and esp. without having seen your ansatz in an explicit form? Solkar

 

[debsankar]

1. I am really sorry to post it so badly... It was my first time and I was in hurry... Thank you all for the suggestions...

2. I assumed people will get the problem... it happens when you are working with a problem long enough...
but it does not work like that... sorry again !
3. Now about the problem- I misguided you all by the comment : 'CODE TO SOLVE DIFFUSION EQUATION
USING SHERMAN-MORRISON METHOD'. It is just solving A*x(n+1)=x(n) in a loop, where x(n) means in
time t0+n*dt value of x... so if you keep on solving and after each solve you put x(n)=x(n+1) you sould be
able to solve time evolution of the quantity x. (A=m x m matrix, x= m x 1 array). Now you can write the 1D
diffusion eq in this format if you consider the implicit finite difference method. That is what I was trying to
solve in addition I had periodic boundary condition.
4. Now what was the problem: the code produces two .txt files 'initial.txt' and 'diffusion.txt' .
'initial.txt' : store the initial configuration for diffusion problem.
'diffusion.txt' : store the end result.
Now if you run the code and plot both files (set the 'time' parameter at Ln20,Col9 to be 0.5 or 1.0, i.e. long
time) you will see the initial config- a rectangle of size 1.0 x 0.2 - area under the curve = 0.2.
As my system size is (x0=0.0, xL=1.0) 1.0, you should expect to get the long time diffusion curve (with PBC)
to be almost a horizontal line of height 0.2 , area under the curve now = 0.2*(xL-x0)=0.2, as there is no loss in the system.
4. I was getting a FLAT LINE AT ZERO! so I guessed it has to be some problem with the periodic BC. That was
my problem with the code.

Hope I have stated my problem with enough details.

 

[debsankar]

Just passing by... so thought I will explain what happened back then >>
When using LAPACK subroutines in your fortran code while passing variable (array or matrix) to LAPACK its better to use a dummy variable to pass, because LAPACK changes your input arrays too. Then when you try to use the input array (after calling LAPACK) and the output from that subroutine, that gives you wrong results. Thats what happened in the above case...
Example >>
---------------------------------------------------------------------------
your input array > A
set dum_A = A
call LAPACK sub( dum_A ,..., OUTPUT)
for later calculations use the actual array A, dum_A has been changed by the subroutine.
----------------------------------------------------------------------------