Hello everyone
I 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 + Ui
This 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!!!

 

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.

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

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.

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



 

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.

Hello everyone
It 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,imax
30      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

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
Instead, you can simply write

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.
 

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

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.

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

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

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.

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" ?  

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.

Hello everyone
salgerman - 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*DT
C      DISTANCE
      Z=Z+DT*(U(I)+U(I+1))

100   continue
      write(UNIT=13,FMT=*) Z,T,U
      STOP
      END


And 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!!!

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 like

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

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.

Thank you, I've done this : open(UNIT=13,STATUS='new',file='dU_by_dT.TXT')
and at the end of the program
write(UNIT=13,FMT=*) 'Z,T,U'
write(unit=13,FMT=*) Z,T,U
And 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?

Please post the entire code, once again...the one above has the write statement outside the LOOP!

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*DT
C      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

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=0
Line 16, you calculate T as function of Tau, Tau=0  =>  T=0
Line 18, U(i) = 0.0; but i has not been initialized, so i=0, U(0)=0
Line 24, for i=1, Re is 0.0 because U(i=1) has not been initialized and so U(1)=0
Line 25, Cd calculation will divide by zero because of 24/Re

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

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

And...

In your do-loop (line 20) i runs to imax which is Nr which is 100
You 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 - 1

Norbert
 

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

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 B2
Then in box B3 you write the equation that updates the initial velocity at T=0.0 to the
velociy 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 Z

Now just copy and paste from boxes A3,B3,C3,D3 into boxes A4,B4,C4,D4 to 101 and you
will have your answer.

You will find it a much easier step from there to FORTRAN

Now 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?

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.