Skip to content

User supplied functions for calculus algorithms (e.g. derivatives) #30

@fluidnumerics-joe

Description

@fluidnumerics-joe

While we can now implement custom functions, I've been thinking about how we might be able to perform calculus operations, e.g. derivatives. The problem is that derivatives require knowledge of dimensionality, which does not seem to fit with the f32 and f64 interfaces, as they are currently defined.

Below is a simple example of what would be nice to accomplish :

program simple_derivative
    use feqparse
    use iso_fortran_env

    implicit none

    type(EquationParser) :: f, dfdx

    f = EquationParser("f = 1.0",["x"])

    call AddFunction("dfdxc2",centered)
    dfdx = EquationParser("g = dfdxc2(f)",["f"])

    ! x = [0,1]
    ! f_array = f % Evaluate(x)
    ! dfdx_array = dfdx % Evaluate(f_array)
    contains

    pure function centered(f) result(dfdx)
    real(real64),intent(in) :: f(:)
    real(real64) :: dfdx(lbound(f,1):ubound(f,1))
    ! Local
    integer :: i

      dfdx(lbound(f,1)) = 0.0_real64
      do i = lbound(f,1)+1, ubound(f,1)-1
        dfdx(i) = f(i+1) - f(i-1)
      enddo
      dfdx(ubound(f,1)) = 0.0_real64

    end function centered

end program simple_derivative

Such a setup would allow someone to provide a handful of pre-canned ways to compute a discrete derivative. Here, we're registering just centered as a new function. However, we could imagine registering multiple functions for computing a derivative and then having a program wherein the end-user could supply a string (through CLI or other input) that evaluates the equation.

Currently, this example fails to compile with the following

src//simple-derivative.f90:11:39:

   11 |     call AddFunction("dfdxc2",centered)
      |                                       1
Error: There is no specific subroutine for the generic ‘addfunction’ at (1)
<ERROR> Compilation failed for object " src_simple-derivative.f90.o "
<ERROR> stopping due to failed compilation
STOP 1

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions