diff options
Diffstat (limited to 'src/EHFinder_Generator_Sources.F90')
| -rw-r--r-- | src/EHFinder_Generator_Sources.F90 | 299 |
1 files changed, 153 insertions, 146 deletions
diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index a89a548..bb8e350 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -15,12 +15,13 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i + CCTK_INT :: i, l CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: gen_interp CCTK_INT :: gen_interp_len character(len=7) :: gen_order + character(len=15) :: vname CCTK_INT, dimension(1) :: lsh CCTK_POINTER, dimension(3) :: interp_coords @@ -70,16 +71,11 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) endif ! Find out how many interpolation points are located on this processor. - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) end if - ! Set the pointers to the points to be interpolated to. - interp_coords(1) = CCTK_PointerTo(xg) - interp_coords(2) = CCTK_PointerTo(yg) - interp_coords(3) = CCTK_PointerTo(zg) - ! Set the pointers to the output arrays. out_arrays(1) = CCTK_PointerTo(alpg) out_arrays(2) = CCTK_PointerTo(betaxg) @@ -96,147 +92,158 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) out_arrays(13) = CCTK_PointerTo(dfzg) out_arrays(14) = CCTK_PointerTo(psig) - ! Set the indices to the input grid functions. - call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) - call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) - call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) - call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) - call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) - call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) - call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) - call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) - call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) - call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) - call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" ) - call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" ) - - ! Check the metric type. At present physical and static_conformal are - ! supported. - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - - ! Set the operand indices table entry, corresponding - ! to interpolation of admbase::alp (1), admbase::shift (3), - ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for - ! a total of 13 output arrays. - call Util_TableSetIntArray ( status, table_handle, 13, & - op_indices(1:13), "operand_indices" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) - endif - - ! Set the corresponding table entry for the operation codes. - call Util_TableSetIntArray ( status, table_handle, 13, & - op_codes(1:13), "operation_codes" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) - endif - - ! Call the interpolator. - call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - lsh(1), CCTK_VARIABLE_REAL, & - interp_coords, 11, in_arrays(1:11), & - 13, out_types(1:13), out_arrays(1:13) ) - - if ( status .lt. 0 ) then - call CCTK_INFO ( 'Interpolation failed.' ) - end if - - ! For each point on this processor calculate the right hand side of the - ! characteristic evolution equation. - do i = 1, lsh(1) - - ! calculate the square of the lapse. - alp2 = alpg(i)**2 - - ! Calculate the inverse of the 3-metric. - guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 - guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) - guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) - - idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz ) - - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz - - guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg - guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg - guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + ! Loop over the level sets + do l = 1, eh_number_level_sets + + ! Set the pointers to the points to be interpolated to. + interp_coords(1) = CCTK_PointerTo(xg(:,l)) + interp_coords(2) = CCTK_PointerTo(yg(:,l)) + interp_coords(3) = CCTK_PointerTo(zg(:,l)) + + ! Set the indices to the input grid functions. + call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) + call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) + call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) + call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) + call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) + call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) + call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) + call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) + call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) + call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) + write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1,']' + call CCTK_VarIndex ( in_arrays(11), vname ) + call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" ) + + ! Check the metric type. At present physical and static_conformal are + ! supported. + if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then + + ! Set the operand indices table entry, corresponding + ! to interpolation of admbase::alp (1), admbase::shift (3), + ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for + ! a total of 13 output arrays. + call Util_TableSetIntArray ( status, table_handle, 13, & + op_indices(1:13), "operand_indices" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + endif + + ! Set the corresponding table entry for the operation codes. + call Util_TableSetIntArray ( status, table_handle, 13, & + op_codes(1:13), "operation_codes" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + endif + + ! Call the interpolator. + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 11, in_arrays(1:11), & + 13, out_types(1:13), out_arrays(1:13) ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + ! For each point on this processor calculate the right hand side of the + ! characteristic evolution equation. + do i = 1, lsh(1) + + ! calculate the square of the lapse. + alp2 = alpg(i)**2 + + ! Calculate the inverse of the 3-metric. + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + + idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz ) + + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg - ! Raise the index of the partial derivatives of f. - dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) - dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) - dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) - - ! Calculate the overall multiplication factor. - factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & - dfuy * dfyg(i) + & - dfuz * dfzg(i) ) ) - - ! Finally obtain dx^i/dt. - dxg(i) = - betaxg(i) + ssign * factor * dfux - dyg(i) = - betayg(i) + ssign * factor * dfuy - dzg(i) = - betazg(i) + ssign * factor * dfuz - end do - else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - call Util_TableSetIntArray ( status, table_handle, 14, & - op_indices, "operand_indices" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) - endif - - call Util_TableSetIntArray ( status, table_handle, 14, & - op_codes, "operation_codes" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) - endif - - call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - lsh(1), CCTK_VARIABLE_REAL, & - interp_coords, 12, in_arrays, & - 14, out_types, out_arrays ) - - if ( status .lt. 0 ) then - call CCTK_INFO ( 'Interpolation failed.' ) + ! Raise the index of the partial derivatives of f. + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + + ! Calculate the overall multiplication factor. + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + + ! Finally obtain dx^i/dt. + dxg(i,l) = - betaxg(i) + ssign * factor * dfux + dyg(i,l) = - betayg(i) + ssign * factor * dfuy + dzg(i,l) = - betazg(i) + ssign * factor * dfuz + end do + else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + call Util_TableSetIntArray ( status, table_handle, 14, & + op_indices, "operand_indices" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + endif + + call Util_TableSetIntArray ( status, table_handle, 14, & + op_codes, "operation_codes" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + endif + + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 12, in_arrays, & + 14, out_types, out_arrays ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + do i = 1, lsh(1) + alp2 = alpg(i)**2 + +! The inverse of psi^4 + psi4 = one / psig(i)**4 + + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + +! The determinant divided by psi^4. + idetg = psi4 / ( gxxg(i) * guxx + & + gxyg(i) * guxy + & + gxzg(i) * guxz ) + +! The inverse metric. Since the determinant is already divided +! by psi^4, this gives the inverse of the physical metric. + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + dxg(i,l) = - betaxg(i) + ssign * factor * dfux + dyg(i,l) = - betayg(i) + ssign * factor * dfuy + dzg(i,l) = - betazg(i) + ssign * factor * dfuz + end do end if + end do - do i = 1, lsh(1) - alp2 = alpg(i)**2 - -! The inverse of psi^4 - psi4 = one / psig(i)**4 - - guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 - guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) - guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) - -! The determinant divided by psi^4. - idetg = psi4 / ( gxxg(i) * guxx + & - gxyg(i) * guxy + & - gxzg(i) * guxz ) - -! The inverse metric. Since the determinant is already divided -! by psi^4, this gives the inverse of the physical metric. - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz - - guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg - guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg - guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg - - dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) - dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) - dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) - factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & - dfuy * dfyg(i) + & - dfuz * dfzg(i) ) ) - dxg(i) = - betaxg(i) + ssign * factor * dfux - dyg(i) = - betayg(i) + ssign * factor * dfuy - dzg(i) = - betazg(i) + ssign * factor * dfuz - end do - end if return end subroutine EHFinder_Generator_Sources |