Cell-integral or cell-average of over grid, using Simpson's 1/3 rule.
f(x) ^
| * * * *
| * * *++++++++++ *
| * *
| ++*++++++ *
| * ++++*+++++
| * *
-|----|----.----|-----.----|----.-----|---->
x_{i-i} x_i x_{i+1}
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
procedure(fx) | :: | fnc |
function to integrate/average over grid cells |
|||
type(grid1), | intent(in) | :: | grid |
|
||
logical, | intent(in), | optional | :: | average |
flag to compute cell-average instead of cell-integral |
pure function quadgrid1(fnc, grid, average) result(res)
!! Cell-integral or cell-average of \( f(x) \) over grid, using Simpson's 1/3 rule.
!! ```
!! f(x) ^
!! | * * * *
!! | * * *++++++++++ *
!! | * *
!! | ++*++++++ *
!! | * ++++*+++++
!! | * *
!! -|----|----.----|-----.----|----.-----|---->
!! x_{i-i} x_i x_{i+1}
!! ```
procedure(fx) :: fnc
!! function \( f(x) \) to integrate/average over grid cells
type(grid1), intent(in) :: grid
!! `grid1` object
logical, intent(in), optional :: average
!! flag to compute cell-average instead of cell-integral
real(rk) :: res(grid%ncells)
real(rk) :: fedges(0:grid%ncells), fcenter(grid%ncells)
integer :: i
associate (nc => grid%ncells)
! Evaluate f(x) at grid edges
do concurrent(i=0:nc)
fedges(i) = fnc(grid%edges(i))
end do
! Evaluate f(x) at grid centers
do concurrent(i=1:nc)
fcenter(i) = fnc(grid%center(i))
end do
! Cell-average using Simpson's 1/3 rule
res = (4*fcenter + fedges(0:nc - 1) + fedges(1:nc))/6
! Convert cell-average to cell-integral
if (.not. (optval(average, .false.))) then
res = res*grid%width
end if
end associate
end function quadgrid1