From 310f0ea48d18866b773136aed11200b6eda6378b Mon Sep 17 00:00:00 2001 From: eschnett <> Date: Thu, 1 Mar 2001 11:40:00 +0000 Subject: Initial revision darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz --- CarpetExtra/WaveToyFO/src/calcrhs.F77 | 98 +++++++++++++++++++++++++++++++++++ 1 file changed, 98 insertions(+) create mode 100644 CarpetExtra/WaveToyFO/src/calcrhs.F77 (limited to 'CarpetExtra/WaveToyFO/src/calcrhs.F77') diff --git a/CarpetExtra/WaveToyFO/src/calcrhs.F77 b/CarpetExtra/WaveToyFO/src/calcrhs.F77 new file mode 100644 index 000000000..55734ed1e --- /dev/null +++ b/CarpetExtra/WaveToyFO/src/calcrhs.F77 @@ -0,0 +1,98 @@ +c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyFO/src/calcrhs.F77,v 1.8 2004/08/28 19:10:33 schnetter Exp $ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + + function WaveToyFO_Deriv (fm2, fm1, f0, fp1, fp2, dx, dir) + implicit none + CCTK_REAL WaveToyFO_Deriv + CCTK_REAL fm2, fm1, f0, fp1, fp2 + CCTK_REAL dx + integer dir + + if (dir .eq. 0) then + WaveToyFO_Deriv = (fp1 - fm1) / (2 * dx) + else if (dir .eq. -1) then +c Use first order one-sided derivatives here. +c These are less dissipative than second order one-sided derivatives. +c (Sic.) + WaveToyFO_Deriv = (f0 - fm1) / dx +c$$$ WaveToyFO_Deriv = (+ 3*f0 - 4*fm1 + fm2) / (2*dx) + else if (dir .eq. +1) then + WaveToyFO_Deriv = (fp1 - f0) / dx +c$$$ WaveToyFO_Deriv = (- 3*f0 + 4*fp1 - fp2) / (2*dx) + else + call CCTK_WARN (0, "internal error") + end if + + end + + + +#define CHOOSEDIR(i,imin,imax) ((min(i+1,imax) - (i+1)) + (max(i-1,imin) - (i-1))) + +#define DIFFX(var,i,j,k) WaveToyFO_Deriv(var((i)-2,j,k), var((i)-1,j,k), var(i,j,k), var((i)+1,j,k), var((i)+2,j,k), dx(1), idir(1)) +#define DIFFY(var,i,j,k) WaveToyFO_Deriv(var(i,(j)-2,k), var(i,(j)-1,k), var(i,j,k), var(i,(j)+1,k), var(i,(j)+2,k), dx(2), idir(2)) +#define DIFFZ(var,i,j,k) WaveToyFO_Deriv(var(i,j,(k)-2), var(i,j,(k)-1), var(i,j,k), var(i,j,(k)+1), var(i,j,(k)+2), dx(3), idir(3)) + + + + subroutine WaveToyFO_CalcRHS (CCTK_ARGUMENTS) + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL dx(3) + integer bndwidth(3) + integer imin(3), imax(3), idir(3) + integer i, j, k + integer d + + external WaveToyFO_Deriv + CCTK_REAL WaveToyFO_Deriv + + do d=1,3 + bndwidth(d) = 0 + end do + + do d=1,3 + if (cctk_bbox(2*d).ne.0) then + imin(d) = 1+bndwidth(d) + else + imin(d) = 1+cctk_nghostzones(d) + end if + if (cctk_bbox(2*d+1).ne.0) then + imax(d) = cctk_lsh(d)-bndwidth(d) + else + imax(d) = cctk_lsh(d)-cctk_nghostzones(d) + end if + end do + + do d=1,3 + dx(d) = CCTK_DELTA_SPACE(d) + end do + + do k=imin(3),imax(3) + idir(3) = CHOOSEDIR(k,imin(3),imax(3)) + do j=imin(2),imax(2) + idir(2) = CHOOSEDIR(j,imin(2),imax(2)) + do i=imin(1),imax(1) + idir(1) = CHOOSEDIR(i,imin(1),imax(1)) + + phidot(i,j,k) = DIFFX (psix,i,j,k) + $ + DIFFY (psiy,i,j,k) + $ + DIFFZ (psiz,i,j,k) + psixdot(i,j,k) = DIFFX (phi,i,j,k) + psiydot(i,j,k) = DIFFY (phi,i,j,k) + psizdot(i,j,k) = DIFFZ (phi,i,j,k) + + end do + end do + end do + + end -- cgit v1.2.3