lundi 21 mai 2018

Program not entering an if branch in Fortran even when condition is satisfied

This is the section of the code that is not running as expected. I had to include some old-style "pause" statements to pin-point the error location.


    iteration_index = 1
    y_refined = y_current + (0.5d0*dx*(dydx_predictor + dydx_refined_corrector))               ! Solution corresponds to i = 2, i.e., the refined Heun's method.

    relative_percent_error = 100.d0*(abs((y_refined - y_next)/y_refined))                         ! Calculation of relative percentage error. This is NOT true error.

    if (relative_percent_error > heun_percent_tolerance) then

        iteration_index = iteration_index + 1

        print*, 'first loop enter', x_next, relative_percent_error, iteration_index
        pause

        if (iteration_index < max_heun_number) then

            y_next = y_refined
            call dydx(x_next, y_next, dydx_refined_corrector)
            y_refined = y_current + (0.5d0*dx*(dydx_predictor + dydx_refined_corrector))
            relative_percent_error = 100.d0*(abs((y_refined - y_next)/y_refined))

            print*, 'second loop enter', x_next, relative_percent_error, iteration_index
            pause

        end if

    end if


Output is seen as follows:

first loop enter   1.0000000000000000        6.7763423346068707          2

PAUSE

To resume execution, type go.  Other input will terminate the job.

go

RESUMED

second loop enter   1.0000000000000000        1.6658644147581094        2

PAUSE 

To resume execution, type go.  Other input will terminate the job.

go

RESUMED

first loop enter   2.0000000000000000        6.6615482639252761         2

PAUSE 

To resume execution, type go.  Other input will terminate the job.


The values of heun_percent_tolerance is 0.01 and max_heun_number is 15. I expect the execution to enter the second if condition for more iterations until the max limit of 15 is reached, but it seems that the code jumps to the next x_next value of 2.

I have even tried to combine the two conditions as If (cond1 .and. cond2), that did not work either.

Aucun commentaire:

Enregistrer un commentaire