convergence_test Program

Uses

  • program~~convergence_test~~UsesGraph program~convergence_test convergence_test IFPORT IFPORT program~convergence_test->IFPORT iso_fortran_env iso_fortran_env program~convergence_test->iso_fortran_env module~bssn_formulation bssn_formulation program~convergence_test->module~bssn_formulation module~cauchy_convergence_test cauchy_convergence_test program~convergence_test->module~cauchy_convergence_test module~id_base id_base program~convergence_test->module~id_base module~sph_particles sph_particles program~convergence_test->module~sph_particles module~sphincs_id_fuka sphincs_id_fuka program~convergence_test->module~sphincs_id_fuka module~sphincs_id_full sphincs_id_full program~convergence_test->module~sphincs_id_full module~sphincs_id_interpolate sphincs_id_interpolate program~convergence_test->module~sphincs_id_interpolate module~sphincs_id_lorene sphincs_id_lorene program~convergence_test->module~sphincs_id_lorene module~standard_tpo_formulation standard_tpo_formulation program~convergence_test->module~standard_tpo_formulation module~utility utility program~convergence_test->module~utility tensor tensor program~convergence_test->tensor timing timing program~convergence_test->timing module~bssn_formulation->module~id_base module~bssn_formulation->module~sph_particles module~bssn_formulation->module~standard_tpo_formulation module~bssn_formulation->module~utility module~bssn_formulation->timing mesh_refinement mesh_refinement module~bssn_formulation->mesh_refinement module~cauchy_convergence_test->module~standard_tpo_formulation module~cauchy_convergence_test->module~utility module~id_base->module~utility module~id_base->timing module~sph_particles->module~id_base module~sph_particles->module~utility module~sph_particles->timing module~sphincs_id_fuka->module~id_base module~bns_fuka bns_fuka module~sphincs_id_fuka->module~bns_fuka module~ejecta_generic ejecta_generic module~sphincs_id_fuka->module~ejecta_generic module~sphincs_id_full->module~id_base module~sphincs_id_full->module~bns_fuka module~bns_lorene bns_lorene module~sphincs_id_full->module~bns_lorene module~diffstar_lorene diffstar_lorene module~sphincs_id_full->module~diffstar_lorene module~sphincs_id_full->module~ejecta_generic module~sphincs_id_interpolate->module~id_base module~sphincs_id_interpolate->module~ejecta_generic module~sphincs_id_lorene->module~id_base module~sphincs_id_lorene->module~bns_lorene module~sphincs_id_lorene->module~diffstar_lorene module~sphincs_id_lorene->module~ejecta_generic module~standard_tpo_formulation->module~id_base module~standard_tpo_formulation->module~sph_particles module~standard_tpo_formulation->module~utility module~standard_tpo_formulation->timing module~standard_tpo_formulation->mesh_refinement constants constants module~utility->constants matrix matrix module~utility->matrix module~bns_fuka->module~id_base module~bns_fuka->module~utility module~bns_fuka->timing module~bns_fuka->mesh_refinement iso_c_binding iso_c_binding module~bns_fuka->iso_c_binding module~bns_base bns_base module~bns_fuka->module~bns_base module~bns_lorene->module~id_base module~bns_lorene->module~utility module~bns_lorene->timing module~bns_lorene->iso_c_binding module~bns_lorene->module~bns_base module~diffstar_lorene->module~id_base module~diffstar_lorene->module~utility module~diffstar_lorene->timing module~diffstar_lorene->iso_c_binding module~diffstar_base diffstar_base module~diffstar_lorene->module~diffstar_base module~ejecta_generic->module~id_base module~ejecta_generic->module~utility module~bns_base->module~id_base module~bns_base->module~utility module~diffstar_base->module~id_base module~diffstar_base->module~utility

Make a Cauchy convergence test to check the validity of the \(\mathrm{ID}\) set up by \(\texttt{SPHINCS_ID}\).

FT 8.12.2020


Calls

program~~convergence_test~~CallsGraph program~convergence_test convergence_test interface~bssn bssn program~convergence_test->interface~bssn interface~find_shared_grid find_shared_grid program~convergence_test->interface~find_shared_grid interface~particles particles program~convergence_test->interface~particles interface~perform_cauchy_convergence_test perform_cauchy_convergence_test program~convergence_test->interface~perform_cauchy_convergence_test proc~allocate_idbase allocate_idbase program~convergence_test->proc~allocate_idbase proc~read_sphincs_id_parameters read_sphincs_id_parameters program~convergence_test->proc~read_sphincs_id_parameters interface~construct_bssn construct_bssn interface~bssn->interface~construct_bssn interface~construct_particles_bin construct_particles_bin interface~particles->interface~construct_particles_bin interface~construct_particles_std construct_particles_std interface~particles->interface~construct_particles_std proc~construct_bssn construct_bssn interface~construct_bssn->proc~construct_bssn proc~construct_particles_bin construct_particles_bin interface~construct_particles_bin->proc~construct_particles_bin proc~construct_particles_std construct_particles_std interface~construct_particles_std->proc~construct_particles_std deallocate_bssn deallocate_bssn proc~construct_bssn->deallocate_bssn initialize_bssn initialize_bssn proc~construct_bssn->initialize_bssn timer timer proc~construct_bssn->timer allocate_sph_memory allocate_sph_memory proc~construct_particles_bin->allocate_sph_memory deallocate_sph_memory deallocate_sph_memory proc~construct_particles_bin->deallocate_sph_memory proc~compute_g4 compute_g4 proc~construct_particles_bin->proc~compute_g4 proc~determinant_sym4x4 determinant_sym4x4 proc~construct_particles_bin->proc~determinant_sym4x4 proc~spacetime_vector_norm_sym4x4 spacetime_vector_norm_sym4x4 proc~construct_particles_bin->proc~spacetime_vector_norm_sym4x4 read_options read_options proc~construct_particles_bin->read_options select_eos_parameters select_eos_parameters proc~construct_particles_bin->select_eos_parameters set_units set_units proc~construct_particles_bin->set_units alive alive proc~construct_particles_std->alive com com proc~construct_particles_std->com interface~check_particle_positions check_particle_positions proc~construct_particles_std->interface~check_particle_positions interface~perform_apm perform_apm proc~construct_particles_std->interface~perform_apm none~check_eos check_eos proc~construct_particles_std->none~check_eos none~place_particles_on_ellipsoidal_surfaces place_particles_on_ellipsoidal_surfaces proc~construct_particles_std->none~place_particles_on_ellipsoidal_surfaces none~place_particles_on_lattices place_particles_on_lattices proc~construct_particles_std->none~place_particles_on_lattices none~read_particles_from_formatted_file read_particles_from_formatted_file proc~construct_particles_std->none~read_particles_from_formatted_file none~read_particles_options read_particles_options proc~construct_particles_std->none~read_particles_options none~reflect_particles_yz_plane reflect_particles_yz_plane proc~construct_particles_std->none~reflect_particles_yz_plane proc~scan_1d_array_for_nans scan_1d_array_for_nans proc~construct_particles_std->proc~scan_1d_array_for_nans proc~spatial_vector_norm_sym3x3 spatial_vector_norm_sym3x3 proc~construct_particles_std->proc~spatial_vector_norm_sym3x3 proc~spherical_from_cartesian spherical_from_cartesian proc~construct_particles_std->proc~spherical_from_cartesian proc~construct_particles_std->timer proc~check_particle_positions check_particle_positions interface~check_particle_positions->proc~check_particle_positions proc~perform_apm perform_apm interface~perform_apm->proc~perform_apm none~place_particles_on_ellipsoidal_surfaces->none~reflect_particles_yz_plane none~place_particles_on_lattices->none~reflect_particles_yz_plane none~read_particles_from_formatted_file->none~reflect_particles_yz_plane interface~impose_equatorial_plane_symmetry impose_equatorial_plane_symmetry none~read_particles_from_formatted_file->interface~impose_equatorial_plane_symmetry none~read_particles_options->read_options none~read_particles_options->set_units ktable ktable none~read_particles_options->ktable proc~is_finite_number is_finite_number proc~scan_1d_array_for_nans->proc~is_finite_number proc~impose_equatorial_plane_symmetry impose_equatorial_plane_symmetry interface~impose_equatorial_plane_symmetry->proc~impose_equatorial_plane_symmetry indexx indexx proc~check_particle_positions->indexx proc~perform_apm->allocate_sph_memory proc~perform_apm->com proc~perform_apm->deallocate_sph_memory proc~perform_apm->interface~check_particle_positions proc~perform_apm->proc~spherical_from_cartesian proc~perform_apm->timer proc~perform_apm->interface~impose_equatorial_plane_symmetry proc~perform_apm->proc~is_finite_number allocate_gradient allocate_gradient proc~perform_apm->allocate_gradient allocate_metric_on_particles allocate_metric_on_particles proc~perform_apm->allocate_metric_on_particles allocate_rcb_tree_memory_3d allocate_rcb_tree_memory_3d proc~perform_apm->allocate_rcb_tree_memory_3d assign_h assign_h proc~perform_apm->assign_h deallocate_gradient deallocate_gradient proc~perform_apm->deallocate_gradient deallocate_metric_on_particles deallocate_metric_on_particles proc~perform_apm->deallocate_metric_on_particles deallocate_rcb_tree_memory_3d deallocate_rcb_tree_memory_3d proc~perform_apm->deallocate_rcb_tree_memory_3d density density proc~perform_apm->density density_loop density_loop proc~perform_apm->density_loop exact_nei_tree_update exact_nei_tree_update proc~perform_apm->exact_nei_tree_update h h proc~perform_apm->h interface~correct_center_of_mass correct_center_of_mass proc~perform_apm->interface~correct_center_of_mass iorig iorig proc~perform_apm->iorig none~allocate_apm_fields allocate_apm_fields proc~perform_apm->none~allocate_apm_fields none~discard_atmosphere discard_atmosphere proc~perform_apm->none~discard_atmosphere none~dump_apm_pos dump_apm_pos proc~perform_apm->none~dump_apm_pos none~get_nstar_id_atm get_nstar_id_atm proc~perform_apm->none~get_nstar_id_atm none~place_and_print_ghost_particles place_and_print_ghost_particles proc~perform_apm->none~place_and_print_ghost_particles none~read_pressure_id read_pressure_id proc~perform_apm->none~read_pressure_id none~reallocate_output_fields reallocate_output_fields proc~perform_apm->none~reallocate_output_fields none~validate_position_final validate_position_final proc~perform_apm->none~validate_position_final nu nu proc~perform_apm->nu position_correction position_correction proc~perform_apm->position_correction proc~cartesian_from_spherical cartesian_from_spherical proc~perform_apm->proc~cartesian_from_spherical proc~find_h_backup find_h_backup proc~perform_apm->proc~find_h_backup proc~correct_center_of_mass correct_center_of_mass interface~correct_center_of_mass->proc~correct_center_of_mass none~discard_atmosphere->h none~discard_atmosphere->nu get_density get_density none~discard_atmosphere->get_density none~dump_apm_pos->get_density none~get_nstar_id_atm->proc~is_finite_number get_nstar_id get_nstar_id none~get_nstar_id_atm->get_nstar_id none~place_and_print_ghost_particles->allocate_sph_memory none~place_and_print_ghost_particles->deallocate_sph_memory none~place_and_print_ghost_particles->proc~spherical_from_cartesian none~place_and_print_ghost_particles->timer none~place_and_print_ghost_particles->proc~is_finite_number none~place_and_print_ghost_particles->allocate_gradient none~place_and_print_ghost_particles->allocate_metric_on_particles none~place_and_print_ghost_particles->allocate_rcb_tree_memory_3d none~place_and_print_ghost_particles->assign_h none~place_and_print_ghost_particles->deallocate_gradient none~place_and_print_ghost_particles->deallocate_metric_on_particles none~place_and_print_ghost_particles->deallocate_rcb_tree_memory_3d none~place_and_print_ghost_particles->density_loop none~place_and_print_ghost_particles->h none~place_and_print_ghost_particles->iorig none~place_and_print_ghost_particles->proc~cartesian_from_spherical none~place_and_print_ghost_particles->proc~find_h_backup bilinear_interpolation bilinear_interpolation none~place_and_print_ghost_particles->bilinear_interpolation center center none~place_and_print_ghost_particles->center compute_pressure compute_pressure none~place_and_print_ghost_particles->compute_pressure none~place_and_print_ghost_particles->get_density none~place_and_print_ghost_particles->get_nstar_id nu_output nu_output none~place_and_print_ghost_particles->nu_output pos_input pos_input none~place_and_print_ghost_particles->pos_input sizes sizes none~place_and_print_ghost_particles->sizes get_pressure_id get_pressure_id none~read_pressure_id->get_pressure_id h_output h_output none~reallocate_output_fields->h_output none~reallocate_output_fields->nu_output none~reallocate_output_fields->pos_input validate_position validate_position none~validate_position_final->validate_position proc~impose_equatorial_plane_symmetry->com interface~find_particles_above_xy_plane find_particles_above_xy_plane proc~impose_equatorial_plane_symmetry->interface~find_particles_above_xy_plane interface~reflect_particles_xy_plane reflect_particles_xy_plane proc~impose_equatorial_plane_symmetry->interface~reflect_particles_xy_plane proc~find_particles_above_xy_plane find_particles_above_xy_plane interface~find_particles_above_xy_plane->proc~find_particles_above_xy_plane proc~reflect_particles_xy_plane reflect_particles_xy_plane interface~reflect_particles_xy_plane->proc~reflect_particles_xy_plane proc~correct_center_of_mass->com proc~correct_center_of_mass->proc~is_finite_number

Contents

Variables

Source Code

convergence_test

Variables

Type Attributes Name Initial
type(bssn), DIMENSION(3) :: bssn_forms

Array storing the bssn objects, containing the BSSN variables on the gravity grid for each idbase object
logical, parameter :: debug = .FALSE.
logical(kind=4) :: dir_out
type(timer) :: execution_timer
logical :: exist
class(idbase), ALLOCATABLE :: idata
integer, parameter :: max_bssn = 3
integer, parameter :: min_bssn = 1
character(len=500) :: name_logfile

String storing the name for the formatted file containing a summary about the \(\mathrm{BSSNOK}\) constraints violations
character(len=500) :: namefile_bssn

String storing the name for the formatted file containing the \(\mathrm{BSSNOK}\) \(\mathrm{ID}\)
character(len=500) :: namefile_bssn_bin

String storing the name for the binary file containing the \(\mathrm{BSSNOK}\) \(\mathrm{ID}\)
character(len=500) :: namefile_parts

String storing the name for the formatted file containing the \(\mathrm{SPH}\) particle \(\mathrm{ID}\)
character(len=500) :: namefile_parts_bin

String storing the name for the binary file containing the \(\mathrm{SPH}\) particle \(\mathrm{ID}\)
character(len=500) :: namefile_sph

String storing the name for ??
double precision :: original_dx
type(particles) :: particles_dist

Array storing the particles objects, containing the particle distributions for each idbase object. Multiple particle objects can contain different particle distributions for the same idbase object.
double precision :: ratio_dx
double precision, DIMENSION(:,:,:,:), ALLOCATABLE :: shared_grid
character(len=:), DIMENSION(:), ALLOCATABLE :: systems

String storing the name of the phyical systems
character(len=:), DIMENSION(:), ALLOCATABLE :: systems_name

String storing the name of the phyical systems
double precision, parameter :: tol_coord = 1.D-10
double precision, parameter :: tol_dx = 1.D-10

Source Code

PROGRAM convergence_test

  !*****************************************************
  !
  !# Make a Cauchy convergence test to check the
  !  validity of the |id| set up by |sphincsid|.
  !
  !  FT 8.12.2020
  !
  !*****************************************************

#ifdef __INTEL_COMPILER

  USE IFPORT,          ONLY: MAKEDIRQQ

#endif

#if flavour == 1

  USE sphincs_id_full,         ONLY: allocate_idbase

#elif flavour == 2

  USE sphincs_id_lorene,       ONLY: allocate_idbase

#elif flavour == 3

  USE sphincs_id_fuka,         ONLY: allocate_idbase

#elif flavour == 4

  USE sphincs_id_interpolate,  ONLY: allocate_idbase

#endif

  USE id_base,                  ONLY: idbase, initialize
  USE sph_particles,            ONLY: particles
  USE bssn_formulation,         ONLY: bssn
  USE standard_tpo_formulation, ONLY: tpo
  USE tensor,                   ONLY: jy, jz
  USE timing,                   ONLY: timer, timers_active
  USE utility,                  ONLY: date, time, zone, values, run_id, &
                                      itr3, hostname, version, &
                                      test_status, show_progress, end_time, &
                                      read_sphincs_id_parameters, one, &
                                      !----------
                                      n_id, common_path, filenames, &
                                      eos_filenames, placer, &
                                      export_bin, export_form, export_form_xy, &
                                      export_form_x, export_constraints_xy, &
                                      export_constraints_x, &
                                      compute_constraints, &
                                      export_constraints, &
                                      export_constraints_details, &
                                      constraints_step, &
                                      compute_parts_constraints, &
                                      numerator_ratio_dx, &
                                      denominator_ratio_dx, &
                                      one_lapse, zero_shift, show_progress, &
                                      run_sph, run_spacetime, sph_path, &
                                      spacetime_path, estimate_length_scale, &
                                      test_int, max_n_parts, ref_lev
  USE cauchy_convergence_test,  ONLY: find_shared_grid, &
                                      perform_cauchy_convergence_test, &
                                      use_constraints_on_mesh, &
                                      use_constraints_with_mapped_hydro
  USE ISO_FORTRAN_ENV,          ONLY: COMPILER_VERSION, COMPILER_OPTIONS

  IMPLICIT NONE

  ! Loop limits for BSSN objects (for debugging; 3 is for production)
  INTEGER, PARAMETER:: min_bssn= 1
  INTEGER, PARAMETER:: max_bssn= 3

  DOUBLE PRECISION, PARAMETER:: tol_dx= 1.D-10
  DOUBLE PRECISION, PARAMETER:: tol_coord= 1.D-10

  ! Grid spacing for the first BSSN object; the other two will have
  ! original_dx/2 and original_dx/4 as grid spacings
  DOUBLE PRECISION:: original_dx
  DOUBLE PRECISION:: ratio_dx

  DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE:: shared_grid

  CHARACTER(LEN=:), DIMENSION(:), ALLOCATABLE:: systems, systems_name
  !! String storing the name of the phyical systems
  CHARACTER(LEN=500):: namefile_parts
  !# String storing the name for the formatted file containing the |sph|
  !  particle |id|
  CHARACTER(LEN=500):: namefile_parts_bin
  !# String storing the name for the binary file containing the |sph|
  !  particle |id|
  CHARACTER(LEN=500):: namefile_sph
  !# String storing the name for ??
  !
  CHARACTER(LEN=500):: namefile_bssn
  !# String storing the name for the formatted file containing the |bssn| |id|
  CHARACTER(LEN=500):: namefile_bssn_bin
  !# String storing the name for the binary file containing the |bssn| |id|
  CHARACTER(LEN=500):: name_logfile
  !# String storing the name for the formatted file containing a summary about
  !  the |bssn| constraints violations

  LOGICAL, PARAMETER:: debug= .FALSE.
  LOGICAL:: exist
  LOGICAL(4):: dir_out

  CLASS(idbase), ALLOCATABLE:: idata

  TYPE(particles):: particles_dist
  !# Array storing the particles objects,
  !  containing the particle distributions for each idbase object.
  !  Multiple particle objects can contain different particle distributions
  !  for the same idbase object.

  TYPE(bssn), DIMENSION(3):: bssn_forms
  !# Array storing the bssn objects,
  !  containing the BSSN variables on the gravity grid for each idbase object

  TYPE(timer):: execution_timer

  !---------------------------!
  !--  End of declarations  --!
  !---------------------------!

  CALL DATE_AND_TIME( date, time, zone, values )
  run_id= date // "-" // time

  !CALL HOSTNM( hostname )
#ifdef host

#ifdef __GFORTRAN__

# define stringize_start(x) "&
# define stringize_end(x) &x"

  hostname= stringize_start(host)
stringize_end(host)

#else

#define stringize(x) tostring(x)
#define tostring(x) #x

  hostname= stringize(host)

#endif

#else

  hostname= "unspecified host."

#endif

#ifdef vers

#ifdef __GFORTRAN__

# define stringize_start(x) "&
# define stringize_end(x) &x"

  version= stringize_start(vers)
stringize_end(vers)

#else

#define stringize(x) tostring(x)
#define tostring(x) #x

  version= stringize(vers)

#endif

#else

  hostname= "unspecified version."

#endif

  PRINT *, "  ________________________________________________________________ "
  PRINT *, "             ____________  ________  __________    __ ___          "
  PRINT *, "            / ___/ _  / /_/ / / __ \/ ___/ ___/   / / __ \         "
  PRINT *, "           (__  ) ___/ __  / / / / / /__(__  )___/ / /_/ /         "
  PRINT *, "          /____/_/  /_/ /_/_/_/ /_/____/____/___/_/_____/          "
  PRINT *
  PRINT *, "  Smoothed Particle Hydrodynamics IN Curved Spacetime              "
  PRINT *, "  Initial Data builder, ", TRIM(version), &
           " - Cauchy convergence test       "
  PRINT *
  PRINT *, "  SPHINCS_ID  Copyright (C) 2020-2023  Francesco Torsello          "
  PRINT *
  PRINT *, "  SPHINCS_ID is free software: you can redistribute it and/or      "
  PRINT *, "  modify it under the terms of the GNU General Public License      "
  PRINT *, "  as published by the Free Software Foundation, either version     "
  PRINT *, "  of the License, or (at your option) any later version.           "
  PRINT *
  PRINT *, "  SPHINCS_ID is distributed in the hope that it will be useful,    "
  PRINT *, "  but WITHOUT ANY WARRANTY; without even the implied warranty of   "
  PRINT *, "  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU "
  PRINT *, "  General Public License for more details.                         "
  PRINT *
  PRINT *, "  You should have received a copy of the GNU General Public License"
  PRINT *, "  along with SPHINCS_ID. If not, see https://www.gnu.org/licenses/."
  PRINT *, "  The copy of the GNU General Public License should be in the file "
  PRINT *, "  'COPYING'.                                                       "
  PRINT *, "  ________________________________________________________________ "
  PRINT *
  PRINT *, "  SPHINCS_ID was compiled with: "
  PRINT *, COMPILER_VERSION()
  PRINT *
  PRINT *, "  using the options: "
  PRINT *, COMPILER_OPTIONS()
  PRINT *
  PRINT *, "  SPHINCS_ID was run on: ", TRIM(hostname)
  PRINT *, "  ________________________________________________________________ "
  PRINT *
  PRINT *, "  Run id: ", run_id
  PRINT *, "  ________________________________________________________________ "
  PRINT *

  timers_active= .TRUE.

  execution_timer= timer( "execution_timer" )
  CALL execution_timer% start_timer()

  CALL read_sphincs_id_parameters()

  ratio_dx= numerator_ratio_dx/denominator_ratio_dx

  ! Check that ratio_dx > 1
  IF( ratio_dx <= one )THEN
    PRINT *, "** ERROR! numerator_ratio_dx has to be larger than ", &
             "denominator_ratio_dx. The current values are ", &
             numerator_ratio_dx, " and ", denominator_ratio_dx, &
             ", respectively."
    PRINT *
    STOP
  ENDIF

#ifdef __INTEL_COMPILER

  INQUIRE( DIRECTORY= TRIM(sph_path), EXIST= exist )
  IF( .NOT.exist )THEN
    dir_out= MAKEDIRQQ( TRIM(sph_path) )
  ELSE
    dir_out= .TRUE.
  ENDIF
  IF( .NOT.dir_out )THEN
    PRINT *, "** ERROR! Failed to create subdirectory ", TRIM(sph_path)
    PRINT *, "Stopping..."
    PRINT *
    STOP
  ENDIF

  INQUIRE( DIRECTORY= TRIM(spacetime_path), EXIST= exist )
  IF( .NOT.exist )THEN
    dir_out= MAKEDIRQQ( TRIM(spacetime_path) )
  ELSE
    dir_out= .TRUE.
  ENDIF
  IF( .NOT.dir_out )THEN
    PRINT *, "** ERROR! Failed to create subdirectory ", TRIM(sph_path)
    PRINT *, "Stopping..."
    PRINT *
    STOP
  ENDIF

#endif

#ifdef __GFORTRAN__

  INQUIRE( FILE= TRIM(sph_path)//"/.", EXIST= exist )
  IF( .NOT.exist )THEN
    CALL EXECUTE_COMMAND_LINE("mkdir "//TRIM(sph_path))
  ENDIF
  INQUIRE( FILE= TRIM(sph_path)//"/.", EXIST= exist )
  IF( .NOT.exist )THEN
    PRINT *, "** ERROR! Failed to create subdirectory ", TRIM(sph_path)
    PRINT *, " * Please create it and re-run the executable."
    PRINT *, " * Stopping..."
    PRINT *
    STOP
  ENDIF

  INQUIRE( FILE= TRIM(spacetime_path)//"/.", EXIST= exist )
  IF( .NOT.exist )THEN
    CALL EXECUTE_COMMAND_LINE("mkdir "//TRIM(spacetime_path))
  ENDIF
  INQUIRE( FILE= TRIM(spacetime_path)//"/.", EXIST= exist )
  IF( .NOT.exist )THEN
    PRINT *, "** ERROR! Failed to create subdirectory ", TRIM(spacetime_path)
    PRINT *, " * Please create it and re-run the executable."
    PRINT *, " * Stopping..."
    PRINT *
    STOP
  ENDIF

#endif

  ALLOCATE( CHARACTER(5):: systems(1) )
  ALLOCATE( CHARACTER(5):: systems_name(1) )

  !
  !-- Construct the idbase objects
  !
  CALL allocate_idbase( idata, TRIM(filenames(1)), systems(1), systems_name(1) )
  PRINT *, "===================================================" &
           // "==============="
  PRINT *, " Constructing idbase object for "//systems(1)
  PRINT *, "===================================================" &
           // "==============="
  PRINT *
  CALL idata% initialize( TRIM(common_path)//TRIM(filenames(1)), eos_filenames )
  ! Set the variables to decide on using the geodesic gauge or not
  ! (lapse=1, shift=0)
  CALL idata% set_one_lapse ( one_lapse )
  CALL idata% set_zero_shift( zero_shift )


  !
  !-- Construct the bssn objects from the idbase object
  !
  construct_bssn_loop: DO itr3 = min_bssn, max_bssn, 1

    PRINT *, "===================================================" &
             // "==============="
    PRINT *, " Setting up BSSN object ", itr3
    PRINT *, "===================================================" &
             // "==============="
    PRINT *

    IF( itr3 == 1 )THEN

      bssn_forms( itr3 )= bssn( idata )
      original_dx= bssn_forms( itr3 )% get_dx(ref_lev)

    ELSE
      IF( itr3 == min_bssn )THEN

        bssn_forms( 1 )= bssn( idata )
        original_dx= bssn_forms( 1 )% get_dx(ref_lev)

      ENDIF
      bssn_forms( itr3 )= bssn( idata, &
                                original_dx/( ratio_dx**( itr3 - 1 ) ), &
                                original_dx/( ratio_dx**( itr3 - 1 ) ), &
                                original_dx/( ratio_dx**( itr3 - 1 ) ) )

      IF( ABS( bssn_forms( itr3 )% get_dx(ref_lev) - &
          original_dx/( ratio_dx**( itr3 - 1 ) ) ) &
          /(original_dx/( ratio_dx**( itr3 - 1 ) )) > tol_dx &
      )THEN
        PRINT *, " ** ERROR! The grid spacing #", itr3, ",", &
                 bssn_forms( itr3 )% get_dx(ref_lev), &
                 " is not equal to dx/", ratio_dx**( itr3 - 1 ), "= ", &
                 original_dx/( ratio_dx**( itr3 - 1 ) )
        STOP
      ENDIF
    ENDIF

    PRINT *, "** The grid spacing is dx=", bssn_forms( itr3 )% get_dx(ref_lev)
    PRINT *, "** The number of grid points for dx is:", &
                            bssn_forms( itr3 )% get_ngrid_x(ref_lev), "**3"
    PRINT *

  ENDDO construct_bssn_loop

  ! Find the grid points shared by the grids
  CALL find_shared_grid( bssn_forms(min_bssn), bssn_forms(min_bssn + 1), &
                         bssn_forms(min_bssn + 2), &
                         numerator_ratio_dx, denominator_ratio_dx, ref_lev, &
                         shared_grid )

  IF( debug )THEN
    PRINT *, "bssn_forms( 1 )% get_ngrid_x=", bssn_forms( 1 )% get_ngrid_x(ref_lev)
    PRINT *, "bssn_forms( 1 )% get_ngrid_y=", bssn_forms( 1 )% get_ngrid_y(ref_lev)
    PRINT *, "bssn_forms( 1 )% get_ngrid_z=", bssn_forms( 1 )% get_ngrid_z(ref_lev)
    PRINT *, "bssn_forms( 2 )% get_ngrid_x=", bssn_forms( 2 )% get_ngrid_x(ref_lev)
    PRINT *, "bssn_forms( 2 )% get_ngrid_y=", bssn_forms( 2 )% get_ngrid_y(ref_lev)
    PRINT *, "bssn_forms( 2 )% get_ngrid_z=", bssn_forms( 2 )% get_ngrid_z(ref_lev)
    PRINT *, "bssn_forms( 3 )% get_ngrid_x=", bssn_forms( 3 )% get_ngrid_x(ref_lev)
    PRINT *, "bssn_forms( 3 )% get_ngrid_y=", bssn_forms( 3 )% get_ngrid_y(ref_lev)
    PRINT *, "bssn_forms( 3 )% get_ngrid_z=", bssn_forms( 3 )% get_ngrid_z(ref_lev)
    PRINT *
    PRINT *, "bssn_forms( 1 )% get_dx ", bssn_forms( 1 )% get_dx(ref_lev)
    PRINT *, "bssn_forms( 2 )% get_dy ", bssn_forms( 2 )% get_dx(ref_lev)
    PRINT *, "bssn_forms( 3 )% get_dz ", bssn_forms( 3 )% get_dx(ref_lev)
    PRINT *
    !STOP
  ENDIF

  !
  !-- Compute the BSSN variables
  !
  compute_export_bssn_loop: DO itr3 = min_bssn, max_bssn, 1
    PRINT *, "===================================================" &
             // "==============="
    PRINT *, " Computing BSSN variables for BSSN formulation", itr3
    PRINT *, "===================================================" &
             // "==============="
    PRINT *
    WRITE( namefile_bssn_bin, "(A10,I1,A4)" ) "BSSN_vars-", itr3, ".bin"
    namefile_bssn_bin= TRIM( spacetime_path ) // TRIM( namefile_bssn_bin )

    bssn_forms( itr3 )% export_bin= export_bin
    bssn_forms( itr3 )% export_form_xy= export_form_xy
    bssn_forms( itr3 )% export_form_x = export_form_x
    CALL bssn_forms( itr3 )% &
                        compute_and_print_tpo_variables( namefile_bssn_bin )
  ENDDO compute_export_bssn_loop

  !
  !-- Print the BSSN initial data to a formatted file
  !
  IF( export_form )THEN
    export_bssn_loop: DO itr3 = min_bssn, max_bssn, 1
      WRITE( namefile_bssn, "(A8,I1,A4)" ) &
                            "bssn-id_", itr3, ".dat"
      namefile_bssn= TRIM( spacetime_path ) // TRIM( namefile_bssn )

      CALL bssn_forms( itr3 )% &
                  print_formatted_id_tpo_variables( namefile= namefile_bssn )
    ENDDO export_bssn_loop
  ENDIF

  !
  !-- Construct the particles object from the idbase object
  !
  IF( compute_parts_constraints )THEN

    PRINT *, "===================================================" &
             // "==============="
    PRINT *, " Placing particles"
    PRINT *, "===================================================" &
             // "==============="
    PRINT *
    particles_dist= particles( idata, placer( 1, 1 ) )

    !
    !-- Compute the SPH variables
    !
    PRINT *, "===================================================" &
             // "====================="
    PRINT *, " Computing SPH variables "
    PRINT *, "===================================================" &
             // "====================="
    PRINT *

    WRITE( namefile_parts_bin, "(A5)" ) systems_name(1)
    namefile_parts_bin= TRIM( sph_path )//TRIM( namefile_parts_bin )

    particles_dist% export_bin    = export_bin
    particles_dist% export_form_xy= export_form_xy
    particles_dist% export_form_x = export_form_x
    CALL particles_dist% compute_and_print_sph_variables( namefile_parts_bin )

    !
    !-- Print the particle initial data to a formatted file
    !
    IF( export_form )THEN
      WRITE( namefile_parts, "(A7,I1,A4)" ) &
                             "sph-id_", itr3, ".dat"
      namefile_parts= TRIM( sph_path ) // TRIM( namefile_parts )
      CALL particles_dist% print_formatted_id_particles( namefile_parts )
    ENDIF

  ENDIF

  !
  !-- Compute the BSSN constraints
  !
  compute_export_bssn_constraints_loop: DO itr3 = min_bssn, max_bssn, 1

    bssn_forms( itr3 )% cons_step= constraints_step
    bssn_forms( itr3 )% export_constraints= export_constraints
    bssn_forms( itr3 )% export_constraints_details= &
                        export_constraints_details
    bssn_forms( itr3 )% export_constraints_xy= export_constraints_xy
    bssn_forms( itr3 )% export_constraints_x = export_constraints_x

    IF( compute_constraints )THEN

      PRINT *, "===================================================" &
               // "==============="
      PRINT *, " Computing BSSN constraints for BSSN formulation", itr3
      PRINT *, "===================================================" &
               // "==============="
      PRINT *

      WRITE( namefile_bssn, "(A17,I1,A4)" ) "bssn-constraints-", itr3, ".dat"
      WRITE( name_logfile, "(A28,I1,A4)" ) &
                          "bssn-constraints-statistics-", itr3

      namefile_bssn= TRIM( spacetime_path ) // TRIM( namefile_bssn )
      name_logfile = TRIM( spacetime_path ) // TRIM( name_logfile )

      CALL bssn_forms( itr3 )% &
                  compute_and_print_tpo_constraints( idata, &
                                                     namefile_bssn, &
                                                     name_logfile, &
                                                     shared_grid )

    ENDIF

    IF( compute_parts_constraints )THEN

      PRINT *, "===================================================" &
               // "==============="
      PRINT *, " Computing BSSN constraints with particle data for BSSN", &
               " formulation", itr3
      PRINT *, "===================================================" &
               // "==============="
      PRINT *

      WRITE( namefile_bssn, "(A23,I1,A4)" ) "bssn-constraints-parts-", &
                                           itr3, ".dat"
      WRITE( namefile_sph, "(A12,I1,A4)" ) "sph-density-", itr3, ".dat"
      WRITE( name_logfile, "(A34,I1,A4)" ) &
                           "bssn-constraints-parts-statistics-", itr3, ".log"

      namefile_bssn= TRIM( spacetime_path ) // TRIM( namefile_bssn )
      namefile_sph = TRIM( sph_path ) // TRIM( namefile_sph )
      name_logfile = TRIM( spacetime_path ) // TRIM( name_logfile )

      CALL bssn_forms( itr3 )% &
                  compute_and_print_tpo_constraints( particles_dist, &
                                                     namefile_bssn, &
                                                     name_logfile, &
                                                     shared_grid )

    ENDIF

  ENDDO compute_export_bssn_constraints_loop

  !
  !-- Perform the convergence test with the appropriate constraints
  !
  IF( compute_constraints )THEN

    PRINT *, "** Performing convergence test with constraints computed ", &
             "without particle data."
    PRINT *

    CALL perform_cauchy_convergence_test &
      ( bssn_forms(1), bssn_forms(2), bssn_forms(3), use_constraints_on_mesh, &
        numerator_ratio_dx, denominator_ratio_dx, ref_lev )

    CALL perform_cauchy_convergence_test &
      ( bssn_forms(2), bssn_forms(3), use_constraints_on_mesh, &
        numerator_ratio_dx, denominator_ratio_dx, ref_lev )

  ENDIF

  IF( compute_parts_constraints )THEN

    PRINT *, "** Performing convergence test with constraints computed ", &
             "with particle data."
    PRINT *

    CALL perform_cauchy_convergence_test &
      ( bssn_forms(1), bssn_forms(2), bssn_forms(3), &
        use_constraints_with_mapped_hydro, &
        numerator_ratio_dx, denominator_ratio_dx, ref_lev )

    CALL perform_cauchy_convergence_test &
      ( bssn_forms(2), bssn_forms(3), &
        use_constraints_with_mapped_hydro, &
        numerator_ratio_dx, denominator_ratio_dx, ref_lev )

  ENDIF

  CALL execution_timer% stop_timer()

  CALL DATE_AND_TIME( date, time, zone, values )
  end_time= date // "-" // time

!STOP

  !
  !-- Print the timers
  !
  PRINT *, "===================================================" &
           // "================================================"
  PRINT *, " Timing and summaries"
  PRINT *, "===================================================" &
           // "================================================"
  PRINT *
  PRINT *
  CALL idata% print_summary()
  !PRINT *, " * BSSN formulation with uniform resolution:", &
  !         bssn_forms( 1 )% get_dx(ref_lev)
  !PRINT *, "    and number of points:", bssn_forms( 1 )% get_ngrid_x(ref_lev), &
  !         "**3"
  !original_dx
  CALL bssn_forms( 1 )% print_summary()
  CALL bssn_forms( 1 )% grid_timer% print_timer( 2 )
  CALL bssn_forms( 1 )% importer_timer% print_timer( 2 )
  CALL bssn_forms( 1 )% bssn_computer_timer% print_timer( 2 )
  PRINT *
  !PRINT *, " * BSSN formulation with uniform resolution:", &
  !         bssn_forms( 2 )% get_dx(ref_lev)
  !PRINT *, "    and number of points:", bssn_forms( 2 )% get_ngrid_x(ref_lev), &
  !         "**3"
  !original_dx/2
  CALL bssn_forms( 2 )% print_summary()
  CALL bssn_forms( 2 )% grid_timer% print_timer( 2 )
  CALL bssn_forms( 2 )% importer_timer% print_timer( 2 )
  CALL bssn_forms( 2 )% bssn_computer_timer% print_timer( 2 )
  PRINT *
  !PRINT *, " * BSSN formulation with uniform resolution:", &
  !         bssn_forms( 3 )% get_dx(ref_lev)
  !PRINT *, "    and number of points:", bssn_forms( 3 )% get_ngrid_x(ref_lev), &
  !         "**3"
  !original_dx/4
  CALL bssn_forms( 3 )% print_summary()
  CALL bssn_forms( 3 )% grid_timer% print_timer( 2 )
  CALL bssn_forms( 3 )% importer_timer% print_timer( 2 )
  CALL bssn_forms( 3 )% bssn_computer_timer% print_timer( 2 )
  PRINT *
  PRINT *, " * Total:"
  CALL execution_timer% print_timer( 2 )
  PRINT *

  PRINT *
  PRINT *, "** Run started on ", run_id, " and ended on ", end_time
  PRINT *

  !
  !-- Destruct the formatted Bin_NS object by hand, since the pointer to it is
  !-- global (because it is bound to C++) and cannot be nullified by the
  !-- destructor of bns. In case of multiple idbase objects, this would lead
  !-- to problems...
  !-- TODO: fix this
  !
  !CALL binary% destruct_binary()


END PROGRAM convergence_test