2 ! { dg-add-options ieee }
3 ! { dg-skip-if "Too big for local store" { spu-*-* } }
5 ! Solve a diffusion problem using an object-oriented approach
7 ! Author: Arjen Markus (comp.lang.fortran)
8 ! This version: pault@gcc.gnu.org
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)
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).
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.
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.
30 ! Module to define the basic objects
32 module base_pde_objects
34 type, abstract :: base_pde_object
36 procedure(process_p), pointer, pass :: process_p
37 procedure(source_p), pointer, pass :: source_p
39 procedure(process), deferred :: process
40 procedure(source), deferred :: source
41 procedure :: initialise
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
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
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
67 function process (obj)
68 import base_pde_object
69 class(base_pde_object), intent(in) :: obj
70 class(base_pde_object), allocatable :: process
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
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
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
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
106 ! Print the concentration field
107 subroutine print (obj)
108 class(base_pde_object) :: obj
112 ! Initialise the concentration field using a specific function
113 subroutine initialise (obj, funcxy)
114 class(base_pde_object) :: obj
116 real function funcxy (coords)
117 real, dimension(:), intent(in) :: coords
121 end subroutine initialise
123 ! Determine the divergence
124 function nabla2 (obj)
125 class(base_pde_object), intent(in) :: obj
126 class(base_pde_object), allocatable :: nabla2
129 end module base_pde_objects
130 ! cartesian_2d_objects --
131 ! PDE object on a 2D cartesian grid
133 module cartesian_2d_objects
136 type, extends(base_pde_object) :: cartesian_2d_object
137 real, dimension(:,:), allocatable :: c
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
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
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)
172 type is (cartesian_2d_object)
173 process_cart2d_p%c = -sign (obj%c, 1.0)*obj%c** 4
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
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
196 end function source_cart2d
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
204 type is (cartesian_2d_object)
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
220 end function source_cart2d_p
223 ! Initialises the grid
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 )
231 type is (cartesian_2d_object)
232 allocate (obj%c(dims(1), dims(2)))
234 obj%dx = sizes(1)/dims(1)
235 obj%dy = sizes(2)/dims(2)
239 end subroutine grid_definition_cart2d
241 ! Print the concentration field to the screen
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
252 subroutine initialise_cart2d (obj, funcxy)
253 class(cartesian_2d_object) :: obj
255 real function funcxy (coords)
256 real, dimension(:), intent(in) :: coords
260 real, dimension(2) :: x
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)
269 end subroutine initialise_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
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
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
300 allocate (cartesian_2d_object :: newobj)
302 type is (cartesian_2d_object)
303 allocate (newobj%c(m,n))
304 newobj%c = factor * obj%c
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
316 allocate (cartesian_2d_object :: newobj)
318 type is (cartesian_2d_object)
319 allocate (newobj%c(m,n))
321 type is (cartesian_2d_object)
322 newobj%c = obj1%c + obj2%c
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
334 type is (cartesian_2d_object)
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
344 module define_pde_objects
346 use cartesian_2d_objects
348 interface grid_definition
349 module procedure grid_definition_general
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
358 case ("cartesian 2d")
359 call grid_definition (obj, sizes, dims)
361 write(*,*) 'Unknown grid type: ', trim (type)
364 end subroutine grid_definition_general
365 end module define_pde_objects
367 ! Module holding the routines specific to the PDE that
373 real function patch (coords)
374 real, dimension(:), intent(in) :: coords
375 if (sum ((coords-[50.0,50.0])**2) < 40.0) then
381 end module pde_specific
383 ! Small test program to demonstrate the usage
385 program test_pde_solver
386 use define_pde_objects
389 class(base_pde_object), allocatable :: solution, deriv
391 real :: time, dtime, diff, chksum(2)
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)
399 type is (cartesian_2d_object)
402 deallocate (solution, deriv)
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
408 subroutine simulation1
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])
415 ! Initialise the concentration field
417 call solution%initialise (patch)
419 ! Set the procedure pointers
421 solution%source_p => source_cart2d_p
422 solution%process_p => process_cart2d_p
424 ! Perform the integration - explicit method
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
436 ! write(*,*) 'Time: ', time, diff
437 ! call solution%print
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
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)
454 subroutine simulation2
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])
461 ! Initialise the concentration field
463 call solution%initialise (patch)
465 ! Set the procedure pointers
467 solution%source_p => source_cart2d_p
468 solution%process_p => process_cart2d_p
470 ! Perform the integration - explicit method
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
482 ! write(*,*) 'Time: ', time, diff
483 ! call solution%print
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
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)
500 end program test_pde_solver