OSDN Git Service

/cp
[pf3gnuchains/gcc-fork.git] / 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" } }