Solving Diffusion Equation with Implicit Method: Problem with DGESV

  • Thread starter Thread starter debsankar
  • Start date Start date
AI Thread Summary
The discussion focuses on solving the diffusion equation using the implicit method with the DGESV subroutine from LAPACK, highlighting issues with periodic boundary conditions. The original poster's code produces unexpected results, specifically a flat line at zero in the output, indicating a potential problem with the boundary conditions. Suggestions include avoiding persistent tabs in code, providing specific details about the issues encountered, and using dummy variables when passing arrays to LAPACK subroutines to prevent unintended modifications. The poster acknowledges the confusion and clarifies that the core problem lies in the implementation of periodic boundary conditions. Overall, the conversation emphasizes the importance of clear coding practices and precise problem descriptions in debugging numerical methods.
debsankar
Messages
11
Reaction score
0
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.
 

Attachments

Technology news on Phys.org
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.
 
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.
 
gsal said:
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

Code:
!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
Code:
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
 
Last edited:
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.
 
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.
----------------------------------------------------------------------------
 
Dear Peeps I have posted a few questions about programing on this sectio of the PF forum. I want to ask you veterans how you folks learn program in assembly and about computer architecture for the x86 family. In addition to finish learning C, I am also reading the book From bits to Gates to C and Beyond. In the book, it uses the mini LC3 assembly language. I also have books on assembly programming and computer architecture. The few famous ones i have are Computer Organization and...
I had a Microsoft Technical interview this past Friday, the question I was asked was this : How do you find the middle value for a dataset that is too big to fit in RAM? I was not able to figure this out during the interview, but I have been look in this all weekend and I read something online that said it can be done at O(N) using something called the counting sort histogram algorithm ( I did not learn that in my advanced data structures and algorithms class). I have watched some youtube...

Similar threads

Replies
22
Views
5K
Replies
1
Views
3K
Replies
8
Views
2K
Replies
6
Views
2K
Replies
4
Views
7K
Replies
2
Views
1K
Back
Top