From ac099068da91d9528088af849caa7aa2e6e1758a Mon Sep 17 00:00:00 2001 From: jthorn Date: Sun, 17 Nov 2002 14:55:50 +0000 Subject: Modified Files: make.code.defn Added Files: Bona_Masso_data.F77 blended_boundary.F77 boundary.F77 decode_pars.F77 gauge.F77 initialize.F77 metric.F77 xyz_blended_boundary.F77 Removed Files: Bona_Masso_data.F blended_boundary.F boundary.F decode_pars.F gauge.F initialize.F metric.F xyz_blended_boundary.F Rename files which really are Fortran 77 from foo.F (= Cactus Fortran 90 fixed form) to foo.F77 (= Fortran 77) This means they're now compiled with a Fortran 77 compiler. This should make no difference to the semantics (they were already Fortran 77 code), but makes it easier to compile this thorn on platforms which have a Fortran 77 compiler but no Fortran 90 compiler. Also small bugfix: blended_boundary.F --> .F77 xyz_blended_boundary.F --> .F77 * change declaration of ierr from CCTK_REAL --> integer git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@135 e296648e-0e4f-0410-bd07-d597d9acff87 --- src/Bona_Masso_data.F | 214 ------------------------------- src/Bona_Masso_data.F77 | 214 +++++++++++++++++++++++++++++++ src/blended_boundary.F | 196 ----------------------------- src/blended_boundary.F77 | 197 +++++++++++++++++++++++++++++ src/boundary.F | 125 ------------------ src/boundary.F77 | 125 ++++++++++++++++++ src/decode_pars.F | 201 ----------------------------- src/decode_pars.F77 | 201 +++++++++++++++++++++++++++++ src/gauge.F | 152 ---------------------- src/gauge.F77 | 152 ++++++++++++++++++++++ src/initialize.F | 78 ------------ src/initialize.F77 | 78 ++++++++++++ src/make.code.defn | 16 +-- src/metric.F | 293 ------------------------------------------- src/metric.F77 | 293 +++++++++++++++++++++++++++++++++++++++++++ src/xyz_blended_boundary.F | 210 ------------------------------- src/xyz_blended_boundary.F77 | 211 +++++++++++++++++++++++++++++++ 17 files changed, 1479 insertions(+), 1477 deletions(-) delete mode 100644 src/Bona_Masso_data.F create mode 100644 src/Bona_Masso_data.F77 delete mode 100644 src/blended_boundary.F create mode 100644 src/blended_boundary.F77 delete mode 100644 src/boundary.F create mode 100644 src/boundary.F77 delete mode 100644 src/decode_pars.F create mode 100644 src/decode_pars.F77 delete mode 100644 src/gauge.F create mode 100644 src/gauge.F77 delete mode 100644 src/initialize.F create mode 100644 src/initialize.F77 delete mode 100644 src/metric.F create mode 100644 src/metric.F77 delete mode 100644 src/xyz_blended_boundary.F create mode 100644 src/xyz_blended_boundary.F77 diff --git a/src/Bona_Masso_data.F b/src/Bona_Masso_data.F deleted file mode 100644 index 9805e0a..0000000 --- a/src/Bona_Masso_data.F +++ /dev/null @@ -1,214 +0,0 @@ -c This routine calculates Bona-Masso initial data, making use of the -c subroutine Exact__metric() to calculate the spacetime metric and its -c inverse. Note that this use of the Bona-Masso variables is independent -c of how (or even if) we are evolving the Einstein equations -- here -c the Bona-Masso variables are "just" used as intermediate variables. - -c $Header$ - -#include "cctk.h" - - subroutine Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x, y, z, t, - $ gxx, gyy, gzz, gxy, gyz, gxz, - $ hxx, hyy, hzz, hxy, hyz, hxz, - $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, - $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, - $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, - $ alp, ax, ay, az, betax, betay, betaz, - $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz) - - implicit none - CCTK_INT decoded_exact_model - CCTK_REAL x, y, z, t, - $ gxx, gyy, gzz, gxy, gyz, gxz, - $ hxx, hyy, hzz, hxy, hyz, hxz, - $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, - $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, - $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, - $ alp, ax, ay, az, betax, betay, betaz, - $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz - -C gxx is g_xx etc. -C hxx is K_xx etc. -C dxgyy is (/2) dg_yy / dx (sic!) -C alp is N, betax is N^x etc. -C bxy is (/2) dN^y / dx (sic and sic!) -C ax is dN / dx / N (sic!) - - CCTK_REAL eps, - $ gdtt, gdtx, gdty, gdtz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, - $ gdtt_p, gdtx_p, gdty_p, gdtz_p, - $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, - $ gutt_p, gutx_p, guty_p, gutz_p, - $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p, - $ gdtt_m, gdtx_m, gdty_m, gdtz_m, - $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, - $ gutt_m, gutx_m, guty_m, gutz_m, - $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m - - parameter (eps=1.d-6) - -C Get the spacetime metric and its inverse at the base point. - - call Exact__metric( - $ decoded_exact_model, - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gxx, gyy, gzz, gxy, gyz, gxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -C Calculate lapse and shift from the upper metric. - - alp = 1.d0 / sqrt(-gutt) - - betax = - gutx / gutt - betay = - guty / gutt - betaz = - gutz / gutt - -C In order to calculate the derivatives of the lapse and shift from -C the contravariant metric, use g^tt = -1/N^2 and g^i0 = N^i / N^2 -C Note that ax is dN/dx / N and that bxy is 1/2 dN^y / dx. - -C Calculate x-derivatives. - - call Exact__metric( - $ decoded_exact_model, - $ x+eps, y, z, t, - $ gdtt_p, gdtx_p, gdty_p, gdtz_p, - $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, - $ gutt_p, gutx_p, guty_p, gutz_p, - $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) - call Exact__metric( - $ decoded_exact_model, - $ x-eps, y, z, t, - $ gdtt_m, gdtx_m, gdty_m, gdtz_m, - $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, - $ gutt_m, gutx_m, guty_m, gutz_m, - $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) - - dxgxx = (gdxx_p - gdxx_m) / 4.d0 / eps - dxgyy = (gdyy_p - gdyy_m) / 4.d0 / eps - dxgzz = (gdzz_p - gdzz_m) / 4.d0 / eps - dxgxy = (gdxy_p - gdxy_m) / 4.d0 / eps - dxgyz = (gdyz_p - gdyz_m) / 4.d0 / eps - dxgxz = (gdxz_p - gdxz_m) / 4.d0 / eps - - ax = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt - - bxx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps - bxy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps - bxz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps - -C Calculate y-derivatives. - - call Exact__metric( - $ decoded_exact_model, - $ x, y+eps, z, t, - $ gdtt_p, gdtx_p, gdty_p, gdtz_p, - $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, - $ gutt_p, gutx_p, guty_p, gutz_p, - $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) - call Exact__metric( - $ decoded_exact_model, - $ x, y-eps, z, t, - $ gdtt_m, gdtx_m, gdty_m, gdtz_m, - $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, - $ gutt_m, gutx_m, guty_m, gutz_m, - $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) - - dygxx = (gdxx_p - gdxx_m) / 4.d0 / eps - dygyy = (gdyy_p - gdyy_m) / 4.d0 / eps - dygzz = (gdzz_p - gdzz_m) / 4.d0 / eps - dygxy = (gdxy_p - gdxy_m) / 4.d0 / eps - dygyz = (gdyz_p - gdyz_m) / 4.d0 / eps - dygxz = (gdxz_p - gdxz_m) / 4.d0 / eps - - ay = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt - - byx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps - byy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps - byz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps - -C Calculate z-derivatives. - - call Exact__metric( - $ decoded_exact_model, - $ x, y, z+eps, t, - $ gdtt_p, gdtx_p, gdty_p, gdtz_p, - $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, - $ gutt_p, gutx_p, guty_p, gutz_p, - $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) - call Exact__metric( - $ decoded_exact_model, - $ x, y, z-eps, t, - $ gdtt_m, gdtx_m, gdty_m, gdtz_m, - $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, - $ gutt_m, gutx_m, guty_m, gutz_m, - $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) - - dzgxx = (gdxx_p - gdxx_m) / 4.d0 / eps - dzgyy = (gdyy_p - gdyy_m) / 4.d0 / eps - dzgzz = (gdzz_p - gdzz_m) / 4.d0 / eps - dzgxy = (gdxy_p - gdxy_m) / 4.d0 / eps - dzgyz = (gdyz_p - gdyz_m) / 4.d0 / eps - dzgxz = (gdxz_p - gdxz_m) / 4.d0 / eps - - az = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt - - bzx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps - bzy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps - bzz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps - -C Calculate t-derivatives, and extrinsic curvature. - - call Exact__metric( - $ decoded_exact_model, - $ x, y, z, t+eps, - $ gdtt_p, gdtx_p, gdty_p, gdtz_p, - $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, - $ gutt_p, gutx_p, guty_p, gutz_p, - $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) - call Exact__metric( - $ decoded_exact_model, - $ x, y, z, t-eps, - $ gdtt_m, gdtx_m, gdty_m, gdtz_m, - $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, - $ gutt_m, gutx_m, guty_m, gutz_m, - $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) - - hxx = - (gdxx_p - gdxx_m) / 4.d0 / eps / alp - $ + (dxgxx * betax + dygxx * betay + dzgxx * betaz - $ + 2.d0 * (bxx * gxx + bxy * gxy + bxz * gxz)) / alp - - hyy = - (gdyy_p - gdyy_m) / 4.d0 / eps / alp - $ + (dxgyy * betax + dygyy * betay + dzgyy * betaz - $ + 2.d0 * (byx * gxy + byy * gyy + byz * gyz)) / alp - - hzz = - (gdzz_p - gdzz_m) / 4.d0 / eps / alp - $ + (dxgzz * betax + dygzz * betay + dzgzz * betaz - $ + 2.d0 * (bzx * gxz + bzy * gyz + bzz * gzz)) / alp - - hxy = - (gdxy_p - gdxy_m) / 4.d0 / eps / alp - $ + (dxgxy * betax + dygxy * betay + dzgxy * betaz - $ + bxx * gxy + bxy * gyy + bxz * gyz - $ + byx * gxx + byy * gxy + byz * gxz) / alp - - hyz = - (gdyz_p - gdyz_m) / 4.d0 / eps / alp - $ + (dxgyz * betax + dygyz * betay + dzgyz * betaz - $ + byx * gxz + byy * gyz + byz * gzz - $ + bzx * gxy + bzy * gyy + bzz * gyz) / alp - - hxz = - (gdxz_p - gdxz_m) / 4.d0 / eps / alp - $ + (dxgxz * betax + dygxz * betay + dzgxz * betaz - $ + bxx * gxz + bxy * gyz + bxz * gzz - $ + bzx * gxx + bzy * gxy + bzz * gxz) / alp - - return - end - - diff --git a/src/Bona_Masso_data.F77 b/src/Bona_Masso_data.F77 new file mode 100644 index 0000000..9805e0a --- /dev/null +++ b/src/Bona_Masso_data.F77 @@ -0,0 +1,214 @@ +c This routine calculates Bona-Masso initial data, making use of the +c subroutine Exact__metric() to calculate the spacetime metric and its +c inverse. Note that this use of the Bona-Masso variables is independent +c of how (or even if) we are evolving the Einstein equations -- here +c the Bona-Masso variables are "just" used as intermediate variables. + +c $Header$ + +#include "cctk.h" + + subroutine Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x, y, z, t, + $ gxx, gyy, gzz, gxy, gyz, gxz, + $ hxx, hyy, hzz, hxy, hyz, hxz, + $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, + $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, + $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, + $ alp, ax, ay, az, betax, betay, betaz, + $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz) + + implicit none + CCTK_INT decoded_exact_model + CCTK_REAL x, y, z, t, + $ gxx, gyy, gzz, gxy, gyz, gxz, + $ hxx, hyy, hzz, hxy, hyz, hxz, + $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, + $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, + $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, + $ alp, ax, ay, az, betax, betay, betaz, + $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz + +C gxx is g_xx etc. +C hxx is K_xx etc. +C dxgyy is (/2) dg_yy / dx (sic!) +C alp is N, betax is N^x etc. +C bxy is (/2) dN^y / dx (sic and sic!) +C ax is dN / dx / N (sic!) + + CCTK_REAL eps, + $ gdtt, gdtx, gdty, gdtz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m + + parameter (eps=1.d-6) + +C Get the spacetime metric and its inverse at the base point. + + call Exact__metric( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gxx, gyy, gzz, gxy, gyz, gxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +C Calculate lapse and shift from the upper metric. + + alp = 1.d0 / sqrt(-gutt) + + betax = - gutx / gutt + betay = - guty / gutt + betaz = - gutz / gutt + +C In order to calculate the derivatives of the lapse and shift from +C the contravariant metric, use g^tt = -1/N^2 and g^i0 = N^i / N^2 +C Note that ax is dN/dx / N and that bxy is 1/2 dN^y / dx. + +C Calculate x-derivatives. + + call Exact__metric( + $ decoded_exact_model, + $ x+eps, y, z, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call Exact__metric( + $ decoded_exact_model, + $ x-eps, y, z, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + + dxgxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dxgyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dxgzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dxgxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dxgyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dxgxz = (gdxz_p - gdxz_m) / 4.d0 / eps + + ax = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + + bxx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + bxy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + bxz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate y-derivatives. + + call Exact__metric( + $ decoded_exact_model, + $ x, y+eps, z, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call Exact__metric( + $ decoded_exact_model, + $ x, y-eps, z, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + + dygxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dygyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dygzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dygxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dygyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dygxz = (gdxz_p - gdxz_m) / 4.d0 / eps + + ay = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + + byx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + byy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + byz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate z-derivatives. + + call Exact__metric( + $ decoded_exact_model, + $ x, y, z+eps, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call Exact__metric( + $ decoded_exact_model, + $ x, y, z-eps, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + + dzgxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dzgyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dzgzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dzgxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dzgyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dzgxz = (gdxz_p - gdxz_m) / 4.d0 / eps + + az = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + + bzx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + bzy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + bzz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate t-derivatives, and extrinsic curvature. + + call Exact__metric( + $ decoded_exact_model, + $ x, y, z, t+eps, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call Exact__metric( + $ decoded_exact_model, + $ x, y, z, t-eps, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + + hxx = - (gdxx_p - gdxx_m) / 4.d0 / eps / alp + $ + (dxgxx * betax + dygxx * betay + dzgxx * betaz + $ + 2.d0 * (bxx * gxx + bxy * gxy + bxz * gxz)) / alp + + hyy = - (gdyy_p - gdyy_m) / 4.d0 / eps / alp + $ + (dxgyy * betax + dygyy * betay + dzgyy * betaz + $ + 2.d0 * (byx * gxy + byy * gyy + byz * gyz)) / alp + + hzz = - (gdzz_p - gdzz_m) / 4.d0 / eps / alp + $ + (dxgzz * betax + dygzz * betay + dzgzz * betaz + $ + 2.d0 * (bzx * gxz + bzy * gyz + bzz * gzz)) / alp + + hxy = - (gdxy_p - gdxy_m) / 4.d0 / eps / alp + $ + (dxgxy * betax + dygxy * betay + dzgxy * betaz + $ + bxx * gxy + bxy * gyy + bxz * gyz + $ + byx * gxx + byy * gxy + byz * gxz) / alp + + hyz = - (gdyz_p - gdyz_m) / 4.d0 / eps / alp + $ + (dxgyz * betax + dygyz * betay + dzgyz * betaz + $ + byx * gxz + byy * gyz + byz * gzz + $ + bzx * gxy + bzy * gyy + bzz * gyz) / alp + + hxz = - (gdxz_p - gdxz_m) / 4.d0 / eps / alp + $ + (dxgxz * betax + dygxz * betay + dzgxz * betaz + $ + bxx * gxz + bxy * gyz + bxz * gzz + $ + bzx * gxx + bzy * gxy + bzz * gxz) / alp + + return + end + + diff --git a/src/blended_boundary.F b/src/blended_boundary.F deleted file mode 100644 index f74beb3..0000000 --- a/src/blended_boundary.F +++ /dev/null @@ -1,196 +0,0 @@ -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - - subroutine Exact__blended_boundary(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - logical doKij, doGij, doLapse, doShift - - integer i,j,k - integer nx,ny,nz - - CCTK_REAL router, rinner, frac, onemfrac - CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze - CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze - CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze - CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze - CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze - - CCTK_REAL alpe, axe, aye, aze - CCTK_REAL betaxe,betaye,betaze - CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze - - CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz - CCTK_REAL vxe,vye,vze,sav - - CCTK_REAL :: dx,dy,dz,time - CCTK_REAL :: ierr,xmx,xmn,ymx,ymn,zmx,zmn,rmx - -C Grid parameters. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - - dx = cctk_delta_space(1) - dy = cctk_delta_space(2) - dz = cctk_delta_space(3) - - time = cctk_time - -C Other parameters. - - doKij = (exblend_Ks.eq.1) - doGij = (exblend_gs.eq.1) - - doLapse = ((exblend_gauge.eq.1).and. - $ (CCTK_Equals(lapse_evolution_method,"exact").ne.0)) - doShift = ((exblend_gauge.eq.1).and. - $ (CCTK_Equals(shift_evolution_method,"exact").ne.0)) - - call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d") - call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d") - call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d") - - rmx = min(xmx,ymx,zmx) - - if (exblend_rout.lt.0) then - router = rmx - 2.0d0*dx - else - router = exblend_rout - endif - - if (exblend_width.lt.0) then - rinner = router + exblend_width*dx - else - rinner = router - exblend_width - endif - - do k=1,nz - do j=1,ny - do i=1,nx - -c We only do anything if r >= rinner so only evaluate exact data -c there. - - if (r(i,j,k) .ge. rinner) then - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), time, - $ gxxe, gyye, gzze, gxye, gyze, gxze, - $ kxxe, kyye, kzze, kxye, kyze, kxze, - $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, - $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, - $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, - $ alpe, axe, aye, aze, betaxe, betaye, betaze, - $ bxxe, bxye, bxze, byxe, - $ byye, byze, bzxe, bzye, bzze) - -c This sucks, but we want the exact vs so we can blend them also. - - det = -(gxze**2*gyye) - & + 2.d0*gxye*gxze*gyze - & - gxxe*gyze**2 - & - gxye**2*gzze - & + gxxe*gyye*gzze - - uxx=(-gyze**2 + gyye*gzze)/det - uxy=(gxze*gyze - gxye*gzze)/det - uyy=(-gxze**2 + gxxe*gzze)/det - uxz=(-gxze*gyye + gxye*gyze)/det - uyz=(gxye*gxze - gxxe*gyze)/det - uzz=(-gxye**2 + gxxe*gyye)/det - -c Outside of router we want to place exact data on our grid - - if (r(i,j,k) .gt. router) then - -c This is one of those things I will invariably screw up if I type -c it in so let the computer do it - -#define exassign(q) q(i,j,k) = q e -#define exassign_grp(p) \ - exassign(p xx) &&\ - exassign(p xy) &&\ - exassign(p xz) &&\ - exassign(p yy) &&\ - exassign(p yz) &&\ - exassign(p zz) - -c Note this plays on the nasty trick that fortran doesnt give a -c damn about spaces so gxx e is the same as gxxe for the parser... -c Grody but effective! - - if (doGij) then - exassign_grp(g) - endif - if (doKij) then - exassign_grp(k) - endif - - if (doLapse) then - exassign(alp) - endif - - if (doShift.and.(shift_state.ne.0)) then - exassign(betax) - exassign(betay) - exassign(betaz) - endif - -c OK so we dont want to blend so use a goto to jump. - else - -c Evaluate the linear weighting fraction. Obvious... - - frac = (r(i,j,k) - rinner) / (router - rinner) - onemfrac = 1.0D0 - frac - -c Once again some c-preprocessor tricks based on the whole fortran -c space thing... - -#define INTPOINT(f,v) f(i,j,k) = frac * v + onemfrac * f(i,j,k) -#define intone(f) INTPOINT(f, f e) -#define int_grp(p) \ - intone(p xx) &&\ - intone(p xy) &&\ - intone(p xz) &&\ - intone(p yy) &&\ - intone(p yz) &&\ - intone(p zz) - - if (doGij) then - int_grp(g) - endif - if (doKij) then - int_grp(k) - endif - - if (doLapse) then - intone(alp) - endif - - if (doShift.and.(shift_state.ne.0)) then - intone(betax) - intone(betay) - intone(betaz) - endif - - endif ! r > router else - endif ! r > rinner - - enddo - enddo - enddo - - return - end diff --git a/src/blended_boundary.F77 b/src/blended_boundary.F77 new file mode 100644 index 0000000..804f679 --- /dev/null +++ b/src/blended_boundary.F77 @@ -0,0 +1,197 @@ +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + subroutine Exact__blended_boundary(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + logical doKij, doGij, doLapse, doShift + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL router, rinner, frac, onemfrac + CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze + CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze + CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze + CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze + CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze + + CCTK_REAL alpe, axe, aye, aze + CCTK_REAL betaxe,betaye,betaze + CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze + + CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL vxe,vye,vze,sav + + CCTK_REAL dx,dy,dz,time + CCTK_REAL xmx,xmn,ymx,ymn,zmx,zmn,rmx + integer ierr + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + time = cctk_time + +C Other parameters. + + doKij = (exblend_Ks.eq.1) + doGij = (exblend_gs.eq.1) + + doLapse = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(lapse_evolution_method,"exact").ne.0)) + doShift = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(shift_evolution_method,"exact").ne.0)) + + call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d") + call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d") + call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d") + + rmx = min(xmx,ymx,zmx) + + if (exblend_rout.lt.0) then + router = rmx - 2.0d0*dx + else + router = exblend_rout + endif + + if (exblend_width.lt.0) then + rinner = router + exblend_width*dx + else + rinner = router - exblend_width + endif + + do k=1,nz + do j=1,ny + do i=1,nx + +c We only do anything if r >= rinner so only evaluate exact data +c there. + + if (r(i,j,k) .ge. rinner) then + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), time, + $ gxxe, gyye, gzze, gxye, gyze, gxze, + $ kxxe, kyye, kzze, kxye, kyze, kxze, + $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, + $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, + $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, + $ alpe, axe, aye, aze, betaxe, betaye, betaze, + $ bxxe, bxye, bxze, byxe, + $ byye, byze, bzxe, bzye, bzze) + +c This sucks, but we want the exact vs so we can blend them also. + + det = -(gxze**2*gyye) + & + 2.d0*gxye*gxze*gyze + & - gxxe*gyze**2 + & - gxye**2*gzze + & + gxxe*gyye*gzze + + uxx=(-gyze**2 + gyye*gzze)/det + uxy=(gxze*gyze - gxye*gzze)/det + uyy=(-gxze**2 + gxxe*gzze)/det + uxz=(-gxze*gyye + gxye*gyze)/det + uyz=(gxye*gxze - gxxe*gyze)/det + uzz=(-gxye**2 + gxxe*gyye)/det + +c Outside of router we want to place exact data on our grid + + if (r(i,j,k) .gt. router) then + +c This is one of those things I will invariably screw up if I type +c it in so let the computer do it + +#define exassign(q) q(i,j,k) = q e +#define exassign_grp(p) \ + exassign(p xx) &&\ + exassign(p xy) &&\ + exassign(p xz) &&\ + exassign(p yy) &&\ + exassign(p yz) &&\ + exassign(p zz) + +c Note this plays on the nasty trick that fortran doesnt give a +c damn about spaces so gxx e is the same as gxxe for the parser... +c Grody but effective! + + if (doGij) then + exassign_grp(g) + endif + if (doKij) then + exassign_grp(k) + endif + + if (doLapse) then + exassign(alp) + endif + + if (doShift.and.(shift_state.ne.0)) then + exassign(betax) + exassign(betay) + exassign(betaz) + endif + +c OK so we dont want to blend so use a goto to jump. + else + +c Evaluate the linear weighting fraction. Obvious... + + frac = (r(i,j,k) - rinner) / (router - rinner) + onemfrac = 1.0D0 - frac + +c Once again some c-preprocessor tricks based on the whole fortran +c space thing... + +#define INTPOINT(f,v) f(i,j,k) = frac * v + onemfrac * f(i,j,k) +#define intone(f) INTPOINT(f, f e) +#define int_grp(p) \ + intone(p xx) &&\ + intone(p xy) &&\ + intone(p xz) &&\ + intone(p yy) &&\ + intone(p yz) &&\ + intone(p zz) + + if (doGij) then + int_grp(g) + endif + if (doKij) then + int_grp(k) + endif + + if (doLapse) then + intone(alp) + endif + + if (doShift.and.(shift_state.ne.0)) then + intone(betax) + intone(betay) + intone(betaz) + endif + + endif ! r > router else + endif ! r > rinner + + enddo + enddo + enddo + + return + end diff --git a/src/boundary.F b/src/boundary.F deleted file mode 100644 index 1862e43..0000000 --- a/src/boundary.F +++ /dev/null @@ -1,125 +0,0 @@ -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - - subroutine Exact__boundary(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - - integer i,j,k - integer nx,ny,nz - - CCTK_REAL tplusone - CCTK_REAL - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ axjunk, ayjunk, azjunk, - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk - -C Grid parameters. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - -C Set all initial data including dijk and vi on all points which -C are on the boundary of the domain if it really is the boundary -C of the complete grid. Treat all six sides of the grid cube this way. - -c Set t = time + dt. This is necessary here because by the time -c we reach this point the geometry has been evolved one time step -c but the variable `time' still hasn't been updated. - - tplusone = cctk_time + cctk_delta_time - -C Note we also always set the lapse and shift at the boundaries at -C time t+1. This is to provide boundary conditions for testing -C elliptic gauge conditions. If they are not used, they will be -C overwritten by Exact__gauge. - -#define EXACTDATAPOINT \ - call Exact__Bona_Masso_data( \ - decoded_exact_model, \ - x(i,j,k), y(i,j,k), z(i,j,k), tplusone, \ - gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), \ - gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), \ - kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), \ - kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), \ - dxgxxjunk, dxgyyjunk, dxgzzjunk, \ - dxgxyjunk, dxgyzjunk, dxgxzjunk, \ - dygxxjunk, dygyyjunk, dygzzjunk, \ - dygxyjunk, dygyzjunk, dygxzjunk, \ - dzgxxjunk, dzgyyjunk, dzgzzjunk, \ - dzgxyjunk, dzgyzjunk, dzgxzjunk, \ - alp(i,j,k), axjunk, ayjunk, azjunk, \ - betax(i,j,k), betay(i,j,k), betaz(i,j,k), \ - bxxjunk, bxyjunk, bxzjunk, \ - byxjunk, byyjunk, byzjunk, \ - bzxjunk, bzyjunk, bzzjunk) - - if (cctk_bbox(1) .eq. 1) then - i=1 - do j=1,ny - do k=1,nz - EXACTDATAPOINT - end do - end do - end if - - if (cctk_bbox(2) .eq. 1) then - i=nx - do j=1,ny - do k=1,nz - EXACTDATAPOINT - end do - end do - end if - - if (cctk_bbox(3) .eq. 1) then - j=1 - do i=1,nx - do k=1,nz - EXACTDATAPOINT - end do - end do - end if - - if (cctk_bbox(4) .eq. 1) then - j=ny - do i=1,nx - do k=1,nz - EXACTDATAPOINT - end do - end do - end if - - if (cctk_bbox(5) .eq. 1) then - k=1 - do j=1,ny - do i=1,nx - EXACTDATAPOINT - end do - end do - end if - - if (cctk_bbox(6) .eq. 1) then - k=nz - do j=1,ny - do i=1,nx - EXACTDATAPOINT - end do - end do - end if - - return - end diff --git a/src/boundary.F77 b/src/boundary.F77 new file mode 100644 index 0000000..1862e43 --- /dev/null +++ b/src/boundary.F77 @@ -0,0 +1,125 @@ +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + subroutine Exact__boundary(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL tplusone + CCTK_REAL + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ axjunk, ayjunk, azjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Set all initial data including dijk and vi on all points which +C are on the boundary of the domain if it really is the boundary +C of the complete grid. Treat all six sides of the grid cube this way. + +c Set t = time + dt. This is necessary here because by the time +c we reach this point the geometry has been evolved one time step +c but the variable `time' still hasn't been updated. + + tplusone = cctk_time + cctk_delta_time + +C Note we also always set the lapse and shift at the boundaries at +C time t+1. This is to provide boundary conditions for testing +C elliptic gauge conditions. If they are not used, they will be +C overwritten by Exact__gauge. + +#define EXACTDATAPOINT \ + call Exact__Bona_Masso_data( \ + decoded_exact_model, \ + x(i,j,k), y(i,j,k), z(i,j,k), tplusone, \ + gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), \ + gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), \ + kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), \ + kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), \ + dxgxxjunk, dxgyyjunk, dxgzzjunk, \ + dxgxyjunk, dxgyzjunk, dxgxzjunk, \ + dygxxjunk, dygyyjunk, dygzzjunk, \ + dygxyjunk, dygyzjunk, dygxzjunk, \ + dzgxxjunk, dzgyyjunk, dzgzzjunk, \ + dzgxyjunk, dzgyzjunk, dzgxzjunk, \ + alp(i,j,k), axjunk, ayjunk, azjunk, \ + betax(i,j,k), betay(i,j,k), betaz(i,j,k), \ + bxxjunk, bxyjunk, bxzjunk, \ + byxjunk, byyjunk, byzjunk, \ + bzxjunk, bzyjunk, bzzjunk) + + if (cctk_bbox(1) .eq. 1) then + i=1 + do j=1,ny + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(2) .eq. 1) then + i=nx + do j=1,ny + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(3) .eq. 1) then + j=1 + do i=1,nx + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(4) .eq. 1) then + j=ny + do i=1,nx + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(5) .eq. 1) then + k=1 + do j=1,ny + do i=1,nx + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(6) .eq. 1) then + k=nz + do j=1,ny + do i=1,nx + EXACTDATAPOINT + end do + end do + end if + + return + end diff --git a/src/decode_pars.F b/src/decode_pars.F deleted file mode 100644 index e81fd8a..0000000 --- a/src/decode_pars.F +++ /dev/null @@ -1,201 +0,0 @@ -c/*@@ -c @file decode_pars.F -c @date Fri Jun 7 19:47:46 CEST 2002 -c @author Jonathan Thornburg -c @desc -c Decode/copy parameters for this thorn into grid scalars -c so we can share this with friends -c so the Calc_Tmunu code in ../include/Scalar_CalcTmunu.inc -c can use them in computing the stress-energy tensor -c -c Actually we only have to copy those parameters which are -c used by the Calc_Tmunu code. For simplicity we decode -c exact_model into the integer decoded_exact_model, then -c copy only the per-model parameters for those models which -c have stress-energy tensor code in ../include/Scalar_CalcTmunu.inc . -c -c @enddesc -c @version $Header$ -c@@*/ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - -#include "param_defs.inc" - -c/*@@ -c @routine Exact__decode_pars -c @date Fri Jun 7 19:47:46 CEST 2002 -c @author Jonathan Thornburg -c @desc -c Decode/copy parameters for this thorn into grid scalars, so -c we can share this with friends for the use of the Calc_Tmunu -c code in ../include/Scalar_CalcTmunu.inc in computing the -c stress-energy tensor. -c @enddesc -c @version $Header$ -c@@*/ - subroutine Exact__decode_pars(CCTK_ARGUMENTS) - implicit none - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c decode exact_model into the integer decoded_exact_model -c - -c Minkowski spacetime - if (CCTK_Equals(exact_model, "Minkowski") .ne. 0) then - decoded_exact_model = EXACT__Minkowski - elseif (CCTK_Equals(exact_model, "Minkowski/shift") .ne. 0) then - decoded_exact_model = EXACT__Minkowski_shift - elseif (CCTK_Equals(exact_model, "Minkowski/funny") .ne. 0) then - decoded_exact_model = EXACT__Minkowski_funny - elseif (CCTK_Equals(exact_model, "Minkowski/gauge wave") .ne. 0) then - decoded_exact_model = EXACT__Minkowski_gauge_wave - -c black hole spacetimes - elseif (CCTK_Equals(exact_model, "Schwarzschild/EF") .ne. 0) then - decoded_exact_model = EXACT__Schwarzschild_EF - elseif (CCTK_Equals(exact_model, "Schwarzschild/PG") .ne. 0) then - decoded_exact_model = EXACT__Schwarzschild_PG - elseif (CCTK_Equals(exact_model, "Schwarzschild/Novikov") .ne. 0) then - decoded_exact_model = EXACT__Schwarzschild_Novikov - elseif (CCTK_Equals(exact_model, "Kerr/Boyer-Lindquist") .ne. 0) then - decoded_exact_model = EXACT__Kerr_BoyerLindquist - elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild") .ne. 0) then - decoded_exact_model = EXACT__Kerr_KerrSchild - elseif (CCTK_Equals(exact_model, "Schwarzschild-Lemaitre") .ne. 0) then - decoded_exact_model = EXACT__Schwarzschild_Lemaitre - elseif (CCTK_Equals(exact_model, "multi-BH") .ne. 0) then - decoded_exact_model = EXACT__multi_BH - elseif (CCTK_Equals(exact_model, "Alvi") .ne. 0) then - decoded_exact_model = EXACT__Alvi - elseif (CCTK_Equals(exact_model, "Thorne-fakebinary") .ne. 0) then - decoded_exact_model = EXACT__Thorne_fakebinary - -c cosmological spacetimes - elseif (CCTK_Equals(exact_model, "Lemaitre") .ne. 0) then - decoded_exact_model = EXACT__Lemaitre - elseif (CCTK_Equals(exact_model, "Robertson-Walker") .ne. 0) then - decoded_exact_model = EXACT__Robertson_Walker - elseif (CCTK_Equals(exact_model, "de Sitter") .ne. 0) then - decoded_exact_model = EXACT__de_Sitter - elseif (CCTK_Equals(exact_model, "de Sitter+Lambda") .ne. 0) then - decoded_exact_model = EXACT__de_Sitter_Lambda - elseif (CCTK_Equals(exact_model, "anti-de Sitter+Lambda") .ne. 0) then - decoded_exact_model = EXACT__anti_de_Sitter_Lambda - elseif (CCTK_Equals(exact_model, "Bianchi I") .ne. 0) then - decoded_exact_model = EXACT__Bianchi_I - elseif (CCTK_Equals(exact_model, "Goedel") .ne. 0) then - decoded_exact_model = EXACT__Goedel - elseif (CCTK_Equals(exact_model, "Bertotti") .ne. 0) then - decoded_exact_model = EXACT__Bertotti - elseif (CCTK_Equals(exact_model, "Kasner-like") .ne. 0) then - decoded_exact_model = EXACT__Kasner_like - elseif (CCTK_Equals(exact_model, "Kasner-axisymmetric") .ne. 0) then - decoded_exact_model = EXACT__Kasner_axisymmetric - elseif (CCTK_Equals(exact_model, "Kasner-generalized") .ne. 0) then - decoded_exact_model = EXACT__Kasner_generalized - elseif (CCTK_Equals(exact_model, "Milne") .ne. 0) then - decoded_exact_model = EXACT__Milne - -c miscellaneous spacetimes - elseif (CCTK_Equals(exact_model, "boost-rotation symmetric") .ne. 0) then - decoded_exact_model = EXACT__boost_rotation_symmetric - elseif (CCTK_Equals(exact_model, "bowl") .ne. 0) then - decoded_exact_model = EXACT__bowl - elseif (CCTK_Equals(exact_model, "constant density star") .ne. 0) then - decoded_exact_model = EXACT__constant_density_star - else - call CCTK_WARN(0, "invalid exact_model parameter!") - endif - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for Schwarzschild-Lemaitre metric -c (Schwarzschiled black hole with cosmological constant) -c - Schwarzschild_Lemaitre___Lambda = Schwarzschild_Lemaitre__Lambda - Schwarzschild_Lemaitre___mass = Schwarzschild_Lemaitre__mass - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for Lemaitre-type spacetime -c - Lemaitre___kappa = Lemaitre__kappa - Lemaitre___Lambda = Lemaitre__Lambda - Lemaitre___epsilon0 = Lemaitre__epsilon0 - Lemaitre___R0 = Lemaitre__R0 - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for Robertson-Walker spacetime -c - Robertson_Walker___R0 = Robertson_Walker__R0 - Robertson_Walker___rho = Robertson_Walker__rho - Robertson_Walker___k = Robertson_Walker__k - Robertson_Walker___pressure = Robertson_Walker__pressure - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for de Sitter spacetime -c - de_Sitter___scale = de_Sitter__scale - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for de Sitter spacetime with cosmological constant -c - de_Sitter_Lambda___scale = de_Sitter_Lambda__scale - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for anti-de Sitter spacetime with cosmological constant -c - anti_de_Sitter_Lambda___scale = anti_de_Sitter_Lambda__scale - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for Bertotti spacetime -c - Bertotti___Lambda = Bertotti__Lambda - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for Kasner-like spacetime -c - Kasner_like___q = Kasner_like__q - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for generalized Kasner spacetime -c - Kasner_generalized___p1 = Kasner_generalized__p1 - Kasner_generalized___p2 = Kasner_generalized__p2 - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c parameters for constant density (Schwarzschild) star -c - constant_density_star___mass = constant_density_star__mass - constant_density_star___radius = constant_density_star__radius - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - return - end diff --git a/src/decode_pars.F77 b/src/decode_pars.F77 new file mode 100644 index 0000000..e81fd8a --- /dev/null +++ b/src/decode_pars.F77 @@ -0,0 +1,201 @@ +c/*@@ +c @file decode_pars.F +c @date Fri Jun 7 19:47:46 CEST 2002 +c @author Jonathan Thornburg +c @desc +c Decode/copy parameters for this thorn into grid scalars +c so we can share this with friends +c so the Calc_Tmunu code in ../include/Scalar_CalcTmunu.inc +c can use them in computing the stress-energy tensor +c +c Actually we only have to copy those parameters which are +c used by the Calc_Tmunu code. For simplicity we decode +c exact_model into the integer decoded_exact_model, then +c copy only the per-model parameters for those models which +c have stress-energy tensor code in ../include/Scalar_CalcTmunu.inc . +c +c @enddesc +c @version $Header$ +c@@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#include "param_defs.inc" + +c/*@@ +c @routine Exact__decode_pars +c @date Fri Jun 7 19:47:46 CEST 2002 +c @author Jonathan Thornburg +c @desc +c Decode/copy parameters for this thorn into grid scalars, so +c we can share this with friends for the use of the Calc_Tmunu +c code in ../include/Scalar_CalcTmunu.inc in computing the +c stress-energy tensor. +c @enddesc +c @version $Header$ +c@@*/ + subroutine Exact__decode_pars(CCTK_ARGUMENTS) + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c decode exact_model into the integer decoded_exact_model +c + +c Minkowski spacetime + if (CCTK_Equals(exact_model, "Minkowski") .ne. 0) then + decoded_exact_model = EXACT__Minkowski + elseif (CCTK_Equals(exact_model, "Minkowski/shift") .ne. 0) then + decoded_exact_model = EXACT__Minkowski_shift + elseif (CCTK_Equals(exact_model, "Minkowski/funny") .ne. 0) then + decoded_exact_model = EXACT__Minkowski_funny + elseif (CCTK_Equals(exact_model, "Minkowski/gauge wave") .ne. 0) then + decoded_exact_model = EXACT__Minkowski_gauge_wave + +c black hole spacetimes + elseif (CCTK_Equals(exact_model, "Schwarzschild/EF") .ne. 0) then + decoded_exact_model = EXACT__Schwarzschild_EF + elseif (CCTK_Equals(exact_model, "Schwarzschild/PG") .ne. 0) then + decoded_exact_model = EXACT__Schwarzschild_PG + elseif (CCTK_Equals(exact_model, "Schwarzschild/Novikov") .ne. 0) then + decoded_exact_model = EXACT__Schwarzschild_Novikov + elseif (CCTK_Equals(exact_model, "Kerr/Boyer-Lindquist") .ne. 0) then + decoded_exact_model = EXACT__Kerr_BoyerLindquist + elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild") .ne. 0) then + decoded_exact_model = EXACT__Kerr_KerrSchild + elseif (CCTK_Equals(exact_model, "Schwarzschild-Lemaitre") .ne. 0) then + decoded_exact_model = EXACT__Schwarzschild_Lemaitre + elseif (CCTK_Equals(exact_model, "multi-BH") .ne. 0) then + decoded_exact_model = EXACT__multi_BH + elseif (CCTK_Equals(exact_model, "Alvi") .ne. 0) then + decoded_exact_model = EXACT__Alvi + elseif (CCTK_Equals(exact_model, "Thorne-fakebinary") .ne. 0) then + decoded_exact_model = EXACT__Thorne_fakebinary + +c cosmological spacetimes + elseif (CCTK_Equals(exact_model, "Lemaitre") .ne. 0) then + decoded_exact_model = EXACT__Lemaitre + elseif (CCTK_Equals(exact_model, "Robertson-Walker") .ne. 0) then + decoded_exact_model = EXACT__Robertson_Walker + elseif (CCTK_Equals(exact_model, "de Sitter") .ne. 0) then + decoded_exact_model = EXACT__de_Sitter + elseif (CCTK_Equals(exact_model, "de Sitter+Lambda") .ne. 0) then + decoded_exact_model = EXACT__de_Sitter_Lambda + elseif (CCTK_Equals(exact_model, "anti-de Sitter+Lambda") .ne. 0) then + decoded_exact_model = EXACT__anti_de_Sitter_Lambda + elseif (CCTK_Equals(exact_model, "Bianchi I") .ne. 0) then + decoded_exact_model = EXACT__Bianchi_I + elseif (CCTK_Equals(exact_model, "Goedel") .ne. 0) then + decoded_exact_model = EXACT__Goedel + elseif (CCTK_Equals(exact_model, "Bertotti") .ne. 0) then + decoded_exact_model = EXACT__Bertotti + elseif (CCTK_Equals(exact_model, "Kasner-like") .ne. 0) then + decoded_exact_model = EXACT__Kasner_like + elseif (CCTK_Equals(exact_model, "Kasner-axisymmetric") .ne. 0) then + decoded_exact_model = EXACT__Kasner_axisymmetric + elseif (CCTK_Equals(exact_model, "Kasner-generalized") .ne. 0) then + decoded_exact_model = EXACT__Kasner_generalized + elseif (CCTK_Equals(exact_model, "Milne") .ne. 0) then + decoded_exact_model = EXACT__Milne + +c miscellaneous spacetimes + elseif (CCTK_Equals(exact_model, "boost-rotation symmetric") .ne. 0) then + decoded_exact_model = EXACT__boost_rotation_symmetric + elseif (CCTK_Equals(exact_model, "bowl") .ne. 0) then + decoded_exact_model = EXACT__bowl + elseif (CCTK_Equals(exact_model, "constant density star") .ne. 0) then + decoded_exact_model = EXACT__constant_density_star + else + call CCTK_WARN(0, "invalid exact_model parameter!") + endif + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for Schwarzschild-Lemaitre metric +c (Schwarzschiled black hole with cosmological constant) +c + Schwarzschild_Lemaitre___Lambda = Schwarzschild_Lemaitre__Lambda + Schwarzschild_Lemaitre___mass = Schwarzschild_Lemaitre__mass + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for Lemaitre-type spacetime +c + Lemaitre___kappa = Lemaitre__kappa + Lemaitre___Lambda = Lemaitre__Lambda + Lemaitre___epsilon0 = Lemaitre__epsilon0 + Lemaitre___R0 = Lemaitre__R0 + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for Robertson-Walker spacetime +c + Robertson_Walker___R0 = Robertson_Walker__R0 + Robertson_Walker___rho = Robertson_Walker__rho + Robertson_Walker___k = Robertson_Walker__k + Robertson_Walker___pressure = Robertson_Walker__pressure + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for de Sitter spacetime +c + de_Sitter___scale = de_Sitter__scale + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for de Sitter spacetime with cosmological constant +c + de_Sitter_Lambda___scale = de_Sitter_Lambda__scale + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for anti-de Sitter spacetime with cosmological constant +c + anti_de_Sitter_Lambda___scale = anti_de_Sitter_Lambda__scale + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for Bertotti spacetime +c + Bertotti___Lambda = Bertotti__Lambda + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for Kasner-like spacetime +c + Kasner_like___q = Kasner_like__q + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for generalized Kasner spacetime +c + Kasner_generalized___p1 = Kasner_generalized__p1 + Kasner_generalized___p2 = Kasner_generalized__p2 + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c parameters for constant density (Schwarzschild) star +c + constant_density_star___mass = constant_density_star__mass + constant_density_star___radius = constant_density_star__radius + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + return + end diff --git a/src/gauge.F b/src/gauge.F deleted file mode 100644 index f5c2e88..0000000 --- a/src/gauge.F +++ /dev/null @@ -1,152 +0,0 @@ -C This routine sets the lapse and/or shift by calling a routine -C that does it pointwise. Note that it could be easily modified -C to set the Bona-Masso variables B_xx etc. -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - - subroutine Exact__gauge(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - integer i,j,k - integer nx,ny,nz - - CCTK_REAL tplushalf - CCTK_REAL gxxjunk, gyyjunk, gzzjunk, - $ gxyjunk, gyzjunk, gxzjunk, - $ hxxjunk, hyyjunk, hzzjunk, - $ hxyjunk, hyzjunk, hxzjunk, - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ alpjunk, axjunk, ayjunk, azjunk, - $ betaxjunk, betayjunk, betazjunk, - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk - -C Grid parameters. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - -C Set t = time + dt/2. This is for consistency with a second -C order numerical scheme. - - tplushalf = cctk_time + 0.5D0*cctk_delta_time - -C Set both lapse and shift. - - if (( (CCTK_Equals(lapse_evolution_method,"exact").ne.0) .and. - $ (CCTK_Equals(shift_evolution_method,"exact").ne.0) ) - $ .or. - $ ( (CCTK_Equals(initial_lapse,"exact").ne.0) .and. - $ (CCTK_Equals(initial_shift,"exact").ne.0) )) then - - do k=1,nz - do j=1,ny - do i=1,nx - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, - $ gxxjunk, gyyjunk, gzzjunk, - $ gxyjunk, gyzjunk, gxzjunk, - $ hxxjunk, hyyjunk, hzzjunk, - $ hxyjunk, hyzjunk, hxzjunk, - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ alp(i,j,k), axjunk, ayjunk, azjunk, - $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk) - end do - end do - end do - -C Set lapse only. - - else if ((CCTK_Equals(lapse_evolution_method,"exact").ne.0) - $ .or. (CCTK_Equals(initial_lapse,"exact").ne.0)) then - - do k=1,nz - do j=1,ny - do i=1,nx - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, - $ gxxjunk, gyyjunk, gzzjunk, - $ gxyjunk, gyzjunk, gxzjunk, - $ hxxjunk, hyyjunk, hzzjunk, - $ hxyjunk, hyzjunk, hxzjunk, - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ alp(i,j,k), axjunk, ayjunk, azjunk, - $ betaxjunk, betayjunk, betazjunk, - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk) - - end do - end do - end do - -C Set shift only. - - else if ((CCTK_Equals(shift_evolution_method,"exact").ne.0) - $ .or. (CCTK_Equals(initial_shift,"exact").ne.0)) then - - do k=1,nz - do j=1,ny - do i=1,nx - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, - $ gxxjunk, gyyjunk, gzzjunk, - $ gxyjunk, gyzjunk, gxzjunk, - $ hxxjunk, hyyjunk, hzzjunk, - $ hxyjunk, hyzjunk, hxzjunk, - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ alpjunk, axjunk, ayjunk, azjunk, - $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk) - - end do - end do - end do - - else - call CCTK_WARN(1,'Exact__gauge has been called without doing anything') - end if - - return - end diff --git a/src/gauge.F77 b/src/gauge.F77 new file mode 100644 index 0000000..f5c2e88 --- /dev/null +++ b/src/gauge.F77 @@ -0,0 +1,152 @@ +C This routine sets the lapse and/or shift by calling a routine +C that does it pointwise. Note that it could be easily modified +C to set the Bona-Masso variables B_xx etc. +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + subroutine Exact__gauge(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL tplushalf + CCTK_REAL gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Set t = time + dt/2. This is for consistency with a second +C order numerical scheme. + + tplushalf = cctk_time + 0.5D0*cctk_delta_time + +C Set both lapse and shift. + + if (( (CCTK_Equals(lapse_evolution_method,"exact").ne.0) .and. + $ (CCTK_Equals(shift_evolution_method,"exact").ne.0) ) + $ .or. + $ ( (CCTK_Equals(initial_lapse,"exact").ne.0) .and. + $ (CCTK_Equals(initial_shift,"exact").ne.0) )) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alp(i,j,k), axjunk, ayjunk, azjunk, + $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + end do + end do + end do + +C Set lapse only. + + else if ((CCTK_Equals(lapse_evolution_method,"exact").ne.0) + $ .or. (CCTK_Equals(initial_lapse,"exact").ne.0)) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alp(i,j,k), axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + +C Set shift only. + + else if ((CCTK_Equals(shift_evolution_method,"exact").ne.0) + $ .or. (CCTK_Equals(initial_shift,"exact").ne.0)) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + + else + call CCTK_WARN(1,'Exact__gauge has been called without doing anything') + end if + + return + end diff --git a/src/initialize.F b/src/initialize.F deleted file mode 100644 index 00fca60..0000000 --- a/src/initialize.F +++ /dev/null @@ -1,78 +0,0 @@ -C Wrapper for boostrotdata. Calls it and vectorini. -C Sets Cauchy data, lapse and shift, and what else is needed -C in the Bona-Masso formalism, at an initial time. -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - - subroutine Exact__initialize(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - integer i,j,k - integer nx,ny,nz - - CCTK_REAL alpjunk, axjunk, ayjunk, azjunk, - $ betaxjunk, betayjunk, betazjunk, - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk - CCTK_REAL - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk - - call CCTK_INFO('setting exact data on slice') - -C Note I assume time has been initialized to physical time. -C Set data pointwise. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - - do k=1,nz - do j=1,ny - do i=1,nx - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, - $ gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), - $ gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), - $ kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), - $ kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), - $ dxgxxjunk, dxgyyjunk, dxgzzjunk, - $ dxgxyjunk, dxgyzjunk, dxgxzjunk, - $ dygxxjunk, dygyyjunk, dygzzjunk, - $ dygxyjunk, dygyzjunk, dygxzjunk, - $ dzgxxjunk, dzgyyjunk, dzgzzjunk, - $ dzgxyjunk, dzgyzjunk, dzgxzjunk, - $ alpjunk, axjunk, ayjunk, azjunk, - $ betaxjunk, betayjunk, betazjunk, - $ bxxjunk, bxyjunk, bxzjunk, - $ byxjunk, byyjunk, byzjunk, - $ bzxjunk, bzyjunk, bzzjunk) - - end do - end do - end do - -C Tell the code there is no need to treat the conformal factor -C as a separate field. That is, we have set the physical metric here. -c Commented out in einstein revamp, now Exact doesnot inherit anything -c about the conformal factor -c conformal_state = 0 -c psi = 1.0D0 - - return - end diff --git a/src/initialize.F77 b/src/initialize.F77 new file mode 100644 index 0000000..00fca60 --- /dev/null +++ b/src/initialize.F77 @@ -0,0 +1,78 @@ +C Wrapper for boostrotdata. Calls it and vectorini. +C Sets Cauchy data, lapse and shift, and what else is needed +C in the Bona-Masso formalism, at an initial time. +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + subroutine Exact__initialize(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + CCTK_REAL + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk + + call CCTK_INFO('setting exact data on slice') + +C Note I assume time has been initialized to physical time. +C Set data pointwise. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do k=1,nz + do j=1,ny + do i=1,nx + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), + $ gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), + $ kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), + $ kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + +C Tell the code there is no need to treat the conformal factor +C as a separate field. That is, we have set the physical metric here. +c Commented out in einstein revamp, now Exact doesnot inherit anything +c about the conformal factor +c conformal_state = 0 +c psi = 1.0D0 + + return + end diff --git a/src/make.code.defn b/src/make.code.defn index 136b470..54d7eef 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -4,20 +4,20 @@ # Source files in this directory SRCS = ParamCheck.c \ Startup.c \ - decode_pars.F \ - initialize.F \ + decode_pars.F77 \ + initialize.F77 \ \ slice_initialize.F \ slice_evolve.F \ slice_data.F \ \ - gauge.F \ - Bona_Masso_data.F \ - metric.F \ + gauge.F77 \ + Bona_Masso_data.F77 \ + metric.F77 \ \ - boundary.F \ - blended_boundary.F \ - xyz_blended_boundary.F \ + boundary.F77 \ + blended_boundary.F77 \ + xyz_blended_boundary.F77 \ linear_extrap_one_bndry.F # Subdirectories containing source files to be compiled diff --git a/src/metric.F b/src/metric.F deleted file mode 100644 index 23797f7..0000000 --- a/src/metric.F +++ /dev/null @@ -1,293 +0,0 @@ -c This subroutine calculates the 4-metric and its inverse at an event, -c by decoding decoded_exact_model and calling the appropriate subroutine -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - -#include "param_defs.inc" - - subroutine Exact__metric( - $ decoded_exact_model, - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz,rama) - - implicit none - DECLARE_CCTK_FUNCTIONS - -c arguments - CCTK_INT decoded_exact_model - CCTK_REAL x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, rama - -c local variables - character*100 warn_buffer - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c Minkowski spacetime -c - - if (decoded_exact_model .eq. EXACT__Minkowski) then - call Exact__Minkowski( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Minkowski_shift) then - call Exact__Minkowski_shift( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Minkowski_funny) then - call Exact__Minkowski_funny( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Minkowski_gauge_wave) then - call Exact__Minkowski_gauge_wave( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c black hole spacetimes -c - - elseif (decoded_exact_model .eq. EXACT__Schwarzschild_EF) then - call Exact__Schwarzschild_EF( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Schwarzschild_PG) then - call Exact__Schwarzschild_PG( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Novikov) then - call Exact__Schwarzschild_Novikov(x,y,z,t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Kerr_BoyerLindquist) then - call Exact__Kerr_BoyerLindquist( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild) then - call Exact__Kerr_KerrSchild( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Lemaitre) then - call Exact__Schwarzschild_Lemaitre( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__multi_BH) then - call Exact__multi_BH( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -c -c not fully implemented yet -- see Nina Jansen for details -c -c elseif (decoded_exact_model .eq. EXACT__Alvi) then -c call Exact__Alvi( -c $ x, y, z, t, -c $ gdtt, gdtx, gdty, gdtz, -c $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, -c $ gutt, gutx, guty, gutz, -c $ guxx, guyy, guzz, guxy, guyz, guxz) -c - - elseif (decoded_exact_model .eq. EXACT__Thorne_fakebinary) then - call Exact__Thorne_fakebinary( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c cosmological spacetimes -c - - elseif (decoded_exact_model .eq. EXACT__Lemaitre) then - call Exact__Lemaitre( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Robertson_Walker) then - call Exact__Robertson_Walker( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz,rama) - - elseif (decoded_exact_model .eq. EXACT__de_Sitter) then - call Exact__de_Sitter( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__de_Sitter_Lambda) then - call Exact__de_Sitter_Lambda( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__anti_de_Sitter_Lambda) then - call Exact__anti_de_Sitter_Lambda( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Bianchi_I) then - call Exact__Bianchi_I( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Goedel) then - call Exact__Goedel( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Bertotti) then - call Exact__Bertotti( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Kasner_like) then - call Exact__Kasner_like( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Kasner_axisymmetric) then - call Exact__Kasner_axisymmetric( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Kasner_generalized) then - call Exact__Kasner_generalized( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__Milne) then - call Exact__Milne( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c -c miscellaneous spacetimes -c - - elseif (decoded_exact_model .eq. EXACT__boost_rotation_symmetric) then - call Exact__boost_rotation_symmetric( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__bowl) then - call Exact__bowl( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - - elseif (decoded_exact_model .eq. EXACT__constant_density_star) then - call Exact__constant_density_star( - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz) - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - else - write (warn_buffer, '(a,i8)') - $ 'Unknown decoded_exact_model = ', decoded_exact_model - call CCTK_WARN(0, warn_buffer) - endif - - return - end diff --git a/src/metric.F77 b/src/metric.F77 new file mode 100644 index 0000000..23797f7 --- /dev/null +++ b/src/metric.F77 @@ -0,0 +1,293 @@ +c This subroutine calculates the 4-metric and its inverse at an event, +c by decoding decoded_exact_model and calling the appropriate subroutine +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#include "param_defs.inc" + + subroutine Exact__metric( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz,rama) + + implicit none + DECLARE_CCTK_FUNCTIONS + +c arguments + CCTK_INT decoded_exact_model + CCTK_REAL x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, rama + +c local variables + character*100 warn_buffer + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c Minkowski spacetime +c + + if (decoded_exact_model .eq. EXACT__Minkowski) then + call Exact__Minkowski( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Minkowski_shift) then + call Exact__Minkowski_shift( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Minkowski_funny) then + call Exact__Minkowski_funny( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Minkowski_gauge_wave) then + call Exact__Minkowski_gauge_wave( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c black hole spacetimes +c + + elseif (decoded_exact_model .eq. EXACT__Schwarzschild_EF) then + call Exact__Schwarzschild_EF( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Schwarzschild_PG) then + call Exact__Schwarzschild_PG( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Novikov) then + call Exact__Schwarzschild_Novikov(x,y,z,t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Kerr_BoyerLindquist) then + call Exact__Kerr_BoyerLindquist( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild) then + call Exact__Kerr_KerrSchild( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Lemaitre) then + call Exact__Schwarzschild_Lemaitre( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__multi_BH) then + call Exact__multi_BH( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +c +c not fully implemented yet -- see Nina Jansen for details +c +c elseif (decoded_exact_model .eq. EXACT__Alvi) then +c call Exact__Alvi( +c $ x, y, z, t, +c $ gdtt, gdtx, gdty, gdtz, +c $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, +c $ gutt, gutx, guty, gutz, +c $ guxx, guyy, guzz, guxy, guyz, guxz) +c + + elseif (decoded_exact_model .eq. EXACT__Thorne_fakebinary) then + call Exact__Thorne_fakebinary( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c cosmological spacetimes +c + + elseif (decoded_exact_model .eq. EXACT__Lemaitre) then + call Exact__Lemaitre( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Robertson_Walker) then + call Exact__Robertson_Walker( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz,rama) + + elseif (decoded_exact_model .eq. EXACT__de_Sitter) then + call Exact__de_Sitter( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__de_Sitter_Lambda) then + call Exact__de_Sitter_Lambda( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__anti_de_Sitter_Lambda) then + call Exact__anti_de_Sitter_Lambda( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Bianchi_I) then + call Exact__Bianchi_I( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Goedel) then + call Exact__Goedel( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Bertotti) then + call Exact__Bertotti( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Kasner_like) then + call Exact__Kasner_like( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Kasner_axisymmetric) then + call Exact__Kasner_axisymmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Kasner_generalized) then + call Exact__Kasner_generalized( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__Milne) then + call Exact__Milne( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c miscellaneous spacetimes +c + + elseif (decoded_exact_model .eq. EXACT__boost_rotation_symmetric) then + call Exact__boost_rotation_symmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__bowl) then + call Exact__bowl( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (decoded_exact_model .eq. EXACT__constant_density_star) then + call Exact__constant_density_star( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + else + write (warn_buffer, '(a,i8)') + $ 'Unknown decoded_exact_model = ', decoded_exact_model + call CCTK_WARN(0, warn_buffer) + endif + + return + end diff --git a/src/xyz_blended_boundary.F b/src/xyz_blended_boundary.F deleted file mode 100644 index d596e79..0000000 --- a/src/xyz_blended_boundary.F +++ /dev/null @@ -1,210 +0,0 @@ -C $Header$ - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" - - subroutine Exact__xyz_blended_boundary(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - logical doKij, doGij, doDs, doVs, doLapse, doShift - - integer i,j,k - integer nx,ny,nz - integer ninterps - - CCTK_REAL xblend, yblend, zblend - CCTK_REAL xmin, xmax, ymin, ymax, zmin, zmax - CCTK_REAL xfrac, yfrac, zfrac, onemxfrac, onemyfrac, onemzfrac - CCTK_REAL oonints - CCTK_REAL sfrac, onemsfrac - CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze - CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze - CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze - CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze - CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze - CCTK_REAL alpe, axe, aye, aze - CCTK_REAL betaxe,betaye,betaze - CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze - CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz - CCTK_REAL vxe,vye,vze,sav - - CCTK_REAL :: dx,dy,dz,time,ierr - -C Grid parameters. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - - dx = cctk_delta_space(1) - dy = cctk_delta_space(2) - dz = cctk_delta_space(3) - - time = cctk_time - -C Other parameters. - - doKij = (exblend_Ks.eq.1) - doGij = (exblend_gs.eq.1) - - doLapse = ((exblend_gauge.eq.1).and. - $ (CCTK_Equals(lapse_evolution_method,"exact").ne.0)) - doShift = ((exblend_gauge.eq.1).and. - $ (CCTK_Equals(shift_evolution_method,"exact").ne.0)) - - if (exblend_width.lt.0) then - xblend = - exblend_width*dx - yblend = - exblend_width*dy - zblend = - exblend_width*dz - else - xblend = exblend_width - yblend = exblend_width - zblend = exblend_width - endif - - call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,-1,"x","cart3d") - call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,-1,"y","cart3d") - call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,-1,"z","cart3d") - - do k=1,nz - do j=1,ny - do i=1,nx - -c We only do anything if in the blending region - - if (x(i,j,k) .ge. xmax - xblend .or. - $ x(i,j,k) .le. xmin + xblend .or. - $ y(i,j,k) .ge. ymax - yblend .or. - $ y(i,j,k) .le. ymin + yblend .or. - $ z(i,j,k) .ge. zmax - zblend .or. - $ z(i,j,k) .le. zmin + zblend) then - - call Exact__Bona_Masso_data( - $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), time, - $ gxxe, gyye, gzze, gxye, gyze, gxze, - $ kxxe, kyye, kzze, kxye, kyze, kxze, - $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, - $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, - $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, - $ alpe, axe, aye, aze, betaxe, betaye, betaze, - $ bxxe, bxye, bxze, byxe, - $ byye, byze, bzxe, bzye, bzze) - -c This sucks, but we want the exact vs so we can blend them also. - - det = -(gxze**2*gyye) - & + 2.d0*gxye*gxze*gyze - & - gxxe*gyze**2 - & - gxye**2*gzze - & + gxxe*gyye*gzze - - uxx=(-gyze**2 + gyye*gzze)/det - uxy=(gxze*gyze - gxye*gzze)/det - uyy=(-gxze**2 + gxxe*gzze)/det - uxz=(-gxze*gyye + gxye*gyze)/det - uyz=(gxye*gxze - gxxe*gyze)/det - uzz=(-gxye**2 + gxxe*gyye)/det - -c OK so 6 blending cases. If frac = 1 we get all exact - - ninterps = 0 - - xfrac = 0.0D0 - onemxfrac = 0.0D0 - yfrac = 0.0D0 - onemyfrac = 0.0D0 - zfrac = 0.0D0 - onemzfrac = 0.0D0 - - if (x(i,j,k) .le. xmin + xblend) then - xfrac = 1.0D0 - (x(i,j,k)-xmin) / xblend - onemxfrac = 1.0D0 - xfrac - ninterps = ninterps + 1 - endif - - if (x(i,j,k) .ge. xmax - xblend) then - xfrac = 1.0D0 - (xmax - x(i,j,k)) / xblend - onemxfrac = 1.0D0 - xfrac - ninterps = ninterps + 1 - endif - - if (y(i,j,k) .le. ymin + yblend) then - yfrac = 1.0D0 - (y(i,j,k)-ymin) / yblend - onemyfrac = 1.0D0 - yfrac - ninterps = ninterps + 1 - endif - - if (y(i,j,k) .ge. ymax - yblend) then - yfrac = 1.0D0 - (ymax - y(i,j,k)) / yblend - onemyfrac = 1.0D0 - yfrac - ninterps = ninterps + 1 - endif - - if (z(i,j,k) .le. zmin + zblend) then - zfrac = 1.0D0 - (z(i,j,k)-zmin) / zblend - onemzfrac = 1.0D0 - zfrac - ninterps = ninterps + 1 - endif - - if (z(i,j,k) .ge. zmax - zblend) then - zfrac = 1.0D0 - (zmax - z(i,j,k)) / zblend - onemzfrac = 1.0D0 - zfrac - ninterps = ninterps + 1 - endif - - oonints = 1.0D0 / ninterps - - if (ninterps .eq. 0 .or. ninterps .gt. 3) then - print *,"NINTERPS error", ninterps - STOP - endif - - sfrac = (xfrac + yfrac + zfrac) * oonints - onemsfrac = 1.0D0 - sfrac - -c Once again some c-preprocessor tricks based on the whole fortran -c space thing... - -#define INTPOINT(f,v) f(i,j,k) = sfrac * v + onemsfrac * f(i,j,k) -#define intone(f) INTPOINT(f, f e) -#define int_grp(p) \ - intone(p xx) &&\ - intone(p xy) &&\ - intone(p xz) &&\ - intone(p yy) &&\ - intone(p yz) &&\ - intone(p zz) - - if (doGij) then - int_grp(g) - endif - - if (doKij) then - int_grp(k) - endif - - if (doLapse) then - intone(alp) - endif - - if (doShift.and.(shift_state.ne.0)) then - intone(betax) - intone(betay) - intone(betaz) - endif - - endif ! r > rinner - - enddo - enddo - enddo - - return - end diff --git a/src/xyz_blended_boundary.F77 b/src/xyz_blended_boundary.F77 new file mode 100644 index 0000000..9c6ef44 --- /dev/null +++ b/src/xyz_blended_boundary.F77 @@ -0,0 +1,211 @@ +C $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + subroutine Exact__xyz_blended_boundary(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + logical doKij, doGij, doDs, doVs, doLapse, doShift + + integer i,j,k + integer nx,ny,nz + integer ninterps + + CCTK_REAL xblend, yblend, zblend + CCTK_REAL xmin, xmax, ymin, ymax, zmin, zmax + CCTK_REAL xfrac, yfrac, zfrac, onemxfrac, onemyfrac, onemzfrac + CCTK_REAL oonints + CCTK_REAL sfrac, onemsfrac + CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze + CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze + CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze + CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze + CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze + CCTK_REAL alpe, axe, aye, aze + CCTK_REAL betaxe,betaye,betaze + CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze + CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL vxe,vye,vze,sav + + CCTK_REAL dx,dy,dz,time + integer ierr + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + time = cctk_time + +C Other parameters. + + doKij = (exblend_Ks.eq.1) + doGij = (exblend_gs.eq.1) + + doLapse = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(lapse_evolution_method,"exact").ne.0)) + doShift = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(shift_evolution_method,"exact").ne.0)) + + if (exblend_width.lt.0) then + xblend = - exblend_width*dx + yblend = - exblend_width*dy + zblend = - exblend_width*dz + else + xblend = exblend_width + yblend = exblend_width + zblend = exblend_width + endif + + call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,-1,"x","cart3d") + call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,-1,"y","cart3d") + call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,-1,"z","cart3d") + + do k=1,nz + do j=1,ny + do i=1,nx + +c We only do anything if in the blending region + + if (x(i,j,k) .ge. xmax - xblend .or. + $ x(i,j,k) .le. xmin + xblend .or. + $ y(i,j,k) .ge. ymax - yblend .or. + $ y(i,j,k) .le. ymin + yblend .or. + $ z(i,j,k) .ge. zmax - zblend .or. + $ z(i,j,k) .le. zmin + zblend) then + + call Exact__Bona_Masso_data( + $ decoded_exact_model, + $ x(i,j,k), y(i,j,k), z(i,j,k), time, + $ gxxe, gyye, gzze, gxye, gyze, gxze, + $ kxxe, kyye, kzze, kxye, kyze, kxze, + $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, + $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, + $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, + $ alpe, axe, aye, aze, betaxe, betaye, betaze, + $ bxxe, bxye, bxze, byxe, + $ byye, byze, bzxe, bzye, bzze) + +c This sucks, but we want the exact vs so we can blend them also. + + det = -(gxze**2*gyye) + & + 2.d0*gxye*gxze*gyze + & - gxxe*gyze**2 + & - gxye**2*gzze + & + gxxe*gyye*gzze + + uxx=(-gyze**2 + gyye*gzze)/det + uxy=(gxze*gyze - gxye*gzze)/det + uyy=(-gxze**2 + gxxe*gzze)/det + uxz=(-gxze*gyye + gxye*gyze)/det + uyz=(gxye*gxze - gxxe*gyze)/det + uzz=(-gxye**2 + gxxe*gyye)/det + +c OK so 6 blending cases. If frac = 1 we get all exact + + ninterps = 0 + + xfrac = 0.0D0 + onemxfrac = 0.0D0 + yfrac = 0.0D0 + onemyfrac = 0.0D0 + zfrac = 0.0D0 + onemzfrac = 0.0D0 + + if (x(i,j,k) .le. xmin + xblend) then + xfrac = 1.0D0 - (x(i,j,k)-xmin) / xblend + onemxfrac = 1.0D0 - xfrac + ninterps = ninterps + 1 + endif + + if (x(i,j,k) .ge. xmax - xblend) then + xfrac = 1.0D0 - (xmax - x(i,j,k)) / xblend + onemxfrac = 1.0D0 - xfrac + ninterps = ninterps + 1 + endif + + if (y(i,j,k) .le. ymin + yblend) then + yfrac = 1.0D0 - (y(i,j,k)-ymin) / yblend + onemyfrac = 1.0D0 - yfrac + ninterps = ninterps + 1 + endif + + if (y(i,j,k) .ge. ymax - yblend) then + yfrac = 1.0D0 - (ymax - y(i,j,k)) / yblend + onemyfrac = 1.0D0 - yfrac + ninterps = ninterps + 1 + endif + + if (z(i,j,k) .le. zmin + zblend) then + zfrac = 1.0D0 - (z(i,j,k)-zmin) / zblend + onemzfrac = 1.0D0 - zfrac + ninterps = ninterps + 1 + endif + + if (z(i,j,k) .ge. zmax - zblend) then + zfrac = 1.0D0 - (zmax - z(i,j,k)) / zblend + onemzfrac = 1.0D0 - zfrac + ninterps = ninterps + 1 + endif + + oonints = 1.0D0 / ninterps + + if (ninterps .eq. 0 .or. ninterps .gt. 3) then + print *,"NINTERPS error", ninterps + STOP + endif + + sfrac = (xfrac + yfrac + zfrac) * oonints + onemsfrac = 1.0D0 - sfrac + +c Once again some c-preprocessor tricks based on the whole fortran +c space thing... + +#define INTPOINT(f,v) f(i,j,k) = sfrac * v + onemsfrac * f(i,j,k) +#define intone(f) INTPOINT(f, f e) +#define int_grp(p) \ + intone(p xx) &&\ + intone(p xy) &&\ + intone(p xz) &&\ + intone(p yy) &&\ + intone(p yz) &&\ + intone(p zz) + + if (doGij) then + int_grp(g) + endif + + if (doKij) then + int_grp(k) + endif + + if (doLapse) then + intone(alp) + endif + + if (doShift.and.(shift_state.ne.0)) then + intone(betax) + intone(betay) + intone(betaz) + endif + + endif ! r > rinner + + enddo + enddo + enddo + + return + end -- cgit v1.2.3