! File: module_standard_tpo_formulation.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'. * !************************************************************************ MODULE standard_tpo_formulation !*********************************************************************** ! !# This module contains the definition of ABSTRACT TYPE [[tpo]] ! !*********************************************************************** USE id_base, ONLY: idbase USE sph_particles, ONLY: particles USE mesh_refinement, ONLY: grid_function_scalar, grid_function, level USE timing, ONLY: timer USE utility, ONLY: ios, err_msg, perc, creturn, run_id, & test_status, show_progress IMPLICIT NONE !******************************************************* ! * ! Definition of abstract TYPE tpo * ! * ! Abstract TYPE for a \(3+1\) formulation of the * ! Einstein equations. It imports the ID on the gravity * ! grid, in the standard 3+1 formulation, and defines * ! DEFERRED PROCEDURES to be implemented in the derived * ! TYPES of the actual 3+1 formulations (for example, * ! the BSSN formulation). * ! * !******************************************************* TYPE, ABSTRACT:: tpo !! ABSTRACT TYPE representing the standard \(3+1\) formulation of the |ee| INTEGER:: tpo_id_number !! Negative integer that identifies the [[tpo]] object INTEGER:: n_matter !! Number of matter objects in the physical system INTEGER, DIMENSION(:), ALLOCATABLE:: npoints_xaxis !# Array containing the number of mesh points of the highest-resolution ! refinement level across the x-axis-size of the matter objects ! !-- Mesh variables ! INTEGER:: nlevels !! Number of refinement levels TYPE(level), DIMENSION(:), ALLOCATABLE :: levels !! Array containing the information on each refinement level TYPE(grid_function):: coords !! Grid function storing the Cartesian coordinates TYPE(grid_function_scalar):: rad_coord !! Grid scalar function storing the radial coordinates of each grid point ! !-- ADM fields ! TYPE(grid_function_scalar):: lapse !! Grid scalar function storing the lapse function TYPE(grid_function):: shift_u !! Grid function storing the shift vector TYPE(grid_function):: g_phys3_ll !! Grid function storing the spatial metric TYPE(grid_function):: K_phys3_ll !! Grid function storing the extrinsic curvature ! !-- Constraint violations ! TYPE(grid_function_scalar):: HC !# Grid scalar function storing the Hamiltonian constraint (violations) ! computed using the ID on the mesh TYPE(grid_function_scalar):: HC_parts !# Grid scalar function storing the Hamiltonian constraint (violations) ! computed using the stress-energy tensor mapped from the particles to the ! mesh TYPE(grid_function_scalar):: rho !# Grid scalar function storing the matter source in the Hamiltonian ! constraint computed using the ID on the mesh, multiplied by \(16\pi\) TYPE(grid_function):: S !# Grid function storing the matter source in the momentum ! constraint computed using the ID on the mesh, multiplied by \(8\pi\) TYPE(grid_function):: MC !# Grid function storing the momentum constraint (violations) ! computed using the ID on the mesh TYPE(grid_function):: MC_parts !# Grid function storing the momentum constraint (violations) ! computed using the stress-energy tensor mapped from the particles to the ! mesh TYPE(grid_function_scalar):: rho_parts !# Grid scalar function storing the matter source in the Hamiltonian ! constraint computed using the stress-energy tensor mapped from the ! particles to the mesh, multiplied by \(16\pi\) TYPE(grid_function):: S_parts !# Grid function storing the matter source in the momentum ! constraint computed using the stress-energy tensor mapped from the ! particles to the mesh, multiplied by \(8\pi\) ! !-- Norms of constraint violations ! DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_l2 !! \(\ell_2\) norm of the Hamiltonian constraint computed on the mesh DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_parts_l2 !# \(\ell_2\) norm of the Hamiltonian constraint computed on the mesh, using ! the stress-energy tensor mapped from the particles DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_loo !# \(\ell_\infty\) norm of the Hamiltonian constraint computed on the mesh ! (i.e., its maximum) DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_parts_loo !# \(\ell_\infty\) norm of the Hamiltonian constraint computed on the mesh, ! using the stress-energy tensor mapped from the particles ! (i.e., its maximum) DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_l2 !! \(\ell_2\) norm of the momentum constraint computed on the mesh DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_parts_l2 !# \(\ell_2\) norm of the Hamiltonian constraint computed on the mesh, using ! the stress-energy tensor mapped from the particles DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_loo !# \(\ell_\infty\) norm of the momentum constraint computed on the mesh ! (i.e., its maximum) DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_parts_loo !# \(\ell_\infty\) norm of the momentum constraint computed on the mesh, ! using the stress-energy tensor mapped from the particles ! (i.e., its maximum) DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_int !# Integral the Hamiltonian constraint computed using the stress-energy ! tensor mapped from the particles DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: HC_parts_int !! Integral the Hamiltonian constraint computed on the mesh DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_int !! Integral of the momentum constraint computed on the mesh DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: MC_parts_int !# Integral of the momentum constraint computed using the stress-energy ! tensor mapped from the particles ! !-- Steering variables ! ! Variables to decide if and how to export the constraints INTEGER, PUBLIC:: cons_step !! Constraint violations are printed to file every cons_step-th grid point ! along each Cartesian direction LOGICAL, PUBLIC:: export_constraints !# `.TRUE.` if the constraint violations are to be printed to file, ! `.FALSE.` otherwise LOGICAL, PUBLIC:: export_constraints_details !# `.TRUE.` if the points at which the constraints violations are within the ! intervals \([0,10^{-7}],[10^3,\infty],[10^n,10^m]\), with ! \(n\in\{-7,2\}\) and \(m=n+1\), are to be printed to file; ! `.FALSE.` otherwise LOGICAL, PUBLIC:: export_constraints_xy !! `.TRUE.` if only the constrain violations on the \(xy\) plane are to be ! printed to file, `.FALSE.` otherwise LOGICAL, PUBLIC:: export_constraints_x !! `.TRUE.` if only the constrain violations on the \(x\) axis are to be ! printed to file, `.FALSE.` otherwise ! !-- Timers ! TYPE(timer):: grid_timer !# Timer that times how long it takes to set up the grid and allocate ! the grid functions TYPE(timer):: importer_timer !# Timer that times how long it takes to import the |lorene| ID ! on the mesh CONTAINS !-----------------------------------! !-- NON_OVERRIDABLE SUBROUTINES --! !-----------------------------------! PROCEDURE, NON_OVERRIDABLE:: setup_standard_tpo_variables !# Set up the refined mesh by reading the `gravity_grid_parameter.dat` ! parameter file, and read the standard 3+1 |id| using the given ! [[idbase]] object, on the refined mesh PROCEDURE, NON_OVERRIDABLE:: deallocate_standard_tpo_variables !! Deallocates memory for the standard 3+1 fields PROCEDURE, NON_OVERRIDABLE:: analyze_constraint !# Analyze a constraint (or an arbitrary scalar grid function) by ! examining its values at the refined mesh. Computes the number of mesh ! points at which the scalar grid function has values lying within ! predefined (hard-coded) intervals ! @todo complete this documentation entries with more details PROCEDURE, PUBLIC, NON_OVERRIDABLE:: test_recovery_m2p !# Computes the conserved variables from the physical ones, and vice versa, ! to test that the recovered physical variables are the same to those ! computed from the |id|. Uses the metric mapped from the grid to the ! particles. @todo add reference for recovery PROCEDURE, PUBLIC, NON_OVERRIDABLE:: compute_adm_momentum_fluid_m2p !# Computes an estimate of the \(\mathrm{ADM}\) linear momentum of the ! fluid using the |sph| hydro fields, and the spacetime metric mapped ! from the mesh ! @todo add reference PROCEDURE, NON_OVERRIDABLE:: print_summary !# Prints a summary about the features of the refined mesh !----------------------------! !-- DEFERRED SUBROUTINES --! !----------------------------! PROCEDURE(define_allocate_fields_interface), & DEFERRED:: define_allocate_fields !# Allocates memory for the fields specific to the formulation identified ! by an EXTENDED TYPE PROCEDURE(deallocate_fields_interface), DEFERRED:: deallocate_fields !# Deallocates memory for the fields specific to the formulation identified ! by an EXTENDED TYPE GENERIC, PUBLIC:: compute_and_print_tpo_constraints => & compute_and_print_tpo_constraints_grid, & compute_and_print_tpo_constraints_particles !# Overloaded PROCEDURE to compute the constraints using only the |id| ! on the refined mesh, or the spacetime |id| on the refined mesh and ! the hydrodynamical |id| mapped from the particles to the refined mesh PROCEDURE(compute_and_print_tpo_constraints_grid_interface), & DEFERRED:: compute_and_print_tpo_constraints_grid !# Computes the constraints specific to the formulation identified by an ! EXTENDED TYPE, using the full |id| on the refined mesh PROCEDURE(compute_and_print_tpo_constraints_particles_interface), & DEFERRED:: compute_and_print_tpo_constraints_particles !# Computes the constraints specific to the formulation identified by an ! EXTENDED TYPE, using the |bssn| |id| on the refined ! mesh and the hydrodynamical |id| mapped from the particles to the mesh PROCEDURE(compute_and_print_tpo_variables_interface), PUBLIC, & DEFERRED:: compute_and_print_tpo_variables !# Compute the fields specific to the formulation identified by an ! EXTENDED TYPE, starting from the standard 3+1 fields PROCEDURE(print_formatted_id_tpo_variables_interface), PUBLIC, & DEFERRED:: print_formatted_id_tpo_variables !! Prints the spacetime |id| to a formatted file !-----------------! !-- FUNCTIONS --! !-----------------! PROCEDURE:: abs_values_in PROCEDURE, PUBLIC:: get_grid_point PROCEDURE, PUBLIC:: get_nlevels PROCEDURE, PUBLIC:: get_levels PROCEDURE, PUBLIC:: get_dx PROCEDURE, PUBLIC:: get_dy PROCEDURE, PUBLIC:: get_dz PROCEDURE, PUBLIC:: get_ngrid_x PROCEDURE, PUBLIC:: get_ngrid_y PROCEDURE, PUBLIC:: get_ngrid_z PROCEDURE, PUBLIC:: get_xR PROCEDURE, PUBLIC:: get_yR PROCEDURE, PUBLIC:: get_zR PROCEDURE, PUBLIC:: get_HC PROCEDURE, PUBLIC:: get_MC PROCEDURE, PUBLIC:: get_HC_parts PROCEDURE, PUBLIC:: get_MC_parts END TYPE tpo ! !-- Interface of the cores of the constructors and destructors of TYPES !-- derived from tpo !-- Their implementations are in submodule tpo_methods.f90 ! INTERFACE MODULE SUBROUTINE setup_standard_tpo_variables( tpof, id, dx, dy, dz ) CLASS(idbase), INTENT(INOUT) :: id CLASS(tpo), INTENT(INOUT) :: tpof DOUBLE PRECISION, INTENT(IN), OPTIONAL:: dx, dy, dz END SUBROUTINE setup_standard_tpo_variables MODULE SUBROUTINE deallocate_standard_tpo_variables( this ) CLASS(tpo), INTENT(INOUT):: this END SUBROUTINE deallocate_standard_tpo_variables END INTERFACE ! !-- Interface of the methods of TYPES derived from tpo !-- Their implementations are in submodule tpo_methods.f90 ! INTERFACE MODULE SUBROUTINE analyze_constraint( this, & l, & constraint, & name_constraint, & unit_logfile, & name_analysis, & l2_norm, & loo_norm, & integral, & source ) CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: constraint CHARACTER( LEN= * ), INTENT(IN) :: name_constraint CHARACTER( LEN= * ), INTENT(IN) :: name_analysis INTEGER, INTENT(IN) :: unit_logfile DOUBLE PRECISION, INTENT(OUT):: l2_norm DOUBLE PRECISION, INTENT(OUT):: loo_norm DOUBLE PRECISION, INTENT(OUT):: integral DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN), OPTIONAL:: source END SUBROUTINE analyze_constraint MODULE SUBROUTINE test_recovery_m2p( this, parts, namefile ) !# Tests the recovery. Computes the conserved variables from the physical ! ones, and then the physical ones from the conserved ones. It then ! compares the variables computed with the recovery PROCEDURES, with ! those computed with |sphincsid|. Uses the mesh-2-particle. ! @todo add reference for recovery CLASS(tpo), INTENT(IN) :: this !! [[tpo]] object which this PROCEDURE is a member of CLASS(particles), INTENT(IN) :: parts !# [[sph_particles:particles]] object used to map the metric ! from the mesh to the particles, and to call the recovery procedures CHARACTER( LEN= * ), INTENT(INOUT):: namefile !! Name of the formatted file where the data is printed END SUBROUTINE test_recovery_m2p MODULE SUBROUTINE compute_adm_momentum_fluid_m2p( this, parts, adm_mom ) !# Computes an estimate of the \(\mathrm{ADM}\) linear momentum of the ! fluid using the |sph| hydro fields, and the spacetime metric mapped ! from the mesh ! @todo add reference CLASS(tpo), INTENT(IN) :: this !! [[tpo]] object which this PROCEDURE is a member of CLASS(particles), INTENT(IN) :: parts !# [[sph_particles:particles]] object used to map the metric ! from the mesh to the particles, and to call the recovery procedures DOUBLE PRECISION, DIMENSION(3), INTENT(OUT):: adm_mom !# ADM linear momentum of the fluid computed using the metric mapped ! with the mesh-to-particle mapping END SUBROUTINE compute_adm_momentum_fluid_m2p MODULE SUBROUTINE print_summary( this )!, filename ) !# Prints a summary of the properties of the refined mesh, ! and optionally, to a formatted file whose name ! is given as the optional argument `filename` CLASS(tpo), INTENT(INOUT):: this !CHARACTER( LEN= * ), INTENT( INOUT ), OPTIONAL:: filename !! Name of the formatted file to print the summary to END SUBROUTINE print_summary MODULE SUBROUTINE abs_values_in( this, lower_bound, upper_bound, & constraint, l, & export, unit_analysis, cnt ) CLASS(tpo), INTENT(INOUT):: this DOUBLE PRECISION, INTENT(IN) :: lower_bound, & upper_bound DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: constraint INTEGER, INTENT(IN) :: l LOGICAL, INTENT(IN) :: export INTEGER, INTENT(IN) :: unit_analysis INTEGER, INTENT(OUT) :: cnt END SUBROUTINE abs_values_in MODULE FUNCTION get_grid_point( this, i, j, k, l ) RESULT( grid_point ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: i, j, k, l ! Result DOUBLE PRECISION, DIMENSION(3) :: grid_point END FUNCTION get_grid_point MODULE FUNCTION get_nlevels( this ) RESULT( nlevels ) ! Arguments CLASS(tpo), INTENT(INOUT):: this ! Result INTEGER:: nlevels END FUNCTION get_nlevels MODULE FUNCTION get_levels( this ) RESULT( levels ) ! Arguments CLASS(tpo), INTENT(INOUT):: this ! Result TYPE(level), DIMENSION(:), ALLOCATABLE:: levels END FUNCTION get_levels MODULE FUNCTION get_dx( this, l ) RESULT( dx ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: dx END FUNCTION get_dx MODULE FUNCTION get_dy( this, l ) RESULT( dy ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: dy END FUNCTION get_dy MODULE FUNCTION get_dz( this, l ) RESULT( dz ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: dz END FUNCTION get_dz MODULE FUNCTION get_ngrid_x( this, l ) RESULT( ngrid_x ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result INTEGER:: ngrid_x END FUNCTION get_ngrid_x MODULE FUNCTION get_ngrid_y( this, l ) RESULT( ngrid_y ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result INTEGER:: ngrid_y END FUNCTION get_ngrid_y MODULE FUNCTION get_ngrid_z( this, l ) RESULT( ngrid_z ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result INTEGER:: ngrid_z END FUNCTION get_ngrid_z MODULE FUNCTION get_xR( this, l ) RESULT( xR ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: xR END FUNCTION get_xR MODULE FUNCTION get_yR( this, l ) RESULT( yR ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: yR END FUNCTION get_yR MODULE FUNCTION get_zR( this, l ) RESULT( zR ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: l ! Result DOUBLE PRECISION:: zR END FUNCTION get_zR MODULE FUNCTION get_HC( this, i, j, k, l ) RESULT( HC_value ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: i, j, k, l ! Result DOUBLE PRECISION :: HC_value END FUNCTION get_HC MODULE FUNCTION get_MC( this, i, j, k, l ) RESULT( MC_value ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: i, j, k, l ! Result DOUBLE PRECISION, DIMENSION(3) :: MC_value END FUNCTION get_MC MODULE FUNCTION get_HC_parts( this, i, j, k, l ) RESULT( HC_value ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: i, j, k, l ! Result DOUBLE PRECISION :: HC_value END FUNCTION get_HC_parts MODULE FUNCTION get_MC_parts( this, i, j, k, l ) RESULT( MC_value ) ! Arguments CLASS(tpo), INTENT(INOUT):: this INTEGER, INTENT(IN) :: i, j, k, l ! Result DOUBLE PRECISION, DIMENSION(3) :: MC_value END FUNCTION get_MC_parts END INTERFACE ! !-- Interfaces of the deferred methods of TYPE tpo !-- Their implementations are deferred to derived TYPES ! ABSTRACT INTERFACE SUBROUTINE define_allocate_fields_interface( this ) IMPORT:: tpo CLASS(tpo), INTENT(INOUT):: this END SUBROUTINE define_allocate_fields_interface SUBROUTINE compute_and_print_tpo_variables_interface( this, namefile ) IMPORT:: tpo CLASS(tpo), INTENT(INOUT) :: this CHARACTER( LEN= * ), INTENT(INOUT), OPTIONAL :: namefile END SUBROUTINE compute_and_print_tpo_variables_interface SUBROUTINE print_formatted_id_tpo_variables_interface & ( this, points, namefile ) IMPORT:: tpo CLASS(tpo), INTENT(INOUT) :: this DOUBLE PRECISION, DIMENSION(:,:,:,:), INTENT(IN), OPTIONAL:: points CHARACTER( LEN= * ), INTENT(INOUT), OPTIONAL :: namefile END SUBROUTINE print_formatted_id_tpo_variables_interface SUBROUTINE compute_and_print_tpo_constraints_grid_interface & ( this, id, namefile, name_logfile, points ) IMPORT:: tpo IMPORT:: idbase CLASS(tpo), INTENT(INOUT), TARGET:: this CLASS(idbase), INTENT(INOUT):: id CHARACTER( LEN= * ), INTENT(INOUT):: namefile CHARACTER( LEN= * ), INTENT(INOUT):: name_logfile DOUBLE PRECISION, DIMENSION(:,:,:,:), INTENT(IN), TARGET, OPTIONAL:: & points END SUBROUTINE compute_and_print_tpo_constraints_grid_interface SUBROUTINE compute_and_print_tpo_constraints_particles_interface & ( this, parts_obj, namefile, name_logfile, points ) IMPORT:: tpo IMPORT:: particles CLASS(tpo), INTENT(INOUT), TARGET:: this CLASS(particles), INTENT(INOUT):: parts_obj CHARACTER( LEN= * ), INTENT(INOUT):: namefile CHARACTER( LEN= * ), INTENT(INOUT):: name_logfile DOUBLE PRECISION, DIMENSION(:,:,:,:), INTENT(IN), TARGET, OPTIONAL:: & points END SUBROUTINE compute_and_print_tpo_constraints_particles_interface SUBROUTINE deallocate_fields_interface( this ) IMPORT:: tpo CLASS(tpo), INTENT(INOUT):: this END SUBROUTINE deallocate_fields_interface END INTERFACE END MODULE standard_tpo_formulation