compute_momentum_rhs Subroutine

private subroutine compute_momentum_rhs(mesh, flow, bc, params, transport, u, p, local_rhs)

Evaluates the advective, diffusive, and pressure terms of the momentum equation.

Supports both Upwind (stable) and Central (high-accuracy) convection schemes. Advection is computed using a flux-form discretization.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
type(flow_mpi_t), intent(in) :: flow
type(bc_set_t), intent(in) :: bc
type(case_params_t), intent(in) :: params
type(transport_properties_t), intent(in) :: transport
real(kind=rk), intent(in) :: u(:,:)
real(kind=rk), intent(in) :: p(:)
real(kind=rk), intent(out) :: local_rhs(:,:)

Calls

proc~~compute_momentum_rhs~~CallsGraph proc~compute_momentum_rhs mod_flow_projection::compute_momentum_rhs proc~boundary_velocity mod_bc::boundary_velocity proc~compute_momentum_rhs->proc~boundary_velocity proc~face_effective_neighbor mod_bc::face_effective_neighbor proc~compute_momentum_rhs->proc~face_effective_neighbor proc~face_linear_vector mod_flow_projection::face_linear_vector proc~compute_momentum_rhs->proc~face_linear_vector proc~face_normal_distance mod_flow_projection::face_normal_distance proc~compute_momentum_rhs->proc~face_normal_distance proc~fatal_error mod_kinds::fatal_error proc~compute_momentum_rhs->proc~fatal_error proc~outward_normal mod_flow_projection::outward_normal proc~compute_momentum_rhs->proc~outward_normal proc~pressure_gradient_cell mod_flow_projection::pressure_gradient_cell proc~compute_momentum_rhs->proc~pressure_gradient_cell proc~is_periodic_face mod_bc::is_periodic_face proc~face_effective_neighbor->proc~is_periodic_face proc~face_neighbor_weight mod_flow_projection::face_neighbor_weight proc~face_linear_vector->proc~face_neighbor_weight proc~face_normal_distance->proc~fatal_error proc~face_normal_distance->proc~outward_normal proc~patch_type_for_face mod_bc::patch_type_for_face proc~face_normal_distance->proc~patch_type_for_face proc~pressure_gradient_cell->proc~face_effective_neighbor proc~pressure_gradient_cell->proc~outward_normal proc~boundary_pressure mod_bc::boundary_pressure proc~pressure_gradient_cell->proc~boundary_pressure proc~boundary_pressure_type mod_bc::boundary_pressure_type proc~pressure_gradient_cell->proc~boundary_pressure_type proc~face_linear_scalar mod_flow_projection::face_linear_scalar proc~pressure_gradient_cell->proc~face_linear_scalar proc~face_linear_scalar->proc~face_neighbor_weight proc~face_neighbor_weight->proc~fatal_error proc~face_neighbor_weight->proc~outward_normal proc~face_neighbor_weight->proc~patch_type_for_face proc~is_periodic_face->proc~patch_type_for_face

Called by

proc~~compute_momentum_rhs~~CalledByGraph proc~compute_momentum_rhs mod_flow_projection::compute_momentum_rhs proc~advance_projection_step mod_flow_projection::advance_projection_step proc~advance_projection_step->proc~compute_momentum_rhs program~lowmach_react_hex lowmach_react_hex program~lowmach_react_hex->proc~advance_projection_step

Source Code

   subroutine compute_momentum_rhs(mesh, flow, bc, params, transport, u, p, local_rhs)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: u(:,:)
      real(rk), intent(in) :: p(:)
      real(rk), intent(out) :: local_rhs(:,:)

      integer :: c, lf, f, nb
      real(rk) :: nvec(3), nhat(3)
      real(rk) :: uf(3), ub(3), advected(3)
      real(rk) :: un_area
      real(rk) :: conv(3), diff(3), gradp(3)
      real(rk) :: dist, diff_face
      logical :: use_central

      use_central = .false.

      select case (trim(params%convection_scheme))
      case ('central', 'central_difference', 'central-difference')
         use_central = .true.
      case ('upwind')
         use_central = .false.
      case default
         call fatal_error('flow', 'unknown convection_scheme '//trim(params%convection_scheme))
      end select

      local_rhs = zero

      do c = flow%first_cell, flow%last_cell
         conv = zero
         diff = zero

         ! Compute local pressure gradient at cell center
         call pressure_gradient_cell(mesh, bc, p, c, gradp)

         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)

            nvec = outward_normal(mesh, f, c)
            nhat = nvec

            nb = face_effective_neighbor(mesh, bc, f, c)

            if (nb > 0) then
               uf = face_linear_vector(mesh, bc, f, c, nb, u(:, c), u(:, nb))
               ub = u(:, nb)
               dist = face_normal_distance(mesh, bc, f, c, nb)
            else
               call boundary_velocity(mesh, bc, f, u(:, c), ub)
               uf = ub
               dist = face_normal_distance(mesh, bc, f, c, 0)
            end if

            un_area = dot_product(uf, nhat) * mesh%faces(f)%area

            ! Scheme selection for advection
            if (use_central) then
               advected = uf
            else
               if (un_area >= zero) then
                  advected = u(:, c)
               else
                  advected = ub
               end if
            end if

            conv = conv + un_area * advected

            ! Viscous diffusion Term
            if (nb > 0) then
               diff_face = half * (transport%nu(c) + transport%nu(nb))
            else
               diff_face = transport%nu(c)
            end if

            diff = diff + diff_face * mesh%faces(f)%area * (ub - u(:, c)) / dist
         end do

         ! Assemble total RHS
         local_rhs(:, c) = -conv / mesh%cells(c)%volume + &
                            diff / mesh%cells(c)%volume + &
                            params%body_force - gradp / params%rho
      end do
   end subroutine compute_momentum_rhs