if statement to determine steady-state

user_4538674 0
Viewed 81

My code below correctly solves a 1D heat equation for a function u(x,t). I now want to find the steady-state solution, the solution that no longer changes in time so it should satisfy u(t+1)-u(t) = 0. What is the most efficient way to find the steady-state solution? I show three different attempts below, but I'm not sure if either are actually doing what I want. The first and third have correct syntax, the second method has a syntax error due to the if statement. Each method is different due to the change in the if structure.

Method 1 :

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: this checks the solutions at every space point  which I don't think is correct.
    do i = 1,n-1
       if ( v(i) - u(i) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if
    end do

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1

Method 2 :

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: (This gives an error message since the if statement doesn't have a logical scalar expression, but I want to compare the full arrays v and u as shown.
       if ( v - u .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1

Method 3 :

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: Perhaps this is the correct expression I want to use
       if( norm2(v) - norm2(u) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1
fortran
1 Answers

Without discussing which method to determine "closeness" is best or correct (not really being a programming problem) we can focus on what the Fortran parts of the methods are doing.

Method 1 and Method 2 are similar ideas (but broken in their execution), while Method 3 is different (and broken in another way).

Note also that in general one wants to compare the magnitude of the difference abs(v-u) rather than the (signed) difference v-u. With non-monotonic changes over iterations these are quite different.

Method 3 uses norm2(v) - norm2(u) to test whether the arrays u and v are similar. This isn't correct. Consider

norm2([1.,0.])-norm2([0.,1.])

instead of the more correct

norm2([1.,0.]-[0.,1.])

Method 2's

if ( v - u .LT. 1.0e-7 ) then

has the problem of being an invalid array expression, but the "are all points close?" can be written appropriately as

if ( ALL( v - u .LT. 1.0e-7 )) then

(You'll find other questions around here about such array reductions).

Method 1 tries something similar, but incorrectly:

    do i = 1,n-1
       if ( v(i) - u(i) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if
    end do

This is incorrect in one big way, and one subtle way.

First, the loop is exited when the condition tests true the first time, with a message saying the steady state is reached. This is incorrect: you need all values close, while this is testing for any value close.

Second, when the condition is met, you exit. But you don't exit the time iteration loop, you exit the closeness testing loop. (exit without a construct name leaves the innermost do construct). You'll be in exactly the same situation, running again immediately after this innermost construct whether the tested condition is ever or never met (if ever met you'll get the message also). You will need to use a construct name on the time loop.

I won't show how to do that (again there are other questions here about that), because you also need to fix the test condition, by which point you'll be better off using if(all(... (corrected Method 2) without that additional do construct.

For Methods 1 and 2 you'll have something like:

if (all(v-u .lt 1e-7)) then
  print *, "Converged"
  exit
end if

And for Method 3:

if (norm2(v-u) .lt. 1e-7) then
  print *, "Converged"
  exit
end if
user_31570760
Related