! File: submodule_bssn_formulation_ricci.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 (bssn_formulation) ricci !************************************************ ! !# Implementation of the method of TYPE [[bssn]] ! that computes the Ricci tensor and scalar ! ! FT 10.02.2022 ! !************************************************ IMPLICIT NONE CONTAINS !-------------------! !-- SUBROUTINES --! !-------------------! MODULE PROCEDURE compute_ricci !************************************************ ! !# Computes the Ricci tensor and the Ricci ! scalar on the mesh ! ! FT 10.02.2022 ! !************************************************ USE utility, ONLY: zero, one, two, three, ten, Msun_geo USE tensor, ONLY: jx, jy, jz, jxx, jxy, jxz, jyy, jyz, jzz USE mesh_refinement, ONLY: allocate_grid_function, levels, nlevels IMPLICIT NONE INTEGER:: l DOUBLE PRECISION:: max_ricci, min_lapse INTEGER, DIMENSION(3) :: imin, imax CHARACTER(LEN= 2):: tpo_id WRITE( tpo_id, "(I2)" ) this% tpo_id_number ALLOCATE ( levels( this% nlevels ), STAT=ios ) IF( ios > 0 )THEN PRINT*,'...allocation error for levels' STOP ENDIF levels = this% levels nlevels= this% nlevels CALL allocate_grid_function( this% Ricci_ll, "Ricci_ll_id"//tpo_id, 6 ) CALL allocate_grid_function & ( this% Ricci_scalar, "Ricci_scalar_id"//tpo_id, 1 ) DO l= 1, this% nlevels, 1 ASSOCIATE( lapse => this% lapse% levels(l)% var, & shift_u => this% shift_u% levels(l)% var, & phi => this% phi% levels(l)% var, & trK => this% trK% levels(l)% var, & g_BSSN3_ll => this% g_BSSN3_ll% levels(l)% var, & A_BSSN3_ll => this% A_BSSN3_ll% levels(l)% var, & Gamma_u => this% Gamma_u% levels(l)% var, & Ricci_ll => this% Ricci_ll% levels(l)% var, & Ricci_scalar => this% Ricci_scalar% levels(l)% var & ) imin(1)= this% levels(l)% nghost_x imin(2)= this% levels(l)% nghost_y imin(3)= this% levels(l)% nghost_z imax(1)= this% get_ngrid_x(l) - this% levels(l)% nghost_x - 1 imax(2)= this% get_ngrid_y(l) - this% levels(l)% nghost_y - 1 imax(3)= this% get_ngrid_z(l) - this% levels(l)% nghost_z - 1 Ricci_ll = zero Ricci_scalar= zero CALL bssn_ricci_interior( & this% get_ngrid_x(l), this% get_ngrid_y(l), this% get_ngrid_z(l), & imin, imax, & this% get_dx(l), this% get_dy(l), this% get_dz(l), & g_BSSN3_ll(:,:,:,jxx), g_BSSN3_ll(:,:,:,jxy), & g_BSSN3_ll(:,:,:,jxz), g_BSSN3_ll(:,:,:,jyy), & g_BSSN3_ll(:,:,:,jyz), g_BSSN3_ll(:,:,:,jzz), & A_BSSN3_ll(:,:,:,jxx), A_BSSN3_ll(:,:,:,jxy), & A_BSSN3_ll(:,:,:,jxz), A_BSSN3_ll(:,:,:,jyy), & A_BSSN3_ll(:,:,:,jyz), A_BSSN3_ll(:,:,:,jzz), & trK(:,:,:), phi(:,:,:), & Gamma_u(:,:,:,jx), & Gamma_u(:,:,:,jy), & Gamma_u(:,:,:,jz), & Ricci_ll(:,:,:,jxx), Ricci_ll(:,:,:,jxy), Ricci_ll(:,:,:,jxz), & Ricci_ll(:,:,:,jyy), Ricci_ll(:,:,:,jyz), Ricci_ll(:,:,:,jzz), & Ricci_scalar(:,:,:) ) END ASSOCIATE ENDDO PRINT *, " * Ricci tensor and scalar computed." PRINT * ! From Wikipedia: ! The sectional curvature of an n-sphere of radius r is K = 1/r^2. ! Hence the scalar curvature is R = n(n − 1)/r^2. ! TODO: Find real reference for this. max_ricci= zero min_lapse= HUGE(one) ASSOCIATE( Ricci_scalar => this% Ricci_scalar% levels(this% nlevels)% var, & lapse => this% lapse% levels(this% nlevels)% var ) !IF( MAXVAL( ABS(Ricci_scalar) ) > max_ricci )THEN max_ricci= MAXVAL( ABS(Ricci_scalar) ) !ENDIF !IF( MINVAL( ABS(lapse) ) < min_lapse )THEN ! ! min_lapse= MINVAL( ABS(lapse) ) ! !ENDIF END ASSOCIATE PRINT *, "* Maximum value of the Ricci scalar over the finest refinement ", & "level= ", max_ricci, "Msun_geo^{-2}" PRINT *, "* Radius of a 3-sphere with the same Ricci scalar=", & SQRT(three*two/max_ricci), "Msun_geo=", & SQRT(three*two/max_ricci)*Msun_geo*ten*ten*ten, "m" !PRINT *, min_lapse END PROCEDURE compute_ricci END SUBMODULE ricci