diff options
author | Paul Thomas <pault@gcc.gnu.org> | 2012-01-05 21:15:52 +0000 |
---|---|---|
committer | Paul Thomas <pault@gcc.gnu.org> | 2012-01-05 21:15:52 +0000 |
commit | 003e0ad60130a4ba700a7b65e58ffcf0f051076c (patch) | |
tree | 31b9d65f981832104905adc6a62d92f7429f99fe | |
parent | f7d6ad0a5c2b3759f1952aa23bf5941013fec280 (diff) | |
download | gcc-003e0ad60130a4ba700a7b65e58ffcf0f051076c.zip gcc-003e0ad60130a4ba700a7b65e58ffcf0f051076c.tar.gz gcc-003e0ad60130a4ba700a7b65e58ffcf0f051076c.tar.bz2 |
PR fortran/PR48946
2012-01-05 Paul Thomas <pault@gcc.gnu.org>
PR fortran/PR48946
* resolve.c (resolve_typebound_static): If the typebound
procedure is 'deferred' try to find the correct specific
procedure in the derived type operator space itself.
2012-01-05 Paul Thomas <pault@gcc.gnu.org>
PR fortran/PR48946
* gfortran.dg/typebound_operator_9.f03: This is now a copy of
the old typebound_operator_8.f03.
* gfortran.dg/typebound_operator_8.f03: New version of
typebound_operator_7.f03 with 'u' a derived type instead of a
class object.
From-SVN: r182929
-rw-r--r-- | gcc/fortran/ChangeLog | 7 | ||||
-rw-r--r-- | gcc/fortran/resolve.c | 33 | ||||
-rw-r--r-- | gcc/testsuite/ChangeLog | 26 | ||||
-rw-r--r-- | gcc/testsuite/gfortran.dg/typebound_operator_8.f03 | 560 | ||||
-rw-r--r-- | gcc/testsuite/gfortran.dg/typebound_operator_9.f03 | 500 |
5 files changed, 628 insertions, 498 deletions
diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index 895d200..879c564 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,10 @@ +2012-01-05 Paul Thomas <pault@gcc.gnu.org> + + PR fortran/PR48946 + * resolve.c (resolve_typebound_static): If the typebound + procedure is 'deferred' try to find the correct specific + procedure in the derived type operator space itself. + 2012-01-04 Mikael Morin <mikael@gcc.gnu.org> PR fortran/50981 diff --git a/gcc/fortran/resolve.c b/gcc/fortran/resolve.c index 82045f8..79245ce 100644 --- a/gcc/fortran/resolve.c +++ b/gcc/fortran/resolve.c @@ -5614,6 +5614,39 @@ resolve_typebound_static (gfc_expr* e, gfc_symtree** target, e->ref = NULL; e->value.compcall.actual = NULL; + /* If we find a deferred typebound procedure, check for derived types + that an over-riding typebound procedure has not been missed. */ + if (e->value.compcall.tbp->deferred + && e->value.compcall.name + && !e->value.compcall.tbp->non_overridable + && e->value.compcall.base_object + && e->value.compcall.base_object->ts.type == BT_DERIVED) + { + gfc_symtree *st; + gfc_symbol *derived; + + /* Use the derived type of the base_object. */ + derived = e->value.compcall.base_object->ts.u.derived; + st = NULL; + + /* If necessary, go throught the inheritance chain. */ + while (!st && derived) + { + /* Look for the typebound procedure 'name'. */ + if (derived->f2k_derived && derived->f2k_derived->tb_sym_root) + st = gfc_find_symtree (derived->f2k_derived->tb_sym_root, + e->value.compcall.name); + if (!st) + derived = gfc_get_derived_super_type (derived); + } + + /* Now find the specific name in the derived type namespace. */ + if (st && st->n.tb && st->n.tb->u.specific) + gfc_find_sym_tree (st->n.tb->u.specific->name, + derived->ns, 1, &st); + if (st) + *target = st; + } return SUCCESS; } diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index 6db7569..3821d7c 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,21 +1,11 @@ -2012-01-05 Jakub Jelinek <jakub@redhat.com> - - PR debug/51762 - * gcc.dg/pr51762.c: New test. - - PR rtl-optimization/51767 - * gcc.c-torture/compile/pr51767.c: New test. - - PR middle-end/51768 - * c-c++-common/pr51768.c: New test. - - PR middle-end/44777 - * gcc.dg/tree-prof/pr44777.c: New test. - -2012-01-05 Jan Hubicka <jh@suse.cz> - - PR middle-end/49710 - * gcc.c-torture/compile/pr49710.c: New file. +2012-01-05 Paul Thomas <pault@gcc.gnu.org> + + PR fortran/PR48946 + * gfortran.dg/typebound_operator_9.f03: This is now a copy of + the old typebound_operator_8.f03. + * gfortran.dg/typebound_operator_8.f03: New version of + typebound_operator_7.f03 with 'u' a derived type instead of a + class object. 2012-01-05 Richard Guenther <rguenther@suse.de> diff --git a/gcc/testsuite/gfortran.dg/typebound_operator_8.f03 b/gcc/testsuite/gfortran.dg/typebound_operator_8.f03 index b27210b..a3726ba 100644 --- a/gcc/testsuite/gfortran.dg/typebound_operator_8.f03 +++ b/gcc/testsuite/gfortran.dg/typebound_operator_8.f03 @@ -1,500 +1,100 @@ ! { dg-do run } -! { dg-add-options ieee } +! PR48946 - complex expressions involving typebound operators of derived types. ! -! 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 +module field_module implicit none - type, abstract :: base_pde_object -! No data - procedure(process_p), pointer, pass :: process_p - procedure(source_p), pointer, pass :: source_p + type ,abstract :: field 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 + procedure(field_op_real) ,deferred :: multiply_real + procedure(field_plus_field) ,deferred :: plus + procedure(assign_field) ,deferred :: assn + generic :: operator(*) => multiply_real + generic :: operator(+) => plus + generic :: ASSIGNMENT(=) => assn 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 + function field_plus_field(lhs,rhs) + import :: field + class(field) ,intent(in) :: lhs + class(field) ,intent(in) :: rhs + class(field) ,allocatable :: field_plus_field + end function 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 + function field_op_real(lhs,rhs) + import :: field + class(field) ,intent(in) :: lhs + real ,intent(in) :: rhs + class(field) ,allocatable :: field_op_real + end function 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 + subroutine assign_field(lhs,rhs) + import :: field + class(field) ,intent(OUT) :: lhs + class(field) ,intent(IN) :: rhs + end subroutine 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 +end module + +module i_field_module + use field_module implicit none - type, extends(base_pde_object) :: cartesian_2d_object - real, dimension(:,:), allocatable :: c - real :: dx - real :: dy + type, extends (field) :: i_field + integer :: i 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 + procedure :: multiply_real => i_multiply_real + procedure :: plus => i_plus_i + procedure :: assn => i_assn + end type 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 - call abort - 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 - call abort + function i_plus_i(lhs,rhs) + class(i_field) ,intent(in) :: lhs + class(field) ,intent(in) :: rhs + class(field) ,allocatable :: i_plus_i + integer :: m = 0 + select type (lhs) + type is (i_field); m = lhs%i 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 - call abort + select type (rhs) + type is (i_field); m = rhs%i + m 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 - call abort - 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 - call abort + allocate (i_plus_i, source = i_field (m)) + end function + function i_multiply_real(lhs,rhs) + class(i_field) ,intent(in) :: lhs + real ,intent(in) :: rhs + class(field) ,allocatable :: i_multiply_real + integer :: m = 0 + select type (lhs) + type is (i_field); m = lhs%i * int (rhs) end select - end function source_cart2d_p + allocate (i_multiply_real, source = i_field (m)) + end function + subroutine i_assn(lhs,rhs) + class(i_field) ,intent(OUT) :: lhs + class(field) ,intent(IN) :: rhs + select type (lhs) + type is (i_field) + select type (rhs) + type is (i_field) + lhs%i = rhs%i + end select + end select + end subroutine +end module -! 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 - call abort - 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 - call abort - 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 - call abort - 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 - call abort - end select - class default - call abort - 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 - call abort - 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 +program main + use i_field_module 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) + type(i_field) ,allocatable :: u + allocate (u, source = i_field (99)) - 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)) call abort - if ((chksum(1) - 0.881868720)**2 > 1e-4) call abort -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 -! { dg-final { cleanup-modules "pde_specific define_pde_objects cartesian_2d_objects base_pde_objects" } } + u = u*2. + u = (u*2.0*4.0) + u*4.0 + u = u%multiply_real (2.0)*4.0 + u = i_multiply_real (u, 2.0) * 4.0 + + if (u%i .ne. 152064) call abort +end program +! { dg-final { cleanup-modules "field_module i_field_module" } } diff --git a/gcc/testsuite/gfortran.dg/typebound_operator_9.f03 b/gcc/testsuite/gfortran.dg/typebound_operator_9.f03 new file mode 100644 index 0000000..b27210b --- /dev/null +++ b/gcc/testsuite/gfortran.dg/typebound_operator_9.f03 @@ -0,0 +1,500 @@ +! { 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + 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 + call abort + end select + class default + call abort + 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 + call abort + 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)) call abort + if ((chksum(1) - 0.881868720)**2 > 1e-4) call abort +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 +! { dg-final { cleanup-modules "pde_specific define_pde_objects cartesian_2d_objects base_pde_objects" } } |