git
/
gcc
mirror of git://gcc.gnu.org/git/gcc.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Paul Thomas 2012-01-05 21:15:52 +00:00
parent f7d6ad0a5c
commit 003e0ad601
5 changed files with 627 additions and 497 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
gcc/fortran/ChangeLog
View File

@ -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

33
gcc/fortran/resolve.c
View File

@ -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;
}

24
gcc/testsuite/ChangeLog
View File

@ -1,21 +1,11 @@
2012-01-05 Jakub Jelinek <jakub@redhat.com>
2012-01-05 Paul Thomas <pault@gcc.gnu.org>
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.
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>

560
gcc/testsuite/gfortran.dg/typebound_operator_8.f03
View File

@ -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
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 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_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 process (obj)
import base_pde_object
class(base_pde_object), intent(in) :: obj
class(base_pde_object), allocatable :: process
end function process
subroutine assign_field(lhs,rhs)
import :: field
class(field) ,intent(OUT) :: lhs
class(field) ,intent(IN) :: rhs
end subroutine
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
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
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
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
select type (rhs)
type is (i_field); m = rhs%i + m
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
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
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
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
program main
use i_field_module
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)
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" } }

500
gcc/testsuite/gfortran.dg/typebound_operator_9.f03 Normal file
View File

@ -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" } }