diff options
author | schnetter <schnetter@b716e942-a2de-43ad-8f52-f3dfe468e4e7> | 2002-07-15 17:15:29 +0000 |
---|---|---|
committer | schnetter <schnetter@b716e942-a2de-43ad-8f52-f3dfe468e4e7> | 2002-07-15 17:15:29 +0000 |
commit | 42a98ca518377949e8cfc902343a4b2d97a25564 (patch) | |
tree | 362e565c054f600cb11a4497593fd2e03d6fd328 /src | |
parent | ae7e90702b85ebba7ab33b40c161815cec69ed5a (diff) |
Added thorn TGRtensor. This is a utility thorn from the TGR code. It
contains generic tensor operations, some useful constants, code for
numeric derivatives, interfaces to Lapack and BLAS, and interfaces for
Cactus routines.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@2 b716e942-a2de-43ad-8f52-f3dfe468e4e7
Diffstat (limited to 'src')
-rw-r--r-- | src/cactus.F90 | 389 | ||||
-rw-r--r-- | src/constants.F90 | 52 | ||||
-rw-r--r-- | src/conversion.F90 | 185 | ||||
-rw-r--r-- | src/covderivs.F90 | 201 | ||||
-rw-r--r-- | src/covderivs2.F90 | 204 | ||||
-rw-r--r-- | src/depend.pl | 46 | ||||
-rw-r--r-- | src/derivs.F90 | 46 | ||||
-rw-r--r-- | src/derivs2.F90 | 40 | ||||
-rw-r--r-- | src/hyperslab.F90 | 33 | ||||
-rw-r--r-- | src/lapack.F90 | 51 | ||||
-rw-r--r-- | src/make.code.defn | 25 | ||||
-rw-r--r-- | src/make.code.deps | 11 | ||||
-rw-r--r-- | src/make.configuration.deps | 6 | ||||
-rw-r--r-- | src/matexp.F90 | 38 | ||||
-rw-r--r-- | src/matinv.F90 | 54 | ||||
-rw-r--r-- | src/pointwise.F90 | 252 | ||||
-rw-r--r-- | src/pointwise2.F90 | 148 | ||||
-rw-r--r-- | src/rotation.F90 | 59 | ||||
-rw-r--r-- | src/tensor.F90 | 127 | ||||
-rw-r--r-- | src/tensor2.F90 | 107 |
20 files changed, 2074 insertions, 0 deletions
diff --git a/src/cactus.F90 b/src/cactus.F90 new file mode 100644 index 0000000..5d0c278 --- /dev/null +++ b/src/cactus.F90 @@ -0,0 +1,389 @@ +! $Header$ + +#include "cctk.h" + +module cactus + implicit none + public + + interface + + ! from Cactus: + + subroutine cctk_barrier (ierr, cctkgh) + implicit none + integer ierr + CCTK_POINTER cctkgh + end subroutine cctk_barrier + + subroutine cctk_coordrange (ierr, cctkgh, lower, upper, dir, name, & + & systemname) + implicit none + integer ierr + CCTK_POINTER cctkgh + CCTK_REAL lower + CCTK_REAL upper + integer dir + character*(*) name + character*(*) systemname + end subroutine cctk_coordrange + + subroutine cctk_coordsystemhandle (handle, coordsystem) + implicit none + integer handle + character*(*) coordsystem + end subroutine cctk_coordsystemhandle + + function cctk_equals (arg1, arg2) + implicit none + integer cctk_equals + CCTK_POINTER arg1 + character*(*) arg2 + end function cctk_equals + + subroutine cctk_fortranstring (clength, cstring, fortstring) + implicit none + CCTK_INT clength ! intent(out) + CCTK_POINTER cstring + character*(*) fortstring + end subroutine cctk_fortranstring + + subroutine cctk_groupbboxgn (ierr, cctkgh, nbbox, bbox, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nbbox + integer bbox(nbbox) + character*(*) groupname + end subroutine cctk_groupbboxgn + + subroutine cctk_groupgshgn (ierr, cctkgh, ngsh, gsh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer ngsh + integer gsh(ngsh) + character*(*) groupname + end subroutine cctk_groupgshgn + + subroutine cctk_grouplbndgn (ierr, cctkgh, nlbnd, lbnd, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nlbnd + integer lbnd(nlbnd) + character*(*) groupname + end subroutine cctk_grouplbndgn + + subroutine cctk_grouplshgn (ierr, cctkgh, nlsh, lsh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nlsh + integer lsh(nlsh) + character*(*) groupname + end subroutine cctk_grouplshgn + + subroutine cctk_groupnghostzonesgn & + (ierr, cctkgh, nnghostzones, nghostzones, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nnghostzones + integer nghostzones(nnghostzones) + character*(*) groupname + end subroutine cctk_groupnghostzonesgn + + subroutine cctk_groupubndgn (ierr, cctkgh, nubnd, ubnd, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nubnd + integer ubnd(nubnd) + character*(*) groupname + end subroutine cctk_groupubndgn + + subroutine cctk_info (thorn, message) + implicit none + character*(*) thorn + character*(*) message + end subroutine cctk_info + + subroutine cctk_interphandle (handle, interp) + implicit none + integer handle + character*(*) interp + end subroutine cctk_interphandle + + function cctk_isthornactive (name) + implicit none + integer cctk_isthornactive + character*(*) name + end function cctk_isthornactive + + function cctk_myproc (cctkgh) + implicit none + integer cctk_myproc + CCTK_POINTER cctkgh + end function cctk_myproc + + function cctk_nprocs (cctkgh) + implicit none + integer cctk_nprocs + CCTK_POINTER cctkgh + end function cctk_nprocs + + subroutine cctk_paramwarn (thorn, message) + implicit none + character*(*) thorn + character*(*) message + end subroutine cctk_paramwarn + + subroutine cctk_reductionarrayhandle (handle, reduction) + implicit none + integer handle + character*(*) reduction + end subroutine cctk_reductionarrayhandle + + subroutine cctk_registerbanner (ierr, banner) + implicit none + integer ierr + character*(*) banner + end subroutine cctk_registerbanner + + subroutine cctk_syncgroup (ierr, cctkgh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) groupname + end subroutine cctk_syncgroup + + subroutine cctk_syncgroupi (ierr, cctkgh, group) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer group + end subroutine cctk_syncgroupi + + subroutine cctk_syncgroupwithvar (ierr, cctkgh, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) varname + end subroutine cctk_syncgroupwithvar + + subroutine cctk_syncgroupwithvari (ierr, cctkgh, var) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer var + end subroutine cctk_syncgroupwithvari + + subroutine cctk_varindex (vi, vname) + implicit none + integer vi + character*(*) vname + end subroutine cctk_varindex + + subroutine cctk_vartypei (vartype, varindex) + implicit none + integer vartype + integer varindex + end subroutine cctk_vartypei + + subroutine cctk_warn (level, line, file, thorn, message) + implicit none + integer level + integer line + character*(*) file + character*(*) thorn + character*(*) message + end subroutine cctk_warn + + + + ! from thorn CactusBase/Boundary: + subroutine bndcopygn (ierr, cctkgh, sw, togroupname, fromgroupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) togroupname + character*(*) fromgroupname + end subroutine bndcopygn + + subroutine bndflatgn (ierr, cctkgh, sw, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) groupname + end subroutine bndflatgn + + subroutine bndradiativegn (ierr, cctkgh, sw, var0, v0, & + & togroupname, fromgroupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + CCTK_REAL v0 + character*(*) togroupname + character*(*) fromgroupname + end subroutine bndradiativegn + + subroutine bndradiativevn (ierr, cctkgh, sw, var0, v0, & + & tovarname, fromvarname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + CCTK_REAL v0 + character*(*) tovarname + character*(*) fromvarname + end subroutine bndradiativevn + + subroutine bndrobinvn (ierr, cctkgh, sw, finf, npow, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL finf + integer npow + character*(*) varname + end subroutine bndrobinvn + + subroutine bndscalargn (ierr, cctkgh, sw, var0, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + character*(*) groupname + end subroutine bndscalargn + + subroutine bndscalarvn (ierr, cctkgh, sw, var0, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + character*(*) varname + end subroutine bndscalarvn + + + + ! from thorn CactusBase/CartGrid3D: + + subroutine cartsymgn (ierr, cctkgh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) groupname + end subroutine cartsymgn + + subroutine setcartsymvn (ierr, cctkgh, sw, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) varname + end subroutine setcartsymvn + + + + ! from thorn AlphaThorns/Cart3D: + + subroutine cart3dsymvi (ierr, cctkgh, vi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer vi + end subroutine cart3dsymvi + + subroutine cart3dsymvn (ierr, cctkgh, vn) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) vn + end subroutine cart3dsymvn + + subroutine cart3dsymgi (ierr, cctkgh, gi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer gi + end subroutine cart3dsymgi + + subroutine cart3dsymgn (ierr, cctkgh, gn) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) gn + end subroutine cart3dsymgn + + subroutine cart3dsettensortypevi (ierr, cctkgh, nvars, vi, typename) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nvars + integer vi(nvars) + character*(*) typename + end subroutine cart3dsettensortypevi + + subroutine cart3dsettensortypevn (ierr, cctkgh, vn, typename) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) vn + character*(*) typename + end subroutine cart3dsettensortypevn + + subroutine cart3dgetsymmetriesvi (ierr, cctkgh, sym, vi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sym(6) + integer vi + end subroutine cart3dgetsymmetriesvi + + subroutine cart3dgetsymmetriesvn (ierr, cctkgh, sym, vn) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sym(6) + character*(*) vn + end subroutine cart3dgetsymmetriesvn + + subroutine cart3dgetsymmetryboundaries (ierr, cctkgh, symbnd) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer symbnd(6) + end subroutine cart3dgetsymmetryboundaries + + subroutine cart3dgetsymmetric (ierr, cctkgh, symmetric) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer symmetric(6) + end subroutine cart3dgetsymmetric + + subroutine cart3dgetstaggered (ierr, cctkgh, staggered) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer staggered(6) + end subroutine cart3dgetstaggered + + subroutine cart3dgetperiodic (ierr, cctkgh, periodic) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer periodic(6) + end subroutine cart3dgetperiodic + +end interface + +end module cactus diff --git a/src/constants.F90 b/src/constants.F90 new file mode 100644 index 0000000..d9b948d --- /dev/null +++ b/src/constants.F90 @@ -0,0 +1,52 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public delta2, delta3, delta4 + public eta4 + public epsilon2, epsilon3 + + public pi + + + + CCTK_REAL, parameter :: zero = 0 + + integer, parameter :: rk = kind(zero) + + + CCTK_REAL, parameter :: delta2(2,2) & + = reshape((/ 1,0, 0,1 /), (/2,2/)) + + CCTK_REAL, parameter :: delta3(3,3) & + = reshape((/ 1,0,0, 0,1,0, 0,0,1 /), (/3,3/)) + + CCTK_REAL, parameter :: delta4(0:3,0:3) & + = reshape((/ 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 /), (/4,4/)) + + + + CCTK_REAL, parameter :: eta4(0:3,0:3) & + = reshape((/ -1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 /), (/4,4/)) + + + + CCTK_REAL, parameter :: epsilon2(2,2) & + = reshape ((/ 0,+1, -1,0 /), (/2,2/)) + + CCTK_REAL, parameter :: epsilon3(3,3,3) & + = reshape ((/ 0,0,0, 0,0,+1, 0,-1,0, & + & 0,0,-1, 0,0,0, +1,0,0, & + & 0,+1,0, -1,0,0, 0,0,0 /), (/3,3,3/)) + + + + CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk + +end module constants diff --git a/src/conversion.F90 b/src/conversion.F90 new file mode 100644 index 0000000..d6e30e9 --- /dev/null +++ b/src/conversion.F90 @@ -0,0 +1,185 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module conversion + implicit none + DECLARE_CCTK_PARAMETERS + private + + public make_cylindrical2cartesian + public make_spherical2cartesian + public make_cartesian2spherical + +contains + + subroutine make_cylindrical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho + + ! cylindrical coordinates: rho phi z + ! Cartesian coordinates: x y z + ! rho^2 = x^2 + y^2 + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[cyl](p)_j + ! where p = p[cart] + + ! Transformation: + ! x = rho cos(phi) + ! y = rho sin(phi) + ! z = z + ! rho = sqrt(x^2+y^2) + ! phi = atan(y/x) + ! z = z + + ! Jacobian: + ! drho/dx = x/rho + ! drho/dy = y/rho + ! drho/dz = 0 + ! dphi/dx = -y/rho^2 + ! dphi/dy = x/rho^2 + ! dphi/dz = 0 + ! dz /dx = 0 + ! dz /dy = 0 + ! dz /dz = 1 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + + tt(1,1) = x/(rho+eps) + tt(2,1) = y/(rho+eps) + tt(3,1) = 0 + tt(1,2) = -y/(rho**2+eps) + tt(2,2) = x/(rho**2+eps) + tt(3,2) = 0 + tt(1,3) = 0 + tt(2,3) = 0 + tt(3,3) = 1 + + end subroutine make_cylindrical2cartesian + + subroutine make_spherical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[spher](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dr /dx = x/r + ! dr /dy = y/r + ! dr /dz = z/r + ! dtheta/dx = xz/r^2rho ! 0 + ! dtheta/dy = yz/r^2rho ! 0 + ! dtheta/dz = -rho/r^2 ! 1/rho + ! dphi /dx = -y/rho^2 + ! dphi /dy = x/rho^2 + ! dphi /dz = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = y/(r+eps) + tt(3,1) = z/(r+eps) + tt(1,2) = x*z/(r**2*rho+eps) + tt(2,2) = y*z/(r**2*rho+eps) + tt(3,2) = -rho**2/(r**2*rho+eps) + tt(1,3) = -y/(rho**2+eps) + tt(2,3) = x/(rho**2+eps) + tt(3,3) = 0 + + end subroutine make_spherical2cartesian + + subroutine make_cartesian2spherical (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[spher](p)_i = TT(p)_i^j A[cart](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dx/dr = sin(theta) cos(phi) = x/r + ! dx/dtheta = r cos(theta) cos(phi) = x z/rho + ! dx/dphi = - r sin(theta) sin(phi) = -y + ! dy/dr = sin(theta) sin(phi) = y/r + ! dy/dtheta = r cos(theta) sin(phi) = y z/rho + ! dy/dphi = r sin(theta) cos(phi) = x + ! dz/dr = cos(theta) = z/r + ! dz/dtheta = - r sin(theta) = -rho + ! dz/dphi = 0 = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = x * z/(rho+eps) + tt(3,1) = -y + tt(1,2) = y/(r+eps) + tt(2,2) = y * z/(rho+eps) + tt(3,2) = x + tt(1,3) = z/(r+eps) + tt(2,3) = -rho + tt(3,3) = 0 + + end subroutine make_cartesian2spherical + +end module conversion diff --git a/src/covderivs.F90 b/src/covderivs.F90 new file mode 100644 index 0000000..50e3624 --- /dev/null +++ b/src/covderivs.F90 @@ -0,0 +1,201 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module covderivs + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_scalargrad + public calc_covectorgrad + public calc_vectorgrad + public calc_tensorgrad + public calc_scalargradgrad + public calc_covectorgradgrad + public calc_longitudinal + public calc_gradlongitudinal + +contains + + subroutine calc_scalargrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3) + integer :: i + ! f;i = f,i + do i=1,3 + gf(i) = df(i) + end do + end subroutine calc_scalargrad + + subroutine calc_covectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3) + integer :: i,j,k + ! f_i;j = f_i,j - Gamma^k_ij f_k + do i=1,3 + do j=1,3 + gf(i,j) = df(i,j) + do k=1,3 + gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k) + end do + end do + end do + end subroutine calc_covectorgrad + + subroutine calc_vectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3) + integer :: i,j,k + ! f^i;j = f^i,j + Gamma^i_kj f^k + do i=1,3 + do j=1,3 + gf(i,j) = df(i,j) + do k=1,3 + gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k) + end do + end do + end do + end subroutine calc_vectorgrad + + subroutine calc_tensorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3,3) + CCTK_REAL, intent(in) :: df(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3,3) + integer :: i,j,k,l + ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il + do i=1,3 + do j=1,3 + do k=1,3 + gf(i,j,k) = df(i,j,k) + do l=1,3 + gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) & + & - gamma(l,j,k) * f(i,l) + end do + end do + end do + end do + end subroutine calc_tensorgrad + + + + subroutine calc_scalargradgrad (f, df, ddf, gamma, ggf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(3) + CCTK_REAL, intent(in) :: ddf(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3) + integer :: i,j,k + ! f;ij = f,ij - Gamma^k_ij f,k + do i=1,3 + do j=1,3 + ggf(i,j) = ddf(i,j) + do k=1,3 + ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k) + end do + end do + end do + end subroutine calc_scalargradgrad + + subroutine calc_covectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: ddf(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3,3) + CCTK_REAL :: gf(3,3) + integer :: i,j,k,l + ! f_i;j = f_i,j - Gamma^l_ij f_l + ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k + ! - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l + call calc_covectorgrad (f, df, gamma, gf) + do i=1,3 + do j=1,3 + do k=1,3 + ggf(i,j,k) = ddf(i,j,k) + do l=1,3 + ggf(i,j,k) = ggf(i,j,k) & + - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) & + - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_covectorgradgrad + + subroutine calc_vectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: ddf(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3,3) + CCTK_REAL :: gf(3,3) + integer :: i,j,k,l + ! f^i;j = f^i,j + Gamma^i_lj f^l + ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k + ! - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l + call calc_vectorgrad (f, df, gamma, gf) + do i=1,3 + do j=1,3 + do k=1,3 + ggf(i,j,k) = ddf(i,j,k) + do l=1,3 + ggf(i,j,k) = ggf(i,j,k) & + + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) & + - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_vectorgradgrad + + + + subroutine calc_longitudinal (gf, gg, gu, lf) + CCTK_REAL, intent(in) :: gf(3,3) + CCTK_REAL, intent(in) :: gg(3,3), gu(3,3) + CCTK_REAL, intent(out) :: lf(3,3) + integer :: i,j,k,l + ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l + do i=1,3 + do j=1,3 + lf(i,j) = gf(i,j) + gf(j,i) + do k=1,3 + do l=1,3 + lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l) + end do + end do + end do + end do + end subroutine calc_longitudinal + + subroutine calc_gradlongitudinal (ggf, gg, gu, glf) + CCTK_REAL, intent(in) :: ggf(3,3,3) + CCTK_REAL, intent(in) :: gg(3,3), gu(3,3) + CCTK_REAL, intent(out) :: glf(3,3,3) + integer :: i,j,k,l,m + ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l) + do i=1,3 + do j=1,3 + do k=1,3 + glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k) + do l=1,3 + do m=1,3 + glf(i,j,k) = glf(i,j,k) & + - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k) + end do + end do + end do + end do + end do + end subroutine calc_gradlongitudinal + +end module covderivs diff --git a/src/covderivs2.F90 b/src/covderivs2.F90 new file mode 100644 index 0000000..15e0143 --- /dev/null +++ b/src/covderivs2.F90 @@ -0,0 +1,204 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module covderivs2 + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_2scalargrad + public calc_2covectorgrad + public calc_2vectorgrad + public calc_2tensorgrad + + public calc_2scalargradgrad + public calc_2covectorgradgrad + public calc_2vectorgradgrad + + public calc_2longitudinal + public calc_2gradlongitudinal + +contains + + subroutine calc_2scalargrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2) + integer :: i + ! f;i = f,i + do i=1,2 + gf(i) = df(i) + end do + end subroutine calc_2scalargrad + + subroutine calc_2covectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2) + integer :: i,j,k + ! f_i;j = f_i,j - Gamma^k_ij f_k + do i=1,2 + do j=1,2 + gf(i,j) = df(i,j) + do k=1,2 + gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k) + end do + end do + end do + end subroutine calc_2covectorgrad + + subroutine calc_2vectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2) + integer :: i,j,k + ! f^i;j = f^i,j + Gamma^i_kj f^k + do i=1,2 + do j=1,2 + gf(i,j) = df(i,j) + do k=1,2 + gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k) + end do + end do + end do + end subroutine calc_2vectorgrad + + subroutine calc_2tensorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2,2) + CCTK_REAL, intent(in) :: df(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2,2) + integer :: i,j,k,l + ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il + do i=1,2 + do j=1,2 + do k=1,2 + gf(i,j,k) = df(i,j,k) + do l=1,2 + gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) & + & - gamma(l,j,k) * f(i,l) + end do + end do + end do + end do + end subroutine calc_2tensorgrad + + + + subroutine calc_2scalargradgrad (f, df, ddf, gamma, ggf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(2) + CCTK_REAL, intent(in) :: ddf(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2) + integer :: i,j,k + ! f;ij = f,ij - Gamma^k_ij f,k + do i=1,2 + do j=1,2 + ggf(i,j) = ddf(i,j) + do k=1,2 + ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k) + end do + end do + end do + end subroutine calc_2scalargradgrad + + subroutine calc_2covectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: ddf(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2), dgamma(2,2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2,2) + CCTK_REAL :: gf(2,2) + integer :: i,j,k,l + ! f_i;j = f_i,j - Gamma^l_ij f_l + ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k + ! - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l + call calc_2covectorgrad (f, df, gamma, gf) + do i=1,2 + do j=1,2 + do k=1,2 + ggf(i,j,k) = ddf(i,j,k) + do l=1,2 + ggf(i,j,k) = ggf(i,j,k) & + - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) & + - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_2covectorgradgrad + + subroutine calc_2vectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: ddf(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2), dgamma(2,2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2,2) + CCTK_REAL :: gf(2,2) + integer :: i,j,k,l + ! f^i;j = f^i,j + Gamma^i_lj f^l + ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k + ! - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l + call calc_2vectorgrad (f, df, gamma, gf) + do i=1,2 + do j=1,2 + do k=1,2 + ggf(i,j,k) = ddf(i,j,k) + do l=1,2 + ggf(i,j,k) = ggf(i,j,k) & + + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) & + - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_2vectorgradgrad + + + + subroutine calc_2longitudinal (gf, gg, gu, lf) + CCTK_REAL, intent(in) :: gf(2,2) + CCTK_REAL, intent(in) :: gg(2,2), gu(2,2) + CCTK_REAL, intent(out) :: lf(2,2) + integer :: i,j,k,l + ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l + do i=1,2 + do j=1,2 + lf(i,j) = gf(i,j) + gf(j,i) + do k=1,2 + do l=1,2 + lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l) + end do + end do + end do + end do + end subroutine calc_2longitudinal + + subroutine calc_2gradlongitudinal (ggf, gg, gu, glf) + CCTK_REAL, intent(in) :: ggf(2,2,2) + CCTK_REAL, intent(in) :: gg(2,2), gu(2,2) + CCTK_REAL, intent(out) :: glf(2,2,2) + integer :: i,j,k,l,m + ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l) + do i=1,2 + do j=1,2 + do k=1,2 + glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k) + do l=1,2 + do m=1,2 + glf(i,j,k) = glf(i,j,k) & + - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k) + end do + end do + end do + end do + end do + end subroutine calc_2gradlongitudinal + +end module covderivs2 diff --git a/src/depend.pl b/src/depend.pl new file mode 100644 index 0000000..b41b8f6 --- /dev/null +++ b/src/depend.pl @@ -0,0 +1,46 @@ +#!/usr/bin/perl -w +# $Header$ + +# Create dependencies for Fortran 90 "use" and "include" statements + +$src = $ARGV[0]; +$srcfile = $ARGV[1]; +$dest = $ARGV[2]; +$srcdir = $ARGV[3]; +$moddir = $ARGV[4]; +$srcsuffix = $ARGV[5]; +$modsuffix = $ARGV[6]; +@otherdirs = @ARGV[7..$#ARGV]; + +print "# $src\n"; +print "$dest:"; + +open (FILE, $src); +while (<FILE>) { + if (/^\s*include\s*['"]([^'"]+)['"]/i) { + print " $srcdir$1"; + } elsif (/^\s*use\s+(\w+)/i) { + $found = 0; + if (! $found) { + if (-e "$srcdir$1$srcsuffix") { + $found = 1; + print " $moddir$1$modsuffix"; + } + } + if (! $found) { + foreach $dir (@otherdirs) { + if (-e "$dir$1$modsuffix") { + $found = 1; + print " $dir$1$modsuffix"; + last; + } + } + } + if (! $found) { + die "\nWhile tracing depencencies:\nFortran module $1 (referenced in file $srcfile) not found.\nAborting.\n\n"; + } + } +} +close (FILE); + +print "\n"; diff --git a/src/derivs.F90 b/src/derivs.F90 new file mode 100644 index 0000000..1d425c8 --- /dev/null +++ b/src/derivs.F90 @@ -0,0 +1,46 @@ +! $Header$ + +#ifndef TGR_INCLUDED + +#include "cctk.h" +#include "cctk_Parameters.h" + +module derivs + implicit none + private + public get_derivs + public get_derivs2 +contains +#endif + subroutine get_derivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + integer :: i + do i=1,3 + f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i)) + end do + end subroutine get_derivs + subroutine get_derivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + integer :: i + do i=1,3 + f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2 + end do + f(1,2) = ( a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) & + & - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2)) + f(2,1) = f(1,2) + f(1,3) = ( a(pos-off(1)-off(3)) - a(pos+off(1)-off(3)) & + & - a(pos-off(1)+off(3)) + a(pos+off(1)+off(3))) / (4*dx(1)*dx(3)) + f(3,1) = f(1,3) + f(2,3) = ( a(pos-off(2)-off(3)) - a(pos+off(2)-off(3)) & + & - a(pos-off(2)+off(3)) + a(pos+off(2)+off(3))) / (4*dx(2)*dx(3)) + f(3,2) = f(2,3) + end subroutine get_derivs2 +#ifndef TGR_INCLUDED +end module derivs +#endif diff --git a/src/derivs2.F90 b/src/derivs2.F90 new file mode 100644 index 0000000..4ca51bf --- /dev/null +++ b/src/derivs2.F90 @@ -0,0 +1,40 @@ +! $Header$ + +#ifndef TGR_INCLUDED + +#include "cctk.h" +#include "cctk_Parameters.h" + +module derivs2 + implicit none + private + public get_2derivs + public get_2derivs2 +contains +#endif + subroutine get_2derivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + integer :: i + do i=1,2 + f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i)) + end do + end subroutine get_2derivs + subroutine get_2derivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + integer :: i + do i=1,2 + f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2 + end do + f(1,2) = ( a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) & + & - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2)) + f(2,1) = f(1,2) + end subroutine get_2derivs2 +#ifndef TGR_INCLUDED +end module derivs2 +#endif diff --git a/src/hyperslab.F90 b/src/hyperslab.F90 new file mode 100644 index 0000000..e540cb2 --- /dev/null +++ b/src/hyperslab.F90 @@ -0,0 +1,33 @@ +! $Header$ + +#include "cctk.h" + +module hyperslab + implicit none + public + + interface + + subroutine hyperslab_fillhyperslab (ierr, cctkgh, target_proc, & + vindex, vtimelevel, & + hdim, global_startpoint, directions, lengths, downsample, & + nhdata, hdata, hsize) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer target_proc + integer vindex + integer vtimelevel + integer hdim + integer global_startpoint(*) ! vdim + integer directions(*) ! vdim + integer lengths(hdim) + integer downsample(hdim) + integer nhdata + CCTK_REAL hdata(*) + integer hsize(hdim) + end subroutine hyperslab_fillhyperslab + + end interface + +end module hyperslab diff --git a/src/lapack.F90 b/src/lapack.F90 new file mode 100644 index 0000000..5bcd70b --- /dev/null +++ b/src/lapack.F90 @@ -0,0 +1,51 @@ +! $Header$ + +#include "cctk.h" + +module lapack + implicit none + public + + interface geev + subroutine sgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, & + & ldvr, work, lwork, info) + character jobvl, jobvr + integer info, lda, ldvl, ldvr, lwork, n + real a(lda,n), vl(ldvl,n), vr(ldvr,n), & + & wi(n), work(lwork), wr(n) + end subroutine sgeev + subroutine dgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, & + & ldvr, work, lwork, info) + character jobvl, jobvr + integer info, lda, ldvl, ldvr, lwork, n + double precision a(lda,n), vl(ldvl,n), vr(ldvr,n), & + & wi(n), work(lwork), wr(n) + end subroutine dgeev + end interface + + interface gesv + subroutine sgesv (n, nrhs, a, lda, ipiv, b, ldb, info) + implicit none + integer n + integer nrhs + integer lda + real a(lda,n) + integer ipiv(n) + integer ldb + real b(ldb,nrhs) + integer info + end subroutine sgesv + subroutine dgesv (n, nrhs, a, lda, ipiv, b, ldb, info) + implicit none + integer n + integer nrhs + integer lda + double precision a(lda,n) + integer ipiv(n) + integer ldb + double precision b(ldb,nrhs) + integer info + end subroutine dgesv + end interface + +end module lapack diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..4b157e4 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,25 @@ +# Main make.code.defn file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +# Source files in this directory +SRCS = cactus.F90 \ + constants.F90 \ + conversion.F90 \ + covderivs.F90 \ + covderivs2.F90 \ + derivs.F90 \ + derivs2.F90 \ + hyperslab.F90 \ + hyperslab.c \ + lapack.F90 \ + matexp.F90 \ + matinv.F90 \ + rotation.F90 \ + tensor.F90 \ + tensor2.F90 \ + pointwise.F90 \ + pointwise2.F90 + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/make.code.deps b/src/make.code.deps new file mode 100644 index 0000000..e8ea1b9 --- /dev/null +++ b/src/make.code.deps @@ -0,0 +1,11 @@ +# Main make.code.deps file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +USESTHORNS = + +# Automatically create dependencies for Fortran modules and Fortran includes +define F90_DEPENDENCIES +$(F_DEPEND) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< $(F_DEPEND_OUT) +$(DEPENDENCY_FIXER) +dir=`pwd`; $(CPP) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< | $(PERL) $(SRCDIR)/depend.pl - $< $(basename $(notdir $<)).F90.o $(SRCDIR)/ ./ .F90 .F90.o $(USESTHORNS:%=$$dir/../%/) >> $@ || { rm $@; exit 1; } +endef diff --git a/src/make.configuration.deps b/src/make.configuration.deps new file mode 100644 index 0000000..6a85f9c --- /dev/null +++ b/src/make.configuration.deps @@ -0,0 +1,6 @@ +# Main make.configuration.deps file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +USESTHORNS = + +$(CCTK_LIBDIR)$(DIRSEP)libTGRtensor.a: $(USESTHORNS:%=$(CCTK_LIBDIR)$(DIRSEP)lib%.a) diff --git a/src/matexp.F90 b/src/matexp.F90 new file mode 100644 index 0000000..d2ecc77 --- /dev/null +++ b/src/matexp.F90 @@ -0,0 +1,38 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module matexp + use constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_exp3 + +contains + + subroutine calc_exp3 (h3, g3) + CCTK_REAL, intent(in) :: h3(3,3) + CCTK_REAL, intent(out) :: g3(3,3) + CCTK_REAL :: nfact + CCTK_REAL :: tmp(3,3) + integer :: n + integer :: i, j + + g3 = delta3 + + nfact = 1 + tmp = delta3 + do n = 1, 18 + nfact = nfact * n + tmp = matmul (h3, tmp) + + ! exp(x) = sum_n x^n / n! + g3 = g3 + tmp / nfact + end do + + end subroutine calc_exp3 + +end module matexp diff --git a/src/matinv.F90 b/src/matinv.F90 new file mode 100644 index 0000000..ab6c515 --- /dev/null +++ b/src/matinv.F90 @@ -0,0 +1,54 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module matinv + use constants + use lapack + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_inv3 + public calc_inv4 + +contains + + subroutine calc_inv3 (g3, gu3) + CCTK_REAL, intent(in) :: g3(3,3) + CCTK_REAL, intent(out) :: gu3(3,3) + CCTK_REAL :: tmp(3,3) + integer :: ipiv(3) + integer :: info + character :: msg*1000 + + tmp = g3 + gu3 = delta3 + call gesv (3, 3, tmp, 3, ipiv, gu3, 3, info) + + if (info /= 0) then + write (msg, '("Error in call to GESV, info=",i2)') info + call CCTK_WARN (0, trim(msg)) + end if + end subroutine calc_inv3 + + subroutine calc_inv4 (g4, gu4) + CCTK_REAL, intent(in) :: g4(0:3,0:3) + CCTK_REAL, intent(out) :: gu4(0:3,0:3) + CCTK_REAL :: tmp(0:3,0:3) + integer :: ipiv(0:3) + integer :: info + character :: msg*1000 + + tmp = g4 + gu4 = delta4 + call gesv (4, 4, tmp, 4, ipiv, gu4, 4, info) + + if (info /= 0) then + write (msg, '("Error in call to GESV, info=",i2)') info + call CCTK_WARN (0, trim(msg)) + end if + end subroutine calc_inv4 + +end module matinv diff --git a/src/pointwise.F90 b/src/pointwise.F90 new file mode 100644 index 0000000..827ad6e --- /dev/null +++ b/src/pointwise.F90 @@ -0,0 +1,252 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module pointwise +!!$ use derivs + implicit none + DECLARE_CCTK_PARAMETERS + private + public calc_position + public calc_offsets + public get_scalar + public get_vector + public get_tensor + public get_connections + public set_scalar + public set_vector + public set_tensor + public set_connections + public get_scalarderivs + public get_vectorderivs + public get_tensorderivs + public get_scalarderivs2 + public get_vectorderivs2 + public get_tensorderivs2 +contains +#define TGR_INCLUDED +#include "derivs.F90" +#undef TGR_INCLUDED + subroutine calc_position (shape, i,j,k, pos) + integer, intent(in) :: shape(3) + integer, intent(in) :: i,j,k + integer, intent(out) :: pos + pos = 1 + i-1 + shape(1) * (j-1 + shape(2) * (k-1)) + end subroutine calc_position + subroutine calc_offsets (shape, off) + integer, intent(in) :: shape(3) + integer, intent(out) :: off(3) + off(1) = 1 + off(2) = shape(1) + off(3) = shape(1) * shape(2) + end subroutine calc_offsets + subroutine get_scalar (a, f, pos) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f + integer, intent(in) :: pos + f = a(pos) + end subroutine get_scalar + subroutine get_vector (ax,ay,az, f, pos) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos + f(1) = ax(pos) + f(2) = ay(pos) + f(3) = az(pos) + end subroutine get_vector + subroutine get_tensor (axx,axy,axz,ayy,ayz,azz, f, pos) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos + f(1,1) = axx(pos) + f(1,2) = axy(pos) + f(1,3) = axz(pos) + f(2,1) = axy(pos) + f(2,2) = ayy(pos) + f(2,3) = ayz(pos) + f(3,1) = axz(pos) + f(3,2) = ayz(pos) + f(3,3) = azz(pos) + end subroutine get_tensor + subroutine get_connections (& + gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, & + gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, & + gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, gamma, pos) + CCTK_REAL, intent(in) :: gammaxxx(*),gammaxxy(*),gammaxxz(*) + CCTK_REAL, intent(in) :: gammaxyy(*),gammaxyz(*),gammaxzz(*) + CCTK_REAL, intent(in) :: gammayxx(*),gammayxy(*),gammayxz(*) + CCTK_REAL, intent(in) :: gammayyy(*),gammayyz(*),gammayzz(*) + CCTK_REAL, intent(in) :: gammazxx(*),gammazxy(*),gammazxz(*) + CCTK_REAL, intent(in) :: gammazyy(*),gammazyz(*),gammazzz(*) + CCTK_REAL, intent(out) :: gamma(3,3,3) + integer, intent(in) :: pos + gamma(1,1,1) = gammaxxx(pos) + gamma(1,1,2) = gammaxxy(pos) + gamma(1,1,3) = gammaxxz(pos) + gamma(1,2,1) = gammaxxy(pos) + gamma(1,2,2) = gammaxyy(pos) + gamma(1,2,3) = gammaxyz(pos) + gamma(1,3,1) = gammaxxz(pos) + gamma(1,3,2) = gammaxyz(pos) + gamma(1,3,3) = gammaxzz(pos) + gamma(2,1,1) = gammayxx(pos) + gamma(2,1,2) = gammayxy(pos) + gamma(2,1,3) = gammayxz(pos) + gamma(2,2,1) = gammayxy(pos) + gamma(2,2,2) = gammayyy(pos) + gamma(2,2,3) = gammayyz(pos) + gamma(2,3,1) = gammayxz(pos) + gamma(2,3,2) = gammayyz(pos) + gamma(2,3,3) = gammayzz(pos) + gamma(3,1,1) = gammazxx(pos) + gamma(3,1,2) = gammazxy(pos) + gamma(3,1,3) = gammazxz(pos) + gamma(3,2,1) = gammazxy(pos) + gamma(3,2,2) = gammazyy(pos) + gamma(3,2,3) = gammazyz(pos) + gamma(3,3,1) = gammazxz(pos) + gamma(3,3,2) = gammazyz(pos) + gamma(3,3,3) = gammazzz(pos) + end subroutine get_connections + subroutine set_scalar (f, a, pos) + CCTK_REAL, intent(in) :: f + CCTK_REAL :: a(*) + integer, intent(in) :: pos + a(pos) = f + end subroutine set_scalar + subroutine set_vector (f, ax,ay,az, pos) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL :: ax(*),ay(*),az(*) + integer, intent(in) :: pos + ax(pos) = f(1) + ay(pos) = f(2) + az(pos) = f(3) + end subroutine set_vector + subroutine set_tensor (f, axx,axy,axz,ayy,ayz,azz, pos) + CCTK_REAL, intent(in) :: f(3,3) + CCTK_REAL :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + integer, intent(in) :: pos + axx(pos) = f(1,1) + axy(pos) = f(1,2) + axz(pos) = f(1,3) + ayy(pos) = f(2,2) + ayz(pos) = f(2,3) + azz(pos) = f(3,3) + end subroutine set_tensor + subroutine set_connections (gamma, & + gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, & + gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, & + gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, pos) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL :: gammaxxx(*),gammaxxy(*),gammaxxz(*) + CCTK_REAL :: gammaxyy(*),gammaxyz(*),gammaxzz(*) + CCTK_REAL :: gammayxx(*),gammayxy(*),gammayxz(*) + CCTK_REAL :: gammayyy(*),gammayyz(*),gammayzz(*) + CCTK_REAL :: gammazxx(*),gammazxy(*),gammazxz(*) + CCTK_REAL :: gammazyy(*),gammazyz(*),gammazzz(*) + integer, intent(in) :: pos + gammaxxx(pos) = gamma(1,1,1) + gammaxxy(pos) = gamma(1,1,2) + gammaxxz(pos) = gamma(1,1,3) + gammaxyy(pos) = gamma(1,2,2) + gammaxyz(pos) = gamma(1,2,3) + gammaxzz(pos) = gamma(1,3,3) + gammayxx(pos) = gamma(2,1,1) + gammayxy(pos) = gamma(2,1,2) + gammayxz(pos) = gamma(2,1,3) + gammayyy(pos) = gamma(2,2,2) + gammayyz(pos) = gamma(2,2,3) + gammayzz(pos) = gamma(2,3,3) + gammazxx(pos) = gamma(3,1,1) + gammazxy(pos) = gamma(3,1,2) + gammazxz(pos) = gamma(3,1,3) + gammazyy(pos) = gamma(3,2,2) + gammazyz(pos) = gamma(3,2,3) + gammazzz(pos) = gamma(3,3,3) + end subroutine set_connections + subroutine get_scalarderivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + call get_derivs (a, f, pos, off, dx) + end subroutine get_scalarderivs + subroutine get_vectorderivs (ax,ay,az, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fx(3),fy(3),fz(3) + call get_derivs (ax, fx, pos, off, dx) + call get_derivs (ay, fy, pos, off, dx) + call get_derivs (az, fz, pos, off, dx) + f(1,:) = fx + f(2,:) = fy + f(3,:) = fz + end subroutine get_vectorderivs + subroutine get_tensorderivs (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fxx(3),fxy(3),fxz(3),fyy(3),fyz(3),fzz(3) + call get_derivs (axx, fxx, pos, off, dx) + call get_derivs (axy, fxy, pos, off, dx) + call get_derivs (axz, fxz, pos, off, dx) + call get_derivs (ayy, fyy, pos, off, dx) + call get_derivs (ayz, fyz, pos, off, dx) + call get_derivs (azz, fzz, pos, off, dx) + f(1,1,:) = fxx + f(1,2,:) = fxy + f(1,3,:) = fxz + f(2,1,:) = fxy + f(2,2,:) = fyy + f(2,3,:) = fyz + f(3,1,:) = fxz + f(3,2,:) = fyz + f(3,3,:) = fzz + end subroutine get_tensorderivs + subroutine get_scalarderivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + call get_derivs2 (a, f, pos, off, dx) + end subroutine get_scalarderivs2 + subroutine get_vectorderivs2 (ax,ay,az, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fx(3,3),fy(3,3),fz(3,3) + call get_derivs2 (ax, fx, pos, off, dx) + call get_derivs2 (ay, fy, pos, off, dx) + call get_derivs2 (az, fz, pos, off, dx) + f(1,:,:) = fx + f(2,:,:) = fy + f(3,:,:) = fz + end subroutine get_vectorderivs2 + subroutine get_tensorderivs2 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fxx(3,3),fxy(3,3),fxz(3,3),fyy(3,3),fyz(3,3),fzz(3,3) + call get_derivs2 (axx, fxx, pos, off, dx) + call get_derivs2 (axy, fxy, pos, off, dx) + call get_derivs2 (axz, fxz, pos, off, dx) + call get_derivs2 (ayy, fyy, pos, off, dx) + call get_derivs2 (ayz, fyz, pos, off, dx) + call get_derivs2 (azz, fzz, pos, off, dx) + f(1,1,:,:) = fxx + f(1,2,:,:) = fxy + f(1,3,:,:) = fxz + f(2,1,:,:) = fxy + f(2,2,:,:) = fyy + f(2,3,:,:) = fyz + f(3,1,:,:) = fxz + f(3,2,:,:) = fyz + f(3,3,:,:) = fzz + end subroutine get_tensorderivs2 +end module pointwise diff --git a/src/pointwise2.F90 b/src/pointwise2.F90 new file mode 100644 index 0000000..1759b28 --- /dev/null +++ b/src/pointwise2.F90 @@ -0,0 +1,148 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module pointwise2 +!!$ use derivs2 + implicit none + DECLARE_CCTK_PARAMETERS + private + public calc_2position + public calc_2offsets + public get_2scalar + public get_2vector + public get_2tensor + public set_2scalar + public set_2vector + public set_2tensor + public get_2scalarderivs + public get_2vectorderivs + public get_2tensorderivs + public get_2scalarderivs2 + public get_2vectorderivs2 + public get_2tensorderivs2 +contains +#define TGR_INCLUDED +#include "derivs2.F90" +#undef TGR_INCLUDED + subroutine calc_2position (shape, i,j, pos) + integer, intent(in) :: shape(2) + integer, intent(in) :: i,j + integer, intent(out) :: pos + pos = 1 + i-1 + shape(1) * (j-1) + end subroutine calc_2position + subroutine calc_2offsets (shape, off) + integer, intent(in) :: shape(2) + integer, intent(out) :: off(2) + off(1) = 1 + off(2) = shape(1) + end subroutine calc_2offsets + subroutine get_2scalar (a, f, pos) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f + integer, intent(in) :: pos + f = a(pos) + end subroutine get_2scalar + subroutine get_2vector (ax,ay, f, pos) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos + f(1) = ax(pos) + f(2) = ay(pos) + end subroutine get_2vector + subroutine get_2tensor (axx,axy,ayy, f, pos) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos + f(1,1) = axx(pos) + f(1,2) = axy(pos) + f(2,1) = axy(pos) + f(2,2) = ayy(pos) + end subroutine get_2tensor + subroutine set_2scalar (f, a, pos) + CCTK_REAL, intent(in) :: f + CCTK_REAL :: a(*) + integer, intent(in) :: pos + a(pos) = f + end subroutine set_2scalar + subroutine set_2vector (f, ax,ay, pos) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL :: ax(*),ay(*) + integer, intent(in) :: pos + ax(pos) = f(1) + ay(pos) = f(2) + end subroutine set_2vector + subroutine set_2tensor (f, axx,axy,ayy, pos) + CCTK_REAL, intent(in) :: f(2,2) + CCTK_REAL :: axx(*),axy(*),ayy(*) + integer, intent(in) :: pos + axx(pos) = f(1,1) + axy(pos) = f(1,2) + ayy(pos) = f(2,2) + end subroutine set_2tensor + subroutine get_2scalarderivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + call get_2derivs (a, f, pos, off, dx) + end subroutine get_2scalarderivs + subroutine get_2vectorderivs (ax,ay, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fx(2),fy(2) + call get_2derivs (ax, fx, pos, off, dx) + call get_2derivs (ay, fy, pos, off, dx) + f(1,:) = fx + f(2,:) = fy + end subroutine get_2vectorderivs + subroutine get_2tensorderivs (axx,axy,ayy, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fxx(2),fxy(2),fyy(2) + call get_2derivs (axx, fxx, pos, off, dx) + call get_2derivs (axy, fxy, pos, off, dx) + call get_2derivs (ayy, fyy, pos, off, dx) + f(1,1,:) = fxx + f(1,2,:) = fxy + f(2,1,:) = fxy + f(2,2,:) = fyy + end subroutine get_2tensorderivs + subroutine get_2scalarderivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + call get_2derivs2 (a, f, pos, off, dx) + end subroutine get_2scalarderivs2 + subroutine get_2vectorderivs2 (ax,ay, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fx(2,2),fy(2,2) + call get_2derivs2 (ax, fx, pos, off, dx) + call get_2derivs2 (ay, fy, pos, off, dx) + f(1,:,:) = fx + f(2,:,:) = fy + end subroutine get_2vectorderivs2 + subroutine get_2tensorderivs2 (axx,axy,ayy, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fxx(2,2),fxy(2,2),fyy(2,2) + call get_2derivs2 (axx, fxx, pos, off, dx) + call get_2derivs2 (axy, fxy, pos, off, dx) + call get_2derivs2 (ayy, fyy, pos, off, dx) + f(1,1,:,:) = fxx + f(1,2,:,:) = fxy + f(2,1,:,:) = fxy + f(2,2,:,:) = fyy + end subroutine get_2tensorderivs2 +end module pointwise2 diff --git a/src/rotation.F90 b/src/rotation.F90 new file mode 100644 index 0000000..d76973e --- /dev/null +++ b/src/rotation.F90 @@ -0,0 +1,59 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module rotation + use constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public make_euler + +contains + + subroutine make_euler1 (phi, tt) + CCTK_REAL, intent(in) :: phi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(1,1) = cos(phi) + tt(1,2) = -sin(phi) + tt(2,1) = sin(phi) + tt(2,2) = cos(phi) + end subroutine make_euler1 + + subroutine make_euler2 (theta, tt) + CCTK_REAL, intent(in) :: theta + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(2,2) = cos(theta) + tt(2,3) = -sin(theta) + tt(3,2) = sin(theta) + tt(3,3) = cos(theta) + end subroutine make_euler2 + + subroutine make_euler3 (psi, tt) + CCTK_REAL, intent(in) :: psi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(1,1) = cos(psi) + tt(1,2) = -sin(psi) + tt(2,1) = sin(psi) + tt(2,2) = cos(psi) + end subroutine make_euler3 + + subroutine make_euler (phi, theta, psi, tt) + CCTK_REAL, intent(in) :: phi, theta, psi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + CCTK_REAL :: tt1(0:3,0:3), tt2(0:3,0:3), tt3(0:3,0:3) + call make_euler1 (phi, tt1) + call make_euler2 (theta, tt2) + call make_euler3 (psi, tt3) + tt = matmul(tt3, matmul(tt2, tt1)) + end subroutine make_euler + +end module rotation diff --git a/src/tensor.F90 b/src/tensor.F90 new file mode 100644 index 0000000..0fef831 --- /dev/null +++ b/src/tensor.F90 @@ -0,0 +1,127 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module tensor + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_trace + + public calc_det + public calc_detderiv + public calc_detdot + + public calc_inv + public calc_invderiv + public calc_invdot + +contains + + subroutine calc_trace (kk, gu, trk) + CCTK_REAL, intent(in) :: kk(3,3) + CCTK_REAL, intent(in) :: gu(3,3) + CCTK_REAL, intent(out) :: trk + integer :: i,j + trk = 0 + do i=1,3 + do j=1,3 + trk = trk + gu(i,j) * kk(i,j) + end do + end do + end subroutine calc_trace + + + + subroutine calc_det (gg, dtg) + CCTK_REAL, intent(in) :: gg(3,3) + CCTK_REAL, intent(out) :: dtg + dtg = gg(1,1) * gg(2,2) * gg(3,3) & + & + 2*gg(1,2) * gg(1,3) * gg(2,3) & + & - gg(1,1) * gg(2,3)**2 & + & - gg(2,2) * gg(1,3)**2 & + & - gg(3,3) * gg(1,2)**2 + end subroutine calc_det + + subroutine calc_pdet (gg, pgg, pdtg) + CCTK_REAL, intent(in) :: gg(3,3), pgg(3,3) + CCTK_REAL, intent(out) :: pdtg + pdtg = pgg(1,1) * gg(2,2) * gg(3,3) & + & + gg(1,1) * pgg(2,2) * gg(3,3) & + & + gg(1,1) * gg(2,2) * pgg(3,3) & + & + 2*pgg(1,2) * gg(1,3) * gg(2,3) & + & + 2* gg(1,2) * pgg(1,3) * gg(2,3) & + & + 2* gg(1,2) * gg(1,3) * pgg(2,3) & + & - pgg(1,1) * gg(2,3)**2 & + & - gg(1,1) * 2*gg(2,3) * pgg(2,3) & + & - pgg(2,2) * gg(1,3)**2 & + & - gg(2,2) * 2*gg(1,3) * pgg(1,3) & + & - pgg(3,3) * gg(1,2)**2 & + & - gg(3,3) * 2*gg(1,2) * pgg(1,2) + end subroutine calc_pdet + + subroutine calc_detderiv (gg, dgg, ddtg) + CCTK_REAL, intent(in) :: gg(3,3), dgg(3,3,3) + CCTK_REAL, intent(out) :: ddtg(3) + integer :: i + do i=1,3 + call calc_pdet (gg, dgg(:,:,i), ddtg(i)) + end do + end subroutine calc_detderiv + + subroutine calc_detdot (gg, gg_dot, dtg_dot) + CCTK_REAL, intent(in) :: gg(3,3), gg_dot(3,3) + CCTK_REAL, intent(out) :: dtg_dot + call calc_pdet (gg, gg_dot, dtg_dot) + end subroutine calc_detdot + + + + subroutine calc_inv (gg, dtg, gu) + CCTK_REAL, intent(in) :: gg(3,3), dtg + CCTK_REAL, intent(out) :: gu(3,3) + gu(1,1) = (gg(2,2) * gg(3,3) - gg(2,3) ** 2 ) / dtg + gu(2,2) = (gg(1,1) * gg(3,3) - gg(1,3) ** 2 ) / dtg + gu(3,3) = (gg(1,1) * gg(2,2) - gg(1,2) ** 2 ) / dtg + gu(1,2) = (gg(1,3) * gg(2,3) - gg(1,2) * gg(3,3)) / dtg + gu(1,3) = (gg(1,2) * gg(2,3) - gg(1,3) * gg(2,2)) / dtg + gu(2,3) = (gg(1,3) * gg(1,2) - gg(2,3) * gg(1,1)) / dtg + gu(2,1) = gu(1,2) + gu(3,1) = gu(1,3) + gu(3,2) = gu(2,3) + end subroutine calc_inv + + subroutine calc_pinv (gu, pgg, pgu) + CCTK_REAL, intent(in) :: gu(3,3), pgg(3,3) + CCTK_REAL, intent(out) :: pgu(3,3) + integer :: i,j,k,l + do i=1,3 + do j=1,3 + pgu(i,j) = 0 + do k=1,3 + do l=1,3 + pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l) + end do + end do + end do + end do + end subroutine calc_pinv + + subroutine calc_invderiv (gu, dgg, dgu) + CCTK_REAL, intent(in) :: gu(3,3), dgg(3,3,3) + CCTK_REAL, intent(out) :: dgu(3,3,3) + integer :: i + do i=1,3 + call calc_pinv (gu, dgg(:,:,i), dgu(:,:,i)) + end do + end subroutine calc_invderiv + + subroutine calc_invdot (gu, gg_dot, gu_dot) + CCTK_REAL, intent(in) :: gu(3,3), gg_dot(3,3) + CCTK_REAL, intent(out) :: gu_dot(3,3) + call calc_pinv (gu, gg_dot, gu_dot) + end subroutine calc_invdot + +end module tensor diff --git a/src/tensor2.F90 b/src/tensor2.F90 new file mode 100644 index 0000000..b02cd06 --- /dev/null +++ b/src/tensor2.F90 @@ -0,0 +1,107 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module tensor2 + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_2trace + + public calc_2det + public calc_2detderiv + public calc_2detdot + + public calc_2inv + public calc_2invderiv + public calc_2invdot + +contains + + subroutine calc_2trace (kk, gu, trk) + CCTK_REAL, intent(in) :: kk(2,2) + CCTK_REAL, intent(in) :: gu(2,2) + CCTK_REAL, intent(out) :: trk + integer :: i,j + trk = 0 + do i=1,2 + do j=1,2 + trk = trk + gu(i,j) * kk(i,j) + end do + end do + end subroutine calc_2trace + + + + subroutine calc_2det (gg, dtg) + CCTK_REAL, intent(in) :: gg(2,2) + CCTK_REAL, intent(out) :: dtg + dtg = gg(1,1) * gg(2,2) - gg(1,2)**2 + end subroutine calc_2det + + subroutine calc_2pdet (gg, pgg, pdtg) + CCTK_REAL, intent(in) :: gg(2,2), pgg(2,2) + CCTK_REAL, intent(out) :: pdtg + pdtg = pgg(1,1) * gg(2,2) + gg(1,1) * pgg(2,2) - 2 * gg(1,2) * pgg(1,2) + end subroutine calc_2pdet + + subroutine calc_2detderiv (gg, dgg, ddtg) + CCTK_REAL, intent(in) :: gg(2,2), dgg(2,2,2) + CCTK_REAL, intent(out) :: ddtg(2) + integer :: i + do i=1,2 + call calc_2pdet (gg, dgg(:,:,i), ddtg(i)) + end do + end subroutine calc_2detderiv + + subroutine calc_2detdot (gg, gg_dot, dtg_dot) + CCTK_REAL, intent(in) :: gg(2,2), gg_dot(2,2) + CCTK_REAL, intent(out) :: dtg_dot + call calc_2pdet (gg, gg_dot, dtg_dot) + end subroutine calc_2detdot + + + + subroutine calc_2inv (gg, dtg, gu) + CCTK_REAL, intent(in) :: gg(2,2), dtg + CCTK_REAL, intent(out) :: gu(2,2) + gu(1,1) = gg(2,2) / dtg + gu(2,2) = gg(1,1) / dtg + gu(1,2) = gg(1,2) / dtg + gu(2,1) = gu(1,2) + end subroutine calc_2inv + + subroutine calc_2pinv (gu, pgg, pgu) + CCTK_REAL, intent(in) :: gu(2,2), pgg(2,2) + CCTK_REAL, intent(out) :: pgu(2,2) + integer :: i,j,k,l + do i=1,2 + do j=1,2 + pgu(i,j) = 0 + do k=1,2 + do l=1,2 + pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l) + end do + end do + end do + end do + end subroutine calc_2pinv + + subroutine calc_2invderiv (gu, dgg, dgu) + CCTK_REAL, intent(in) :: gu(2,2), dgg(2,2,2) + CCTK_REAL, intent(out) :: dgu(2,2,2) + integer :: i + do i=1,2 + call calc_2pinv (gu, dgg(:,:,i), dgu(:,:,i)) + end do + end subroutine calc_2invderiv + + subroutine calc_2invdot (gu, gg_dot, gu_dot) + CCTK_REAL, intent(in) :: gu(2,2), gg_dot(2,2) + CCTK_REAL, intent(out) :: gu_dot(2,2) + call calc_2pinv (gu, gg_dot, gu_dot) + end subroutine calc_2invdot + +end module tensor2 |