diff options
Diffstat (limited to 'src/qlm_import_surface.F90')
| -rw-r--r-- | src/qlm_import_surface.F90 | 190 |
1 files changed, 190 insertions, 0 deletions
diff --git a/src/qlm_import_surface.F90 b/src/qlm_import_surface.F90 new file mode 100644 index 0000000..8eec453 --- /dev/null +++ b/src/qlm_import_surface.F90 @@ -0,0 +1,190 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_import_surface (CCTK_ARGUMENTS, hn) + use cctk + use qlm_boundary + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + integer :: i, j, sn + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Importing surface shape") + end if + + if (surface_index(hn) < 0 .or. surface_index(hn) >= nsurfaces) then + call CCTK_WARN (0, "Illegal spherical surface index specified") + end if + + if (verbose/=0 .or. veryverbose/=0) then + write (msg, '("Importing from spherical surface ",i4)') surface_index(hn) + call CCTK_INFO (msg) + end if + + sn = surface_index(hn) + 1 + + + + if (qlm_calc_error(hn) == 0 .and. cctk_iteration > qlm_iteration(hn)) then + + ! Cycle time levels + qlm_have_valid_data_p_p(hn) = qlm_have_valid_data_p(hn) + qlm_have_valid_data_p (hn) = qlm_have_valid_data (hn) + qlm_have_killing_vector_p_p(hn) = qlm_have_killing_vector_p(hn) + qlm_have_killing_vector_p (hn) = qlm_have_killing_vector (hn) + qlm_iteration(hn) = cctk_iteration + + qlm_time_p_p(hn) = qlm_time_p(hn) + qlm_time_p (hn) = qlm_time (hn) + qlm_radius_p_p(hn) = qlm_radius_p(hn) + qlm_radius_p (hn) = qlm_radius (hn) + + qlm_origin_x_p_p(hn) = qlm_origin_x_p(hn) + qlm_origin_x_p (hn) = qlm_origin_x (hn) + qlm_origin_y_p_p(hn) = qlm_origin_y_p(hn) + qlm_origin_y_p (hn) = qlm_origin_y (hn) + qlm_origin_z_p_p(hn) = qlm_origin_z_p(hn) + qlm_origin_z_p (hn) = qlm_origin_z (hn) + + qlm_shape_p_p(:,:,hn) = qlm_shape_p(:,:,hn) + qlm_shape_p (:,:,hn) = qlm_shape (:,:,hn) + + qlm_l0_p_p(:,:,hn) = qlm_l0_p(:,:,hn) + qlm_l0_p (:,:,hn) = qlm_l0 (:,:,hn) + qlm_l1_p_p(:,:,hn) = qlm_l1_p(:,:,hn) + qlm_l1_p (:,:,hn) = qlm_l1 (:,:,hn) + qlm_l2_p_p(:,:,hn) = qlm_l2_p(:,:,hn) + qlm_l2_p (:,:,hn) = qlm_l2 (:,:,hn) + qlm_l3_p_p(:,:,hn) = qlm_l3_p(:,:,hn) + qlm_l3_p (:,:,hn) = qlm_l3 (:,:,hn) + + qlm_n0_p_p(:,:,hn) = qlm_n0_p(:,:,hn) + qlm_n0_p (:,:,hn) = qlm_n0 (:,:,hn) + qlm_n1_p_p(:,:,hn) = qlm_n1_p(:,:,hn) + qlm_n1_p (:,:,hn) = qlm_n1 (:,:,hn) + qlm_n2_p_p(:,:,hn) = qlm_n2_p(:,:,hn) + qlm_n2_p (:,:,hn) = qlm_n2 (:,:,hn) + qlm_n3_p_p(:,:,hn) = qlm_n3_p(:,:,hn) + qlm_n3_p (:,:,hn) = qlm_n3 (:,:,hn) + + qlm_m0_p_p(:,:,hn) = qlm_m0_p(:,:,hn) + qlm_m0_p (:,:,hn) = qlm_m0 (:,:,hn) + qlm_m1_p_p(:,:,hn) = qlm_m1_p(:,:,hn) + qlm_m1_p (:,:,hn) = qlm_m1 (:,:,hn) + qlm_m2_p_p(:,:,hn) = qlm_m2_p(:,:,hn) + qlm_m2_p (:,:,hn) = qlm_m2 (:,:,hn) + qlm_m3_p_p(:,:,hn) = qlm_m3_p(:,:,hn) + qlm_m3_p (:,:,hn) = qlm_m3 (:,:,hn) + + qlm_xi_t_p_p(:,:,hn) = qlm_xi_t_p(:,:,hn) + qlm_xi_t_p (:,:,hn) = qlm_xi_t (:,:,hn) + qlm_xi_p_p_p(:,:,hn) = qlm_xi_p_p(:,:,hn) + qlm_xi_p_p (:,:,hn) = qlm_xi_p (:,:,hn) + + end if + + + + ! Check for valid horizon data + if (sf_valid(sn) <= 0) then + if (verbose/=0) then + call CCTK_INFO ("No valid horizon data found") + end if + qlm_calc_error(hn) = 1 + qlm_have_valid_data(hn) = 0 + qlm_have_killing_vector(hn) = 0 + return + end if + + + + ! Import surface description + qlm_calc_error(hn) = 0 + qlm_have_valid_data(hn) = 0 + qlm_have_killing_vector(hn) = 1 + + if (qlm_have_valid_data_p(hn) == 0) then + qlm_timederiv_order(hn) = 0 + else if (qlm_have_valid_data_p_p(hn) == 0) then + qlm_timederiv_order(hn) = 1 + else + qlm_timederiv_order(hn) = 2 + end if + + qlm_time(hn) = cctk_time + +#warning "TODO: Ensure that the surface parameters don't change" + qlm_nghoststheta(hn) = sf_nghoststheta(sn) + qlm_nghostsphi (hn) = sf_nghostsphi(sn) + qlm_ntheta (hn) = sf_ntheta(sn) + qlm_nphi (hn) = sf_nphi(sn) + + qlm_origin_x (hn) = sf_origin_x(sn) + qlm_origin_y (hn) = sf_origin_y(sn) + qlm_origin_z (hn) = sf_origin_z(sn) + qlm_origin_theta(hn) = sf_origin_theta(sn) + qlm_origin_phi (hn) = sf_origin_phi(sn) + qlm_delta_theta (hn) = sf_delta_theta(sn) + qlm_delta_phi (hn) = sf_delta_phi(sn) + + if (veryverbose /= 0) then + write (msg, '("calc error : ",i6)') qlm_calc_error(hn) + call CCTK_INFO (msg) + write (msg, '("time deriv order: ",i6)') qlm_timederiv_order(hn) + call CCTK_INFO (msg) + write (msg, '("time : ",g16.6)') qlm_time(hn) + call CCTK_INFO (msg) + write (msg, '("nghosts : ",2i6)') qlm_nghoststheta(hn), qlm_nghostsphi(hn) + call CCTK_INFO (msg) + write (msg, '("n : ",2i6)') qlm_ntheta(hn), qlm_nphi(hn) + call CCTK_INFO (msg) + write (msg, '("origin : ",3g16.6)') & + qlm_origin_x(hn), qlm_origin_y(hn), qlm_origin_z(hn) + call CCTK_INFO (msg) + write (msg, '("origin : ",2g16.6)') & + qlm_origin_theta(hn), qlm_origin_phi(hn) + call CCTK_INFO (msg) + write (msg, '("delta : ",2g16.6)') & + qlm_delta_theta(hn), qlm_delta_phi(hn) + call CCTK_INFO (msg) + end if + + if (qlm_ntheta(hn) > maxntheta .or. qlm_nphi(hn) > maxnphi) then + call CCTK_WARN (0, "Surface is too large") + end if + + if (qlm_nghoststheta(hn)<1 .or. qlm_nghostsphi(hn)<1) then + call CCTK_WARN (0, "Not enough ghost zones for the horizon surface -- need at least 1") + end if + + + ! Import the surface + ! Calculate the coordinates + do j = 1, qlm_nphi(hn) + do i = 1, qlm_ntheta(hn) + qlm_shape(i,j,hn) = sf_radius(i,j,sn) + end do + end do + + + + if (mod(int(qlm_ntheta(hn) - 2*qlm_nghoststheta(hn)),2) /= 1) then + ! We need a grid point on the equator + call CCTK_WARN (0, "The number of interior grid points in theta direction must be odd") + end if + + if (mod(int(qlm_nphi(hn) - 2*qlm_nghostsphi(hn)),4) /= 0) then + ! We need grid points on the four major meridians + call CCTK_WARN (0, "The number of interior grid points in phi direction must be a multiple of four") + end if + +end subroutine qlm_import_surface |