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! 32module 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 104contains 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 129end module base_pde_objects 130! cartesian_2d_objects -- 131! PDE object on a 2D cartesian grid 132! 133module 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 153contains 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 call abort 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 call abort 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 call abort 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 call abort 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 call abort 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 call abort 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 call abort 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 call abort 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 call abort 325 end select 326 class default 327 call abort 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 call abort 338 end select 339 end subroutine obj_assign_cart2d 340end module cartesian_2d_objects 341! define_pde_objects -- 342! Module to bring all the PDE object types together 343! 344module 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 351contains 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 365end module define_pde_objects 366! pde_specific -- 367! Module holding the routines specific to the PDE that 368! we are solving 369! 370module pde_specific 371 implicit none 372contains 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 381end module pde_specific 382! test_pde_solver -- 383! Small test program to demonstrate the usage 384! 385program 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)) call abort 406 if ((chksum(1) - 0.881868720)**2 > 1e-4) call abort 407contains 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 500end program test_pde_solver 501