gcc/testsuite/gfortran.dg/typebound_operator_9.f03

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