| ! { dg-do run } |
| ! { dg-add-options ieee } |
| ! |
| ! Solve a diffusion problem using an object-oriented approach |
| ! |
| ! Author: Arjen Markus (comp.lang.fortran) |
| ! This version: pault@gcc.gnu.org |
| ! |
| ! Note: |
| ! (i) This could be turned into a more sophisticated program |
| ! using the techniques described in the chapter on |
| ! mathematical abstractions. |
| ! (That would allow the selection of the time integration |
| ! method in a transparent way) |
| ! |
| ! (ii) The target procedures for process_p and source_p are |
| ! different to the typebound procedures for dynamic types |
| ! because the passed argument is not type(base_pde_object). |
| ! |
| ! (iii) Two solutions are calculated, one with the procedure |
| ! pointers and the other with typebound procedures. The sums |
| ! of the solutions are compared. |
| |
| ! (iv) The source is a delta function in the middle of the |
| ! mesh, whilst the process is quartic in the local value, |
| ! when it is positive. |
| ! |
| ! base_pde_objects -- |
| ! Module to define the basic objects |
| ! |
| module base_pde_objects |
| implicit none |
| type, abstract :: base_pde_object |
| ! No data |
| procedure(process_p), pointer, pass :: process_p |
| procedure(source_p), pointer, pass :: source_p |
| contains |
| procedure(process), deferred :: process |
| procedure(source), deferred :: source |
| procedure :: initialise |
| procedure :: nabla2 |
| procedure :: print |
| procedure(real_times_obj), pass(obj), deferred :: real_times_obj |
| procedure(obj_plus_obj), deferred :: obj_plus_obj |
| procedure(obj_assign_obj), deferred :: obj_assign_obj |
| generic :: operator(*) => real_times_obj |
| generic :: operator(+) => obj_plus_obj |
| generic :: assignment(=) => obj_assign_obj |
| end type |
| abstract interface |
| function process_p (obj) |
| import base_pde_object |
| class(base_pde_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: process_p |
| end function process_p |
| end interface |
| abstract interface |
| function source_p (obj, time) |
| import base_pde_object |
| class(base_pde_object), intent(in) :: obj |
| real, intent(in) :: time |
| class(base_pde_object), allocatable :: source_p |
| end function source_p |
| end interface |
| abstract interface |
| function process (obj) |
| import base_pde_object |
| class(base_pde_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: process |
| end function process |
| end interface |
| abstract interface |
| function source (obj, time) |
| import base_pde_object |
| class(base_pde_object), intent(in) :: obj |
| real, intent(in) :: time |
| class(base_pde_object), allocatable :: source |
| end function source |
| end interface |
| abstract interface |
| function real_times_obj (factor, obj) result(newobj) |
| import base_pde_object |
| real, intent(in) :: factor |
| class(base_pde_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: newobj |
| end function real_times_obj |
| end interface |
| abstract interface |
| function obj_plus_obj (obj1, obj2) result(newobj) |
| import base_pde_object |
| class(base_pde_object), intent(in) :: obj1 |
| class(base_pde_object), intent(in) :: obj2 |
| class(base_pde_object), allocatable :: newobj |
| end function obj_plus_obj |
| end interface |
| abstract interface |
| subroutine obj_assign_obj (obj1, obj2) |
| import base_pde_object |
| class(base_pde_object), intent(inout) :: obj1 |
| class(base_pde_object), intent(in) :: obj2 |
| end subroutine obj_assign_obj |
| end interface |
| contains |
| ! print -- |
| ! Print the concentration field |
| subroutine print (obj) |
| class(base_pde_object) :: obj |
| ! Dummy |
| end subroutine print |
| ! initialise -- |
| ! Initialise the concentration field using a specific function |
| subroutine initialise (obj, funcxy) |
| class(base_pde_object) :: obj |
| interface |
| real function funcxy (coords) |
| real, dimension(:), intent(in) :: coords |
| end function funcxy |
| end interface |
| ! Dummy |
| end subroutine initialise |
| ! nabla2 -- |
| ! Determine the divergence |
| function nabla2 (obj) |
| class(base_pde_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: nabla2 |
| ! Dummy |
| end function nabla2 |
| end module base_pde_objects |
| ! cartesian_2d_objects -- |
| ! PDE object on a 2D cartesian grid |
| ! |
| module cartesian_2d_objects |
| use base_pde_objects |
| implicit none |
| type, extends(base_pde_object) :: cartesian_2d_object |
| real, dimension(:,:), allocatable :: c |
| real :: dx |
| real :: dy |
| contains |
| procedure :: process => process_cart2d |
| procedure :: source => source_cart2d |
| procedure :: initialise => initialise_cart2d |
| procedure :: nabla2 => nabla2_cart2d |
| procedure :: print => print_cart2d |
| procedure, pass(obj) :: real_times_obj => real_times_cart2d |
| procedure :: obj_plus_obj => obj_plus_cart2d |
| procedure :: obj_assign_obj => obj_assign_cart2d |
| end type cartesian_2d_object |
| interface grid_definition |
| module procedure grid_definition_cart2d |
| end interface |
| contains |
| function process_cart2d (obj) |
| class(cartesian_2d_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: process_cart2d |
| allocate (process_cart2d,source = obj) |
| select type (process_cart2d) |
| type is (cartesian_2d_object) |
| process_cart2d%c = -sign (obj%c, 1.0)*obj%c** 4 |
| class default |
| STOP 1 |
| end select |
| end function process_cart2d |
| function process_cart2d_p (obj) |
| class(base_pde_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: process_cart2d_p |
| allocate (process_cart2d_p,source = obj) |
| select type (process_cart2d_p) |
| type is (cartesian_2d_object) |
| select type (obj) |
| type is (cartesian_2d_object) |
| process_cart2d_p%c = -sign (obj%c, 1.0)*obj%c** 4 |
| end select |
| class default |
| STOP 2 |
| end select |
| end function process_cart2d_p |
| function source_cart2d (obj, time) |
| class(cartesian_2d_object), intent(in) :: obj |
| real, intent(in) :: time |
| class(base_pde_object), allocatable :: source_cart2d |
| integer :: m, n |
| m = size (obj%c, 1) |
| n = size (obj%c, 2) |
| allocate (source_cart2d, source = obj) |
| select type (source_cart2d) |
| type is (cartesian_2d_object) |
| if (allocated (source_cart2d%c)) deallocate (source_cart2d%c) |
| allocate (source_cart2d%c(m, n)) |
| source_cart2d%c = 0.0 |
| if (time .lt. 5.0) source_cart2d%c(m/2, n/2) = 0.1 |
| class default |
| STOP 3 |
| end select |
| end function source_cart2d |
| |
| function source_cart2d_p (obj, time) |
| class(base_pde_object), intent(in) :: obj |
| real, intent(in) :: time |
| class(base_pde_object), allocatable :: source_cart2d_p |
| integer :: m, n |
| select type (obj) |
| type is (cartesian_2d_object) |
| m = size (obj%c, 1) |
| n = size (obj%c, 2) |
| class default |
| STOP 4 |
| end select |
| allocate (source_cart2d_p,source = obj) |
| select type (source_cart2d_p) |
| type is (cartesian_2d_object) |
| if (allocated (source_cart2d_p%c)) deallocate (source_cart2d_p%c) |
| allocate (source_cart2d_p%c(m,n)) |
| source_cart2d_p%c = 0.0 |
| if (time .lt. 5.0) source_cart2d_p%c(m/2, n/2) = 0.1 |
| class default |
| STOP 5 |
| end select |
| end function source_cart2d_p |
| |
| ! grid_definition -- |
| ! Initialises the grid |
| ! |
| subroutine grid_definition_cart2d (obj, sizes, dims) |
| class(base_pde_object), allocatable :: obj |
| real, dimension(:) :: sizes |
| integer, dimension(:) :: dims |
| allocate( cartesian_2d_object :: obj ) |
| select type (obj) |
| type is (cartesian_2d_object) |
| allocate (obj%c(dims(1), dims(2))) |
| obj%c = 0.0 |
| obj%dx = sizes(1)/dims(1) |
| obj%dy = sizes(2)/dims(2) |
| class default |
| STOP 6 |
| end select |
| end subroutine grid_definition_cart2d |
| ! print_cart2d -- |
| ! Print the concentration field to the screen |
| ! |
| subroutine print_cart2d (obj) |
| class(cartesian_2d_object) :: obj |
| character(len=20) :: format |
| write( format, '(a,i0,a)' ) '(', size(obj%c,1), 'f6.3)' |
| write( *, format ) obj%c |
| end subroutine print_cart2d |
| ! initialise_cart2d -- |
| ! Initialise the concentration field using a specific function |
| ! |
| subroutine initialise_cart2d (obj, funcxy) |
| class(cartesian_2d_object) :: obj |
| interface |
| real function funcxy (coords) |
| real, dimension(:), intent(in) :: coords |
| end function funcxy |
| end interface |
| integer :: i, j |
| real, dimension(2) :: x |
| obj%c = 0.0 |
| do j = 2,size (obj%c, 2)-1 |
| x(2) = obj%dy * (j-1) |
| do i = 2,size (obj%c, 1)-1 |
| x(1) = obj%dx * (i-1) |
| obj%c(i,j) = funcxy (x) |
| enddo |
| enddo |
| end subroutine initialise_cart2d |
| ! nabla2_cart2d |
| ! Determine the divergence |
| function nabla2_cart2d (obj) |
| class(cartesian_2d_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: nabla2_cart2d |
| integer :: m, n |
| real :: dx, dy |
| m = size (obj%c, 1) |
| n = size (obj%c, 2) |
| dx = obj%dx |
| dy = obj%dy |
| allocate (cartesian_2d_object :: nabla2_cart2d) |
| select type (nabla2_cart2d) |
| type is (cartesian_2d_object) |
| allocate (nabla2_cart2d%c(m,n)) |
| nabla2_cart2d%c = 0.0 |
| nabla2_cart2d%c(2:m-1,2:n-1) = & |
| -(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 & |
| -(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 |
| class default |
| STOP 7 |
| end select |
| end function nabla2_cart2d |
| function real_times_cart2d (factor, obj) result(newobj) |
| real, intent(in) :: factor |
| class(cartesian_2d_object), intent(in) :: obj |
| class(base_pde_object), allocatable :: newobj |
| integer :: m, n |
| m = size (obj%c, 1) |
| n = size (obj%c, 2) |
| allocate (cartesian_2d_object :: newobj) |
| select type (newobj) |
| type is (cartesian_2d_object) |
| allocate (newobj%c(m,n)) |
| newobj%c = factor * obj%c |
| class default |
| STOP 8 |
| end select |
| end function real_times_cart2d |
| function obj_plus_cart2d (obj1, obj2) result( newobj ) |
| class(cartesian_2d_object), intent(in) :: obj1 |
| class(base_pde_object), intent(in) :: obj2 |
| class(base_pde_object), allocatable :: newobj |
| integer :: m, n |
| m = size (obj1%c, 1) |
| n = size (obj1%c, 2) |
| allocate (cartesian_2d_object :: newobj) |
| select type (newobj) |
| type is (cartesian_2d_object) |
| allocate (newobj%c(m,n)) |
| select type (obj2) |
| type is (cartesian_2d_object) |
| newobj%c = obj1%c + obj2%c |
| class default |
| STOP 9 |
| end select |
| class default |
| STOP 10 |
| end select |
| end function obj_plus_cart2d |
| subroutine obj_assign_cart2d (obj1, obj2) |
| class(cartesian_2d_object), intent(inout) :: obj1 |
| class(base_pde_object), intent(in) :: obj2 |
| select type (obj2) |
| type is (cartesian_2d_object) |
| obj1%c = obj2%c |
| class default |
| STOP 11 |
| end select |
| end subroutine obj_assign_cart2d |
| end module cartesian_2d_objects |
| ! define_pde_objects -- |
| ! Module to bring all the PDE object types together |
| ! |
| module define_pde_objects |
| use base_pde_objects |
| use cartesian_2d_objects |
| implicit none |
| interface grid_definition |
| module procedure grid_definition_general |
| end interface |
| contains |
| subroutine grid_definition_general (obj, type, sizes, dims) |
| class(base_pde_object), allocatable :: obj |
| character(len=*) :: type |
| real, dimension(:) :: sizes |
| integer, dimension(:) :: dims |
| select case (type) |
| case ("cartesian 2d") |
| call grid_definition (obj, sizes, dims) |
| case default |
| write(*,*) 'Unknown grid type: ', trim (type) |
| stop |
| end select |
| end subroutine grid_definition_general |
| end module define_pde_objects |
| ! pde_specific -- |
| ! Module holding the routines specific to the PDE that |
| ! we are solving |
| ! |
| module pde_specific |
| implicit none |
| contains |
| real function patch (coords) |
| real, dimension(:), intent(in) :: coords |
| if (sum ((coords-[50.0,50.0])**2) < 40.0) then |
| patch = 1.0 |
| else |
| patch = 0.0 |
| endif |
| end function patch |
| end module pde_specific |
| ! test_pde_solver -- |
| ! Small test program to demonstrate the usage |
| ! |
| program test_pde_solver |
| use define_pde_objects |
| use pde_specific |
| implicit none |
| class(base_pde_object), allocatable :: solution, deriv |
| integer :: i |
| real :: time, dtime, diff, chksum(2) |
| |
| call simulation1 ! Use proc pointers for source and process define_pde_objects |
| select type (solution) |
| type is (cartesian_2d_object) |
| deallocate (solution%c) |
| end select |
| select type (deriv) |
| type is (cartesian_2d_object) |
| deallocate (deriv%c) |
| end select |
| deallocate (solution, deriv) |
| |
| call simulation2 ! Use typebound procedures for source and process |
| if (chksum(1) .ne. chksum(2)) STOP 12 |
| if ((chksum(1) - 0.881868720)**2 > 1e-4) STOP 13 |
| contains |
| subroutine simulation1 |
| ! |
| ! Create the grid |
| ! |
| call grid_definition (solution, "cartesian 2d", [100.0, 100.0], [16, 16]) |
| call grid_definition (deriv, "cartesian 2d", [100.0, 100.0], [16, 16]) |
| ! |
| ! Initialise the concentration field |
| ! |
| call solution%initialise (patch) |
| ! |
| ! Set the procedure pointers |
| ! |
| solution%source_p => source_cart2d_p |
| solution%process_p => process_cart2d_p |
| ! |
| ! Perform the integration - explicit method |
| ! |
| time = 0.0 |
| dtime = 0.1 |
| diff = 5.0e-3 |
| |
| ! Give the diffusion coefficient correct dimensions. |
| select type (solution) |
| type is (cartesian_2d_object) |
| diff = diff * solution%dx * solution%dy / dtime |
| end select |
| |
| ! write(*,*) 'Time: ', time, diff |
| ! call solution%print |
| do i = 1,100 |
| deriv = solution%nabla2 () |
| solution = solution + diff * dtime * deriv + solution%source_p (time) + solution%process_p () |
| ! if ( mod(i, 25) == 0 ) then |
| ! write(*,*)'Time: ', time |
| ! call solution%print |
| ! endif |
| time = time + dtime |
| enddo |
| ! write(*,*) 'End result 1: ' |
| ! call solution%print |
| select type (solution) |
| type is (cartesian_2d_object) |
| chksum(1) = sum (solution%c) |
| end select |
| end subroutine |
| subroutine simulation2 |
| ! |
| ! Create the grid |
| ! |
| call grid_definition (solution, "cartesian 2d", [100.0, 100.0], [16, 16]) |
| call grid_definition (deriv, "cartesian 2d", [100.0, 100.0], [16, 16]) |
| ! |
| ! Initialise the concentration field |
| ! |
| call solution%initialise (patch) |
| ! |
| ! Set the procedure pointers |
| ! |
| solution%source_p => source_cart2d_p |
| solution%process_p => process_cart2d_p |
| ! |
| ! Perform the integration - explicit method |
| ! |
| time = 0.0 |
| dtime = 0.1 |
| diff = 5.0e-3 |
| |
| ! Give the diffusion coefficient correct dimensions. |
| select type (solution) |
| type is (cartesian_2d_object) |
| diff = diff * solution%dx * solution%dy / dtime |
| end select |
| |
| ! write(*,*) 'Time: ', time, diff |
| ! call solution%print |
| do i = 1,100 |
| deriv = solution%nabla2 () |
| solution = solution + diff * dtime * deriv + solution%source (time) + solution%process () |
| ! if ( mod(i, 25) == 0 ) then |
| ! write(*,*)'Time: ', time |
| ! call solution%print |
| ! endif |
| time = time + dtime |
| enddo |
| ! write(*,*) 'End result 2: ' |
| ! call solution%print |
| select type (solution) |
| type is (cartesian_2d_object) |
| chksum(2) = sum (solution%c) |
| end select |
| end subroutine |
| end program test_pde_solver |
