gcc/testsuite/gfortran.dg/typebound_operator_9.f03

   1 ! { dg-do run }
   2 ! { dg-add-options ieee }
   3 !
   4 !     Solve a diffusion problem using an object-oriented approach
   5 !
   6 !     Author: Arjen Markus (comp.lang.fortran)
   7 !     This version: pault@gcc.gnu.org
   8 !
   9 !     Note:
  10 !     (i) This could be turned into a more sophisticated program
  11 !     using the techniques described in the chapter on
  12 !     mathematical abstractions.
  13 !     (That would allow the selection of the time integration
  14 !     method in a transparent way)
  15 !
  16 !     (ii) The target procedures for process_p and source_p are
  17 !     different to the typebound procedures for dynamic types
  18 !     because the passed argument is not type(base_pde_object).
  19 !
  20 !     (iii) Two solutions are calculated, one with the procedure
  21 !     pointers and the other with typebound procedures. The sums
  22 !     of the solutions are compared.
  23
  24 !     (iv) The source is a delta function in the middle of the
  25 !     mesh, whilst the process is quartic in the local value,
  26 !     when it is positive.
  27 !
  28 ! base_pde_objects --
  29 !     Module to define the basic objects
  30 !
  31 module base_pde_objects
  32   implicit none
  33   type, abstract :: base_pde_object
  34 ! No data
  35     procedure(process_p), pointer, pass :: process_p
  36     procedure(source_p), pointer, pass  :: source_p
  37   contains
  38     procedure(process), deferred :: process
  39     procedure(source), deferred :: source
  40     procedure :: initialise
  41     procedure :: nabla2
  42     procedure :: print
  43     procedure(real_times_obj), pass(obj), deferred :: real_times_obj
  44     procedure(obj_plus_obj),              deferred :: obj_plus_obj
  45     procedure(obj_assign_obj),            deferred :: obj_assign_obj
  46     generic :: operator(*)    => real_times_obj
  47     generic :: operator(+)    => obj_plus_obj
  48     generic :: assignment(=)  => obj_assign_obj
  49   end type
  50   abstract interface
  51     function process_p (obj)
  52       import base_pde_object
  53       class(base_pde_object), intent(in)  :: obj
  54       class(base_pde_object), allocatable :: process_p
  55     end function process_p
  56   end interface
  57   abstract interface
  58     function source_p (obj, time)
  59       import base_pde_object
  60       class(base_pde_object), intent(in)  :: obj
  61       real, intent(in)                    :: time
  62       class(base_pde_object), allocatable :: source_p
  63     end function source_p
  64   end interface
  65   abstract interface
  66     function process (obj)
  67       import base_pde_object
  68       class(base_pde_object), intent(in)  :: obj
  69       class(base_pde_object), allocatable :: process
  70     end function process
  71   end interface
  72   abstract interface
  73     function source (obj, time)
  74       import base_pde_object
  75       class(base_pde_object), intent(in)  :: obj
  76       real, intent(in)                    :: time
  77       class(base_pde_object), allocatable :: source
  78     end function source
  79   end interface
  80   abstract interface
  81     function real_times_obj (factor, obj) result(newobj)
  82       import base_pde_object
  83       real, intent(in)                    :: factor
  84       class(base_pde_object), intent(in)  :: obj
  85       class(base_pde_object), allocatable :: newobj
  86     end function real_times_obj
  87   end interface
  88   abstract interface
  89     function obj_plus_obj (obj1, obj2) result(newobj)
  90       import base_pde_object
  91       class(base_pde_object), intent(in)  :: obj1
  92       class(base_pde_object), intent(in)  :: obj2
  93       class(base_pde_object), allocatable :: newobj
  94     end function obj_plus_obj
  95   end interface
  96   abstract interface
  97     subroutine obj_assign_obj (obj1, obj2)
  98       import base_pde_object
  99       class(base_pde_object), intent(inout)  :: obj1
 100       class(base_pde_object), intent(in)     :: obj2
 101     end subroutine obj_assign_obj
 102   end interface
 103 contains
 104 ! print --
 105 !     Print the concentration field
 106   subroutine print (obj)
 107     class(base_pde_object) :: obj
 108     ! Dummy
 109   end subroutine print
 110 ! initialise --
 111 !     Initialise the concentration field using a specific function
 112   subroutine initialise (obj, funcxy)
 113     class(base_pde_object) :: obj
 114     interface
 115       real function funcxy (coords)
 116         real, dimension(:), intent(in) :: coords
 117       end function funcxy
 118     end interface
 119     ! Dummy
 120   end subroutine initialise
 121 ! nabla2 --
 122 !     Determine the divergence
 123   function nabla2 (obj)
 124     class(base_pde_object), intent(in)  :: obj
 125     class(base_pde_object), allocatable :: nabla2
 126     ! Dummy
 127   end function nabla2
 128 end module base_pde_objects
 129 ! cartesian_2d_objects --
 130 !     PDE object on a 2D cartesian grid
 131 !
 132 module cartesian_2d_objects
 133   use base_pde_objects
 134   implicit none
 135   type, extends(base_pde_object) :: cartesian_2d_object
 136     real, dimension(:,:), allocatable :: c
 137     real                              :: dx
 138     real                              :: dy
 139   contains
 140     procedure            :: process       => process_cart2d
 141     procedure            :: source         => source_cart2d
 142     procedure            :: initialise     => initialise_cart2d
 143     procedure            :: nabla2         => nabla2_cart2d
 144     procedure            :: print          => print_cart2d
 145     procedure, pass(obj) :: real_times_obj => real_times_cart2d
 146     procedure            :: obj_plus_obj   => obj_plus_cart2d
 147     procedure            :: obj_assign_obj => obj_assign_cart2d
 148   end type cartesian_2d_object
 149   interface grid_definition
 150     module procedure grid_definition_cart2d
 151   end interface
 152 contains
 153   function process_cart2d (obj)
 154     class(cartesian_2d_object), intent(in)  :: obj
 155     class(base_pde_object), allocatable :: process_cart2d
 156     allocate (process_cart2d,source = obj)
 157     select type (process_cart2d)
 158       type is (cartesian_2d_object)
 159         process_cart2d%c = -sign (obj%c, 1.0)*obj%c** 4
 160       class default
 161         STOP 1
 162     end select
 163   end function process_cart2d
 164   function process_cart2d_p (obj)
 165     class(base_pde_object), intent(in)  :: obj
 166     class(base_pde_object), allocatable :: process_cart2d_p
 167     allocate (process_cart2d_p,source = obj)
 168     select type (process_cart2d_p)
 169       type is (cartesian_2d_object)
 170         select type (obj)
 171           type is (cartesian_2d_object)
 172             process_cart2d_p%c = -sign (obj%c, 1.0)*obj%c** 4
 173         end select
 174       class default
 175         STOP 2
 176     end select
 177   end function process_cart2d_p
 178   function source_cart2d (obj, time)
 179     class(cartesian_2d_object), intent(in)  :: obj
 180     real, intent(in)                    :: time
 181     class(base_pde_object), allocatable :: source_cart2d
 182     integer :: m, n
 183     m = size (obj%c, 1)
 184     n = size (obj%c, 2)
 185     allocate (source_cart2d, source = obj)
 186     select type (source_cart2d)
 187       type is (cartesian_2d_object)
 188         if (allocated (source_cart2d%c)) deallocate (source_cart2d%c)
 189         allocate (source_cart2d%c(m, n))
 190         source_cart2d%c = 0.0
 191         if (time .lt. 5.0) source_cart2d%c(m/2, n/2) = 0.1
 192       class default
 193         STOP 3
 194     end select
 195   end function source_cart2d
 196
 197   function source_cart2d_p (obj, time)
 198     class(base_pde_object), intent(in)  :: obj
 199     real, intent(in)                    :: time
 200     class(base_pde_object), allocatable :: source_cart2d_p
 201     integer :: m, n
 202     select type (obj)
 203       type is (cartesian_2d_object)
 204         m = size (obj%c, 1)
 205         n = size (obj%c, 2)
 206       class default
 207        STOP 4
 208     end select
 209     allocate (source_cart2d_p,source = obj)
 210     select type (source_cart2d_p)
 211       type is (cartesian_2d_object)
 212         if (allocated (source_cart2d_p%c)) deallocate (source_cart2d_p%c)
 213         allocate (source_cart2d_p%c(m,n))
 214         source_cart2d_p%c = 0.0
 215         if (time .lt. 5.0) source_cart2d_p%c(m/2, n/2) = 0.1
 216       class default
 217         STOP 5
 218     end select
 219   end function source_cart2d_p
 220
 221 ! grid_definition --
 222 !     Initialises the grid
 223 !
 224   subroutine grid_definition_cart2d (obj, sizes, dims)
 225     class(base_pde_object), allocatable :: obj
 226     real, dimension(:)                  :: sizes
 227     integer, dimension(:)               :: dims
 228     allocate( cartesian_2d_object :: obj )
 229     select type (obj)
 230       type is (cartesian_2d_object)
 231         allocate (obj%c(dims(1), dims(2)))
 232         obj%c  = 0.0
 233         obj%dx = sizes(1)/dims(1)
 234         obj%dy = sizes(2)/dims(2)
 235       class default
 236         STOP 6
 237     end select
 238   end subroutine grid_definition_cart2d
 239 ! print_cart2d --
 240 !     Print the concentration field to the screen
 241 !
 242   subroutine print_cart2d (obj)
 243     class(cartesian_2d_object) :: obj
 244     character(len=20)          :: format
 245     write( format, '(a,i0,a)' ) '(', size(obj%c,1), 'f6.3)'
 246     write( *, format ) obj%c
 247   end subroutine print_cart2d
 248 ! initialise_cart2d --
 249 !     Initialise the concentration field using a specific function
 250 !
 251   subroutine initialise_cart2d (obj, funcxy)
 252     class(cartesian_2d_object) :: obj
 253     interface
 254       real function funcxy (coords)
 255         real, dimension(:), intent(in) :: coords
 256       end function funcxy
 257     end interface
 258     integer                    :: i, j
 259     real, dimension(2)         :: x
 260     obj%c = 0.0
 261     do j = 2,size (obj%c, 2)-1
 262       x(2) = obj%dy * (j-1)
 263       do i = 2,size (obj%c, 1)-1
 264         x(1) = obj%dx * (i-1)
 265         obj%c(i,j) = funcxy (x)
 266       enddo
 267     enddo
 268   end subroutine initialise_cart2d
 269 ! nabla2_cart2d
 270 !     Determine the divergence
 271   function nabla2_cart2d (obj)
 272     class(cartesian_2d_object), intent(in)  :: obj
 273     class(base_pde_object), allocatable     :: nabla2_cart2d
 274     integer                                 :: m, n
 275     real                                    :: dx, dy
 276     m = size (obj%c, 1)
 277     n = size (obj%c, 2)
 278     dx = obj%dx
 279     dy = obj%dy
 280     allocate (cartesian_2d_object :: nabla2_cart2d)
 281     select type (nabla2_cart2d)
 282       type is (cartesian_2d_object)
 283         allocate (nabla2_cart2d%c(m,n))
 284         nabla2_cart2d%c = 0.0
 285         nabla2_cart2d%c(2:m-1,2:n-1) = &
 286           -(2.0 * obj%c(2:m-1,2:n-1) - obj%c(1:m-2,2:n-1) - obj%c(3:m,2:n-1)) / dx**2 &
 287           -(2.0 * obj%c(2:m-1,2:n-1) - obj%c(2:m-1,1:n-2) - obj%c(2:m-1,3:n)) / dy**2
 288       class default
 289         STOP 7
 290     end select
 291   end function nabla2_cart2d
 292   function real_times_cart2d (factor, obj) result(newobj)
 293     real, intent(in)                        :: factor
 294     class(cartesian_2d_object), intent(in)  :: obj
 295     class(base_pde_object), allocatable     :: newobj
 296     integer                                 :: m, n
 297     m = size (obj%c, 1)
 298     n = size (obj%c, 2)
 299     allocate (cartesian_2d_object :: newobj)
 300     select type (newobj)
 301       type is (cartesian_2d_object)
 302         allocate (newobj%c(m,n))
 303         newobj%c = factor * obj%c
 304       class default
 305         STOP 8
 306     end select
 307   end function real_times_cart2d
 308   function obj_plus_cart2d (obj1, obj2) result( newobj )
 309     class(cartesian_2d_object), intent(in)  :: obj1
 310     class(base_pde_object), intent(in)      :: obj2
 311     class(base_pde_object), allocatable     :: newobj
 312     integer                                 :: m, n
 313     m = size (obj1%c, 1)
 314     n = size (obj1%c, 2)
 315     allocate (cartesian_2d_object :: newobj)
 316     select type (newobj)
 317       type is (cartesian_2d_object)
 318         allocate (newobj%c(m,n))
 319           select type (obj2)
 320             type is (cartesian_2d_object)
 321               newobj%c = obj1%c + obj2%c
 322             class default
 323               STOP 9
 324           end select
 325       class default
 326         STOP 10
 327     end select
 328   end function obj_plus_cart2d
 329   subroutine obj_assign_cart2d (obj1, obj2)
 330     class(cartesian_2d_object), intent(inout) :: obj1
 331     class(base_pde_object), intent(in)        :: obj2
 332     select type (obj2)
 333       type is (cartesian_2d_object)
 334         obj1%c = obj2%c
 335       class default
 336         STOP 11
 337     end select
 338   end subroutine obj_assign_cart2d
 339 end module cartesian_2d_objects
 340 ! define_pde_objects --
 341 !     Module to bring all the PDE object types together
 342 !
 343 module define_pde_objects
 344   use base_pde_objects
 345   use cartesian_2d_objects
 346   implicit none
 347   interface grid_definition
 348     module procedure grid_definition_general
 349   end interface
 350 contains
 351   subroutine grid_definition_general (obj, type, sizes, dims)
 352     class(base_pde_object), allocatable :: obj
 353     character(len=*)                    :: type
 354     real, dimension(:)                  :: sizes
 355     integer, dimension(:)               :: dims
 356     select case (type)
 357       case ("cartesian 2d")
 358         call grid_definition (obj, sizes, dims)
 359       case default
 360         write(*,*) 'Unknown grid type: ', trim (type)
 361         stop
 362     end select
 363   end subroutine grid_definition_general
 364 end module define_pde_objects
 365 ! pde_specific --
 366 !     Module holding the routines specific to the PDE that
 367 !     we are solving
 368 !
 369 module pde_specific
 370   implicit none
 371 contains
 372   real function patch (coords)
 373     real, dimension(:), intent(in) :: coords
 374     if (sum ((coords-[50.0,50.0])**2) < 40.0) then
 375       patch = 1.0
 376     else
 377       patch = 0.0
 378     endif
 379   end function patch
 380 end module pde_specific
 381 ! test_pde_solver --
 382 !     Small test program to demonstrate the usage
 383 !
 384 program test_pde_solver
 385   use define_pde_objects
 386   use pde_specific
 387   implicit none
 388   class(base_pde_object), allocatable :: solution, deriv
 389   integer                             :: i
 390   real                                :: time, dtime, diff, chksum(2)
 391
 392   call simulation1     ! Use proc pointers for source and process define_pde_objects
 393   select type (solution)
 394     type is (cartesian_2d_object)
 395       deallocate (solution%c)
 396   end select
 397   select type (deriv)
 398     type is (cartesian_2d_object)
 399       deallocate (deriv%c)
 400   end select
 401   deallocate (solution, deriv)
 402
 403   call simulation2     ! Use typebound procedures for source and process
 404   if (chksum(1) .ne. chksum(2)) STOP 12
 405   if ((chksum(1) - 0.881868720)**2 > 1e-4) STOP 13
 406 contains
 407   subroutine simulation1
 408 !
 409 ! Create the grid
 410 !
 411     call grid_definition (solution, "cartesian 2d", [100.0, 100.0], [16, 16])
 412     call grid_definition (deriv,    "cartesian 2d", [100.0, 100.0], [16, 16])
 413 !
 414 ! Initialise the concentration field
 415 !
 416     call solution%initialise (patch)
 417 !
 418 ! Set the procedure pointers
 419 !
 420     solution%source_p => source_cart2d_p
 421     solution%process_p => process_cart2d_p
 422 !
 423 ! Perform the integration - explicit method
 424 !
 425     time  = 0.0
 426     dtime = 0.1
 427     diff =  5.0e-3
 428
 429 ! Give the diffusion coefficient correct dimensions.
 430     select type (solution)
 431       type is (cartesian_2d_object)
 432         diff  = diff * solution%dx * solution%dy / dtime
 433     end select
 434
 435 !     write(*,*) 'Time: ', time, diff
 436 !     call solution%print
 437     do i = 1,100
 438       deriv    =  solution%nabla2 ()
 439       solution = solution + diff * dtime * deriv + solution%source_p (time) + solution%process_p ()
 440 !         if ( mod(i, 25) == 0 ) then
 441 !             write(*,*)'Time: ', time
 442 !             call solution%print
 443 !         endif
 444     time = time + dtime
 445     enddo
 446 !    write(*,*) 'End result 1: '
 447 !    call solution%print
 448     select type (solution)
 449       type is (cartesian_2d_object)
 450         chksum(1) = sum (solution%c)
 451     end select
 452   end subroutine
 453   subroutine simulation2
 454 !
 455 ! Create the grid
 456 !
 457     call grid_definition (solution, "cartesian 2d", [100.0, 100.0], [16, 16])
 458     call grid_definition (deriv,    "cartesian 2d", [100.0, 100.0], [16, 16])
 459 !
 460 ! Initialise the concentration field
 461 !
 462     call solution%initialise (patch)
 463 !
 464 ! Set the procedure pointers
 465 !
 466     solution%source_p => source_cart2d_p
 467     solution%process_p => process_cart2d_p
 468 !
 469 ! Perform the integration - explicit method
 470 !
 471     time  = 0.0
 472     dtime = 0.1
 473     diff =  5.0e-3
 474
 475 ! Give the diffusion coefficient correct dimensions.
 476     select type (solution)
 477       type is (cartesian_2d_object)
 478         diff  = diff * solution%dx * solution%dy / dtime
 479     end select
 480
 481 !     write(*,*) 'Time: ', time, diff
 482 !     call solution%print
 483     do i = 1,100
 484       deriv    =  solution%nabla2 ()
 485       solution = solution + diff * dtime * deriv + solution%source (time) + solution%process ()
 486 !         if ( mod(i, 25) == 0 ) then
 487 !             write(*,*)'Time: ', time
 488 !             call solution%print
 489 !         endif
 490       time = time + dtime
 491     enddo
 492 !    write(*,*) 'End result 2: '
 493 !    call solution%print
 494     select type (solution)
 495       type is (cartesian_2d_object)
 496         chksum(2) = sum (solution%c)
 497     end select
 498   end subroutine
 499 end program test_pde_solver