CalculateTendency_ECDGModel2D_t Subroutine

  subroutine CalculateTendency_ECDGModel2D_t(this)
    !! Computes du/dt = source - EC-DG flux divergence.
    !!
    !! Volume term : MappedDivergence of twoPointFlux (EC split-form sum, /J)
    !! Surface term: (1/J) * M^{-1} B^T f_Riemann  (weak-form Riemann flux)
    implicit none
    class(ECDGModel2D_t),intent(inout) :: this
    ! Local
    integer :: i,j,iEl,iVar
    real(prec) :: bM1,bM2,qwi,qwj,jac

    call this%solution%BoundaryInterp()
    call this%solution%SideExchange(this%mesh)

    call this%PreTendency()
    call this%SetBoundaryCondition()

    if(this%gradient_enabled) then
      call this%CalculateSolutionGradient()
      call this%SetGradientBoundaryCondition()
      call this%solutionGradient%AverageSides()
    endif

    call this%SourceMethod()
    call this%BoundaryFlux() ! fills flux%boundaryNormal with Riemann fluxes

    ! EC volume divergence (includes J weighting via MappedDivergence)
    call this%TwoPointFluxMethod()
    call this%twoPointFlux%UpdateDevice()
    call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior)

    ! Add weak-form DG surface term: (1/J) * M^{-1} B^T f_Riemann
    ! Boundary side ordering: 1=South, 2=East, 3=North, 4=West
    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)

      bM1 = this%solution%interp%bMatrix(i,2) ! right/north boundary value at node i
      bM2 = this%solution%interp%bMatrix(i,1) ! left/south boundary value at node i
      qwi = this%solution%interp%qWeights(i)
      qwj = this%solution%interp%qWeights(j)
      jac = this%geometry%J%interior(i,j,iEl,1)

      ! East (side 2) and West (side 4) contributions — vary in xi^1 (i-direction)
      this%fluxDivergence%interior(i,j,iEl,iVar) = &
        this%fluxDivergence%interior(i,j,iEl,iVar)+ &
        (bM1*this%flux%boundaryNormal(j,2,iEl,iVar)+ &
         bM2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ &
        (this%solution%interp%bMatrix(j,2)*this%flux%boundaryNormal(i,3,iEl,iVar)+ &
         this%solution%interp%bMatrix(j,1)*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac)

    enddo

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)

      this%dSdt%interior(i,j,iEl,iVar) = &
        this%source%interior(i,j,iEl,iVar)- &
        this%fluxDivergence%interior(i,j,iEl,iVar)

    enddo

  endsubroutine CalculateTendency_ECDGModel2D_t