2 ! { dg-add-options ieee }
4 ! Solve a diffusion problem using an object-oriented approach
6 ! Author: Arjen Markus (comp.lang.fortran)
7 ! This version: pault@gcc.gnu.org
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)
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).
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.
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.
29 ! Module to define the basic objects
31 module base_pde_objects
33 type, abstract :: base_pde_object
35 procedure(process_p), pointer, pass :: process_p
36 procedure(source_p), pointer, pass :: source_p
38 procedure(process), deferred :: process
39 procedure(source), deferred :: source
40 procedure :: initialise
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
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
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
66 function process (obj)
67 import base_pde_object
68 class(base_pde_object), intent(in) :: obj
69 class(base_pde_object), allocatable :: process
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
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
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
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
105 ! Print the concentration field
106 subroutine print (obj)
107 class(base_pde_object) :: obj
111 ! Initialise the concentration field using a specific function
112 subroutine initialise (obj, funcxy)
113 class(base_pde_object) :: obj
115 real function funcxy (coords)
116 real, dimension(:), intent(in) :: coords
120 end subroutine initialise
122 ! Determine the divergence
123 function nabla2 (obj)
124 class(base_pde_object), intent(in) :: obj
125 class(base_pde_object), allocatable :: nabla2
128 end module base_pde_objects
129 ! cartesian_2d_objects --
130 ! PDE object on a 2D cartesian grid
132 module cartesian_2d_objects
135 type, extends(base_pde_object) :: cartesian_2d_object
136 real, dimension(:,:), allocatable :: c
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
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
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)
171 type is (cartesian_2d_object)
172 process_cart2d_p%c = -sign (obj%c, 1.0)*obj%c** 4
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
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
195 end function source_cart2d
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
203 type is (cartesian_2d_object)
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
219 end function source_cart2d_p
222 ! Initialises the grid
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 )
230 type is (cartesian_2d_object)
231 allocate (obj%c(dims(1), dims(2)))
233 obj%dx = sizes(1)/dims(1)
234 obj%dy = sizes(2)/dims(2)
238 end subroutine grid_definition_cart2d
240 ! Print the concentration field to the screen
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
251 subroutine initialise_cart2d (obj, funcxy)
252 class(cartesian_2d_object) :: obj
254 real function funcxy (coords)
255 real, dimension(:), intent(in) :: coords
259 real, dimension(2) :: x
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)
268 end subroutine initialise_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
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
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
299 allocate (cartesian_2d_object :: newobj)
301 type is (cartesian_2d_object)
302 allocate (newobj%c(m,n))
303 newobj%c = factor * obj%c
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
315 allocate (cartesian_2d_object :: newobj)
317 type is (cartesian_2d_object)
318 allocate (newobj%c(m,n))
320 type is (cartesian_2d_object)
321 newobj%c = obj1%c + obj2%c
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
333 type is (cartesian_2d_object)
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
343 module define_pde_objects
345 use cartesian_2d_objects
347 interface grid_definition
348 module procedure grid_definition_general
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
357 case ("cartesian 2d")
358 call grid_definition (obj, sizes, dims)
360 write(*,*) 'Unknown grid type: ', trim (type)
363 end subroutine grid_definition_general
364 end module define_pde_objects
366 ! Module holding the routines specific to the PDE that
372 real function patch (coords)
373 real, dimension(:), intent(in) :: coords
374 if (sum ((coords-[50.0,50.0])**2) < 40.0) then
380 end module pde_specific
382 ! Small test program to demonstrate the usage
384 program test_pde_solver
385 use define_pde_objects
388 class(base_pde_object), allocatable :: solution, deriv
390 real :: time, dtime, diff, chksum(2)
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)
398 type is (cartesian_2d_object)
401 deallocate (solution, deriv)
403 call simulation2 ! Use typebound procedures for source and process
404 if (chksum(1) .ne. chksum(2)) STOP 12
405 if ((chksum(1) - 0.881868720)**2 > 1e-4) STOP 13
407 subroutine simulation1
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])
414 ! Initialise the concentration field
416 call solution%initialise (patch)
418 ! Set the procedure pointers
420 solution%source_p => source_cart2d_p
421 solution%process_p => process_cart2d_p
423 ! Perform the integration - explicit method
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
435 ! write(*,*) 'Time: ', time, diff
436 ! call solution%print
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
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)
453 subroutine simulation2
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])
460 ! Initialise the concentration field
462 call solution%initialise (patch)
464 ! Set the procedure pointers
466 solution%source_p => source_cart2d_p
467 solution%process_p => process_cart2d_p
469 ! Perform the integration - explicit method
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
481 ! write(*,*) 'Time: ', time, diff
482 ! call solution%print
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
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)
499 end program test_pde_solver