! { 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" } }
--- /dev/null
+! { 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" } }
projects / platform / upstream / gcc.git / commitdiff