method hho_navier_stokes {
  parameter pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198

  // Viscosity
  parameter nu = 0.025

  // Reynolds number
  parameter reynolds = pow(nu, -1.0)

  // Kovasznay solution
  parameter lambda = 0.5 * reynolds - sqrt(0.25 * pow(reynolds, 2.0) + 4.0 * pow(pi, 2.0))

  parameter exact_pressure_average =
    0.5 * lambda * (exp(4.0 * lambda) - 1.0)
    - (exp(3.0 * lambda) - exp(-lambda)) * pow(8.0 * lambda, -1.0)

  function exact_velocity(vector X) -> vector
  = vector(
            1.0 - exp(lambda * X[0]) * cos(2.0 * pi * X[1]),
            lambda * pow(2.0 * pi, -1.0) * exp(lambda * X[0]) * sin(2.0 * pi * X[1])
          )

  function exact_pressure(vector X) -> scalar =
    -0.5 * exp(2.0 * lambda * X[0])
    + 0.5 * lambda * (exp(4.0 * lambda) - 1.0)

  function source(vector X) -> vector = vector(0.0, 0.0)

  function pressure_average_source(vector X) -> scalar = pow(nu, -1.0) * exact_pressure_average

  // Velocity space
  space Uh {
    element Poly(k, vector)
    edge Poly(k, vector)
  }

  // Pressure space
  space Ph {
    element Poly(k, scalar)
  }

  // Velocity-pressure Cartesian product space
  product space Xh = Uh times Ph

  // Velocity interpolant
  interpolant interpolant_exact_velocity on Uh {
    on element T:
      dof(T) = l2_projection(exact_velocity, Poly(k, vector))
    on edge E:
      dof(E) = l2_projection(exact_velocity, Poly(k, vector))
  }

  // Pressure interpolant
  interpolant interpolant_exact_pressure on Ph {
    on element T:
      dof(T) = l2_projection(exact_pressure, Poly(k, scalar))
  }

  // Operators
  operator velocity_gradient : Uh(u) -> Poly(k, matrix) on element T {
    forall tau in Poly(k, matrix):
      int(T) velocity_gradient(u) dot tau =
      - int(T) dof(u, T) dot div(tau)
      + int(dT) dof(u, E) dot (tau dot normal)
  }

  operator velocity_potential : Uh(u) -> Poly(k+1, vector) on element T {
    forall w in Poly(k+1, vector):
      int(T) grad(velocity_potential(u)) dot grad(w) =
      int(T) velocity_gradient(u) dot grad(w)
    constraint int(T) velocity_potential(u) = int(T) dof(u, T)
  }

  operator velocity_divergence : Uh(u) -> Poly(k, scalar) on element T {
    forall q in Poly(k, scalar):
      int(T) velocity_divergence(u) * q =
      - int(T) dof(u, T) dot grad(q)
      + int(dT) (dof(u, E) dot normal) * q
  }

  operator element_difference : Uh(u) -> Poly(k, vector) on element T {
    forall w in Poly(k, vector):
      int(T) element_difference(u) dot w =
      int(T) (velocity_potential(u) - dof(u, T)) dot w
  }

  operator edge_difference : Uh(u) -> Poly(k, vector) on edge E of element T {
    forall w in Poly(k, vector):
      int(E) edge_difference(u) dot w =
      int(E) (velocity_potential(u) - dof(u, E) - element_difference(u)) dot w
  }

  operator edge_jump : Uh(u) -> Poly(k, vector) on edge E of element T {
    forall w in Poly(k, vector):
      int(E) edge_jump(u) dot w =
      int(E) (dof(u, E) - dof(u, T)) dot w
  }

  operator pressure_field : Ph(p) -> Poly(k) on element T {
    pressure_field(p) = dof(p, T)
  }

  // Forms
  bilinear form viscous_term : Uh(trial u) times Uh(test v) {
    sum_elements(
                  nu * int(T) velocity_gradient(u) dot velocity_gradient(v)
                  + nu * pow(diameter(T), -1.0) * int(dT) edge_difference(u) dot edge_difference(v)
                )
  }

  bilinear form pressure_velocity : Ph(trial p) times Uh(test v) {
    -sum_elements(int(T) dof(p, T) * velocity_divergence(v))
  }

  bilinear form velocity_pressure : Uh(trial u) times Ph(test q) {
    sum_elements(int(T) velocity_divergence(u) * dof(q, T))
  }

  bilinear form pressure_pressure : Ph(trial p) times Ph(test q) {
    sum_elements(pow(nu, -1.0) * int(T) dof(p, T) * global_average(q))
  }

  linear form momentum_balance_load : Uh(test v) {
    sum_elements(int(T) source dot dof(v, T))
  }

  linear form mass_balance_load : Ph(test q) {
    sum_elements(int(T) pressure_average_source * dof(q, T))
  }

  // Residuals
  residual form convective_residual : Uh(state current) times Uh(test v) on element T {
    0.5 * (
            int(T) (grad(dof(current, T)) dot dof(current, T)) dot dof(v, T)
            - int(T) dof(current, T) dot (grad(dof(v, T)) dot dof(current, T))
            + sum_edges(
                         int(E) (dof(current, E) dot normal)
                         * (
                             dof(current, E) dot dof(v, T)
                             + dof(current, E) dot dof(v, E)
                             - dof(v, E) dot dof(current, T)
                           )
                       )
          )
  }

  // Newton derivative of the convective term
  jacobian form convective_jacobian : Uh(state current) times Uh(trial du) times Uh(test v) on element T {
    0.5 * (
            int(T) (grad(dof(current, T)) dot dof(du, T)) dot dof(v, T)
            - int(T) dof(current, T) dot (grad(dof(v, T)) dot dof(du, T))
            + sum_edges(
                         int(E) (dof(du, E) dot normal)
                         * (
                             dof(current, E) dot dof(v, T)
                             + dof(current, E) dot dof(v, E)
                             - dof(v, E) dot dof(current, T)
                           )
                       )
            + int(T) (grad(dof(du, T)) dot dof(current, T)) dot dof(v, T)
            - int(T) dof(du, T) dot (grad(dof(v, T)) dot dof(current, T))
            + sum_edges(
                         int(E) (dof(current, E) dot normal)
                         * (
                             dof(du, E) dot dof(v, T)
                             + dof(du, E) dot dof(v, E)
                             - dof(v, E) dot dof(du, T)
                           )
                       )
          )
    }

  // Picard matrix t_h(current, du, v). Substituting du = current recovers
  // convective_residual
  jacobian form convective_fixed_point : Uh(state current) times Uh(trial du) times Uh(test v) on element T {
      0.5 * (
              int(T) (grad(dof(du, T)) dot dof(current, T)) dot dof(v, T)
              - int(T) dof(du, T) dot (grad(dof(v, T)) dot dof(current, T))
              + sum_edges(
                           int(E) (dof(current, E) dot normal)
                           * (
                               dof(du, E) dot dof(v, T)
                               + dof(du, E) dot dof(v, E)
                               - dof(v, E) dot dof(du, T)
                             )
                         )
            )
    }

  // Boundary conditions
  boundary conditions dirichlet_boundary_conditions on Uh {
    on edge E:
      dof(E) = l2_project(exact_velocity, Poly(k, vector))
  }

  // Complete mixed residual and its two linearisations
  nonlinear problem navier_stokes_problem on Xh {
    residual {
      viscous_term
      + pressure_velocity
      + velocity_pressure
      + pressure_pressure
      + convective_residual
      - momentum_balance_load
      - mass_balance_load
    }
    jacobian {
      viscous_term
      + pressure_velocity
      + velocity_pressure
      + pressure_pressure
      + convective_jacobian
    }
    fixed_point {
      viscous_term
      + pressure_velocity
      + velocity_pressure
      + pressure_pressure
      + convective_fixed_point
    }
    boundary conditions dirichlet_boundary_conditions
    export { velocity_potential, pressure_field }
    compute errors using interpolant_exact_velocity {
      velocity_h1_norm, velocity_l2_norm // , divergence_l2_norm
    }
    compute errors using interpolant_exact_pressure {
      pressure_l2_norm // , pressure_average_norm
    }
  }

  // Error functionals
  function velocity_h1_norm : Uh(v) -> scalar {
    sqrt(
          sum_elements(
                        int(T) squared_norm(grad(dof(v, T)))
                        + pow(diameter(T), -1.0) * int(dT) squared_norm(edge_jump(v))
                      )
        )
  }

  function velocity_l2_norm : Uh(v) -> scalar {
    sqrt(
          sum_elements(
                        int(T) squared_norm(dof(v, T))
                        + diameter(T) * int(dT) squared_norm(edge_jump(v))
                      )
        )
  }

  function divergence_l2_norm : Uh(v) -> scalar {
    sqrt(sum_elements(int(T) pow(velocity_divergence(v), 2.0)))
  }

  function pressure_l2_norm : Ph(q) -> scalar {
    sqrt(sum_elements(int(T) squared_norm(dof(q, T))))
  }

  function pressure_average_norm : Ph(q) -> scalar {
    sum_elements(int(T) dof(q, T))
  }
}