Smart questions
Smart people
 Find A ForumFind An Expert
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Remember Me

Are you a
Computer / IT professional?
Join Tek-Tips now!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Just copy and paste the

#### Feedback

"...It is good to know that there are groups such as this willing to share knowledge in this money driven economy..."

#### Geography

Where in the world do Tek-Tips members come from?

# Some help with a Fortran source code

 Forum Search FAQs Links Jobs Whitepapers MVPs
 bigduke88 (TechnicalUser) 3 May 12 6:27
 Hello everyoneI would like to ask you for some help with a Fortran source code.I have to make a program, which solves numerically an ordinary differential equation and after that an integral.The equation, is going to be solved using the finite difference approximation.This is equation "developed" in finite differences:Ui+1 = 0.998*h - 0.7366*(Ui)^2*Cd*h + UiThis is an equation, which describes the velocity of a moving spherical particle."Cd" is a drag coefficient and is function of U:And this is the equation, which describes the drag coefficient:Cd = (0.087108/U) + (0.622126/U^0.28)I'm quite new in any use of an algorythmic "language", but I have to solve the equation and then an intergral in Fortran.Could you help with some advice, how to do this?Thank you!!!
 gummibaer (Programmer) 3 May 12 7:01
 And you expect us to deliver the code ??Norbert  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 bigduke88 (TechnicalUser) 3 May 12 7:13
 No, I already have created a code, for some pde-s and an itegral, but if I include this equation in the "original" source code for the pde-s I think, that I'm going to mess things, because I'm working with a dimensionless time, but for the ODE, have to transfer from dimensionless to a "normal" one... I thought about a subroutine - to solve the ode with a different program and then "call" the result in the original one
 gummibaer (Programmer) 3 May 12 8:18
 Why don't you post what you have ?I do not have any idea what a PDE or ODE might be (partial differential equation, ordinary differential equation ??) and how this is introduced into your code and what you want to do with it in terms of coding.With some code we would have a baseline for discussion.Norbert  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 salgerman (Programmer) 3 May 12 12:35
 As you just realized, in order to program something, you really, REALLY need to know what you are doing.Contrary to popular believe, computer programs are not magic...other than artificial intelligence, neural networks and other stuff that I don't know...for the rest of us, computers suffer of artificial stupidity and we need to tell them exactly what to do.But don't let programming scare you...at the minimum, think of programming as a way of automating your hand calculations...so, how would you go about solving this problem by hand? do it, at least once; or, if it is tooo long, sketch the steps up, take a look at the steps and program them.Lastly, do not expect any of us be experts on the field that you are and so we may not be able to help you go from "dimensionless" time to a "normal" one...if we are experts on something is fortran and whatever other field we use fortran for.So...think and ask the question in a way we can understand...try to make it a fortran question
 fanta2 (TechnicalUser) 4 May 12 0:27
 Your equation seems to be the final product (explicit solution) after you have substituted the finite difference form into the PDE. If that is the case you have to just calculate Ui+1 from an initial condition of Ui.
 bigduke88 (TechnicalUser) 8 May 12 8:17
 Hello everyoneIt took a long time, but I had to do some other things, about my project. The project describes the process of diffusion of an acid gas, with a freely falling drop of an alkalyne solution. One of the basics in the transport phenomena.That's why I was talking about dimensionless times etc :)Here is source, that I've written about the equation describing the velocity of the falling drop in function of time.At the end, I've implemented a "Z", this is the distance passed by the drop.The code has some errors and now I'm trying to fix them. When I do it, will be posted here.Best regards!!!P.S.The source code      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      open(unit=3,file='du/dt.txt')      NIT=200      DT1=1D-8      DT=DT1      DTau=4.D0*DT      Dp=0.00271D0      Ro=998.2D0      Mu=0.001005D0      Diff=2.97271D-9      Tau=(4.D0*Diff*T)/(Dp**2.D0)      T=(Tau*(Dp**2.D0))/(Diff*4.D0)      imax=100      DO 30 i=1,imax30      U(i)=0.      T=0.      DO i=1,imax      T=T+DT      Tau=Tau+DTau      Re=(Dp*U(i)*Ro)/(Mu)      Cd=(24/Re)*(1+0.125*(Re**0.72))      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT      DO i=1,imax      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT      Z=Z+DT*(U(i+1)+U(i))100   continue      write(18,*) Z,T,Tau      STOP      END
 salgerman (Programmer) 8 May 12 9:37
Yeah, I can see a couple of errors...

First, you have not declared the array U().  When you do, you can simply declare it of fix length like U(100), if you are going to stick to imax=100.

Second, you have "100 continue", but the "100" label is missing from your "DO" statements.

I would recommend you stop using numeric labels, altogether, and start using the "DO..END DO" construct.  As it is, you have 2 do-loops ended using the same exact continue...not cool, anymore.

Also, in Fortran 90, you can deal with arrays a-la-matlab; in other words, you do no need to initialize U() one item at a time as in

#### CODE

DO 30 i=1,imax
30      U(i)=0.

#### CODE

U = 0.

Lastly, I don't see why you have to re-calculate Cd every time inside the loop if it is a constant.

my 2 cents.

 salgerman (Programmer) 8 May 12 9:40
 oops...sorry, I guess Cd is changing.Also, I am confused about your nested DO-loops and you re-calculating the entire U() array in the inner one...without know anything, something does not look right...
 gummibaer (Programmer) 8 May 12 11:41
 Well,tracking down bugs is sure a tedious business - but the best way to learn programming.But one (or two) hint(s) as to good programming practice:Never use implicit type declaration - unless you want to train your eyesight and attention by tracking down typos. It is more work at first - but I'd recommend a data dictionary anyway. Or are you really sure you still know what you thought today when you have to edit your code in six months ?Norbert.  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 bigduke88 (TechnicalUser) 8 May 12 13:12
 Really thank you all for the responses, they are really usefull!!!I've done some corrections on the source, now it works, but when I open the folder, where should be stored the results from the calculations, there are only F.files and application... I thought, that when I enter the statement"open(unit=3,file='du/dt.DAT') i could find a .dat file with the results.The work will continue:)Best regards to all!!!
 bigduke88 (TechnicalUser) 8 May 12 14:03
 Excuse me for the repost, but I've forgotten to post the new source code. If it is interessting, I could explain, what is the exact purpose of my Fortran "battle".The source code:      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      DIMENSION U(100)      open(unit=3,file='du/dt.DAT')      NIT=2000      jprint=500      DT1=1D-8      DT=DT1      DTau=4.D0*DT      Dp=0.00271D0      Ro=998.2D0      Mu=0.001005D0      Diff=2.97271D-9      Tau=(4.D0*Diff*T)/(Dp**2.D0)      T=(Tau*(Dp**2.D0))/(Diff*4.D0)      imax=100      U(i)=0.      T=0.      DO 100 i=1,imax      T=T+DT      Tau=Tau+DTau      Re=(Dp*U(i)*Ro)/(Mu)      Cd=(24/Re)*(1+0.125*(Re**0.72))      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT      Z=Z+DT*(U(i+1)+U(i))100   continue      write(18,*) Z,T,Tau      STOP      END
 gummibaer (Programmer) 8 May 12 15:23
 The filename you selected seems somewhat inappropriate ('/' is not legal in filenames). I wonder why you do not receive a runtime error on this.Norbert  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 salgerman (Programmer) 8 May 12 17:01
 Nothing wrong with "/" in file name...just think of it as being in a sub-directory.bigduke:  is the *.dat file in a subfolder named "du" ?
 gummibaer (Programmer) 8 May 12 17:16
 ohh, that's new to me. I allways used '\' for subdirectories instead of '/'Learning never ends....Norbert The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 bigduke88 (TechnicalUser) 9 May 12 11:33
 Hello everyonesalgerman - yes I've created a subfolder named "dU_dT", because windows does not accept the "/". At this folder I can just find some .dat files, an application file named source and an F file. I was wondering, why there is no text file with the reuslts of the calculations...This is the newer version of the source code:      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      PARAMETER (NR=100)      DIMENSION U(NR)      open(UNIT=13,STATUS='new',file='dU\dT.TXT')      NIT=2000      jprint=500      KPRINT=1      DT1=1D-8      DT=DT1      DTau=4.D0*DT      Dp=0.00271D0      Ro=998.2D0      Mu=0.001005D0      Diff=2.97271D-9      Tau=(4.D0*Diff*T)/(Dp**2.D0)      T=(Tau*(Dp**2.D0))/(Diff*4.D0)      imax=NR      U(i)=0.      T=0.      DO 100 i=1,imax      T=T+DT      Tau=Tau+DTau      Re=(Dp*U(i)*Ro)/(Mu)      Cd=(24/Re)*(1+0.125*(Re**0.72))      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DTC      DISTANCE      Z=Z+DT*(U(I)+U(I+1))100   continue      write(UNIT=13,FMT=*) Z,T,U      STOP      ENDAnd there are no problems about running it, but there are no results and I don't know why.The results should be separated in three colums: one for the distance "Z", one for the time "T" and one for the velocity "U".Now I'm going to work, because there was plenty of other work today.Best regrads!!!
 gummibaer (Programmer) 9 May 12 11:55
 You still have the '\' in the open-statement. If '\' creates a subfolder, then your program creates the subfolder. In your current directory you should have a directory called dU and in there you may find a file dt.txt. At least that is what I would expect from your code (with what I have learned ).Rewrite your open to something likeopen (unit = 13, file = 'dU_by_dT.txt')and after your prog executed you will find a textfile of that name in your current directory. Note, I skipped the 'new' clause, because you will receive a runtime error for every run of your program when you did not delete your outputfile from the run before first. Could be annoying while debugging.Norbert  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 salgerman (Programmer) 9 May 12 12:32
 bigduke:Hhhmmm, there might be some kind of confusion here...I am not sure if YOU mean for "\" to be part of the filename alone...which it wouldn't...I think Windows is interpreting "\" as part of a relative full path name, which includes a subdirectory "dU"  and a file named "dT.TXT".So, I will ask again...do you have a directory named "dU" ?!! Not "dU_dT", but "dU"?Is the program actually running? Or is it crashing because it cannot find the file "dU\dT.TXT"?Are you running the program from the command line to see feedback or are you simply double-clicking on the executable?Also, follow Norbert's suggestion about the name change...that should most probably clear up your problem.
 bigduke88 (TechnicalUser) 9 May 12 14:12
 Thank you, I've done this : open(UNIT=13,STATUS='new',file='dU_by_dT.TXT')and at the end of the programwrite(UNIT=13,FMT=*) 'Z,T,U'write(unit=13,FMT=*) Z,T,UAnd it creates a text document named dU_by_dT, but there are only Z,T,U as head of column, but there are no data. When I start the program, it initiates a window, where normally I colud see the results of the iterations, but this "window" just appears for a moment and then, nothing happens.Maybe the mistake is in the print statements?
 salgerman (Programmer) 9 May 12 14:28
 Please post the entire code, once again...the one above has the write statement outside the LOOP!
 bigduke88 (TechnicalUser) 9 May 12 16:54
 Here is the code, I've tried to fix the things, but no success :(The text file is still empty      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      PARAMETER (NR=100)      DIMENSION U(NR)      open(UNIT=13,STATUS='new',file='dU_by_dT.txt')      NIT=2500      jprint=500      KPRINT=1      DT1=1D-8      DT=DT1      DTau=4.D0*DT      Dp=0.00271D0      Ro=998.2D0      Mu=0.001005D0      Diff=2.97271D-9      Tau=(4.D0*Diff*T)/(Dp**2.D0)      T=(Tau*(Dp**2.D0))/(Diff*4.D0)      imax=NR      U(i)=0.      T=0.      DO 100 i=1,imax      T=T+DT      Tau=Tau+DTau      Re=(Dp*U(i)*Ro)/(Mu)      Cd=(24/Re)*(1+0.125*(Re**0.72))      U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DTC      DISTANCE      Z=Z+DT*(U(I)+U(I+1))100   continue      write(UNIT=13,FMT=*) 'Z,T,U'      write(unit=13,FMT=*) Z,T,U      STOP      END
 salgerman (Programmer) 9 May 12 20:59
 bigduke:You program has a few problems.Line 15, you calculate Tau as funtion of T, which has not been initialized, so T=0  => Tau=0Line 16, you calculate T as function of Tau, Tau=0  =>  T=0Line 18, U(i) = 0.0; but i has not been initialized, so i=0, U(0)=0Line 24, for i=1, Re is 0.0 because U(i=1) has not been initialized and so U(1)=0Line 25, Cd calculation will divide by zero because of 24/ReAlso, Z and T are scalars, if you care to see their intermediate values, you may want to move the write statement into the loop and write Z,T,U(i)You may want to move the header line to before the loop.If you just wants to see the final Z and T and the entire U vector, you can leave the write statements where they are.....and I am afraid I am running out of good mood, here.
 bigduke88 (TechnicalUser) 10 May 12 2:34
 Thank you for the help!!!Please don't be running out fo good mood, obviously the "art of programming" is beyond my skills, at this moment, but hpefully I'll advance in it.... but the progress is really slow...Best regards!!!
 gummibaer (Programmer) 10 May 12 6:45
 And...In your do-loop (line 20) i runs to imax which is Nr which is 100You reference u(i+1) in line 21, whereas u is diemnsioned to NR = 100. This will give an 'array bounds exceeded' error. You should set iMax = NR - 1Norbert  The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 PaulWDent (TechnicalUser) 18 May 12 12:55
 If you have never programmed before, it is good to start with solving the problem with an EXCEL spreadsheet. Your problem is very appropriate to solve on an Excel spread sheet. You jsut have columns labelled T, U, Cd, Z etc along the top, (columsn A, B, C, Din axecl) and then you start with 0.0 in box A2, put the initial velocity in box B2 and I guess the initial integral value Z is zero in box D2. Then in boc A3 you write "=A1 + 1.0E-8" which adds dT (do you really want 10 nanoseconds!)Then copy that box A3 to A4 to A101 and you have your T's.Then in box C2 write the formula that calculates the drag coeffcient Cd from the velocity in box B2Then in box B3 you write the equation that updates the initial velocity at T=0.0 to thevelociy at T+dT based on the drag coefficient in box C2. Then copy and paste from boc C2 to C3.In box D3 you write ="D2+1.0E-8*(C2+C3)/2.0 which is doing trapezoidal rule integration for ZNow just copy and paste from boxes A3,B3,C3,D3 into boxes A4,B4,C4,D4 to 101 and youwill have your answer.You will find it a much easier step from there to FORTRANNow one of the big problems with your thinking is this: Drag reduces velocity, so the particle is slowing. But from what intial velocity? You have set the initial velocity to zero!!You must know the initial velocity to solve this problem.You must write, before the Do-loop, U(1) = something positive and probably large.Alternatively, you could start with a very small velocity (not zero) and work in time-reversed order to find out what the initial velocity would have to be 10uS ago to have a terminal velocity of say 1m/sec 10uS later.Another thing you must do for physical problems is to work in real engineering units, i.e. all quantities should be in the meter-kilo-second units, as otherwise you cannot write Acceleration = Force/Mass without some fudge factors. Also, with a Do-loop of 100 and dT = 1e-8, you do realize that your simulation is only simulating 10uS of real time, do you?
 gummibaer (Programmer) 18 May 12 14:08

#### Quote:

Another thing you must do for physical problems is to work in real engineering units

Depends on the point of view.
If you are dealing with physical scientific problems you should work with dimension-free numbers to do your calculations. I am not a specialist for two phase flow but I would guess Reynolds-Number, Grashof-Number and maybe a few others would be present in this problem when you did it.

If you are dealing with an engineering problem or this is a side-aspect of your task only, then you most probably will do your computations in engineering units.

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!