Program to numerically solve 3 ODE's

Help Homework
1

Hi all,
I am tasked with writing a program to solve Stagnation point flow (Hiemenz) flow. The formulation of the solution requires us to guess at a value for f’’ and watch for boundary conditions to converge.

I’ve already received so much help here and have enough snips of code running that I’ve proved to myself it can be done - but now I need to write my own code and feel ready to do it. If you guys don’t mind and are interested I’d like to post my progress here and if you can comment I’d be most thankful.

I haven’t really programmed since Fortran 77 and some of the variable format codes look intimidating. So I’d like to ask for help along the way.

As part of the solution I have chosen to use a “Bisection” method to hone in on the correct value of f’’ and that must be part of the program. To get my thoughts in order I wrote a flow chart, for those curious I’l post it below.

Much Thanks in advance and thank you for following.
Kind regards,
Rick

Screen Shot 2021-11-15 at 11.54.21 AM
Screen Shot 2021-11-15 at 11.54.21 AM545×707 139 KB

2 Likes
2

Go ahead and post your progress here. Yes, the “shooting method” (using bisection to converge the BC) is a common method to solve boundary value problems.

3

Hey Thanks a lot - Shooting method = Bisection okay!

Some Code I looked at had the following format descriptor. I’m kind of scratching my head here. Isn’t “1d0” just an integer?

integer, parameter :: dp = kind(1d0)
1 Like
4

No, it’s double precision 1.0. It can also be written integer, parameter :: dp = kind(1.0d0), which is what I do.

2 Likes
5

Okay, so then I see later on in the program a variable defined later with the formatting, is this not redundant or not required?

integer, parameter :: dp = kind(1d0)
    real(dp) :: tstart
    tstart = 0.0d0

tstart is 0.0d0? What is going on here?
RP

6

It is not redundant. Declaring a variable does not initialize it, so it is necessary to initialize tstart before using it.

7

I am really struggling with the syntax of the formatting. 1d0 says to use one digit with zero digits to the right of the decimal? then what does tstart = 0.0d0 mean?
RP

8

No, 1d0 is just an example of a double precision constant. Any double precision constant could have used as the argument of kind. Note that

print*,kind(1d0),kind(1.0d0),kind(1.000d0)

will print the same integer (8 for gfortran and Intel Fortran) thrice. Similarly,

print*,kind(1e0),kind(1.0e0),kind(1.000e0)

will print the kind corresponding to single precision (4 for gfortran and Intel Fortran) thrice.

1 Like
9

The d in 1.0d0 is a form of e-notation used to represent “times ten raised to the power of”.

For single precision real literal constants, the letter used is e, while for double precision reals, d is used instead.

1 Like
10

Oh, excellent, Thank you Beliavsky and Ivanpribec!

Next topic, subroutine structure, This is a rather simple program, I would like to avoid having to call any external subroutines.

Can I have an internal function and subroutine in the same program? I’ll have both I think.

The Contains function seems really neat, but it also appears a little advanced. I’d like to write something on the simple side.

What structure do you like?
RP

11

Yes you can nest several subroutines and/or functions in the contains section.

1 Like
12

Yes, you can, but internal functions and subroutines inherit all the variables defined in the main program. To increase modularity it is generally better to define procedures in a separate module and use the module in the main program to get access to those procedures.

2 Likes
13

The test “g = 0” in your flow chart is too strict for a finite-precision calculation. Replace it with something reasonable such as |g| < 10-8.

1 Like
14

Hi Ivan,
I am struggling to understand the “Contains” section. Unlike a Subroutine or a Function which I have to interface in the beginning of a program, a Contains section I do not?
Can you recommend a good reference for program branching so I can study the options?
Much Thanks
Rick

15

One resource is the fortran-lang tutorial: Organising code structure - Fortran Programming Language

If available to you, I would recommend the latest edition of “Modern Fortran explained” by Metcalf, Reid & Cohen, chapter 5 - Program units and procedures.

A second option would be “Guide to Fortran 2008 Programming” by Walter Brainerd: Guide to Fortran 2008 Programming | SpringerLink. Perhaps you can access the digital version with your university account.

The keyword contains is simply used to collect internal subprograms (subroutines and functions) in a main program or module:

program main
  implicit none
  
  ! ... declarations ...

  ! ... main program unit commands ...

contains

  ! ... internal subprograms ...

  subroutine runge(...)

  end subroutine

  function dydt(t,y)

  end function

end program

Instead of placing the procedures in the main program, you can place them in a module:

module flow_solver_routines
  implicit none

  ! ... declarations ...

contains

  ! ... internal subprograms ...

  subroutine runge(...)

  end subroutine

  function dydt(t,y)

  end function

end module

that you then use in the main program as follows:

program main
  use flow_solver_routines, only: runge, dydt
  implicit none

  ! ...

  call runge(dydt,tstart,tend,y)

end program

In fact even subroutines and functions can have procedures nested in their own contains section.

On the other hand, your old program used so-called external subprograms, i.e. subprograms that are not located within a program unit or module. External procedures were the norm in F77, but since F90 modules are the preferred approach to organize your (sub)programs.

When using internal subprograms the compiler can check you are calling the function or subroutine correctly since the interface is visible. Internal subprograms also have access to the objects in their host program, which can be useful for many things.

2 Likes
16

The literature of existing software to solve boundary-value problems is extensive. If this is a homework problem, then I guess you need to attack the problem yourself. Shooting is one way. Multiple shooting, where the code gives up and restarts when errors look too large, combined with a least-squares method to minimize the “jumps” at restarts, and collocation, are other methods. If your problem is a Sturm-Liouville problem, and you need its eigenvalues, I recommend sleign (or dleign).

But if it’s a “real” problem you should just get an existing code. Start at netlib.org.

1 Like
17

vsnyder,
This is a programming exercise, so I’ll have to write something in my own code. But thank you for the advice, you are right on.
Rick

1 Like
18

Hi Ivan,
Okay, Can I do the following, it doesn’t seem to work. I want to put the use the function within the subroutine. Is there something fundamentally wrong with this?
Rick

program main
  implicit none
  
  ! ... declarations ...

  ! ... main program unit commands ...

	Call runge
	
contains

  subroutine runge(...)
	
	dydt(t,y)
	
  end subroutine

  function dydt(t,y)

  end function

end program
19

Assuming the calling sequences are correct, it looks okay to me. Are you passing the function dydt as a dummy variable to the subroutine, or is it just found directly since it’s located in the same program scope? (The function call should probably be assigned to a left-hand side.)

20

@RickP330 , you may find it very helpful to go through this book by Chapman first and then follow-up on online forums such as this one with targeted questions that help clarify aspects related to your particular exercise.

1 Like