expression.f90

Go to the documentation of this file.

! Copyright (c) 2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module expression
  use num_types, only : rp
  use utils, only : neko_error
  use registry, only : neko_const_registry
  use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
  implicit none
  private
 
  public :: expression_eval_static, expression_check_finite
 
  integer, public, parameter :: neko_expr_len = 256
 
  integer, parameter :: expr_max_ident = 64
 
  integer, parameter :: expr_chunk = 1024
 
  !
  ! Opcodes. The evaluator itself does not care about the numbering, but the
  ! parser derives the arity of an opcode from the ranges below, so nullary
  ! opcodes must stay below OP_UNARY and binary ones at or above OP_BINARY.
  !
  integer, parameter :: op_lit = 1, op_x = 2, op_y = 3, op_z = 4, &
       op_t = 5, op_dt = 6
 
  integer, parameter :: op_unary = 20
  integer, parameter :: op_neg = 20, op_sin = 21, op_cos = 22, op_tan = 23, &
       op_asin = 24, op_acos = 25, op_atan = 26, op_sinh = 27, &
       op_cosh = 28, op_tanh = 29, op_exp = 30, op_log = 31, &
       op_log10 = 32, op_sqrt = 33, op_abs = 34, op_erf = 35, &
       op_erfc = 36, op_step = 37, op_ipow = 38
 
  integer, parameter :: op_binary = 60
  integer, parameter :: op_add = 60, op_sub = 61, op_mul = 62, op_div = 63, &
       op_pow = 64, op_atan2 = 65, op_min = 66, op_max = 67, op_mod = 68
 
  !
  ! Token kinds produced by the tokenizer.
  !
  integer, parameter :: tk_end = 0, tk_num = 1, tk_ident = 2, tk_plus = 3, &
       tk_minus = 4, tk_star = 5, tk_slash = 6, tk_pow = 7, tk_lpar = 8, &
       tk_rpar = 9, tk_comma = 10
 
  type :: token_t
     integer :: kind = tk_end
     real(kind=rp) :: num = 0.0_rp
     character(len=EXPR_MAX_IDENT) :: name = ''
     integer :: pos = 0
  end type token_t
 
  type :: parser_t
     type(token_t), allocatable :: tok(:)
     integer :: n_tok = 0
     integer :: cur = 1
     character(len=:), allocatable :: src
     integer, allocatable :: op(:)
     real(kind=rp), allocatable :: lit(:)
     integer :: n_op = 0
     integer :: depth = 0
     integer :: max_depth = 0
     logical :: time_dependent = .false.
  end type parser_t
 
  type, public :: expression_t
     integer, allocatable :: op(:)
     real(kind=rp), allocatable :: lit(:)
     integer :: n_op = 0
     integer :: depth = 0
     logical :: time_dependent = .false.
     character(len=:), allocatable :: src
     real(kind=rp), allocatable, private :: stk(:,:)
   contains
     procedure, pass(this) :: init => expression_init
     procedure, pass(this) :: eval => expression_eval
     procedure, pass(this) :: free => expression_free
  end type expression_t
 
contains
 
  subroutine expression_init(this, str)
    class(expression_t), intent(inout) :: this
    character(len=*), intent(in) :: str
    type(parser_t) :: p
 
    call this%free()
 
    p%src = trim(str)
    if (len(p%src) .eq. 0) call neko_error("Empty expression")
 
    call expr_tokenize(p)
 
    allocate(p%op(p%n_tok + 1))
    allocate(p%lit(p%n_tok + 1))
 
    call parse_expr(p)
    if (p%tok(p%cur)%kind .ne. tk_end) then
       call expr_error(p, "unexpected trailing input", p%tok(p%cur)%pos)
    end if
 
    this%n_op = p%n_op
    this%depth = p%max_depth
    this%time_dependent = p%time_dependent
    this%src = p%src
 
    allocate(this%op(this%n_op), this%lit(this%n_op))
    this%op = p%op(1:this%n_op)
    this%lit = p%lit(1:this%n_op)
 
    allocate(this%stk(expr_chunk, max(this%depth, 1)))
 
    call parser_free(p)
 
  end subroutine expression_init
 
  subroutine parser_free(p)
    type(parser_t), intent(inout) :: p
 
    if (allocated(p%tok)) deallocate(p%tok)
    if (allocated(p%src)) deallocate(p%src)
    if (allocated(p%op)) deallocate(p%op)
    if (allocated(p%lit)) deallocate(p%lit)
 
    p%n_tok = 0
    p%cur = 1
    p%n_op = 0
    p%depth = 0
    p%max_depth = 0
    p%time_dependent = .false.
 
  end subroutine parser_free
 
  subroutine expression_free(this)
    class(expression_t), intent(inout) :: this
 
    if (allocated(this%op)) deallocate(this%op)
    if (allocated(this%lit)) deallocate(this%lit)
    if (allocated(this%src)) deallocate(this%src)
    if (allocated(this%stk)) deallocate(this%stk)
    this%n_op = 0
    this%depth = 0
    this%time_dependent = .false.
 
  end subroutine expression_free
 
  subroutine expression_eval(this, res, n, x, y, z, t, dt)
    class(expression_t), intent(inout) :: this
    integer, intent(in) :: n
    real(kind=rp), intent(inout) :: res(n)
    real(kind=rp), intent(in) :: x(n)
    real(kind=rp), intent(in) :: y(n)
    real(kind=rp), intent(in) :: z(n)
    real(kind=rp), intent(in), optional :: t
    real(kind=rp), intent(in), optional :: dt
    real(kind=rp) :: t_, dt_
    integer :: i0, nb, ip, sp, i, k
 
    if (this%n_op .eq. 0) call neko_error("Evaluating an uncompiled expression")
 
    if (this%time_dependent .and. .not. present(t)) then
       call neko_error("The expression '" // this%src // "' depends on " // &
            "time, but no time is available where it is used")
    end if
 
    t_ = 0.0_rp
    dt_ = 0.0_rp
    if (present(t)) t_ = t
    if (present(dt)) dt_ = dt
 
    do i0 = 1, n, expr_chunk
       nb = min(expr_chunk, n - i0 + 1)
       sp = 0
 
       do ip = 1, this%n_op
          select case (this%op(ip))
 
             !
             ! Operand pushes
             !
          case (op_lit)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = this%lit(ip)
             end do
          case (op_x)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = x(i0 + i - 1)
             end do
          case (op_y)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = y(i0 + i - 1)
             end do
          case (op_z)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = z(i0 + i - 1)
             end do
          case (op_t)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = t_
             end do
          case (op_dt)
             sp = sp + 1
             do i = 1, nb
                this%stk(i, sp) = dt_
             end do
 
             !
             ! Unary operators
             !
          case (op_neg)
             do i = 1, nb
                this%stk(i, sp) = -this%stk(i, sp)
             end do
          case (op_sin)
             do i = 1, nb
                this%stk(i, sp) = sin(this%stk(i, sp))
             end do
          case (op_cos)
             do i = 1, nb
                this%stk(i, sp) = cos(this%stk(i, sp))
             end do
          case (op_tan)
             do i = 1, nb
                this%stk(i, sp) = tan(this%stk(i, sp))
             end do
          case (op_asin)
             do i = 1, nb
                this%stk(i, sp) = asin(this%stk(i, sp))
             end do
          case (op_acos)
             do i = 1, nb
                this%stk(i, sp) = acos(this%stk(i, sp))
             end do
          case (op_atan)
             do i = 1, nb
                this%stk(i, sp) = atan(this%stk(i, sp))
             end do
          case (op_sinh)
             do i = 1, nb
                this%stk(i, sp) = sinh(this%stk(i, sp))
             end do
          case (op_cosh)
             do i = 1, nb
                this%stk(i, sp) = cosh(this%stk(i, sp))
             end do
          case (op_tanh)
             do i = 1, nb
                this%stk(i, sp) = tanh(this%stk(i, sp))
             end do
          case (op_exp)
             do i = 1, nb
                this%stk(i, sp) = exp(this%stk(i, sp))
             end do
          case (op_log)
             do i = 1, nb
                this%stk(i, sp) = log(this%stk(i, sp))
             end do
          case (op_log10)
             do i = 1, nb
                this%stk(i, sp) = log10(this%stk(i, sp))
             end do
          case (op_sqrt)
             do i = 1, nb
                this%stk(i, sp) = sqrt(this%stk(i, sp))
             end do
          case (op_abs)
             do i = 1, nb
                this%stk(i, sp) = abs(this%stk(i, sp))
             end do
          case (op_erf)
             do i = 1, nb
                this%stk(i, sp) = erf(this%stk(i, sp))
             end do
          case (op_erfc)
             do i = 1, nb
                this%stk(i, sp) = erfc(this%stk(i, sp))
             end do
          case (op_step)
             do i = 1, nb
                if (this%stk(i, sp) .lt. 0.0_rp) then
                   this%stk(i, sp) = 0.0_rp
                else
                   this%stk(i, sp) = 1.0_rp
                end if
             end do
          case (op_ipow)
             ! Integer exponent, folded by the parser so that a negative base
             ! raised to a whole power stays well defined.
             k = nint(this%lit(ip))
             do i = 1, nb
                this%stk(i, sp) = this%stk(i, sp) ** k
             end do
 
             !
             ! Binary operators
             !
          case (op_add)
             do i = 1, nb
                this%stk(i, sp-1) = this%stk(i, sp-1) + this%stk(i, sp)
             end do
             sp = sp - 1
          case (op_sub)
             do i = 1, nb
                this%stk(i, sp-1) = this%stk(i, sp-1) - this%stk(i, sp)
             end do
             sp = sp - 1
          case (op_mul)
             do i = 1, nb
                this%stk(i, sp-1) = this%stk(i, sp-1) * this%stk(i, sp)
             end do
             sp = sp - 1
          case (op_div)
             do i = 1, nb
                this%stk(i, sp-1) = this%stk(i, sp-1) / this%stk(i, sp)
             end do
             sp = sp - 1
          case (op_pow)
             do i = 1, nb
                this%stk(i, sp-1) = this%stk(i, sp-1) ** this%stk(i, sp)
             end do
             sp = sp - 1
          case (op_atan2)
             do i = 1, nb
                this%stk(i, sp-1) = atan2(this%stk(i, sp-1), this%stk(i, sp))
             end do
             sp = sp - 1
          case (op_min)
             do i = 1, nb
                this%stk(i, sp-1) = min(this%stk(i, sp-1), this%stk(i, sp))
             end do
             sp = sp - 1
          case (op_max)
             do i = 1, nb
                this%stk(i, sp-1) = max(this%stk(i, sp-1), this%stk(i, sp))
             end do
             sp = sp - 1
          case (op_mod)
             do i = 1, nb
                this%stk(i, sp-1) = mod(this%stk(i, sp-1), this%stk(i, sp))
             end do
             sp = sp - 1
 
          case default
             call neko_error("Corrupt expression opcode stream")
          end select
       end do
 
       if (sp .ne. 1) call neko_error("Corrupt expression opcode stream")
 
       do i = 1, nb
          res(i0 + i - 1) = this%stk(i, 1)
       end do
    end do
 
  end subroutine expression_eval
 
  subroutine expression_eval_static(str, res, n, x, y, z, usage)
    character(len=*), intent(in) :: str
    integer, intent(in) :: n
    real(kind=rp), intent(inout) :: res(n)
    real(kind=rp), intent(in) :: x(n)
    real(kind=rp), intent(in) :: y(n)
    real(kind=rp), intent(in) :: z(n)
    character(len=*), intent(in) :: usage
    type(expression_t) :: e
 
    if (len_trim(str) .eq. 0) then
       call neko_error("The " // usage // " needs a non-empty expression")
    end if
 
    call e%init(str)
 
    if (e%time_dependent) then
       call neko_error("The " // usage // " expression '" // trim(str) // &
            "' depends on time, which is not available where it is used")
    end if
 
    call e%eval(res, n, x, y, z)
    call e%free()
 
    call expression_check_finite(str, res, n, usage)
 
  end subroutine expression_eval_static
 
  subroutine expression_check_finite(str, res, n, usage)
    character(len=*), intent(in) :: str
    integer, intent(in) :: n
    real(kind=rp), intent(in) :: res(n)
    character(len=*), intent(in) :: usage
    integer :: i
 
    do i = 1, n
       if (.not. ieee_is_finite(res(i))) then
          call neko_error("The " // usage // " expression '" // trim(str) // &
               "' does not evaluate to a finite value everywhere")
       end if
    end do
 
  end subroutine expression_check_finite
 
  !
  ! Tokenizer
  !
 
  subroutine expr_tokenize(p)
    type(parser_t), intent(inout) :: p
    integer :: i, j, k, n, ios
    logical :: seen_dot
    character(len=1) :: c
    character(len=EXPR_MAX_IDENT) :: buf
 
    n = len(p%src)
    allocate(p%tok(n + 1))
    p%n_tok = 0
 
    i = 1
    do while (i .le. n)
       c = p%src(i:i)
 
       if (c .eq. ' ' .or. c .eq. achar(9)) then
          i = i + 1
          cycle
       end if
 
       if (is_digit(c) .or. (c .eq. '.' .and. i .lt. n .and. &
            is_digit(p%src(min(i+1, n):min(i+1, n))))) then
          ! Mantissa
          j = i
          seen_dot = .false.
          do while (j .le. n)
             if (is_digit(p%src(j:j))) then
                j = j + 1
             else if (p%src(j:j) .eq. '.' .and. .not. seen_dot) then
                seen_dot = .true.
                j = j + 1
             else
                exit
             end if
          end do
          ! Exponent, but only if it is actually followed by digits, so that
          ! e.g. `2*d` is not mistaken for a malformed number.
          if (j .le. n) then
             if (scan(p%src(j:j), 'eEdD') .gt. 0) then
                k = j + 1
                if (k .le. n) then
                   if (scan(p%src(k:k), '+-') .gt. 0) k = k + 1
                end if
                if (k .le. n) then
                   if (is_digit(p%src(k:k))) then
                      do while (k .le. n)
                         if (.not. is_digit(p%src(k:k))) exit
                         k = k + 1
                      end do
                      j = k
                   end if
                end if
             end if
          end if
 
          if (j - i .gt. expr_max_ident) then
             call expr_error(p, "numeric literal is too long", i)
          end if
          buf = p%src(i:j-1)
          ! Fortran list-directed input does not accept a `d` exponent in
          ! every implementation, so normalise it.
          do k = 1, len_trim(buf)
             if (buf(k:k) .eq. 'd' .or. buf(k:k) .eq. 'D') buf(k:k) = 'e'
          end do
          call push_token(p, tk_num, i)
          read (buf, *, iostat = ios) p%tok(p%n_tok)%num
          if (ios .ne. 0) call expr_error(p, "malformed number", i)
          i = j
          cycle
       end if
 
       if (is_alpha(c)) then
          j = i
          do while (j .le. n)
             if (is_alpha(p%src(j:j)) .or. is_digit(p%src(j:j)) .or. &
                  p%src(j:j) .eq. '_') then
                j = j + 1
             else
                exit
             end if
          end do
          if (j - i .gt. expr_max_ident) then
             call expr_error(p, "identifier is too long", i)
          end if
          call push_token(p, tk_ident, i)
          p%tok(p%n_tok)%name = p%src(i:j-1)
          i = j
          cycle
       end if
 
       select case (c)
       case ('+')
          call push_token(p, tk_plus, i)
       case ('-')
          call push_token(p, tk_minus, i)
       case ('*')
          if (i .lt. n) then
             if (p%src(i+1:i+1) .eq. '*') then
                call push_token(p, tk_pow, i)
                i = i + 2
                cycle
             end if
          end if
          call push_token(p, tk_star, i)
       case ('/')
          call push_token(p, tk_slash, i)
       case ('^')
          call push_token(p, tk_pow, i)
       case ('(')
          call push_token(p, tk_lpar, i)
       case (')')
          call push_token(p, tk_rpar, i)
       case (',')
          call push_token(p, tk_comma, i)
       case default
          call expr_error(p, "unexpected character '" // c // "'", i)
       end select
       i = i + 1
    end do
 
    call push_token(p, tk_end, n + 1)
 
  end subroutine expr_tokenize
 
  subroutine push_token(p, kind, pos)
    type(parser_t), intent(inout) :: p
    integer, intent(in) :: kind
    integer, intent(in) :: pos
 
    p%n_tok = p%n_tok + 1
    p%tok(p%n_tok)%kind = kind
    p%tok(p%n_tok)%pos = pos
 
  end subroutine push_token
 
  !
  ! Recursive descent parser, emitting postfix opcodes as it goes.
  !
  !   expr    := term (('+' | '-') term)*
  !   term    := factor (('*' | '/') factor)*
  !   factor  := ('+' | '-') factor | power
  !   power   := primary ['^' factor]
  !   primary := NUM | IDENT | IDENT '(' expr (',' expr)* ')' | '(' expr ')'
  !
 
  recursive subroutine parse_expr(p)
    type(parser_t), intent(inout) :: p
 
    call parse_term(p)
    do
       if (p%tok(p%cur)%kind .eq. tk_plus) then
          p%cur = p%cur + 1
          call parse_term(p)
          call emit(p, op_add)
       else if (p%tok(p%cur)%kind .eq. tk_minus) then
          p%cur = p%cur + 1
          call parse_term(p)
          call emit(p, op_sub)
       else
          exit
       end if
    end do
 
  end subroutine parse_expr
 
  recursive subroutine parse_term(p)
    type(parser_t), intent(inout) :: p
 
    call parse_factor(p)
    do
       if (p%tok(p%cur)%kind .eq. tk_star) then
          p%cur = p%cur + 1
          call parse_factor(p)
          call emit(p, op_mul)
       else if (p%tok(p%cur)%kind .eq. tk_slash) then
          p%cur = p%cur + 1
          call parse_factor(p)
          call emit(p, op_div)
       else
          exit
       end if
    end do
 
  end subroutine parse_term
 
  recursive subroutine parse_factor(p)
    type(parser_t), intent(inout) :: p
    integer :: before
 
    if (p%tok(p%cur)%kind .eq. tk_minus) then
       p%cur = p%cur + 1
       before = p%n_op
       call parse_factor(p)
       if (p%n_op .eq. before + 1 .and. p%op(p%n_op) .eq. op_lit) then
          ! Fold the sign into the literal, so that e.g. `y^-2` can use an
          ! integer power rather than the real valued one.
          p%lit(p%n_op) = -p%lit(p%n_op)
       else
          call emit(p, op_neg)
       end if
    else if (p%tok(p%cur)%kind .eq. tk_plus) then
       p%cur = p%cur + 1
       call parse_factor(p)
    else
       call parse_power(p)
    end if
 
  end subroutine parse_factor
 
  recursive subroutine parse_power(p)
    type(parser_t), intent(inout) :: p
    integer :: before
 
    call parse_primary(p)
 
    if (p%tok(p%cur)%kind .eq. tk_pow) then
       p%cur = p%cur + 1
       before = p%n_op
       call parse_factor(p)
 
       if (p%n_op .eq. before + 1 .and. p%op(p%n_op) .eq. op_lit) then
          ! Tested exactly, and not within a tolerance: an exponent such as
          ! `2.0000005` is not a whole number, and folding it would silently
          ! turn an undefined negative base power into a finite value.
          if (p%lit(p%n_op) .eq. anint(p%lit(p%n_op)) .and. &
               abs(p%lit(p%n_op)) .lt. 1.0e4_rp) then
             ! Whole exponent. Rewrite the literal push into a unary integer
             ! power, which unlike `**` with a real exponent is also defined
             ! for a negative base.
             p%op(p%n_op) = op_ipow
             p%depth = p%depth - 1
             return
          end if
       end if
 
       call emit(p, op_pow)
    end if
 
  end subroutine parse_power
 
  recursive subroutine parse_primary(p)
    type(parser_t), intent(inout) :: p
    integer :: k, op, nargs
 
    select case (p%tok(p%cur)%kind)
    case (tk_num)
       call emit(p, op_lit, p%tok(p%cur)%num)
       p%cur = p%cur + 1
 
    case (tk_lpar)
       p%cur = p%cur + 1
       call parse_expr(p)
       if (p%tok(p%cur)%kind .ne. tk_rpar) then
          call expr_error(p, "expected ')'", p%tok(p%cur)%pos)
       end if
       p%cur = p%cur + 1
 
    case (tk_ident)
       k = p%cur
       p%cur = p%cur + 1
       if (p%tok(p%cur)%kind .eq. tk_lpar) then
          p%cur = p%cur + 1
          call parse_expr(p)
          nargs = 1
          do while (p%tok(p%cur)%kind .eq. tk_comma)
             p%cur = p%cur + 1
             call parse_expr(p)
             nargs = nargs + 1
          end do
          if (p%tok(p%cur)%kind .ne. tk_rpar) then
             call expr_error(p, "expected ')' closing the argument list of '" &
                  // trim(p%tok(k)%name) // "'", p%tok(p%cur)%pos)
          end if
          p%cur = p%cur + 1
 
          op = func_op(trim(p%tok(k)%name), nargs)
          if (op .eq. -1) then
             call expr_error(p, "unknown function '" // &
                  trim(p%tok(k)%name) // "'", p%tok(k)%pos)
          else if (op .eq. -2) then
             call expr_error(p, "wrong number of arguments to '" // &
                  trim(p%tok(k)%name) // "'", p%tok(k)%pos)
          end if
          call emit(p, op)
       else
          call parse_symbol(p, k)
       end if
 
    case default
       call expr_error(p, "expected a value", p%tok(p%cur)%pos)
    end select
 
  end subroutine parse_primary
 
  subroutine parse_symbol(p, k)
    type(parser_t), intent(inout) :: p
    integer, intent(in) :: k
    character(len=:), allocatable :: name
    real(kind=rp), pointer :: rval
    integer, pointer :: ival
 
    name = trim(p%tok(k)%name)
 
    select case (name)
    case ('x')
       call emit(p, op_x)
    case ('y')
       call emit(p, op_y)
    case ('z')
       call emit(p, op_z)
    case ('t')
       call emit(p, op_t)
       p%time_dependent = .true.
    case ('dt')
       call emit(p, op_dt)
       p%time_dependent = .true.
    case ('pi')
       call emit(p, op_lit, 4.0_rp * atan(1.0_rp))
    case default
       if (neko_const_registry%real_scalar_exists(name)) then
          rval => neko_const_registry%get_real_scalar(name)
          call emit(p, op_lit, rval)
       else if (neko_const_registry%integer_scalar_exists(name)) then
          ival => neko_const_registry%get_integer_scalar(name)
          call emit(p, op_lit, real(ival, kind=rp))
       else
          call expr_error(p, "unknown symbol '" // name // "', it is " // &
               "neither a coordinate nor a constant declared under " // &
               "case.constants", p%tok(k)%pos)
       end if
    end select
 
    if (allocated(name)) deallocate(name)
 
  end subroutine parse_symbol
 
  function func_op(name, nargs) result(op)
    character(len=*), intent(in) :: name
    integer, intent(in) :: nargs
    integer :: op, arity
 
    arity = 1
    select case (name)
    case ('sin')
       op = op_sin
    case ('cos')
       op = op_cos
    case ('tan')
       op = op_tan
    case ('asin')
       op = op_asin
    case ('acos')
       op = op_acos
    case ('atan')
       op = op_atan
    case ('sinh')
       op = op_sinh
    case ('cosh')
       op = op_cosh
    case ('tanh')
       op = op_tanh
    case ('exp')
       op = op_exp
    case ('log')
       op = op_log
    case ('log10')
       op = op_log10
    case ('sqrt')
       op = op_sqrt
    case ('abs')
       op = op_abs
    case ('erf')
       op = op_erf
    case ('erfc')
       op = op_erfc
    case ('step')
       op = op_step
    case ('atan2')
       op = op_atan2
       arity = 2
    case ('min')
       op = op_min
       arity = 2
    case ('max')
       op = op_max
       arity = 2
    case ('mod')
       op = op_mod
       arity = 2
    case default
       op = -1
       return
    end select
 
    if (nargs .ne. arity) op = -2
 
  end function func_op
 
  subroutine emit(p, op, val)
    type(parser_t), intent(inout) :: p
    integer, intent(in) :: op
    real(kind=rp), intent(in), optional :: val
 
    if (p%n_op .ge. size(p%op)) then
       call neko_error("Expression '" // p%src // "' is too complex")
    end if
 
    p%n_op = p%n_op + 1
    p%op(p%n_op) = op
    if (present(val)) then
       p%lit(p%n_op) = val
    else
       p%lit(p%n_op) = 0.0_rp
    end if
 
    if (op .lt. op_unary) then
       p%depth = p%depth + 1
    else if (op .ge. op_binary) then
       p%depth = p%depth - 1
    end if
    p%max_depth = max(p%max_depth, p%depth)
 
  end subroutine emit
 
  subroutine expr_error(p, msg, pos)
    type(parser_t), intent(in) :: p
    character(len=*), intent(in) :: msg
    integer, intent(in) :: pos
    character(len=16) :: buf
 
    write (buf, '(I0)') pos
    call neko_error("Invalid expression '" // p%src // "', at character " // &
         trim(buf) // ": " // msg)
 
  end subroutine expr_error
 
  pure function is_digit(c) result(res)
    character(len=1), intent(in) :: c
    logical :: res
 
    res = (c .ge. '0' .and. c .le. '9')
 
  end function is_digit
 
  pure function is_alpha(c) result(res)
    character(len=1), intent(in) :: c
    logical :: res
 
    res = (c .ge. 'a' .and. c .le. 'z') .or. (c .ge. 'A' .and. c .le. 'Z')
 
  end function is_alpha
 
end module expression