Purity refactor

Purity refactor:

> all math-operators return real array, all math-operators are **pure**
functions.

Why:

+ completely avoid memory leaks;
+ boost performance;
+ parallel-ready.

This change addresses the need by:

+ eliminate all polymorphic, non pure function results.

Side effects:

It could be not fully backward compatible, however it seems quite so.

Author: szaghi

src/lib/foodie_integrand_object.f90

@@ -11,13 +11,10 @@ module foodie_integrand_object
 
 type, abstract :: integrand_object
   !< Abstract type for building FOODIE ODE integrators.
-  integer(I_P) :: n=0_I_P !< ODE *space* dimension.
 #ifdef CAF
   class(*), allocatable :: dummy_to_allow_extensions[:] !< Dummy member to allow concrete extensions with coarray members.
 #endif
   contains
-    ! public methods
-    procedure, pass(self), public :: set_dimension !< Set ODE *space* dimension.
     ! public deferred procedures that concrete integrand-field must implement
     procedure(time_derivative), pass(self), deferred, public :: t !< Time derivative, residuals.
     ! operators
@@ -60,9 +57,9 @@ module foodie_integrand_object
   function time_derivative(self, t) result(dState_dt)
   !< Time derivative function of integrand class, i.e. the residuals function.
   import :: integrand_object, R_P
-  class(integrand_object), intent(in) :: self                !< Integrand field.
-  real(R_P),    optional,  intent(in) :: t                   !< Time.
-  real(R_P)                           :: dState_dt(1:self%n) !< Result of the time derivative function of integrand field.
+  class(integrand_object), intent(in) :: self         !< Integrand field.
+  real(R_P),    optional,  intent(in) :: t            !< Time.
+  real(R_P), allocatable              :: dState_dt(:) !< Result of the time derivative function of integrand field.
   endfunction time_derivative
 
   ! operators
@@ -77,41 +74,41 @@ function local_error_operator(lhs, rhs) result(error)
   pure function integrand_op_real(lhs, rhs) result(operator_result)
   !< Asymmetric type operator `integrand.op.real`.
   import :: integrand_object, R_P
-  class(integrand_object), intent(in) :: lhs                      !< Left hand side.
-  real(R_P),               intent(in) :: rhs(1:lhs%n)             !< Right hand side.
-  real(R_P)                           :: operator_result(1:lhs%n) !< Operator result.
+  class(integrand_object), intent(in) :: lhs                !< Left hand side.
+  real(R_P),               intent(in) :: rhs(1:)            !< Right hand side.
+  real(R_P), allocatable              :: operator_result(:) !< Operator result.
   endfunction integrand_op_real
 
   pure function real_op_integrand(lhs, rhs) result(operator_result)
   !< Asymmetric type operator `real.op.integrand`.
   import :: integrand_object, R_P
-  class(integrand_object), intent(in) :: rhs                      !< Right hand side.
-  real(R_P),               intent(in) :: lhs(1:rhs%n)             !< Left hand side.
-  real(R_P)                           :: operator_result(1:rhs%n) !< Operator result.
+  class(integrand_object), intent(in) :: rhs                !< Right hand side.
+  real(R_P),               intent(in) :: lhs(1:)            !< Left hand side.
+  real(R_P), allocatable              :: operator_result(:) !< Operator result.
   endfunction real_op_integrand
 
   pure function integrand_op_real_scalar(lhs, rhs) result(operator_result)
   !< Asymmetric type operator `integrand.op.real`.
   import :: integrand_object, R_P
-  class(integrand_object), intent(in) :: lhs                      !< Left hand side.
-  real(R_P),               intent(in) :: rhs                      !< Right hand side.
-  real(R_P)                           :: operator_result(1:lhs%n) !< Operator result.
+  class(integrand_object), intent(in) :: lhs                !< Left hand side.
+  real(R_P),               intent(in) :: rhs                !< Right hand side.
+  real(R_P), allocatable              :: operator_result(:) !< Operator result.
   endfunction integrand_op_real_scalar
 
   pure function real_scalar_op_integrand(lhs, rhs) result(operator_result)
   !< Asymmetric type operator `real.op.integrand`.
   import :: integrand_object, R_P
-  real(R_P),               intent(in) :: lhs                      !< Left hand side.
-  class(integrand_object), intent(in) :: rhs                      !< Right hand side.
-  real(R_P)                           :: operator_result(1:rhs%n) !< Operator result.
+  real(R_P),               intent(in) :: lhs                !< Left hand side.
+  class(integrand_object), intent(in) :: rhs                !< Right hand side.
+  real(R_P), allocatable              :: operator_result(:) !< Operator result.
   endfunction real_scalar_op_integrand
 
   pure function symmetric_operator(lhs, rhs) result(operator_result)
   !< Symmetric type operator integrand.op.integrand.
   import :: integrand_object, R_P
-  class(integrand_object), intent(in) :: lhs                      !< Left hand side.
-  class(integrand_object), intent(in) :: rhs                      !< Right hand side.
-  real(R_P)                           :: operator_result(1:lhs%n) !< Operator result.
+  class(integrand_object), intent(in) :: lhs                !< Left hand side.
+  class(integrand_object), intent(in) :: rhs                !< Right hand side.
+  real(R_P), allocatable              :: operator_result(:) !< Operator result.
   endfunction symmetric_operator
 
   pure subroutine assignment_integrand(lhs, rhs)
@@ -124,17 +121,8 @@ pure subroutine assignment_integrand(lhs, rhs)
   pure subroutine assignment_real(lhs, rhs)
   !< Symmetric assignment integrand = integrand.
   import :: integrand_object, R_P
-  class(integrand_object), intent(inout) :: lhs          !< Left hand side.
-  real(R_P),               intent(in)    :: rhs(1:lhs%n) !< Right hand side.
+  class(integrand_object), intent(inout) :: lhs     !< Left hand side.
+  real(R_P),               intent(in)    :: rhs(1:) !< Right hand side.
   endsubroutine assignment_real
 endinterface
-
-contains
-  pure subroutine set_dimension(self, n)
-  !< Set ODE *space* dimension.
-  class(integrand_object), intent(inout) :: self !< Ingetrand.
-  integer(I_P),            intent(in)    :: n    !< ODE *space* dimension.
-
-  self%n = n
-  endsubroutine set_dimension
 endmodule foodie_integrand_object

src/tests/accuracy/oscillation/oscillation_oscillator.f90

@@ -79,7 +79,6 @@ pure subroutine init(self, initial_state, frequency)
   real(R_P),         intent(in)    :: initial_state(1:2) !< Initial state of the Oscillation field vector.
   real(R_P),         intent(in)    :: frequency          !< Frequency of oscillation.
 
-  self%n = 2
   self%U = initial_state
   self%f = frequency
   endsubroutine init
@@ -95,12 +94,12 @@ pure function output(self) result(state)
   ! deferred methods
   pure function doscillator_dt(self, t) result(dState_dt)
   !< Time derivative of oscillator field.
-  class(oscillator), intent(in)           :: self                !< Oscillation field.
-  real(R_P),         intent(in), optional :: t                   !< Time.
-  real(R_P)                               :: dState_dt(1:self%n) !< Oscillation field time derivative.
+  class(oscillator), intent(in)           :: self         !< Oscillation field.
+  real(R_P),         intent(in), optional :: t            !< Time.
+  real(R_P), allocatable                  :: dState_dt(:) !< Oscillation field time derivative.
 
-  dState_dt(1) = -self%f * self%U(2)
-  dState_dt(2) =  self%f * self%U(1)
+  dState_dt = [-self%f * self%U(2), &
+                self%f * self%U(1)]
   endfunction doscillator_dt
 
   pure function local_error(lhs, rhs) result(error)
@@ -117,7 +116,7 @@ pure function local_error(lhs, rhs) result(error)
   select type(rhs)
   class is(oscillator)
     error = 0._R_P
-    do i=1, lhs%n
+    do i=1, size(lhs%U, dim=1)
       error = error + (lhs%U(i) - rhs%U(i)) ** 2 / lhs%U(i) ** 2
     enddo
     error = sqrt(error)
@@ -127,9 +126,9 @@ pure function local_error(lhs, rhs) result(error)
   ! +
   pure function integrand_add_integrand(lhs, rhs) result(opr)
   !< `+` operator.
-  class(oscillator),       intent(in) :: lhs          !< Left hand side.
-  class(integrand_object), intent(in) :: rhs          !< Right hand side.
-  real(R_P)                           :: opr(1:lhs%n) !< Operator result.
+  class(oscillator),       intent(in) :: lhs    !< Left hand side.
+  class(integrand_object), intent(in) :: rhs    !< Right hand side.
+  real(R_P), allocatable              :: opr(:) !< Operator result.
 
   select type(rhs)
   class is(oscillator)
@@ -139,28 +138,28 @@ pure function integrand_add_integrand(lhs, rhs) result(opr)
 
   pure function integrand_add_real(lhs, rhs) result(opr)
   !< `+ real` operator.
-  class(oscillator), intent(in) :: lhs          !< Left hand side.
-  real(R_P),         intent(in) :: rhs(1:lhs%n) !< Right hand side.
-  real(R_P)                     :: opr(1:lhs%n) !< Operator result.
+  class(oscillator), intent(in) :: lhs     !< Left hand side.
+  real(R_P),         intent(in) :: rhs(1:) !< Right hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs%U + rhs
   endfunction integrand_add_real
 
   pure function real_add_integrand(lhs, rhs) result(opr)
   !< `real +` operator.
-  class(oscillator), intent(in) :: rhs          !< Left hand side.
-  real(R_P),         intent(in) :: lhs(1:rhs%n) !< Left hand side.
-  real(R_P)                     :: opr(1:rhs%n) !< Operator result.
+  real(R_P),         intent(in) :: lhs(1:) !< Left hand side.
+  class(oscillator), intent(in) :: rhs     !< Left hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs + rhs%U
   endfunction real_add_integrand
 
   ! *
   pure function integrand_multiply_integrand(lhs, rhs) result(opr)
   !< `*` operator.
-  class(oscillator),       intent(in) :: lhs          !< Left hand side.
-  class(integrand_object), intent(in) :: rhs          !< Right hand side.
-  real(R_P)                           :: opr(1:lhs%n) !< Operator result.
+  class(oscillator),       intent(in) :: lhs    !< Left hand side.
+  class(integrand_object), intent(in) :: rhs    !< Right hand side.
+  real(R_P), allocatable              :: opr(:) !< Operator result.
 
   select type(rhs)
   class is(oscillator)
@@ -170,45 +169,46 @@ pure function integrand_multiply_integrand(lhs, rhs) result(opr)
 
   pure function integrand_multiply_real(lhs, rhs) result(opr)
   !< `* real_scalar` operator.
-  class(oscillator), intent(in) :: lhs          !< Left hand side.
-  real(R_P),         intent(in) :: rhs(1:lhs%n) !< Right hand side.
-  real(R_P)                     :: opr(1:lhs%n) !< Operator result.
+  class(oscillator), intent(in) :: lhs     !< Left hand side.
+  real(R_P),         intent(in) :: rhs(1:) !< Right hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs%U * rhs
   endfunction integrand_multiply_real
 
   pure function real_multiply_integrand(lhs, rhs) result(opr)
   !< `real_scalar *` operator.
-  class(oscillator), intent(in) :: rhs          !< Right hand side.
-  real(R_P),         intent(in) :: lhs(1:rhs%n) !< Left hand side.
-  real(R_P)                     :: opr(1:rhs%n) !< Operator result.
+  class(oscillator), intent(in) :: rhs     !< Right hand side.
+  real(R_P),         intent(in) :: lhs(1:) !< Left hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs * rhs%U
   endfunction real_multiply_integrand
 
   pure function integrand_multiply_real_scalar(lhs, rhs) result(opr)
   !< `* real_scalar` operator.
-  class(oscillator), intent(in) :: lhs          !< Left hand side.
-  real(R_P),         intent(in) :: rhs          !< Right hand side.
-  real(R_P)                     :: opr(1:lhs%n) !< Operator result.
+  class(oscillator), intent(in) :: lhs    !< Left hand side.
+  real(R_P),         intent(in) :: rhs    !< Right hand side.
+  real(R_P), allocatable        :: opr(:) !< Operator result.
 
   opr = lhs%U * rhs
   endfunction integrand_multiply_real_scalar
 
   pure function real_scalar_multiply_integrand(lhs, rhs) result(opr)
   !< `real_scalar *` operator.
-  real(R_P),         intent(in) :: lhs          !< Left hand side.
-  class(oscillator), intent(in) :: rhs          !< Right hand side.
-  real(R_P)                     :: opr(1:rhs%n) !< Operator result.
+  real(R_P),         intent(in) :: lhs    !< Left hand side.
+  class(oscillator), intent(in) :: rhs    !< Right hand side.
+  real(R_P), allocatable        :: opr(:) !< Operator result.
 
   opr = lhs * rhs%U
   endfunction real_scalar_multiply_integrand
+
    ! -
   pure function integrand_sub_integrand(lhs, rhs) result(opr)
   !< `-` operator.
-  class(oscillator),       intent(in) :: lhs          !< Left hand side.
-  class(integrand_object), intent(in) :: rhs          !< Right hand side.
-  real(R_P)                           :: opr(1:lhs%n) !< Operator result.
+  class(oscillator),       intent(in) :: lhs    !< Left hand side.
+  class(integrand_object), intent(in) :: rhs    !< Right hand side.
+  real(R_P), allocatable              :: opr(:) !< Operator result.
 
   select type(rhs)
   class is(oscillator)
@@ -218,18 +218,18 @@ pure function integrand_sub_integrand(lhs, rhs) result(opr)
 
   pure function integrand_sub_real(lhs, rhs) result(opr)
   !< `- real` operator.
-  class(oscillator), intent(in) :: lhs          !< Left hand side.
-  real(R_P),         intent(in) :: rhs(1:lhs%n) !< Right hand side.
-  real(R_P)                     :: opr(1:lhs%n) !< Operator result.
+  class(oscillator), intent(in) :: lhs     !< Left hand side.
+  real(R_P),         intent(in) :: rhs(1:) !< Right hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs%U - rhs
   endfunction integrand_sub_real
 
   pure function real_sub_integrand(lhs, rhs) result(opr)
   !< `real -` operator.
-  class(oscillator), intent(in) :: rhs          !< Left hand side.
-  real(R_P),         intent(in) :: lhs(1:rhs%n) !< Left hand side.
-  real(R_P)                     :: opr(1:rhs%n) !< Operator result.
+  class(oscillator), intent(in) :: rhs     !< Left hand side.
+  real(R_P),         intent(in) :: lhs(1:) !< Left hand side.
+  real(R_P), allocatable        :: opr(:)  !< Operator result.
 
   opr = lhs - rhs%U
   endfunction real_sub_integrand
@@ -242,16 +242,15 @@ pure subroutine assign_integrand(lhs, rhs)
 
   select type(rhs)
   class is(oscillator)
-    lhs%n = rhs%n
     lhs%U = rhs%U
     lhs%f = rhs%f
   endselect
   endsubroutine assign_integrand
 
   pure subroutine assign_real(lhs, rhs)
   !< `= real` operator.
-  class(oscillator), intent(inout) :: lhs          !< Left hand side.
-  real(R_P),         intent(in)    :: rhs(1:lhs%n) !< Right hand side.
+  class(oscillator), intent(inout) :: lhs     !< Left hand side.
+  real(R_P),         intent(in)    :: rhs(1:) !< Right hand side.
 
   lhs%U = rhs
   endsubroutine assign_real