gcc/testsuite/gfortran.dg/typebound_operator_8.f03

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