OSDN Git Service

* gfortran.dg/typebound_operator_8.f03: Use dg-add-options ieee.
[pf3gnuchains/gcc-fork.git] / gcc / testsuite / gfortran.dg / typebound_operator_8.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         call abort
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         call abort
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         call abort
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        call abort
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         call abort
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         call abort
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         call abort
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         call abort
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               call abort
324           end select
325       class default
326         call abort
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         call abort
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)) call abort
405   if ((chksum(1) - 0.881868720)**2 > 1e-4) call abort
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
500 ! { dg-final { cleanup-modules "pde_specific define_pde_objects cartesian_2d_objects base_pde_objects" } }