From 9f9d618e54f91e0eb54c9de197b597ff99ae4749 Mon Sep 17 00:00:00 2001 From: schnetter Date: Sat, 9 Apr 2005 13:54:39 +0000 Subject: Add some functionality: Conversion between 4 and 3+1 dimensions Check for nans (isnan for Fortran) Gram-Schmidt orthonormalisation Calculate Riemann, Ricci, and Weyl tensors in 2, 3, and 4 dimensions Calculate time derivatives git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@21 b716e942-a2de-43ad-8f52-f3dfe468e4e7 --- src/adm_metric.F90 | 277 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/depend.pl | 46 --------- src/gram_schmidt.F90 | 121 ++++++++++++++++++++++ src/isnan.c | 32 ++++++ src/make.code.defn | 9 +- src/ricci.F90 | 85 ++++++++++++++++ src/ricci2.F90 | 74 ++++++++++++++ src/ricci4.F90 | 181 +++++++++++++++++++++++++++++++++ src/timederivs.F90 | 225 +++++++++++++++++++++++++++++++++++++++++ 9 files changed, 1003 insertions(+), 47 deletions(-) create mode 100644 src/adm_metric.F90 delete mode 100644 src/depend.pl create mode 100644 src/gram_schmidt.F90 create mode 100644 src/isnan.c create mode 100644 src/ricci.F90 create mode 100644 src/ricci2.F90 create mode 100644 src/ricci4.F90 create mode 100644 src/timederivs.F90 diff --git a/src/adm_metric.F90 b/src/adm_metric.F90 new file mode 100644 index 0000000..3df077e --- /dev/null +++ b/src/adm_metric.F90 @@ -0,0 +1,277 @@ +! $Header$ + +#include "cctk.h" + +module adm_metric + use tensor + implicit none + private + public calc_4metric + public calc_4metricderivs + public calc_4metricderivs2 + public calc_3metric + public calc_3metricderivs + + public calc_3metricdot + public calc_extcurv +contains + + subroutine calc_4metric (gg, alfa, beta, g4) + CCTK_REAL, intent(in) :: gg(3,3), alfa, beta(3) + CCTK_REAL, intent(out) :: g4(0:3,0:3) + CCTK_REAL :: betal(3) + + ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt) + betal = matmul(gg, beta) + g4(0 ,0 ) = -alfa**2 + sum(betal*beta) + g4(1:3,0 ) = betal + g4(0 ,1:3) = betal + g4(1:3,1:3) = gg + end subroutine calc_4metric + + subroutine calc_4metricderivs (gg,alfa,beta, dgg,dalfa,dbeta, & + gg_dot,alfa_dot,beta_dot, g4,dg4) + CCTK_REAL, intent(in) :: gg(3,3),alfa,beta(3) + CCTK_REAL, intent(in) :: dgg(3,3,3),dalfa(3),dbeta(3,3) + CCTK_REAL, intent(in) :: gg_dot(3,3),alfa_dot,beta_dot(3) + CCTK_REAL, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3) + CCTK_REAL :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3) + CCTK_REAL :: betal(3),d4betal(3,0:3) + integer :: i,j + + ! 4-metric + forall (i=1:3) + betal(i) = sum(gg(i,:) * beta(:)) + end forall + g4(0 ,0 ) = -alfa**2 + sum(betal(:) * beta(:)) + g4(1:3,0 ) = betal + g4(0 ,1:3) = betal + g4(1:3,1:3) = gg + + ! derivatives + d4gg (:,:,0 ) = gg_dot(:,:) + d4gg (:,:,1:3) = dgg(:,:,:) + d4alfa( 0 ) = alfa_dot + d4alfa( 1:3) = dalfa(:) + d4beta(:, 0 ) = beta_dot(:) + d4beta(:, 1:3) = dbeta(:,:) + + forall (i=1:3, j=0:3) + d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j)) + end forall + + forall (i=0:3) + dg4(0 ,0 ,i) = - 2 * alfa * d4alfa(i) & + & + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i)) + dg4(1:3,0 ,i) = d4betal(:,i) + dg4(0 ,1:3,i) = d4betal(:,i) + dg4(1:3,1:3,i) = d4gg(:,:,i) + end forall + end subroutine calc_4metricderivs + + subroutine calc_4metricderivs2 (gg,alfa,beta, dgg,dalfa,dbeta, & + ddgg,ddalfa,ddbeta, gg_dot,alfa_dot,beta_dot, & + gg_dot2,alfa_dot2,beta_dot2, dgg_dot,dalfa_dot,dbeta_dot, g4,dg4,ddg4) + CCTK_REAL, intent(in) :: gg(3,3),alfa,beta(3) + CCTK_REAL, intent(in) :: dgg(3,3,3),dalfa(3),dbeta(3,3) + CCTK_REAL, intent(in) :: ddgg(3,3,3,3),ddalfa(3,3),ddbeta(3,3,3) + CCTK_REAL, intent(in) :: gg_dot(3,3),alfa_dot,beta_dot(3) + CCTK_REAL, intent(in) :: gg_dot2(3,3),alfa_dot2,beta_dot2(3) + CCTK_REAL, intent(in) :: dgg_dot(3,3,3),dalfa_dot(3),dbeta_dot(3,3) + CCTK_REAL, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3) + CCTK_REAL, intent(out) :: ddg4(0:3,0:3,0:3,0:3) + CCTK_REAL :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3) + CCTK_REAL :: dd4gg(3,3,0:3,0:3),dd4alfa(0:3,0:3),dd4beta(3,0:3,0:3) + CCTK_REAL :: betal(3),d4betal(3,0:3),dd4betal(3,0:3,0:3) + integer :: i,j,k + + ! 4-metric + forall (i=1:3) + betal(i) = sum(gg(i,:) * beta(:)) + end forall + g4(0 ,0 ) = -alfa**2 + sum(betal(:) * beta(:)) + g4(1:3,0 ) = betal + g4(0 ,1:3) = betal + g4(1:3,1:3) = gg + + ! first derivatives + d4gg (:,:,0 ) = gg_dot(:,:) + d4gg (:,:,1:3) = dgg(:,:,:) + d4alfa( 0 ) = alfa_dot + d4alfa( 1:3) = dalfa(:) + d4beta(:, 0 ) = beta_dot(:) + d4beta(:, 1:3) = dbeta(:,:) + + forall (i=1:3, j=0:3) + d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j)) + end forall + + forall (i=0:3) + dg4(0 ,0 ,i) = - 2 * alfa * d4alfa(i) & + & + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i)) + dg4(1:3,0 ,i) = d4betal(:,i) + dg4(0 ,1:3,i) = d4betal(:,i) + dg4(1:3,1:3,i) = d4gg(:,:,i) + end forall + + ! second derivatives + dd4gg (:,:,0 ,0 ) = gg_dot2(:,:) + dd4gg (:,:,1:3,0 ) = dgg_dot(:,:,:) + dd4gg (:,:,0 ,1:3) = dgg_dot(:,:,:) + dd4gg (:,:,1:3,1:3) = ddgg(:,:,:,:) + dd4alfa( 0 ,0 ) = alfa_dot2 + dd4alfa( 1:3,0 ) = dalfa_dot(:) + dd4alfa( 0 ,1:3) = dalfa_dot(:) + dd4alfa( 1:3,1:3) = ddalfa(:,:) + dd4beta(:, 0 ,0 ) = beta_dot2(:) + dd4beta(:, 1:3,0 ) = dbeta_dot(:,:) + dd4beta(:, 0 ,1:3) = dbeta_dot(:,:) + dd4beta(:, 1:3,1:3) = ddbeta(:,:,:) + + ! betal(i) = gg(i,m) * beta(m) + ! d4betal(i,j) = d4gg(i,m,j) * beta(m) + gg(i,m) * d4beta(m,j) + ! dd4betal(i,j,k) = dd4gg(i,m,j,k) * beta(m) + d4gg(i,m,j) * d4beta(m,k) + ! + d4gg(i,m,k) * d4beta(m,j) + gg(i,m) * dd4beta(m,j,k) + forall (i=1:3, j=0:3, k=0:3) + dd4betal(i,j,k) = sum(+ dd4gg(i,:,j,k) * beta(:) & + & + d4gg(i,:,j) * d4beta(:,k) & + & + d4gg(i,:,k) * d4beta(:,j) & + & + gg(i,:) * dd4beta(:,j,k)) + end forall + + ! g4(0 ,0 ) = -alfa**2 + sum(betal*beta) + ! g4(1:3,0 ) = betal + ! g4(0 ,1:3) = betal + ! g4(1:3,1:3) = gg + + ! dg4(0 ,0 ,i) = -2*alfa*d4alfa(i) & + ! & + d4betal(m,i)*beta(m) + betal(m)*d4beta(m,i) + ! dg4(1:3,0 ,i) = d4betal(:,i) + ! dg4(0 ,1:3,i) = d4betal(:,i) + ! dg4(1:3,1:3,i) = d4gg(:,:,i) + forall (i=0:3, j=0:3) + ddg4(0 ,0 ,i,j) = - 2 * d4alfa(j) * d4alfa(i) & + & - 2 * alfa * dd4alfa(i,j) & + & + sum(+ dd4betal(:,i,j) * beta(:) & + & + d4betal(:,i) * d4beta(:,j) & + & + d4betal(:,j) * d4beta(:,i) & + & + betal(:) * dd4beta(:,i,j)) + ddg4(1:3,0 ,i,j) = dd4betal(:,i,j) + ddg4(0 ,1:3,i,j) = dd4betal(:,i,j) + ddg4(1:3,1:3,i,j) = dd4gg(:,:,i,j) + end forall + end subroutine calc_4metricderivs2 + + + + subroutine calc_3metric (g4, gg, alfa, beta) + CCTK_REAL, intent(in) :: g4(0:3,0:3) + CCTK_REAL, intent(out) :: gg(3,3), alfa, beta(3) + CCTK_REAL :: betal(3) + CCTK_REAL :: dtg, gu(3,3) + + ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt) + betal = g4(1:3,0) + gg = g4(1:3,1:3) + call calc_det (gg, dtg) + call calc_inv (gg, dtg, gu) + beta = matmul(gu, betal) + alfa = sqrt(sum(betal*beta) - g4(0,0)) + end subroutine calc_3metric + + subroutine calc_3metricderivs (g4,dg4, gg,alfa,beta, dgg,dalfa,dbeta, & + gg_dot,alfa_dot,beta_dot) + CCTK_REAL, intent(in) :: g4(0:3,0:3),dg4(0:3,0:3,0:3) + CCTK_REAL, intent(out) :: gg(3,3),alfa,beta(3) + CCTK_REAL, intent(out) :: dgg(3,3,3),dalfa(3),dbeta(3,3) + CCTK_REAL, intent(out) :: gg_dot(3,3),alfa_dot,beta_dot(3) + CCTK_REAL :: betal(3),d4betal(3,0:3) + CCTK_REAL :: dtg,gu(3,3),dgu(3,3,3),gu_dot(3,3) + CCTK_REAL :: d4gg(3,3,0:3),d4gu(3,3,0:3) + CCTK_REAL :: d4alfa(0:3),d4beta(3,0:3) + integer :: i,j + + ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt) + gg = g4(1:3,1:3) + call calc_det (gg, dtg) + call calc_inv (gg, dtg, gu) + betal = g4(1:3,0) + beta = matmul(gu, betal) + alfa = sqrt(sum(betal*beta) - g4(0,0)) + + forall (i=0:3) + d4gg(:,:,i) = dg4(1:3,1:3,i) + end forall + gg_dot = d4gg(:,:,0) + dgg = d4gg(:,:,1:3) + + call calc_invderiv (gu, dgg, dgu) + call calc_invdot (gu, gg_dot, gu_dot) + d4gu(:,:,0) = gu_dot + d4gu(:,:,1:3) = dgu + + forall (i=0:3) + d4betal(:,i) = dg4(1:3,0,i) + end forall + forall (i=0:3, j=1:3) + d4beta(j,i) = sum(d4gu(j,:,i)*betal) + sum(gu(j,:)*d4betal(:,i)) + end forall + forall (i=0:3) + d4alfa(i) = 1/(2*alfa) & + * (sum(d4betal(:,i)*beta) + sum(betal*d4beta(:,i)) - dg4(0,0,i)) + end forall + + alfa_dot = d4alfa(0) + dalfa = d4alfa(1:3) + beta_dot = d4beta(:,0) + dbeta = d4beta(:,1:3) + end subroutine calc_3metricderivs + + + + subroutine calc_3metricdot (gg, dgg, kk, alfa, beta, dbeta, gg_dot) + CCTK_REAL, intent(in) :: gg(3,3), dgg(3,3,3) + CCTK_REAL, intent(in) :: kk(3,3) + CCTK_REAL, intent(in) :: alfa + CCTK_REAL, intent(in) :: beta(3), dbeta(3,3) + CCTK_REAL, intent(out) :: gg_dot(3,3) + integer :: i,j,k + + ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k + + do i=1,3 + do j=1,3 + gg_dot(i,j) = -2*alfa * kk(i,j) + do k=1,3 + gg_dot(i,j) = gg_dot(i,j) & + + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) & + + beta(k) * dgg(i,j,k) + end do + end do + end do + end subroutine calc_3metricdot + + subroutine calc_extcurv (gg, dgg, gg_dot, alfa, beta, dbeta, kk) + CCTK_REAL, intent(in) :: gg(3,3), dgg(3,3,3), gg_dot(3,3) + CCTK_REAL, intent(in) :: alfa + CCTK_REAL, intent(in) :: beta(3), dbeta(3,3) + CCTK_REAL, intent(out) :: kk(3,3) + integer :: i,j,k + + ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k + + do i=1,3 + do j=1,3 + kk(i,j) = - gg_dot(i,j) + do k=1,3 + kk(i,j) = kk(i,j) & + + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) & + + beta(k) * dgg(i,j,k) + end do + kk(i,j) = kk(i,j) / (2*alfa) + end do + end do + + end subroutine calc_extcurv + +end module adm_metric diff --git a/src/depend.pl b/src/depend.pl deleted file mode 100644 index b41b8f6..0000000 --- a/src/depend.pl +++ /dev/null @@ -1,46 +0,0 @@ -#!/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 () { - 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/gram_schmidt.F90 b/src/gram_schmidt.F90 new file mode 100644 index 0000000..61e865a --- /dev/null +++ b/src/gram_schmidt.F90 @@ -0,0 +1,121 @@ +! $Header$ + +#include "cctk.h" + +module gram_schmidt + implicit none + private + public gram_schmidt_project4 + public gram_schmidt_normalise4 + +contains + + subroutine gram_schmidt_project4 (g, y, dy, ddy, y2, x, dx, ddx) + CCTK_REAL, intent(in) :: g(0:3,0:3) + CCTK_REAL, intent(in) :: y(0:3), dy(0:3,0:3), ddy(0:3,0:3,0:3), y2 + CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3) + CCTK_REAL :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3) + + integer :: a, b, c, d, e + + do a=0,3 + z(a) = x(a) + do d=0,3 + do e=0,3 + z(a) = z(a) - g(d,e) * x(d) * y(e) * y(a) / y2 + end do + end do + end do + + do a=0,3 + do b=0,3 + dz(a,b) = dx(a,b) + do d=0,3 + do e=0,3 + dz(a,b) = dz(a,b) - g(d,e) * ((dx(d,b) * y(e) + x(d) * dy(e,b)) * y(a) + x(d) * y(e) * dy(a,b)) / y2 + end do + end do + end do + end do + + do a=0,3 + do b=0,3 + do c=0,3 + ddz(a,b,c) = ddx(a,b,c) + do d=0,3 + do e=0,3 + ddz(a,b,c) = ddz(a,b,c) - g(d,e) * ( (ddx(d,b,c) * y(e) + dx(d,b) * dy(e,c) + dx(d,c) * dy(e,b) + x(d) * ddy(e,b,c)) * y(a) + (dx(d,b) * y(e) + x(d) * dy(e,b)) * dy(a,c) & + & + dx(d,c) * y(e) * dy(a,b) + x(d) * dy(e,c) * dy(a,b) + x(d) * y(e) * ddy(a,b,c)) / y2 + end do + end do + end do + end do + end do + + x = z + dx = dz + ddx = ddz + end subroutine gram_schmidt_project4 + + subroutine gram_schmidt_normalise4 (g, x, dx, ddx, x2) + CCTK_REAL, intent(in) :: g(0:3,0:3) + CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3) + CCTK_REAL, intent(in) :: x2 + CCTK_REAL, parameter :: two=2, half=1/two + CCTK_REAL :: z2, dz2(0:3), ddz2(0:3,0:3) + CCTK_REAL :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3) + + integer :: a, b, c, d, e + + z2 = 0 + do d=0,3 + do e=0,3 + z2 = z2 + g(d,e) * x(d) * x(e) + end do + end do + + do a=0,3 + dz2(a) = 0 + do d=0,3 + do e=0,3 + dz2(a) = dz2(a) + 2 * g(d,e) * dx(d,a) * x(e) + end do + end do + end do + + do a=0,3 + do b=0,3 + ddz2(a,b) = 0 + do d=0,3 + do e=0,3 + ddz2(a,b) = ddz2(a,b) + 2 * g(d,e) * ddx(d,a,b) * x(e) + 2 * g(d,e) * dx(d,a) * dx(e,b) + end do + end do + end do + end do + + do a=0,3 + z(a) = x(a) / sqrt(z2 / x2) + end do + + do a=0,3 + do b=0,3 + dz(a,b) = (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) / sqrt(z2 / x2)**3 + end do + end do + + do a=0,3 + do b=0,3 + do c=0,3 + ddz(a,b,c) = ((ddx(a,b,c) * z2 + dx(a,b) * dz2(c) / x2 - half * (dx(a,c) * dz2(b) / x2 + x(a) * ddz2(b,c) / x2**2)) * z2 & + & - 3*half * (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) * dz2(c) / x2) / sqrt(z2 / x2)**5 + end do + end do + end do + + x = z + dx = dz + ddx = ddz + end subroutine gram_schmidt_normalise4 + +end module gram_schmidt diff --git a/src/isnan.c b/src/isnan.c new file mode 100644 index 0000000..a2feeb6 --- /dev/null +++ b/src/isnan.c @@ -0,0 +1,32 @@ +/* $Header$ */ + +#include + +#include "cctk.h" + +int CCTK_FCALL CCTK_FNAME(TAT_isnan) (const CCTK_REAL * restrict const x) +{ +#ifdef HAVE_ISNAN + return isnan(*x); +#else + return 0; +#endif +} + +int CCTK_FCALL CCTK_FNAME(TAT_isinf) (const CCTK_REAL * restrict const x) +{ +#ifdef HAVE_ISINF + return isinf(*x); +#else + return 0; +#endif +} + +int CCTK_FCALL CCTK_FNAME(TAT_finite) (const CCTK_REAL * restrict const x) +{ +#ifdef HAVE_FINITE + return finite(*x); +#else + return 1; +#endif +} diff --git a/src/make.code.defn b/src/make.code.defn index 69e3078..462985e 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,21 +2,28 @@ # $Header$ # Source files in this directory -SRCS = cactus.F90 \ +SRCS = adm_metric.F90 \ + cactus.F90 \ constants.F90 \ conversion.F90 \ covderivs.F90 \ covderivs2.F90 \ derivs.F90 \ derivs2.F90 \ + gram_schmidt.F90 \ + isnan.c \ lapack.F90 \ matdet.F90 \ matexp.F90 \ matinv.F90 \ + ricci.F90 \ + ricci2.F90 \ + ricci4.F90 \ rotation.F90 \ tensor.F90 \ tensor2.F90 \ tensor4.F90 \ + timederivs.F90 \ pointwise.F90 \ pointwise2.F90 diff --git a/src/ricci.F90 b/src/ricci.F90 new file mode 100644 index 0000000..ebe01fc --- /dev/null +++ b/src/ricci.F90 @@ -0,0 +1,85 @@ +! $Header$ + +#include "cctk.h" + +module ricci + implicit none + private + public calc_connections + public calc_connectionderivs + public calc_ricci + +contains + + subroutine calc_connections (gu, dg, gamma) + CCTK_REAL, intent(in) :: gu(3,3), dg(3,3,3) + CCTK_REAL, intent(out) :: gamma(3,3,3) + integer :: i,j,k,l + ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l) + do i=1,3 + do j=1,3 + do k=1,3 + gamma(i,j,k) = 0 + do l=1,3 + gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) & + * (dg(l,j,k) + dg(l,k,j) - dg(j,k,l)) + end do + end do + end do + end do + end subroutine calc_connections + + subroutine calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma) + CCTK_REAL, intent(in) :: gu(3,3), dgg(3,3,3), dgu(3,3,3), ddgg(3,3,3,3) + CCTK_REAL, intent(out) :: dgamma(3,3,3,3) + integer :: i,j,k,l,m + ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m) + ! + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml) + do i=1,3 + do j=1,3 + do k=1,3 + do l=1,3 + dgamma(i,j,k,l) = 0 + do m=1,3 + dgamma(i,j,k,l) = dgamma(i,j,k,l) & + + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) & + + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l)) + end do + end do + end do + end do + end do + end subroutine calc_connectionderivs + + subroutine calc_ricci (gamma, dgamma, ri) + CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3) + CCTK_REAL, intent(out) :: ri(3,3) + CCTK_REAL :: sum, cnt + integer :: i,j,k,l + ! R_ij = Gamma^k_ij,k - Gamma^k_ik,j + ! + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki + do i=1,3 + do j=1,3 + ri(i,j) = 0 + do k=1,3 + ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j) + do l=1,3 + ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) & + & - gamma(k,l,j) * gamma(l,k,i) + end do + end do + end do + end do + ! check symmetries + sum = 0 + cnt = 0 + do i=1,3 + do j=1,3 + sum = sum + (ri(i,j) - ri(j,i))**2 + cnt = cnt + ri(i,j)**2 + end do + end do + if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "Ricci tensor is not symmetric") + end subroutine calc_ricci + +end module ricci diff --git a/src/ricci2.F90 b/src/ricci2.F90 new file mode 100644 index 0000000..dd2c429 --- /dev/null +++ b/src/ricci2.F90 @@ -0,0 +1,74 @@ +! $Header$ + +#include "cctk.h" + +module ricci2 + implicit none + private + public calc_2connections + public calc_2connectionderivs + public calc_2ricci + +contains + + subroutine calc_2connections (gu, dgg, gamma) + CCTK_REAL, intent(in) :: gu(2,2), dgg(2,2,2) + CCTK_REAL, intent(out) :: gamma(2,2,2) + integer :: i,j,k,l + ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l) + do i=1,2 + do j=1,2 + do k=1,2 + gamma(i,j,k) = 0 + do l=1,2 + gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) & + * (dgg(l,j,k) + dgg(l,k,j) - dgg(j,k,l)) + end do + end do + end do + end do + end subroutine calc_2connections + + subroutine calc_2connectionderivs (gu, dgg, dgu, ddgg, dgamma) + CCTK_REAL, intent(in) :: gu(2,2), dgg(2,2,2), dgu(2,2,2), ddgg(2,2,2,2) + CCTK_REAL, intent(out) :: dgamma(2,2,2,2) + integer :: i,j,k,l,m + ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m) + ! + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml) + do i=1,2 + do j=1,2 + do k=1,2 + do l=1,2 + dgamma(i,j,k,l) = 0 + do m=1,2 + dgamma(i,j,k,l) = dgamma(i,j,k,l) & + + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) & + + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l)) + end do + end do + end do + end do + end do + end subroutine calc_2connectionderivs + + subroutine calc_2ricci (gamma, dgamma, ri) + CCTK_REAL, intent(in) :: gamma(2,2,2), dgamma(2,2,2,2) + CCTK_REAL, intent(out) :: ri(2,2) + integer :: i,j,k,l + ! R_ij = Gamma^k_ij,k - Gamma^k_ik,j + ! + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki + do i=1,2 + do j=1,2 + ri(i,j) = 0 + do k=1,2 + ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j) + do l=1,2 + ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) & + & - gamma(k,l,j) * gamma(l,k,i) + end do + end do + end do + end do + end subroutine calc_2ricci + +end module ricci2 diff --git a/src/ricci4.F90 b/src/ricci4.F90 new file mode 100644 index 0000000..ada7733 --- /dev/null +++ b/src/ricci4.F90 @@ -0,0 +1,181 @@ +! $Header$ + +#include "cctk.h" + +module ricci4 + implicit none + private + public calc_4connections + public calc_4connectionderivs + public calc_4ricci + public calc_4riemann + public calc_4weyl + +contains + + subroutine calc_4connections (gu, dg, gamma) + CCTK_REAL, intent(in) :: gu(4,4), dg(4,4,4) + CCTK_REAL, intent(out) :: gamma(4,4,4) + integer :: i,j,k,l + ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l) + do i=1,4 + do j=1,4 + do k=1,4 + gamma(i,j,k) = 0 + do l=1,4 + gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) & + * (dg(l,j,k) + dg(l,k,j) - dg(j,k,l)) + end do + end do + end do + end do + end subroutine calc_4connections + + subroutine calc_4connectionderivs (gu, dgg, dgu, ddgg, dgamma) + CCTK_REAL, intent(in) :: gu(4,4), dgg(4,4,4), dgu(4,4,4), ddgg(4,4,4,4) + CCTK_REAL, intent(out) :: dgamma(4,4,4,4) + integer :: i,j,k,l,m + ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m) + ! + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml) + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + dgamma(i,j,k,l) = 0 + do m=1,4 + dgamma(i,j,k,l) = dgamma(i,j,k,l) & + + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) & + + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l)) + end do + end do + end do + end do + end do + end subroutine calc_4connectionderivs + + subroutine calc_4ricci (gamma, dgamma, ri) + CCTK_REAL, intent(in) :: gamma(4,4,4), dgamma(4,4,4,4) + CCTK_REAL, intent(out) :: ri(4,4) + CCTK_REAL :: sum, cnt + integer :: i,j,k,l + ! R_ij = Gamma^k_ij,k - Gamma^k_ik,j + ! + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki + do i=1,4 + do j=1,4 + ri(i,j) = 0 + do k=1,4 + ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j) + do l=1,4 + ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) & + & - gamma(k,l,j) * gamma(l,k,i) + end do + end do + end do + end do + ! check symmetries + sum = 0 + cnt = 0 + do i=1,4 + do j=1,4 + sum = sum + (ri(i,j) - ri(j,i))**2 + cnt = cnt + ri(i,j)**2 + end do + end do + if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Ricci tensor is not symmetric") + end subroutine calc_4ricci + + subroutine calc_4riemann (gg, gamma, dgamma, rm) + CCTK_REAL, intent(in) :: gg(4,4), gamma(4,4,4), dgamma(4,4,4,4) + CCTK_REAL, intent(out) :: rm(4,4,4,4) + CCTK_REAL :: rmu(4,4,4,4) + CCTK_REAL :: sum, cnt + integer :: i,j,k,l,m + ! R^i_jkl = Gamma^i_jl,k - Gamma^i_jk,l + ! + Gamma^m_jl Gamma^i_mk - Gamma^m_jk Gamma^i_ml + ! (Wald, 3.4, eqn. (3.4.4), p. 48) + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + rmu(i,j,k,l) = dgamma(i,j,l,k) - dgamma(i,j,k,l) + do m=1,4 + rmu(i,j,k,l) = rmu(i,j,k,l) & + + gamma(m,j,l) * gamma(i,m,k) & + - gamma(m,j,k) * gamma(i,m,l) + end do + end do + end do + end do + end do + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + rm(i,j,k,l) = 0 + do m=1,4 + rm(i,j,k,l) = rm(i,j,k,l) + gg(i,m) * rmu(m,j,k,l) + end do + end do + end do + end do + end do + ! check symmetries + sum = 0 + cnt = 0 + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + sum = sum + (rm(i,j,k,l) + rm(i,j,l,k))**2 + sum = sum + (rm(i,j,k,l) - rm(k,l,i,j))**2 + sum = sum + (rm(i,j,k,l) + rm(j,i,k,l))**2 + sum = sum + (rm(i,j,k,l) + rm(j,k,i,l) + rm(k,i,j,l))**2 + cnt = cnt + rm(i,j,k,l)**2 + end do + end do + end do + end do + if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Riemann tensor is not symmetric") + end subroutine calc_4riemann + + subroutine calc_4weyl (gg, rm, ri, rsc, we) + CCTK_REAL, parameter :: one=1, half=one/2, third=one/3 + CCTK_REAL, intent(in) :: gg(4,4), rm(4,4,4,4), ri(4,4), rsc + CCTK_REAL, intent(out) :: we(4,4,4,4) + CCTK_REAL :: sum, cnt + integer :: i,j,k,l + ! R_ijkl = C_ijkl + 2/(n-2) (g_i[k R_l]j - g_j[k R_l]i) + ! - 2/(n-1)(n-2) R g_i[k g_l]j + ! (Wald, 3.2, eqn. (3.2.28), p. 40) + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + we(i,j,k,l) = rm(i,j,k,l) & + - ( (gg(i,k) * ri(l,j) - gg(i,l) * ri(k,j)) & + & - (gg(j,k) * ri(l,i) - gg(j,l) * ri(k,i))) & + + third * rsc * (gg(i,k) * gg(l,j) - gg(i,l) * gg(k,j)) + end do + end do + end do + end do + ! check symmetries + sum = 0 + cnt = 0 + do i=1,4 + do j=1,4 + do k=1,4 + do l=1,4 + sum = sum + (we(i,j,k,l) + we(i,j,l,k))**2 + sum = sum + (we(i,j,k,l) - we(k,l,i,j))**2 + sum = sum + (we(i,j,k,l) + we(j,i,k,l))**2 + sum = sum + (we(i,j,k,l) + we(j,k,i,l) + we(k,i,j,l))**2 + cnt = cnt + we(i,j,k,l)**2 + end do + end do + end do + end do + if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Weyl tensor is not symmetric") + end subroutine calc_4weyl + +end module ricci4 diff --git a/src/timederivs.F90 b/src/timederivs.F90 new file mode 100644 index 0000000..42e107d --- /dev/null +++ b/src/timederivs.F90 @@ -0,0 +1,225 @@ +! $Header$ + +#include "cctk.h" + +module timederivs + implicit none + private + public abs2 + public timederiv + public timederiv2 + public timederiv_uneven + public timederiv2_uneven + public operator(.outer.) + public operator(.dot.) + + interface timederiv + module procedure rtimederiv + end interface + + interface timederiv2 + module procedure rtimederiv2 + end interface + + interface timederiv_uneven + module procedure rtimederiv_uneven + module procedure ctimederiv_uneven + end interface + + interface timederiv2_uneven + module procedure rtimederiv2_uneven + module procedure ctimederiv2_uneven + end interface + + interface operator(.outer.) + module procedure router + module procedure couter + end interface + + interface operator(.dot.) + module procedure rdot + module procedure cdot + end interface + +contains + + ! abs(c)**2 for complex c without a square root + elemental function abs2 (a) + CCTK_COMPLEX, intent(in) :: a + CCTK_REAL :: abs2 + abs2 = a * conjg(a) + end function abs2 + + + + ! Calculate a time derivate from several time levels + elemental function rtimederiv (f0, f1, f2, dt) result (fdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: dt + CCTK_REAL :: fdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdot1, fdot2 + + dt1 = dt + dt2 = 2*dt + + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f'(0) = [f(dt1) - f(0)] / dt1 - dt1/2 f''(0) + ! f'(0) = [f(dt2) - f(0)] / dt2 - dt2/2 f''(0) + fdot1 = (f0 - f1) / dt1 + fdot2 = (f0 - f2) / dt2 + fdot = (fdot1 * dt2 - fdot2 * dt1) / (dt2 - dt1) + end function rtimederiv + + elemental function rtimederiv2 (f0, f1, f2, dt) result (fdotdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: dt + CCTK_REAL :: fdotdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdotdot1, fdotdot2 + + dt1 = dt + dt2 = 2*dt + + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f''(0) = [f(dt1) - f(0)] / [dt1^2/2] - f'(0) / [dt1/2] + ! f''(0) = [f(dt2) - f(0)] / [dt2^2/2] - f'(0) / [dt2/2] + fdotdot1 = (f1 - f0) / (dt1**2/2) + fdotdot2 = (f2 - f0) / (dt2**2/2) + fdotdot = (fdotdot1 * dt1 - fdotdot2 * dt2) / (dt1 - dt2) + end function rtimederiv2 + + + + ! Calculate a time derivate from several time levels with uneven spacing + elemental function rtimederiv_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + CCTK_INT, intent(in) :: ce0, ce1, ce2 + CCTK_REAL :: fdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdot1, fdot2 + +!!$ dt1 = ih_time(hn) - ih_time_p(hn) +!!$ dt2 = ih_time(hn) - ih_time_p_p(hn) + dt1 = t0 - t1 + dt2 = t0 - t2 + + if (ce1/=0) then + fdot = 0 + else if (ce2/=0) then + fdot = (f0 - f1) / dt1 + else + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f'(0) = [f(dt1) - f(0)] / dt1 - dt1/2 f''(0) + ! f'(0) = [f(dt2) - f(0)] / dt2 - dt2/2 f''(0) + fdot1 = (f0 - f1) / dt1 + fdot2 = (f0 - f2) / dt2 + fdot = (fdot1 * dt2 - fdot2 * dt1) / (dt2 - dt1) + end if + end function rtimederiv_uneven + + elemental function ctimederiv_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot) + CCTK_COMPLEX, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + CCTK_INT, intent(in) :: ce0, ce1, ce2 + CCTK_COMPLEX :: fdot + + fdot = cmplx(timederiv_uneven(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), & + & timederiv_uneven(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), & + & kind(fdot)) + end function ctimederiv_uneven + + elemental function rtimederiv2_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + CCTK_INT, intent(in) :: ce0, ce1, ce2 + CCTK_REAL :: fdotdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdotdot1, fdotdot2 + +!!$ dt1 = ih_time(hn) - ih_time_p(hn) +!!$ dt2 = ih_time(hn) - ih_time_p_p(hn) + dt1 = t0 - t1 + dt2 = t0 - t2 + + if (ce1/=0) then + fdotdot = 0 + else if (ce2/=0) then + fdotdot = 0 + else + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f''(0) = [f(dt1) - f(0)] / [dt1^2/2] - f'(0) / [dt1/2] + ! f''(0) = [f(dt2) - f(0)] / [dt2^2/2] - f'(0) / [dt2/2] + fdotdot1 = (f1 - f0) / (dt1**2/2) + fdotdot2 = (f2 - f0) / (dt2**2/2) + fdotdot = (fdotdot1 * dt1 - fdotdot2 * dt2) / (dt1 - dt2) + end if + end function rtimederiv2_uneven + + elemental function ctimederiv2_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot) + CCTK_COMPLEX, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + CCTK_INT, intent(in) :: ce0, ce1, ce2 + CCTK_COMPLEX :: fdotdot + + fdotdot = cmplx(timederiv2_uneven(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), & + & timederiv2_uneven(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), & + & kind(fdotdot)) + end function ctimederiv2_uneven + + + + function router (left, right) result (result) + CCTK_REAL, intent(in) :: left(:), right(:) + CCTK_REAL :: result(size(left,1),size(right,1)) + integer :: i, j + forall (i=1:size(left,1), j=1:size(right,1)) + result(i,j) = left(i) * right(j) + end forall + end function router + + function couter (left, right) result (result) + CCTK_COMPLEX, intent(in) :: left(:), right(:) + CCTK_COMPLEX :: result(size(left,1),size(right,1)) + integer :: i, j + forall (i=1:size(left,1), j=1:size(right,1)) + result(i,j) = left(i) * right(j) + end forall + end function couter + + + + function rdot (left, right) result (result) + CCTK_REAL, intent(in) :: left(:), right(:) + CCTK_REAL :: result + integer :: i + if (size(left,1) /= size(right,1)) then + call CCTK_WARN (0, "Array sizes must have the same sizes") + end if +!!$ result = sum((/( left(i) * right(i), i=1,size(left,1) )/)) + result = 0 + do i=1,size(left,1) + result = result + left(i) * right(i) + end do + end function rdot + + function cdot (left, right) result (result) + CCTK_COMPLEX, intent(in) :: left(:), right(:) + CCTK_COMPLEX :: result + integer :: i + if (size(left,1) /= size(right,1)) then + call CCTK_WARN (0, "Array sizes must have the same sizes") + end if +!!$ result = sum((/( left(i) * right(i), i=1,size(left,1) )/)) + result = 0 + do i=1,size(left,1) + result = result + left(i) * right(i) + end do + end function cdot + +end module timederivs -- cgit v1.2.3