Changeset 727


Ignore:
Timestamp:
05/09/12 15:10:03 (13 months ago)
Author:
micael
Message:
  • Renamed some functions and variables to improve clarity.
  • Added some documentation.
Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/ode_integrator.F90

    r723 r727  
    8989contains 
    9090 
    91   subroutine ode_integrator_null(integrator) 
    92     !-----------------------------------------------------------------------! 
    93     !-----------------------------------------------------------------------! 
    94     type(ode_integrator_t), intent(out) :: integrator 
    95  
    96     integrator%ode = 0 
    97     integrator%dim = 0 
    98     integrator%nstepmax = 0 
    99     integrator%qn = QN_NULL 
    100     integrator%e = M_ZERO 
    101     call potential_null(integrator%potential) 
     91  subroutine ode_integrator_null(odeint) 
     92    !-----------------------------------------------------------------------! 
     93    ! Nullifies and sets to zero all the components of the ODE integrator.  ! 
     94    !-----------------------------------------------------------------------! 
     95    type(ode_integrator_t), intent(out) :: odeint 
     96 
     97    odeint%ode = 0 
     98    odeint%dim = 0 
     99    odeint%nstepmax = 0 
     100    odeint%qn = QN_NULL 
     101    odeint%e = M_ZERO 
     102    call potential_null(odeint%potential) 
    102103 
    103104  end subroutine ode_integrator_null 
    104105 
    105   subroutine ode_integrator_init(integrator, ode, stepping_function, nstepmax, tol, qn, ev, potential) 
    106     !-----------------------------------------------------------------------! 
    107     !  qn         - set of quantum numbers                                  ! 
    108     !  ev         - energy                                                  ! 
    109     !  potential  - potential object                                        ! 
    110     !-----------------------------------------------------------------------! 
    111     type(ode_integrator_t), intent(inout) :: integrator 
     106  subroutine ode_integrator_init(odeint, ode, stepping_function, nstepmax, tol, qn, ev, potential) 
     107    !-----------------------------------------------------------------------! 
     108    ! Initialize the ODE integrator.                                        ! 
     109    !                                                                       ! 
     110    !  odeint          - ODE integrator                                     ! 
     111    !  ode             - equation to be integrated                          ! 
     112    !  stepping_function - stepping function used during integration        ! 
     113    !  nstepmax        - maximum number of steps allowed during integration ! 
     114    !  tol             - tolerance                                          ! 
     115    !  qn              - set of quantum numbers                             ! 
     116    !  ev              - energy                                             ! 
     117    !  potential       - potential object                                   ! 
     118    !-----------------------------------------------------------------------! 
     119    type(ode_integrator_t), intent(inout) :: odeint 
    112120    integer,                intent(in)    :: ode 
    113121    integer,                intent(in)    :: stepping_function 
    … …  
    120128    call push_sub("ode_integrator_init") 
    121129 
    122     integrator%ode = ode 
    123     integrator%dim = 2 
    124     if (ode == ODE_DIRAC_POL2) integrator%dim = 8 
    125     integrator%nstepmax = nstepmax 
    126     integrator%qn = qn 
    127     integrator%e = ev 
    128     integrator%potential = potential 
     130    odeint%ode = ode 
     131    odeint%dim = 2 
     132    if (ode == ODE_DIRAC_POL2) odeint%dim = 8 
     133    odeint%nstepmax = nstepmax 
     134    odeint%qn = qn 
     135    odeint%e = ev 
     136    odeint%potential = potential 
    129137 
    130138    !Initialize GSL objects 
    131     call gsl_odeiv_step_alloc(stepping_function, integrator%dim, integrator%stp) 
    132     call gsl_odeiv_evolve_alloc(integrator%dim, integrator%evl) 
    133     call gsl_odeiv_control_standart_new(integrator%ctrl, 1.0E-22_r8, tol, M_ONE, M_ONE) 
     139    call gsl_odeiv_step_alloc(stepping_function, odeint%dim, odeint%stp) 
     140    call gsl_odeiv_evolve_alloc(odeint%dim, odeint%evl) 
     141    call gsl_odeiv_control_standart_new(odeint%ctrl, 1.0E-22_r8, tol, M_ONE, M_ONE) 
    134142 
    135143    call pop_sub() 
    136144  end subroutine ode_integrator_init 
    137145 
    138   subroutine ode_integrator_end(integrator) 
     146  subroutine ode_integrator_end(odeint) 
    139147    !-----------------------------------------------------------------------! 
    140148    ! Frees all the memory associated with the GSL objects used by the      ! 
    141149    ! ODE solver.                                                           ! 
    142150    !-----------------------------------------------------------------------! 
    143     type(ode_integrator_t), intent(inout) :: integrator 
     151    type(ode_integrator_t), intent(inout) :: odeint 
    144152 
    145153    call push_sub("ode_integrator_end") 
    146154 
    147     if (integrator%nstepmax > 0) then 
     155    if (odeint%nstepmax > 0) then 
    148156      !Free GSL objects for simple precision integrations 
    149       call gsl_odeiv_step_free(integrator%stp) 
    150       call gsl_odeiv_evolve_free(integrator%evl) 
    151       call gsl_odeiv_control_free(integrator%ctrl) 
     157      call gsl_odeiv_step_free(odeint%stp) 
     158      call gsl_odeiv_evolve_free(odeint%evl) 
     159      call gsl_odeiv_control_free(odeint%ctrl) 
    152160    end if 
    153161 
    154     integrator%ode = 0 
    155     integrator%dim = 0 
    156     integrator%nstepmax = 0 
    157     integrator%qn = QN_NULL 
    158     integrator%e = 0 
    159     call potential_end(integrator%potential) 
     162    odeint%ode = 0 
     163    odeint%dim = 0 
     164    odeint%nstepmax = 0 
     165    odeint%qn = QN_NULL 
     166    odeint%e = 0 
     167    call potential_end(odeint%potential) 
    160168 
    161169    call pop_sub() 
    162170  end subroutine ode_integrator_end 
    163171 
    164   subroutine ode_integration(integrator, r1, r2, nstep, r, f) 
     172  subroutine ode_integration(odeint, r1, r2, nstep, r, f) 
    165173    !-----------------------------------------------------------------------! 
    166174    ! Integrates the radial wave-equation from r1 to r2.                    ! 
    167175    !                                                                       ! 
    168     !  integrator - ode integrator object                                   ! 
    169     !  r1         - starting point for integration                          ! 
    170     !  r2         - final point for integration                             ! 
    171     !  nstep      - number of steps taken by the solver                     ! 
    172     !  r          - mesh used by the ODE solver                             ! 
    173     !  f          - functions                                               ! 
    174     !-----------------------------------------------------------------------! 
    175     type(ode_integrator_t), intent(in)  :: integrator 
     176    !  odeint - ode integrator object                                       ! 
     177    !  r1     - starting point for integration                              ! 
     178    !  r2     - final point for integration                                 ! 
     179    !  nstep  - number of steps taken by the solver                         ! 
     180    !  r      - mesh used by the ODE solver                                 ! 
     181    !  f      - functions                                                   ! 
     182    !-----------------------------------------------------------------------! 
     183    type(ode_integrator_t), intent(in)  :: odeint 
    176184    real(R8),               intent(in)  :: r1, r2 
    177185    integer,                intent(out) :: nstep 
    … …  
    184192 
    185193    !Allocate work arrays 
    186     allocate(rmax(integrator%nstepmax), fmax(integrator%nstepmax, integrator%dim)) 
     194    allocate(rmax(odeint%nstepmax), fmax(odeint%nstepmax, odeint%dim)) 
    187195 
    188196    !Set initial points for integration 
    … …  
    190198    if (r1 < r2) then 
    191199      !Outward integration 
    192       call ode_bc_origin(integrator, rmax(1), fmax(1,:)) 
     200      call ode_bc_origin(odeint, rmax(1), fmax(1,:)) 
    193201    else 
    194202      !Inward integration 
    195       call ode_bc_infinity(integrator, rmax(1), fmax(1,:)) 
     203      call ode_bc_infinity(odeint, rmax(1), fmax(1,:)) 
    196204    end if 
    197205 
    198206    !Integrate the equation from r1 to r2 
    199     ierr = ode_solve(r2, integrator%stp, integrator%evl, integrator%ctrl, & 
    200                      (r1+r2)/M_TWO, integrator%nstepmax, nstep, rmax, & 
    201                      integrator%dim, fmax, integrator, ode_derivatives) 
     207    ierr = ode_solve(r2, odeint%stp, odeint%evl, odeint%ctrl, & 
     208                     (r1+r2)/M_TWO, odeint%nstepmax, nstep, rmax, & 
     209                     odeint%dim, fmax, odeint, ode_derivatives) 
    202210    if (ierr /= 0) then 
    203211      if (in_debug_mode) then 
    204         call ode_integrator_debug(integrator, r1, r2, nstep, rmax, fmax) 
     212        call ode_integrator_debug(odeint, r1, r2, nstep, rmax, fmax) 
    205213      end if 
    206214      message(1) = "Error in subtoutine ode_integration. Error message:" 
    … …  
    210218 
    211219    !Copy the functions to new arrays 
    212     allocate(r(nstep), f(nstep, integrator%dim)) 
     220    allocate(r(nstep), f(nstep, odeint%dim)) 
    213221    r(1:nstep) = rmax(1:nstep) 
    214     f(1:nstep,1:integrator%dim) = fmax(1:nstep,1:integrator%dim) 
     222    f(1:nstep,1:odeint%dim) = fmax(1:nstep,1:odeint%dim) 
    215223 
    216224    !Deallocate arrays 
    … …  
    220228  end subroutine ode_integration 
    221229 
    222   subroutine ode_function_to_wavefunction(integrator, r, f, wf, wfp) 
    223     !-----------------------------------------------------------------------! 
    224     !-----------------------------------------------------------------------! 
    225     type(ode_integrator_t), intent(in)  :: integrator 
     230  subroutine ode_function_to_wavefunction(odeint, r, f, wf, wfp) 
     231    !-----------------------------------------------------------------------! 
     232    ! Given a set of solutions to the ODE, returns the corresponding        ! 
     233    ! wave-functions and wave-functions derivatives.                        ! 
     234    !                                                                       ! 
     235    !  odeint - ode integrator object                                       ! 
     236    !  r      - r coodinate                                                 ! 
     237    !  f      - ODE solution at r                                           ! 
     238    !  wf     - wave-functions at r                                         ! 
     239    !  wfp    - wave-functions derivatives at r                             ! 
     240    !-----------------------------------------------------------------------! 
     241    type(ode_integrator_t), intent(in)  :: odeint 
    226242    real(R8),               intent(in)  :: r 
    227     real(R8),               intent(in)  :: f(integrator%dim) 
     243    real(R8),               intent(in)  :: f(odeint%dim) 
    228244    real(R8),               intent(out) :: wf(:), wfp(:) 
    229245 
    230     real(R8) :: fp(integrator%dim) 
     246    real(R8) :: fp(odeint%dim) 
    231247 
    232248    call push_sub("ode_function_to_wavefunction") 
    233249 
    234     call ode_derivatives(integrator, r, f, fp) 
    235  
    236     select case (integrator%ode) 
     250    call ode_derivatives(odeint, r, f, fp) 
     251 
     252    select case (odeint%ode) 
    237253    case (ODE_SCHRODINGER, ODE_SCALAR_REL) 
    238254      wf(1) = f(1) 
    … …  
    256272  end subroutine ode_function_to_wavefunction 
    257273 
    258   subroutine ode_match_functions(integrator, nstep_out, nstep_in, f_out, f_in) 
    259     !-----------------------------------------------------------------------! 
    260     !-----------------------------------------------------------------------! 
    261     type(ode_integrator_t), intent(in)    :: integrator 
     274  subroutine ode_match_functions(odeint, nstep_out, nstep_in, f_out, f_in) 
     275    !-----------------------------------------------------------------------! 
     276    ! Given the outward and inward solutions to the ODE, match all          ! 
     277    ! functions minus one.                                                  ! 
     278    !                                                                       ! 
     279    !  odeint    - ode integrator object                                    ! 
     280    !  nstep_out - number of points of f_out                                ! 
     281    !  nstep_in  - number of points of f_in                                 ! 
     282    !  f_out     - outward solutions                                        ! 
     283    !  f_in      - inward solutions                                         ! 
     284    !-----------------------------------------------------------------------! 
     285    type(ode_integrator_t), intent(in)    :: odeint 
    262286    integer,                intent(in)    :: nstep_out, nstep_in 
    263     real(R8),               intent(inout) :: f_out(nstep_out, integrator%dim) 
    264     real(R8),               intent(inout) :: f_in (nstep_in,  integrator%dim) 
     287    real(R8),               intent(inout) :: f_out(nstep_out, odeint%dim) 
     288    real(R8),               intent(inout) :: f_in (nstep_in,  odeint%dim) 
    265289 
    266290    real(R8) :: a, b, c 
    … …  
    272296    call push_sub("ode_match_functions") 
    273297 
    274     select case (integrator%ode) 
     298    select case (odeint%ode) 
    275299    case (ODE_SCHRODINGER, ODE_SCALAR_REL, ODE_DIRAC, ODE_DIRAC_POL1) 
    276300      a = f_in(nstep_in,1)/f_out(nstep_out,1) 
    … …  
    323347  end subroutine ode_match_functions 
    324348 
    325   function ode_function_mismatch(integrator, r, f_out, f_in) 
    326     !-----------------------------------------------------------------------! 
    327     !-----------------------------------------------------------------------! 
    328     type(ode_integrator_t), intent(in) :: integrator 
     349  function ode_function_mismatch(odeint, r, f_out, f_in) 
     350    !-----------------------------------------------------------------------! 
     351    ! Given the outward and inward solutions to the ODE, match all          ! 
     352    ! functions minus one and return the mismatch of the remaining function.! 
     353    !                                                                       ! 
     354    !  odeint    - ode integrator object                                    ! 
     355    !  r         - coordinate where to evaluate the mismatch                ! 
     356    !  f_out     - outward solutions                                        ! 
     357    !  f_in      - inward solutions                                         ! 
     358    !-----------------------------------------------------------------------! 
     359    type(ode_integrator_t), intent(in) :: odeint 
    329360    real(R8),               intent(in) :: r 
    330     real(R8),               intent(in) :: f_out(integrator%dim) 
    331     real(R8),               intent(in) :: f_in(integrator%dim) 
     361    real(R8),               intent(in) :: f_out(odeint%dim) 
     362    real(R8),               intent(in) :: f_in(odeint%dim) 
    332363    real(R8) :: ode_function_mismatch 
    333364 
    334     real(R8) :: fo(1, integrator%dim), fi(1, integrator%dim) 
    335     real(R8) :: fpo(integrator%dim), fpi(integrator%dim) 
     365    real(R8) :: fo(1, odeint%dim), fi(1, odeint%dim) 
     366    real(R8) :: fpo(odeint%dim), fpi(odeint%dim) 
    336367 
    337368    call push_sub("ode_function_mismatch") 
    … …  
    340371    fo(1, :) = f_out(:) 
    341372    fi(1, :) = f_in(:) 
    342     call ode_match_functions(integrator, 1, 1, fo, fi) 
     373    call ode_match_functions(odeint, 1, 1, fo, fi) 
    343374 
    344375    !Calculate the function mismatch 
    345     select case (integrator%ode) 
     376    select case (odeint%ode) 
    346377    case (ODE_SCHRODINGER, ODE_SCALAR_REL, ODE_DIRAC, ODE_DIRAC_POL1) 
    347       call ode_derivatives(integrator, r, fo(1,:), fpo) 
    348       call ode_derivatives(integrator, r, fi(1,:), fpi) 
     378      call ode_derivatives(odeint, r, fo(1,:), fpo) 
     379      call ode_derivatives(odeint, r, fi(1,:), fpi) 
    349380      ode_function_mismatch = fpo(1)/fo(1,1) - fpi(1)/fi(1,1) 
    350381 
    … …  
    361392  end function ode_function_mismatch 
    362393 
    363   subroutine ode_bc_origin(integrator, r0, f) 
     394  subroutine ode_bc_origin(odeint, r0, f) 
    364395    !-----------------------------------------------------------------------! 
    365396    ! Set the value of the functions at a point close to the origin.        ! 
    366397    !                                                                       !     
    367     !  integrator - integrator object                                       ! 
    368     !  r0         - point close to origin                                   ! 
    369     !  f          - function at r0                                          ! 
    370     !-----------------------------------------------------------------------! 
    371     type(ode_integrator_t), intent(in)  :: integrator 
     398    !  odeint - integrator object                                           ! 
     399    !  r0     - point close to origin                                       ! 
     400    !  f      - function at r0                                              ! 
     401    !-----------------------------------------------------------------------! 
     402    type(ode_integrator_t), intent(in)  :: odeint 
    372403    real(R8),               intent(in)  :: r0 
    373     real(R8),               intent(out) :: f(integrator%dim) 
     404    real(R8),               intent(out) :: f(odeint%dim) 
    374405 
    375406    real(R8) :: e, z, vp0, s, a0, a1, a2, a3, b0, b1, b2, w 
    … …  
    431462    !       equation Z=0 is a special case. 
    432463 
    433     z = potential_nuclear_charge(integrator%potential) 
    434     vp0 = v(integrator%potential, r0, integrator%qn) + z/r0 
    435     e = integrator%e 
    436     qn = integrator%qn 
    437  
    438     select case (integrator%ode) 
     464    z = potential_nuclear_charge(odeint%potential) 
     465    vp0 = v(odeint%potential, r0, odeint%qn) + z/r0 
     466    e = odeint%e 
     467    qn = odeint%qn 
     468 
     469    select case (odeint%ode) 
    439470    case (ODE_SCHRODINGER) 
    440471      s = qn%l 
    … …  
    508539    case (ODE_DIRAC_POL2) 
    509540      clm = -sqrt( (M_TWO*qn%l + M_ONE)**2 - M_FOUR*qn%m**2  )/(M_TWO*qn%l + M_ONE) 
    510       bxc_0 = bxc(integrator%potential, r0, qn) 
     541      bxc_0 = bxc(odeint%potential, r0, qn) 
    511542 
    512543      !Set 1 
    … …  
    540571    end select 
    541572 
    542     if (integrator%ode /= ODE_DIRAC_POL2) then 
     573    if (odeint%ode /= ODE_DIRAC_POL2) then 
    543574      f(1) = r0**(s - M_ONE)*(a0 + a1*r0 + a2*r0**2 + a3*r0**3) 
    544575      f(2) = r0**(s - M_ONE)*(b0 + b1*r0 + b2*r0**2) 
    … …  
    546577 
    547578    ! MGGA term 
    548     if (integrator%ode == ODE_SCHRODINGER) then 
    549       f(2) = f(2)*(M_ONE + M_TWO*vtau(integrator%potential, r0, qn)) 
     579    if (odeint%ode == ODE_SCHRODINGER) then 
     580      f(2) = f(2)*(M_ONE + M_TWO*vtau(odeint%potential, r0, qn)) 
    550581    end if 
    551582 
    … …  
    553584  end subroutine ode_bc_origin 
    554585 
    555   subroutine ode_bc_infinity(integrator, rinf, f) 
     586  subroutine ode_bc_infinity(odeint, rinf, f) 
    556587    !-----------------------------------------------------------------------! 
    557588    ! Set the value of the functions at a point far awy from the nucleus.   ! 
    558589    !                                                                       ! 
    559     !  integrator - integrator object                                       ! 
    560     !  rinf       - point far away                                          ! 
    561     !  f          - function at rinf                                        ! 
    562     !-----------------------------------------------------------------------! 
    563     type(ode_integrator_t), intent(in)  :: integrator 
     590    !  odeint - integrator object                                           ! 
     591    !  rinf   - point far away                                              ! 
     592    !  f      - function at rinf                                            ! 
     593    !-----------------------------------------------------------------------! 
     594    type(ode_integrator_t), intent(in)  :: odeint 
    564595    real(R8),               intent(in)  :: rinf 
    565     real(R8),               intent(out) :: f(integrator%dim) 
     596    real(R8),               intent(out) :: f(odeint%dim) 
    566597 
    567598    real(R8) :: e, vinf, minf 
    … …  
    570601    call push_sub("ode_bc_infinity") 
    571602 
    572     e = integrator%e 
    573     qn = integrator%qn 
     603    e = odeint%e 
     604    qn = odeint%qn 
    574605 
    575606    !Values of the functions at infinity 
    576     select case (integrator%ode) 
     607    select case (odeint%ode) 
    577608    case (ODE_SCHRODINGER) 
    578609      f(1) = ODE_FINF 
    579       f(2) = -sqrt(-M_TWO*e)*f(1)*(M_ONE + M_TWO*vtau(integrator%potential, rinf, qn)) 
     610      f(2) = -sqrt(-M_TWO*e)*f(1)*(M_ONE + M_TWO*vtau(odeint%potential, rinf, qn)) 
    580611 
    581612    case (ODE_SCALAR_REL) 
    582613      f(1) = ODE_FINF 
    583       vinf = v(integrator%potential, rinf, qn) 
     614      vinf = v(odeint%potential, rinf, qn) 
    584615      minf = M_ONE + (e - vinf)/M_TWO/M_C2 
    585616      f(2) = -sqrt(-M_TWO*minf*e)/(M_TWO*minf*M_C)*f(1) 
    … …  
    618649  end subroutine ode_bc_infinity 
    619650 
    620   function ode_practical_infinity(integrator, ri, tol) 
     651  function ode_practical_infinity(odeint, ri, tol) 
    621652    !-----------------------------------------------------------------------! 
    622653    ! Returns the value of the practical infinity. The practical infinity   ! 
    … …  
    624655    ! magnitude of ODE_FINF.                                                ! 
    625656    !                                                                       ! 
    626     !  e         - energy                                                   ! 
    627     !  ode       - equation to solve                                        ! 
    628     !  tol       - tolerance                                                !  
    629     !  qn        -                                                          !   
    630     !  ri        -                                                          ! 
    631     !  potential -                                                          ! 
    632     !-----------------------------------------------------------------------! 
    633     type(ode_integrator_t), intent(in) :: integrator 
     657    !  odeint - integrator object                                           ! 
     658    !  tol    - tolerance                                                   !  
     659    !-----------------------------------------------------------------------! 
     660    type(ode_integrator_t), intent(in) :: odeint 
    634661    real(R8),               intent(in) :: ri 
    635662    real(R8),               intent(in) :: tol 
    … …  
    640667    call push_sub("ode_practical_infinity") 
    641668 
    642     e_max = min(integrator%e, -1.0e-12_r8) 
    643  
    644     select case (integrator%ode) 
     669    e_max = min(odeint%e, -1.0e-12_r8) 
     670 
     671    select case (odeint%ode) 
    645672    case (ODE_SCHRODINGER) 
    646673      ode_practical_infinity = -log(ODE_FINF)/sqrt(-M_TWO*e_max) 
    … …  
    652679      k = sqrt(-e_max)*sqrt(M_TWO + e_max/M_C2) 
    653680      x = ri 
    654       vinf = v(integrator%potential, x, integrator%qn)*x 
     681      vinf = v(odeint%potential, x, odeint%qn)*x 
    655682      h = -vinf*(e_max + M_C2)/(k*M_C2) - M_ONE 
    656683      f = h*log(x) - k*x - log(ODE_FINF) 
    … …  
    659686        d = d*M_HALF 
    660687        xm = x + d 
    661         vinf = v(integrator%potential, xm, integrator%qn)*xm 
     688        vinf = v(odeint%potential, xm, odeint%qn)*xm 
    662689        h = -vinf*(e_max+M_C2)/(k*M_C2) - 1.0  
    663690        if (f*(h*log(xm) - k*xm - log(ODE_FINF)) > M_ZERO) x = xm 
    … …  
    675702  end function ode_practical_infinity 
    676703 
    677   subroutine ode_derivatives(integrator, r, y, f) 
     704  subroutine ode_derivatives(odeint, r, y, f) 
    678705    !-----------------------------------------------------------------------! 
    679706    ! This subroutine is supposed to be called from the GSL ODE solver. It  ! 
    680     ! basically defines the system of equations to solve by returning the   ! 
    681     ! values of the derivatives of the functions to solve.                  ! 
     707    ! basically defines the system of equations to be solved by returning   ! 
     708    ! the values of the derivatives of the functions.                       ! 
    682709    !                                                                       ! 
    683710    ! The system of differential equations should be written in the         ! 
    … …  
    688715    !    d r                                                                ! 
    689716    !-----------------------------------------------------------------------! 
    690     type(ode_integrator_t), intent(in) :: integrator 
     717    type(ode_integrator_t), intent(in) :: odeint 
    691718    real(R8),               intent(in)  :: r 
    692     real(R8),               intent(in)  :: y(integrator%dim) 
    693     real(R8),               intent(out) :: f(integrator%dim) 
     719    real(R8),               intent(in)  :: y(odeint%dim) 
     720    real(R8),               intent(out) :: f(odeint%dim) 
    694721 
    695722    real(R8) :: k, v_r, dvdr_r, m, dmdr_r, mi, ci, c2i, ri, r2i, optvtau, llp1 
    … …  
    697724    real(R8) :: c11, c12, c33, c34, c21, c23, c22, c41, c43, c44 
    698725 
    699     k = real(integrator%qn%k, R8) 
    700     v_r = v(integrator%potential, r, integrator%qn) 
     726    k = real(odeint%qn%k, R8) 
     727    v_r = v(odeint%potential, r, odeint%qn) 
    701728    ri = M_ONE/r 
    702729    r2i = ri*ri 
    703730 
    704     select case (integrator%ode) 
     731    select case (odeint%ode) 
    705732    case (ODE_SCHRODINGER) 
    706733      ! Schrodinger equation 
    … …  
    715742      ! f_2 = - - y_2 + [(1 + 2 v_t)------ + 2(v_ks - e)] y_1 
    716743      !         r                    r^2       
    717       optvtau = M_ONE + M_TWO*vtau(integrator%potential, r, integrator%qn) 
    718       llp1 = real(integrator%qn%l, R8)*(real(integrator%qn%l, R8) + M_ONE) 
     744      optvtau = M_ONE + M_TWO*vtau(odeint%potential, r, odeint%qn) 
     745      llp1 = real(odeint%qn%l, R8)*(real(odeint%qn%l, R8) + M_ONE) 
    719746 
    720747      f(1) = y(2)/optvtau 
    721       f(2) = (M_TWO*(v_r - integrator%e) + optvtau*llp1*r2i)*y(1) - M_TWO*ri*y(2) 
     748      f(2) = (M_TWO*(v_r - odeint%e) + optvtau*llp1*r2i)*y(1) - M_TWO*ri*y(2) 
    722749 
    723750    case (ODE_SCALAR_REL)  
    724751      ! Scalar-relativistic equation 
    725       dvdr_r = dvdr(integrator%potential, r, integrator%qn) 
     752      dvdr_r = dvdr(odeint%potential, r, odeint%qn) 
    726753      ci = M_ONE/M_C  
    727754      c2i = ci*ci 
    728       m = M_ONE + (integrator%e - v_r)*M_HALF*c2i 
     755      m = M_ONE + (odeint%e - v_r)*M_HALF*c2i 
    729756      mi = M_ONE/m 
    730757      dmdr_r = -dvdr_r*M_HALF*c2i 
    731       llp1 = real(integrator%qn%l, R8)*(real(integrator%qn%l, R8) + M_ONE) 
     758      llp1 = real(odeint%qn%l, R8)*(real(odeint%qn%l, R8) + M_ONE) 
    732759 
    733760      f(1) = M_TWO*m*M_C*y(2) 
    734       f(2) = (llp1*M_HALF*mi*r2i + v_r - integrator%e)*ci*y(1) - M_TWO*ri*y(2) 
     761      f(2) = (llp1*M_HALF*mi*r2i + v_r - odeint%e)*ci*y(1) - M_TWO*ri*y(2) 
    735762 
    736763    case (ODE_DIRAC) 
    737764      ! Spin-unpolarized Dirac equation 
    738765      ci = M_ONE/M_C 
    739       f(1) = -(k + M_ONE)*ri*y(1) + (integrator%e + M_TWO*M_C2 - v_r)*ci*y(2) 
    740       f(2) = (v_r - integrator%e)*ci*y(1) + (k - M_ONE)*ri*y(2) 
     766      f(1) = -(k + M_ONE)*ri*y(1) + (odeint%e + M_TWO*M_C2 - v_r)*ci*y(2) 
     767      f(2) = (v_r - odeint%e)*ci*y(1) + (k - M_ONE)*ri*y(2) 
    741768 
    742769    case (ODE_DIRAC_POL1) 
    743770      ! Spin-polarized Dirac equation for m == 2l+1 
    744771      ci = M_ONE/M_C 
    745       bxc_r = bxc(integrator%potential, r, integrator%qn) 
    746  
    747       c11 = integrator%qn%l*ri 
    748       c12 = ci*(M_TWO*M_C2 - (v_r - integrator%e + M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l + M_THREE)*bxc_r)) 
    749       c21 = ci*(v_r - integrator%e + M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l + M_ONE)*bxc_r) 
    750       c22 = -(integrator%qn%l + M_TWO)*ri 
     772      bxc_r = bxc(odeint%potential, r, odeint%qn) 
     773 
     774      c11 = odeint%qn%l*ri 
     775      c12 = ci*(M_TWO*M_C2 - (v_r - odeint%e + M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l + M_THREE)*bxc_r)) 
     776      c21 = ci*(v_r - odeint%e + M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l + M_ONE)*bxc_r) 
     777      c22 = -(odeint%qn%l + M_TWO)*ri 
    751778 
    752779      f(1) = c11*y(1) + c12*y(2) 
    … …  
    756783      ! Spin-polarized Dirac equation for m /= 2l+1 
    757784      ci = M_ONE/M_C 
    758       bxc_r = bxc(integrator%potential, r, integrator%qn) 
    759       clm = -sqrt( (M_TWO*integrator%qn%l + M_ONE)**2 - M_FOUR*integrator%qn%m**2  )/(M_TWO*integrator%qn%l + M_ONE) 
    760  
    761       c11 = integrator%qn%l*ri 
    762       c12 = ci*(M_TWO*M_C2 - (v_r - integrator%e + M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l + M_THREE)*bxc_r)) 
    763       c21 = ci*(v_r - integrator%e + M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l + M_ONE)*bxc_r) 
    764       c22 = -(integrator%qn%l + M_TWO)*ri 
     785      bxc_r = bxc(odeint%potential, r, odeint%qn) 
     786      clm = -sqrt( (M_TWO*odeint%qn%l + M_ONE)**2 - M_FOUR*odeint%qn%m**2  )/(M_TWO*odeint%qn%l + M_ONE) 
     787 
     788      c11 = odeint%qn%l*ri 
     789      c12 = ci*(M_TWO*M_C2 - (v_r - odeint%e + M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l + M_THREE)*bxc_r)) 
     790      c21 = ci*(v_r - odeint%e + M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l + M_ONE)*bxc_r) 
     791      c22 = -(odeint%qn%l + M_TWO)*ri 
    765792      c23 = ci*clm*bxc_r 
    766       c33 = -(integrator%qn%l + M_ONE)*ri 
    767       c34 = ci*(M_TWO*M_C2 - (v_r - integrator%e - M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l - M_ONE) *bxc_r)) 
     793      c33 = -(odeint%qn%l + M_ONE)*ri 
     794      c34 = ci*(M_TWO*M_C2 - (v_r - odeint%e - M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l - M_ONE) *bxc_r)) 
    768795      c41 = c23 
    769       c43 = ci*(v_r - integrator%e - M_TWO*integrator%qn%m/(M_TWO*integrator%qn%l + M_ONE)*bxc_r) 
    770       c44 = (integrator%qn%l - M_ONE)*ri 
     796      c43 = ci*(v_r - odeint%e - M_TWO*odeint%qn%m/(M_TWO*odeint%qn%l + M_ONE)*bxc_r) 
     797      c44 = (odeint%qn%l - M_ONE)*ri 
    771798 
    772799      f(1) = c11*y(1) + c12*y(2) 
    … …  
    783810  end subroutine ode_derivatives 
    784811 
    785   subroutine ode_integrator_debug(integrator, ri, rf, nstep, r, f) 
    786     !-----------------------------------------------------------------------! 
    787     ! Prints debug information to the "debug_info/integrator" file.         ! 
    788     !                                                                       ! 
    789     !  integrator - integrator object                                       ! 
    790     !  ri         - initial integration point                               ! 
    791     !  rf         - final integration point                                 ! 
    792     !  nstep      - number of steps taken by the solver                     ! 
    793     !  r          - mesh used by the ODE solver                             ! 
    794     !  f          - integrated function                                     ! 
    795     !-----------------------------------------------------------------------! 
    796     type(ode_integrator_t), intent(in) :: integrator 
     812  subroutine ode_integrator_debug(odeint, ri, rf, nstep, r, f) 
     813    !-----------------------------------------------------------------------! 
     814    ! Prints debug information to the "debug_info/ode_integrator" file.     ! 
     815    !                                                                       ! 
     816    !  odeint - integrator object                                           ! 
     817    !  ri     - initial integration point                                   ! 
     818    !  rf     - final integration point                                     ! 
     819    !  nstep  - number of steps taken by the solver                         ! 
     820    !  r      - mesh used by the ODE solver                                 ! 
     821    !  f      - integrated function                                         ! 
     822    !-----------------------------------------------------------------------! 
     823    type(ode_integrator_t), intent(in) :: odeint 
    797824    real(R8),               intent(in) :: ri, rf 
    798825    integer,                intent(in) :: nstep 
    799     real(R8),               intent(in) :: r(integrator%nstepmax) 
    800     real(R8),               intent(in) :: f(integrator%nstepmax, integrator%dim) 
     826    real(R8),               intent(in) :: r(odeint%nstepmax) 
     827    real(R8),               intent(in) :: f(odeint%nstepmax, odeint%dim) 
    801828     
    802829    character(len=80) :: fmt 
    … …  
    806833 
    807834    call io_open(unit, file='debug_info/ode_integrator') 
    808     write(unit,'("ODE Integrator maximum number of steps: ",I6)') integrator%nstepmax 
     835    write(unit,'("ODE Integrator maximum number of steps: ",I6)') odeint%nstepmax 
    809836    write(unit,'("Quantum numbers:")') 
    810     write(unit,'("  l =  ",I2)') integrator%qn%l 
    811     write(unit,'("  k =  ",I2)') integrator%qn%k 
    812     write(unit,'("  s =  ",F4.1)') integrator%qn%s 
    813     write(unit,'("Energy:",F12.5)') integrator%e 
     837    write(unit,'("  l =  ",I2)') odeint%qn%l 
     838    write(unit,'("  k =  ",I2)') odeint%qn%k 
     839    write(unit,'("  s =  ",F4.1)') odeint%qn%s 
     840    write(unit,'("Energy:",F12.5)') odeint%e 
    814841    close(unit) 
    815842 
    … …  
    817844    write(unit,'("# Integration Starting Point: ",ES10.3E2)') ri 
    818845    write(unit,'("# Integration Ending Point:   ",ES10.3E2)') rf 
    819     write(fmt,'(A,I1,A)') "(ES14.8E2,1X,", integrator%dim,"(ES14.8E2,1X))" 
     846    write(fmt,'(A,I1,A)') "(ES14.8E2,1X,", odeint%dim,"(ES14.8E2,1X))" 
    820847    do i = 1, nstep 
    821       write(unit,fmt) r(i), f(i,1:integrator%dim) 
     848      write(unit,fmt) r(i), f(i,1:odeint%dim) 
    822849    end do 
    823850    close(unit) 

  • trunk/src/potentials.F90

    r725 r727  
    9494            potential_min, & 
    9595            potential_max, & 
    96             smallest_safe_r, & 
     96            potential_rmin, & 
     97            potential_rmax, & 
    9798            potential_kb_energy, & 
    9899            potential_output, & 
    … …  
    916917  end function potential_max 
    917918 
    918   function smallest_safe_r(potential) 
     919  function potential_rmin(potential) 
    919920    !-----------------------------------------------------------------------! 
    920921    ! Closest point to the origin that can be safely used.                  ! 
    921922    !-----------------------------------------------------------------------! 
    922923    type(potential_t), intent(in) :: potential 
    923     real(R8) :: smallest_safe_r 
    924  
    925     call push_sub("smallest_safe_r") 
    926  
    927     ASSERT(potential%type /= 0) 
    928  
    929     smallest_safe_r = potential%m%r(1) 
    930  
    931     call pop_sub() 
    932   end function smallest_safe_r 
     924    real(R8) :: potential_rmin 
     925 
     926    call push_sub("potential_rmin") 
     927 
     928    ASSERT(potential%type /= 0) 
     929 
     930    potential_rmin = potential%m%r(1) 
     931 
     932    call pop_sub() 
     933  end function potential_rmin 
     934 
     935  function potential_rmax(potential) 
     936    !-----------------------------------------------------------------------! 
     937    ! Larger point that can be safely used.                                 ! 
     938    !-----------------------------------------------------------------------! 
     939    type(potential_t), intent(in) :: potential 
     940    real(R8) :: potential_rmax 
     941 
     942    call push_sub("potential_rmax") 
     943 
     944    ASSERT(potential%type /= 0) 
     945 
     946    potential_rmax = potential%m%r(potential%m%np) 
     947 
     948    call pop_sub() 
     949  end function potential_rmax 
    933950 
    934951  function potential_kb_energy(potential, qn) 

  • trunk/src/wave_equations.F90

    r721 r727  
    255255    integer  :: i, nstep_out, nstep_in, ode 
    256256    real(R8) :: r0, ri, rinf 
    257     type(ode_integrator_t) :: ode_integrator 
     257    type(ode_integrator_t) :: odeint 
    258258    real(R8), pointer :: r_out(:), f_out(:,:) 
    259259    real(R8), pointer :: r_in(:),  f_in(:,:) 
    … …  
    263263    !Set the parameters needed to compute the derivatives of the functions 
    264264    ode = wave_equation_to_ode(wave_eq, qn) 
    265     call ode_integrator_null(ode_integrator) 
    266     call ode_integrator_init(ode_integrator, ode, integrator%stepping_function, & 
     265    call ode_integrator_null(odeint) 
     266    call ode_integrator_init(odeint, ode, integrator%stepping_function, & 
    267267                             integrator%nstepmax, integrator%tol, qn, e, potential) 
    268268 
    269269    !Set initial, final, and intermidiate points 
    270     r0 = smallest_safe_r(potential) 
     270    r0 = potential_rmin(potential) 
    271271    ri = classical_turning_point(potential, e, qn) 
    272     rinf = ode_practical_infinity(ode_integrator, ri, integrator%tol) 
     272    rinf = ode_practical_infinity(odeint, ri, integrator%tol) 
    273273 
    274274    !Outward and inward integrations 
    275     call ode_integration(ode_integrator, r0,   ri, nstep_out, r_out, f_out) 
    276     call ode_integration(ode_integrator, rinf, ri, nstep_in,  r_in,  f_in) 
     275    call ode_integration(odeint, r0,   ri, nstep_out, r_out, f_out) 
     276    call ode_integration(odeint, rinf, ri, nstep_in,  r_in,  f_in) 
    277277 
    278278    !Calculate functions mismatch 
    279279    wave_equation_ld_diff = & 
    280          ode_function_mismatch(ode_integrator, ri, f_out(nstep_out,:), f_in(nstep_in,:)) 
     280         ode_function_mismatch(odeint, ri, f_out(nstep_out,:), f_in(nstep_in,:)) 
    281281 
    282282    !Unset the derivatives parameters 
    283     call ode_integrator_end(ode_integrator) 
     283    call ode_integrator_end(odeint) 
    284284 
    285285    !Calculate the number of nodes 
    … …  
    320320    integer  :: ode, wf_dim, i, nstep_out, nstep_in, nstep_tot, r_tot_max 
    321321    real(R8) :: r0, ri, rinf, factor 
    322     type(ode_integrator_t) :: ode_integrator 
     322    type(ode_integrator_t) :: odeint 
    323323    real(R8), pointer     :: r_out(:), f_out(:,:) 
    324324    real(R8), pointer     :: r_in(:),  f_in(:,:) 
    … …  
    329329    !Set the parameters needed to compute the derivatives of the functions 
    330330    ode = wave_equation_to_ode(wave_eq, qn) 
    331     call ode_integrator_null(ode_integrator) 
    332     call ode_integrator_init(ode_integrator, ode, integrator%stepping_function, & 
     331    call ode_integrator_null(odeint) 
     332    call ode_integrator_init(odeint, ode, integrator%stepping_function, & 
    333333                             integrator%nstepmax, integrator%tol, qn, e, potential) 
    334334 
    335335    !Set initial, final, and intermidiate points 
    336     r0 = smallest_safe_r(potential) 
     336    r0 = potential_rmin(potential) 
    337337    ri = classical_turning_point(potential, e, qn) 
    338     rinf = ode_practical_infinity(ode_integrator, ri, integrator%tol) 
     338    rinf = ode_practical_infinity(odeint, ri, integrator%tol) 
    339339 
    340340    !Outward and inward integrations 
    341     call ode_integration(ode_integrator, r0,   ri, nstep_out, r_out, f_out) 
    342     call ode_integration(ode_integrator, rinf, ri, nstep_in,  r_in,  f_in) 
     341    call ode_integration(odeint, r0,   ri, nstep_out, r_out, f_out) 
     342    call ode_integration(odeint, rinf, ri, nstep_in,  r_in,  f_in) 
    343343 
    344344    !Match inward and outward functions at the classical turning point 
    345     call ode_match_functions(ode_integrator, nstep_out, nstep_in, f_out, f_in) 
     345    call ode_match_functions(odeint, nstep_out, nstep_in, f_out, f_in) 
    346346     
    347347    !Get the inward and outward wavefunctions on a single array 
    … …  
    351351    do i = 1, nstep_out 
    352352      r_tot(i) = r_out(i) 
    353       call ode_function_to_wavefunction(ode_integrator, & 
     353      call ode_function_to_wavefunction(odeint, & 
    354354                             r_tot(i), f_out(i,:), wf_tot(i,:), wfp_tot(i,:)) 
    355355    end do 
    356356    do i = 1,nstep_in-1 
    357357      r_tot(nstep_out+i) = r_in(nstep_in-i) 
    358       call ode_function_to_wavefunction(ode_integrator, & 
     358      call ode_function_to_wavefunction(odeint, & 
    359359                             r_tot(nstep_out+i), f_in(nstep_in-i,:), & 
    360360                             wf_tot(nstep_out+i,:), wfp_tot(nstep_out+i,:)) 
    … …  
    362362 
    363363    !Unset the derivatives parameters 
    364     call ode_integrator_end(ode_integrator) 
     364    call ode_integrator_end(odeint) 
    365365 
    366366    !Get the wavefunctions in the correct mesh 
    … …  
    427427    integer :: ode, i, nstep, nnodes 
    428428    real(R8) :: e 
    429     type(ode_integrator_t) :: ode_integrator 
     429    type(ode_integrator_t) :: odeint 
    430430    real(R8), pointer :: r(:), f(:,:) 
    431431 
    … …  
    436436    !Set the parameters needed to compute the derivatives of the functions 
    437437    ode = wave_equation_to_ode(wave_eq, qn) 
    438     call ode_integrator_null(ode_integrator) 
    439     call ode_integrator_init(ode_integrator, ode, integrator%stepping_function, & 
     438    call ode_integrator_null(odeint) 
     439    call ode_integrator_init(odeint, ode, integrator%stepping_function, & 
    440440                             integrator%nstepmax, integrator%tol, qn, e, potential) 
    441441 
    442442    !Outward integration 
    443     call ode_integration(ode_integrator, m%r(1), m%r(m%np), nstep, r, f) 
     443    call ode_integration(odeint, m%r(1), m%r(m%np), nstep, r, f) 
    444444 
    445445    !Unset the derivatives parameters 
    446     call ode_integrator_end(ode_integrator) 
     446    call ode_integrator_end(odeint) 
    447447 
    448448    !Calculate the number of nodes 
    … …  
    484484    integer :: ode, i, nstep, wf_dim 
    485485    real(R8) :: rc 
    486     type(ode_integrator_t) :: ode_integrator 
     486    type(ode_integrator_t) :: odeint 
    487487    real(R8), allocatable :: wf_rc(:), wfp_rc(:) 
    488488    real(R8), pointer :: r(:), f(:,:) 
    … …  
    492492    !Set the parameters needed to compute the derivatives of the functions 
    493493    ode = wave_equation_to_ode(wave_eq, qn) 
    494     call ode_integrator_null(ode_integrator) 
    495     call ode_integrator_init(ode_integrator, ode, integrator%stepping_function, & 
     494    call ode_integrator_null(odeint) 
     495    call ode_integrator_init(odeint, ode, integrator%stepping_function, & 
    496496                             integrator%nstepmax, integrator%tol, qn, e, potential) 
    497497 
    498498    rc = m%r(m%np) 
    499     call ode_integration(ode_integrator, smallest_safe_r(potential), rc, nstep, r, f) 
     499    call ode_integration(odeint, potential_rmin(potential), rc, nstep, r, f) 
    500500 
    501501    !Calculate the logaritmic derivative at rc 
    502502    wf_dim = qn_wf_dim(qn) 
    503503    allocate(wf_rc(wf_dim), wfp_rc(wf_dim)) 
    504     call ode_function_to_wavefunction(ode_integrator, rc, f(nstep,:), wf_rc, wfp_rc) 
     504    call ode_function_to_wavefunction(odeint, rc, f(nstep,:), wf_rc, wfp_rc) 
    505505    hamann_ld = wfp_rc(1)/wf_rc(1) 
    506506    deallocate(wf_rc, wfp_rc) 
    507507 
    508508    !Unset the derivatives parameters 
    509     call ode_integrator_end(ode_integrator) 
     509    call ode_integrator_end(odeint) 
    510510 
    511511    !Calculate the number of nodes 
    … …  
    546546    integer :: ode, i, nstep_out, wf_dim 
    547547    real(R8) :: factor 
    548     type(ode_integrator_t) :: ode_integrator 
     548    type(ode_integrator_t) :: odeint 
    549549    real(R8), allocatable :: wf_out(:,:), wfp_out(:,:) 
    550550    real(R8), pointer :: r_out(:), f_out(:,:) 
    … …  
    554554    !Set the parameters needed to compute the derivatives of the functions 
    555555    ode = wave_equation_to_ode(wave_eq, qn) 
    556     call ode_integrator_null(ode_integrator) 
    557     call ode_integrator_init(ode_integrator, ode, integrator%stepping_function, & 
     556    call ode_integrator_null(odeint) 
     557    call ode_integrator_init(odeint, ode, integrator%stepping_function, & 
    558558                             integrator%nstepmax, integrator%tol, qn, e, potential) 
    559559 
    560560    !Outward integration 
    561     call ode_integration(ode_integrator, smallest_safe_r(potential), m%r(m%np), nstep_out, r_out, f_out) 
     561    call ode_integration(odeint, potential_rmin(potential), m%r(m%np), nstep_out, r_out, f_out) 
    562562 
    563563    wf_dim = qn_wf_dim(qn) 
    564564    allocate(wf_out(nstep_out, wf_dim), wfp_out(nstep_out, wf_dim)) 
    565565    do i = 1, nstep_out 
    566       call ode_function_to_wavefunction(ode_integrator, r_out(i), & 
     566      call ode_function_to_wavefunction(odeint, r_out(i), & 
    567567                                  f_out(i,:), wf_out(i,:), wfp_out(i,:))       
    568568    end do 
    569569 
    570570    !Unset the derivatives parameters 
    571     call ode_integrator_end(ode_integrator) 
     571    call ode_integrator_end(odeint) 
    572572 
    573573    ! 
    … …  
    629629  function wave_function_cutoff(integrator) result(cutoff) 
    630630    !-----------------------------------------------------------------------! 
     631    ! We consider the wave-function to be zero if it is smaller than cutoff ! 
    631632    !-----------------------------------------------------------------------! 
    632633    type(integrator_t), intent(in) :: integrator 
    … …  
    639640  function wave_equation_to_ode(wave_eq, qn) 
    640641    !-----------------------------------------------------------------------! 
     642    ! Returns what is the set of ODE's that should be solve for a given     ! 
     643    ! wave-equation.                                                        ! 
    641644    !-----------------------------------------------------------------------! 
    642645    integer,    intent(in) :: wave_eq 
Note: See TracChangeset for help on using the changeset viewer.