diff options
Diffstat (limited to 'CarpetExtra/CarpetProlongateTest/src/init.F90')
| -rw-r--r-- | CarpetExtra/CarpetProlongateTest/src/init.F90 | 216 |
1 files changed, 197 insertions, 19 deletions
diff --git a/CarpetExtra/CarpetProlongateTest/src/init.F90 b/CarpetExtra/CarpetProlongateTest/src/init.F90 index 3d12dc3f2..9784dd8e8 100644 --- a/CarpetExtra/CarpetProlongateTest/src/init.F90 +++ b/CarpetExtra/CarpetProlongateTest/src/init.F90 @@ -1,31 +1,45 @@ #include "cctk.h" #include "cctk_Arguments.h" +#include "cctk_Functions.h" #include "cctk_Parameters.h" -subroutine CarpetProlongateTest_Init (CCTK_ARGUMENTS) +module CarpetProlongateTest implicit none - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - CCTK_INT, parameter :: izero = 0 - - ! Add +1 to coordinates so that the domain is not symmetric about - ! the origin (which may accidendally cancel out some errors) - u = 0 + & - 1 * density(x+1,y+1,z+1, power_x,izero ,izero ) + & - 2 * density(x+1,y+1,z+1, izero ,power_y,izero ) + & - 3 * density(x+1,y+1,z+1, izero ,izero ,power_z) + & - 4 * density(x+1,y+1,z+1, power_x,power_y,izero ) + & - 5 * density(x+1,y+1,z+1, power_x,izero ,power_z) + & - 6 * density(x+1,y+1,z+1, izero ,power_y,power_z) + & - 7 * density(x+1,y+1,z+1, power_x,power_y,power_z) - u0 = u +contains - uscaled = u * product(cctk_delta_space(:)) + ! Set the time scaling factor to a sum of powers, ensuring that the + ! scaling factor is not zero for t=0. Also, for power_t=0, the + ! scaling factor is 1 and has no effect. + function density_time_scale (tt) result (tscale) + DECLARE_CCTK_PARAMETERS + CCTK_REAL, intent(in) :: tt + CCTK_REAL :: tscale + integer :: n + + tscale = 1 + do n = 1, power_t + tscale = tscale + tt ** n + end do + end function density_time_scale -contains + elemental function density_sum (xx, yy, zz) result (res) + DECLARE_CCTK_PARAMETERS + CCTK_INT, parameter :: izero = 0 + CCTK_REAL, intent(in) :: xx, yy, zz + CCTK_REAL :: res + + res = 0 + & + 1 * density(xx,yy,zz, power_x,izero ,izero ) + & + 2 * density(xx,yy,zz, izero ,power_y,izero ) + & + 3 * density(xx,yy,zz, izero ,izero ,power_z) + & + 4 * density(xx,yy,zz, power_x,power_y,izero ) + & + 5 * density(xx,yy,zz, power_x,izero ,power_z) + & + 6 * density(xx,yy,zz, izero ,power_y,power_z) + & + 7 * density(xx,yy,zz, power_x,power_y,power_z) + end function density_sum ! Calculate an integrand that is polynomial in the coordinates ! (xx,yy,zz), and which is zero if integrated over [0,2]^3 @@ -50,27 +64,182 @@ contains volume = (xmax - xmin) * (ymax - ymin) * (zmax - zmin) res = xx**px * yy**py * zz**pz - integral / volume - end function density +end module CarpetProlongateTest + + + +subroutine CarpetProlongateTest_Init (CCTK_ARGUMENTS) + use CarpetProlongateTest + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + print '("CarpetProlongateTest_Init it=",i0," t=",g0.6, " rl=",i0)', & + cctk_iteration, cctk_time, nint(log(dble(cctk_levfac(1)))/log(2.0d0)) + + ! Add +1 to coordinates so that the domain is not symmetric about + ! the origin (which may accidendally cancel out some errors) + u = density_time_scale(cctk_time) * density_sum(x+1,y+1,z+1) + + uscaled = u * product(cctk_delta_space(:)) end subroutine CarpetProlongateTest_Init subroutine CarpetProlongateTest_Diff (CCTK_ARGUMENTS) + use CarpetProlongateTest implicit none DECLARE_CCTK_ARGUMENTS + u0 = density_time_scale(cctk_time) * density_sum(x+1,y+1,z+1) du = u - u0 end subroutine CarpetProlongateTest_Diff +subroutine CarpetProlongateTest_InterpInit (CCTK_ARGUMENTS) + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL :: dx, dy, dz + integer :: lbnd(3), lsh(3), gsh(3) + integer :: i,j,k + integer :: ierr + + print '("CarpetProlongateTest_InterpInit it=",i0," t=",g0.6, " rl=",i0)', & + cctk_iteration, cctk_time, nint(log(dble(cctk_levfac(1)))/log(2.0d0)) + + call CCTK_GrouplbndGN (ierr, cctkGH, 3, lbnd, "CarpetProlongateTest::interp_difference") + if (ierr/=0) call CCTK_WARN (0, "internal error") + call CCTK_GrouplshGN (ierr, cctkGH, 3, lsh , "CarpetProlongateTest::interp_difference") + if (ierr/=0) call CCTK_WARN (0, "internal error") + call CCTK_GroupgshGN (ierr, cctkGH, 3, gsh , "CarpetProlongateTest::interp_difference") + if (ierr/=0) call CCTK_WARN (0, "internal error") + + dx = (interp_xmax - interp_xmin) / (gsh(1)-1) + dy = (interp_ymax - interp_ymin) / (gsh(2)-1) + dz = (interp_zmax - interp_zmin) / (gsh(3)-1) + + do k=1,lsh(3) + do j=1,lsh(2) + do i=1,lsh(1) + interp_x(i,j,k) = interp_xmin + (lbnd(1)+i-1) * dx + interp_y(i,j,k) = interp_ymin + (lbnd(2)+j-1) * dy + interp_z(i,j,k) = interp_zmin + (lbnd(3)+k-1) * dz + end do + end do + end do +end subroutine CarpetProlongateTest_InterpInit + +subroutine CarpetProlongateTest_Interp (CCTK_ARGUMENTS) + use cctk + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + integer :: lsh(3) + + integer :: len_interpolator + integer :: len_interpolator_options + character :: fort_interpolator*100 + character :: fort_interpolator_options*1000 + integer :: coord_handle + integer :: interp_handle + integer :: options_table + + integer :: coord_type + CCTK_POINTER :: coords(3) + integer :: ninputs + CCTK_INT :: inputs(1) + + integer :: noutputs + CCTK_INT :: output_types(1) + CCTK_POINTER :: outputs(1) + integer :: npoints + + integer :: ierr + + print '("CarpetProlongateTest_Interp it=",i0," t=",g0.6, " rl=",i0)', & + cctk_iteration, cctk_time, nint(log(dble(cctk_levfac(1)))/log(2.0d0)) + + call CCTK_GrouplshGN (ierr, cctkGH, 3, lsh, "CarpetProlongateTest::interp_difference") + if (ierr/=0) call CCTK_WARN (0, "internal error") + + ! Get coordinate system + call CCTK_CoordSystemHandle (coord_handle, "cart3d") + if (coord_handle<0) then + call CCTK_WARN (0, "Could not obtain coordinate system handle") + end if + + ! Get interpolator + call CCTK_FortranString & + (len_interpolator, interpolator, fort_interpolator) + call CCTK_InterpHandle (interp_handle, fort_interpolator) + if (interp_handle<0) then + call CCTK_WARN (0, "Could not obtain interpolator handle") + end if + + ! Get interpolator options + call CCTK_FortranString & + (len_interpolator_options, interpolator_options, fort_interpolator_options) + call Util_TableCreateFromString (options_table, fort_interpolator_options) + if (options_table<0) then + call CCTK_WARN (0, "Could not create interpolator option table") + end if + + npoints = product(lsh) + + coord_type = CCTK_VARIABLE_REAL + coords = (/ CCTK_PointerTo(interp_x), CCTK_PointerTo(interp_y), CCTK_PointerTo(interp_z) /) + + ninputs = 1 + call CCTK_VarIndex(inputs(1), "CarpetProlongateTest::u") + if (any(inputs<0)) then + call CCTK_WARN (0, "Could not obtain interpolation variable index") + end if + + noutputs = 1 + output_types = (/ CCTK_VARIABLE_REAL /) + outputs = (/ CCTK_PointerTo(interp_u) /) + + call CCTK_InterpGridArrays & + (ierr, cctkGH, 3, & + interp_handle, options_table, coord_handle, & + npoints, coord_type, coords, & + ninputs, inputs, & + noutputs, output_types, outputs) + + if (ierr/=0) then + call CCTK_WARN (0, "Interpolator error") + end if +end subroutine CarpetProlongateTest_Interp + +subroutine CarpetProlongateTest_InterpDiff (CCTK_ARGUMENTS) + use CarpetProlongateTest + implicit none + DECLARE_CCTK_ARGUMENTS + + print '("CarpetProlongateTest_InterpDiff it=",i0," t=",g0.6, " rl=",i0)', & + cctk_iteration, cctk_time, nint(log(dble(cctk_levfac(1)))/log(2.0d0)) + + interp_u0 = & + density_time_scale(cctk_time) * & + density_sum(interp_x+1,interp_y+1,interp_z+1) + interp_du = interp_u - interp_u0 +end subroutine CarpetProlongateTest_InterpDiff + + + subroutine CarpetProlongateTest_NormInit (CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS errornorm = 0 + interp_errornorm = 0 end subroutine CarpetProlongateTest_NormInit subroutine CarpetProlongateTest_NormCalc (CCTK_ARGUMENTS) @@ -78,6 +247,7 @@ subroutine CarpetProlongateTest_NormCalc (CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS errornorm = max (errornorm, maxval(abs(du))) + interp_errornorm = max (interp_errornorm, maxval(abs(interp_du))) end subroutine CarpetProlongateTest_NormCalc subroutine CarpetProlongateTest_NormReduce (CCTK_ARGUMENTS) @@ -92,10 +262,18 @@ subroutine CarpetProlongateTest_NormReduce (CCTK_ARGUMENTS) if (max_handle < 0) then call CCTK_WARN (CCTK_WARN_ABORT, "Could not obtain norm handle") end if + call CCTK_ReduceLocScalar & (ierr, cctkGH, -1, max_handle, errornorm, tmp, CCTK_VARIABLE_REAL) if (ierr/=0) then call CCTK_WARN (CCTK_WARN_ABORT, "Could not evaluate norm") end if errornorm = tmp + + call CCTK_ReduceLocScalar & + (ierr, cctkGH, -1, max_handle, interp_errornorm, tmp, CCTK_VARIABLE_REAL) + if (ierr/=0) then + call CCTK_WARN (CCTK_WARN_ABORT, "Could not evaluate norm") + end if + interp_errornorm = tmp end subroutine CarpetProlongateTest_NormReduce |