! File:         submodule_standard_tpo_formulation_sph_adm_variables.f90
! Authors:      Francesco Torsello (FT)
!************************************************************************
! Copyright (C) 2020-2023 Francesco Torsello                            *
!                                                                       *
! This file is part of SPHINCS_ID                                       *
!                                                                       *
! SPHINCS_ID is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by  *
! the Free Software Foundation, either version 3 of the License, or     *
! (at your option) any later version.                                   *
!                                                                       *
! SPHINCS_ID is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of        *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
! GNU General Public License for more details.                          *
!                                                                       *
! You should have received a copy of the GNU General Public License     *
! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>.   *
! The copy of the GNU General Public License should be in the file      *
! 'COPYING'.                                                            *
!************************************************************************

SUBMODULE (standard_tpo_formulation) sph_adm_variables

  !************************************************
  !
  !# This SUBMODULE contains the implementation
  !  of the methods that compute an estimate of
  !  the ADM variables on the |sph| fluid, using
  !  the metric mapped from the mesh to the particles
  !
  !  FT 12.04.2020
  !
  !************************************************


  IMPLICIT NONE


  CONTAINS


  MODULE PROCEDURE compute_adm_momentum_fluid_m2p

    !************************************************
    !
    !# Computes an estimate of the \(\mathrm{ADM}\)
    !  linear momentum using the |sph| fields
    !  on the particles
    !  @todo add reference
    !
    !  FT 12.04.2022
    !
    !************************************************

    USE tensor,               ONLY: jx, jy, jz, n_sym4x4, lower_index_4vector
    USE constants,            ONLY: amu, MSun
    USE utility,              ONLY: compute_g4, determinant_sym4x4, &
                                    spatial_vector_norm_sym3x3, &
                                    compute_tpo_metric, is_finite_number, &
                                    zero, one, two

    USE metric_on_particles,  ONLY: allocate_metric_on_particles, &
                                    deallocate_metric_on_particles, &
                                    get_metric_on_particles, g4_ll

    USE ADM_refine,           ONLY: lapse, shift_u, g_phys3_ll, &
                                    allocate_ADM, deallocate_ADM
    USE BSSN_refine,          ONLY: allocate_BSSN, deallocate_BSSN
    USE mesh_refinement,      ONLY: nlevels, levels, rad_coord, coords, &
                                    allocate_grid_function, &
                                    deallocate_grid_function
    USE gradient,             ONLY: allocate_gradient, deallocate_gradient
    USE McLachlan_refine,     ONLY: allocate_Ztmp, deallocate_Ztmp
    USE alive_flag,           ONLY: alive
    USE input_output,         ONLY: read_options
    USE units,                ONLY: set_units
    USE sph_variables,        ONLY: npart, allocate_SPH_memory, &
                                    deallocate_SPH_memory

    IMPLICIT NONE

    INTEGER, PARAMETER:: unit_recovery= 34956

    INTEGER:: a, j, l, i_matter

    DOUBLE PRECISION:: det, shift_norm2, adm_mom_element

    !DOUBLE PRECISION, DIMENSION(n_sym4x4,npart)  :: g4
    DOUBLE PRECISION, DIMENSION(0:3)             :: v_u
    DOUBLE PRECISION, DIMENSION(0:3,parts% get_npart()):: v_l

    DOUBLE PRECISION              :: lapse_loc
    DOUBLE PRECISION, DIMENSION(3):: shift
    DOUBLE PRECISION, DIMENSION(6):: g3

    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: pos
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: vel_loc
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nu_loc
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nlrf_loc
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: u_loc
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: pr_loc
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: theta_loc

    DOUBLE PRECISION, DIMENSION(parts% get_n_matter(),3):: adm_mom_i

    LOGICAL, PARAMETER:: debug= .FALSE.

    npart    = parts% get_npart()
    pos      = parts% get_pos()
    vel_loc  = parts% get_vel()
    nu_loc   = parts% get_nu()
    nlrf_loc = parts% get_nlrf_sph()
    u_loc    = parts% get_u_sph()
    pr_loc   = parts% get_pressure_sph()
    theta_loc= parts% get_theta()

    ALLOCATE ( levels( this% nlevels ), STAT=ios )
    IF( ios > 0 )THEN
     PRINT*,'...allocation error for levels'
     STOP
    ENDIF
    nlevels= this% nlevels
    levels = this% levels
    coords = this% coords

    CALL allocate_ADM()
    CALL allocate_BSSN()

    ! Allocate temporary memory for time integration
    CALL allocate_Ztmp()

    ! Allocate memory for the derivatives of the ADM variables
    ! CALL allocate_GravityAcceleration()

    CALL allocate_grid_function( rad_coord, 'rad_coord', 1 )

    ! Initialize the stress-energy tensor to 0
    DO l= 1, this% nlevels, 1
      rad_coord%  levels(l)% var= this% rad_coord%  levels(l)% var
      g_phys3_ll% levels(l)% var= this% g_phys3_ll% levels(l)% var
      shift_u%    levels(l)% var= this% shift_u%    levels(l)% var
      lapse%      levels(l)% var= this% lapse%      levels(l)% var
    ENDDO

    CALL set_units('NSM')
    CALL read_options

    CALL allocate_SPH_memory

    ! flag that particles are 'alive'
    IF( .NOT.ALLOCATED( alive ) ) ALLOCATE( alive( npart ) )
    alive( 1:npart )= 1

    CALL allocate_gradient( npart )

    IF( ALLOCATED(g4_ll) )THEN
      DEALLOCATE(g4_ll)
    ENDIF
    CALL allocate_metric_on_particles(npart)

    PRINT *, " * Mapping metric from the grid to the particles..."
    PRINT *
    CALL get_metric_on_particles( npart, pos )

 !   adm_mom= zero
 !   !$OMP PARALLEL DO DEFAULT( NONE ) &
 !   !$OMP             SHARED( npart, parts, g4_ll, vel_loc, v_l, nu_loc, &
 !   !$OMP                     theta_loc, pr_loc, nlrf_loc, u_loc ) &
 !   !$OMP             PRIVATE( a, det, v_u, shift_norm2, j, g3, lapse_loc, &
 !   !$OMP                      shift ) &
 !   !$OMP             REDUCTION( +: adm_mom )
 !   DO a= 1, npart, 1
 !
 !     CALL compute_tpo_metric( g4_ll(:,a), lapse_loc, shift, g3 )
 !
 !     v_u= [one, vel_loc(1:3,a)]
 !     CALL lower_index_4vector( v_u, g4_ll(:,a), v_l(:,a) )
 !
 !     CALL spatial_vector_norm_sym3x3( g3, shift, shift_norm2 )
 !
 !     DO j= jx, jz, 1
 !
 !       adm_mom(j)= adm_mom(j) &
 !           - ( nu_loc(a)*amu/Msun )*( shift_norm2/(lapse_loc**two) - one ) &
 !             *theta_loc(a)*( one + u_loc(a) + pr_loc(a)/nlrf_loc(a) )*v_l(j,a)
 !
 !     ENDDO
 !
 !   ENDDO
 !   !$OMP END PARALLEL DO

    adm_mom= zero
    matter_objects_loop: DO i_matter= 1, parts% get_n_matter(), 1

      adm_mom_i(i_matter,:)= zero
      !$OMP PARALLEL DO DEFAULT( NONE ) &
      !$OMP             SHARED( parts, g4_ll, vel_loc, &
      !$OMP                     v_l, nu_loc, theta_loc, pr_loc, nlrf_loc, &
      !$OMP                     u_loc, i_matter ) &
      !$OMP             PRIVATE( a, det, v_u, shift_norm2, j, g3, lapse_loc, &
      !$OMP                      shift, adm_mom_element ) &
      !$OMP             REDUCTION( +: adm_mom_i )
      DO a= parts% get_npart_i(i_matter-1) + 1, &
            parts% get_npart_i(i_matter-1) + parts% get_npart_i(i_matter), 1

        CALL compute_tpo_metric( g4_ll(:,a), lapse_loc, shift, g3 )

        DO j= 1, 6, 1
          IF( .NOT.is_finite_number(g3(j)) )THEN
            PRINT *, "** ERROR! g3(", j, ") is not a finite number on ", &
                     "particle ", a, ", in SUBROUTINE ", &
                     "compute_adm_momentum_fluid_m2p!"
            PRINT *, " * Stopping..."
            STOP
          ENDIF
        ENDDO
        IF( .NOT.is_finite_number(lapse_loc) )THEN
          PRINT *, "** ERROR! lapse_loc is not a finite number on ", &
                   "particle ", a, ", in SUBROUTINE ", &
                   "compute_adm_momentum_fluid_m2p!"
          PRINT *, " * Stopping..."
          STOP
        ENDIF
        DO j= 1, 3, 1
          IF( .NOT.is_finite_number(shift(j)) )THEN
            PRINT *, "** ERROR! shift(", j, ") is not a finite number on ", &
                     "particle ", a, ", in SUBROUTINE ", &
                     "compute_adm_momentum_fluid_m2p!"
            PRINT *, " * Stopping..."
            STOP
          ENDIF
        ENDDO

        v_u= [one, vel_loc(1:3,a)]
        CALL lower_index_4vector( v_u, g4_ll(:,a), v_l(:,a) )
        DO j= 0, 3, 1
          IF( .NOT.is_finite_number(v_l(j,a)) )THEN
            PRINT *, "** ERROR! v_l(", j, ") is not a finite number on ", &
                     "particle ", a, ", in SUBROUTINE ", &
                     "compute_adm_momentum_fluid_m2p!"
            PRINT *, " * Stopping..."
            STOP
          ENDIF
        ENDDO

        CALL spatial_vector_norm_sym3x3( g3, shift, shift_norm2 )
        IF( .NOT.is_finite_number(shift_norm2) )THEN
          PRINT *, "** ERROR! shift_norm2 is not a finite number on ", &
                   "particle ", a, ", in SUBROUTINE ", &
                   "compute_adm_momentum_fluid_m2p!"
          PRINT *, " * Stopping..."
          STOP
        ENDIF

        DO j= jx, jz, 1

          adm_mom_element= ( nu_loc(a)*amu/Msun ) &
            *theta_loc(a)*( one + u_loc(a) + pr_loc(a)/nlrf_loc(a) )*v_l(j,a)

          IF( is_finite_number(adm_mom_element) )THEN
            adm_mom_i(i_matter,j)= adm_mom_i(i_matter,j) + adm_mom_element
          ELSE
            CYCLE
          ENDIF

        ENDDO

      ENDDO
      !$OMP END PARALLEL DO

      PRINT *, "   Estimate of the ADM momentum of matter object", i_matter, &
               " computed using the SPH hydro fields and", &
               " the metric mapped with mesh-to-particle mapping= "
      PRINT *, "   (", adm_mom_i(i_matter,1), ","
      PRINT *, "    ", adm_mom_i(i_matter,2), ","
      PRINT *, "    ", adm_mom_i(i_matter,3), ") Msun*c"
      PRINT *

      adm_mom= adm_mom + adm_mom_i(i_matter,:)

    ENDDO matter_objects_loop

    PRINT *, "   Estimate of the ADM momentum of the fluid computed using ",&
             "the SPH hydro fields and the metric mapped with ", &
             "mesh-to-particle mapping= "
    PRINT *, "   (", adm_mom(1), ","
    PRINT *, "    ", adm_mom(2), ","
    PRINT *, "    ", adm_mom(3), ") Msun*c"
    PRINT *

    CALL deallocate_metric_on_particles
    CALL deallocate_gradient
    DEALLOCATE(alive)
    CALL deallocate_sph_memory

    CALL deallocate_grid_function ( rad_coord, 'rad_coord' )
    CALL deallocate_ADM()
    CALL deallocate_Ztmp()
    !CALL deallocate_GravityAcceleration()
    CALL deallocate_BSSN()
    DEALLOCATE( levels )
    DEALLOCATE( pos     )
    DEALLOCATE( vel_loc )
    DEALLOCATE( nu_loc    )
    DEALLOCATE( nlrf_loc  )
    DEALLOCATE( u_loc     )
    DEALLOCATE( pr_loc    )
    DEALLOCATE( theta_loc )

  END PROCEDURE compute_adm_momentum_fluid_m2p


END SUBMODULE sph_adm_variables