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/IDFOScalarWave/src/InitialData.F77 | 279 +++++++++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100644 CarpetExtra/IDFOScalarWave/src/InitialData.F77 (limited to 'CarpetExtra/IDFOScalarWave/src/InitialData.F77') diff --git a/CarpetExtra/IDFOScalarWave/src/InitialData.F77 b/CarpetExtra/IDFOScalarWave/src/InitialData.F77 new file mode 100644 index 000000000..aa361f8a1 --- /dev/null +++ b/CarpetExtra/IDFOScalarWave/src/InitialData.F77 @@ -0,0 +1,279 @@ +c -*-Fortran-*- +c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/InitialData.F77,v 1.8 2003/11/05 16:18:40 schnetter Exp $ + + /*@@ + @file InitialData.F77 + @date + @author Scott Hawley, using code from Tom Goodale, Erik Schnetter + @desc + Initial data for the 3D Wave Equation + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + /*@@ + @routine IDFOScalarWave_InitialData + @date + @author Scott Hawley, using code from Tom Goodale, Erik Schnetter + @desc + Set up initial data for the wave equation + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + subroutine IDFOScalarWave_InitialData (CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + INTEGER i,j,k + CCTK_REAL dt,omega, cpi + CCTK_REAL ri3 + +c call CCTK_INFO ("IDFOScalarWave_InitialData") + + cpi = 4.d0*atan(1.d0) + + dt = CCTK_DELTA_TIME + + omega = sqrt(kx**2+ky**2+kz**2) + + if (CCTK_EQUALS(initial_data,"plane")) then + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + phi(i,j,k) = amplitude + $ * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) + $ + omega*cctk_time) * cpi) + + phi_p(i,j,k) = amplitude + $ * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) + $ + omega*(cctk_time - dt)) * cpi) + + phi_p_p(i,j,k) = amplitude + $ * cos((kx*x(i,j,k) + ky*y(i,j,k) + kz*z(i,j,k) + $ + omega*(cctk_time - 2*dt)) * cpi) + + end do + end do + end do + + else if (CCTK_EQUALS(initial_data,"gaussian")) then + + do k=1, cctk_lsh(3) + do j=1, cctk_lsh(2) + do i=1, cctk_lsh(1) + + phi(i,j,k) = amplitude + $ * exp(- (x(i,j,k) - radius)**2 / sigma**2) + $ * exp(- (y(i,j,k) - radius)**2 / sigma**2) + $ * exp(- (z(i,j,k) - radius)**2 / sigma**2) + + pi(i,j,k) = 0.0 + phix(i,j,k) = phi(i,j,k)* (-2) * (x(i,j,k) - radius) + $ / sigma**2 + phiy(i,j,k) = phi(i,j,k)* (-2) * (y(i,j,k) - radius) + $ / sigma**2 + phiz(i,j,k) = phi(i,j,k)* (-2) * (z(i,j,k) - radius) + $ / sigma**2 + + pi_p(i,j,k) = amplitude * + & exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 )* + & ((x(i,j,k)-radius-dt) + (y(i,j,k)-radius-dt) + $ + (z(i,j,k)-radius-dt))/sigma**2 + & - amplitude * + & exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 )* + & ((x(i,j,k)-radius+dt) + (y(i,j,k)-radius+dt) + & + (z(i,j,k)-radius+dt))/sigma**2 + + pi_p_p(i,j,k) = amplitude * + & exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & ((x(i,j,k)-radius-2*dt) + (y(i,j,k)-radius-2*dt) + & + (z(i,j,k)-radius-2*dt))/sigma**2 + & - amplitude * + & exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & ((x(i,j,k)-radius+2*dt) + (y(i,j,k)-radius+2*dt) + & + (z(i,j,k)-radius+2*dt))/sigma**2 + + phix_p(i,j,k) = - amplitude * (x(i,j,k) - radius - dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) + & - amplitude * (x(i,j,k) - radius + dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) + + phix_p_p(i,j,k) = - amplitude * (x(i,j,k) - radius - 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (x(i,j,k) - radius + 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 ) + + phiy_p(i,j,k) = - amplitude * (y(i,j,k) - radius - dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) + & - amplitude * (y(i,j,k) - radius + dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) + + phiy_p_p(i,j,k) = - amplitude * (y(i,j,k) - radius - 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (y(i,j,k) - radius + 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 ) + + phiz_p(i,j,k) = - amplitude * (z(i,j,k) - radius - dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) + & - amplitude * (z(i,j,k) - radius + dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) + + phiz_p_p(i,j,k) = - amplitude * (z(i,j,k) - radius - 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (z(i,j,k) - radius + 2*dt) + & / sigma**2 + & * exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 ) + + phi_p(i,j,k) = amplitude / 2 * + & exp( -(x(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - dt)**2 / sigma**2 ) + & + amplitude / 2 * + & exp( -(x(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + dt)**2 / sigma**2 ) + + phi_p_p(i,j,k) = amplitude / 2 * + & exp( -(x(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius - 2*dt)**2 / sigma**2 ) + & + amplitude / 2 * + & exp( -(x(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y(i,j,k) - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z(i,j,k) - radius + 2*dt)**2 / sigma**2 ) + + end do + end do + end do + + else if (CCTK_EQUALS(initial_data, "box")) then + +c Use kx,ky,kz as number of modes in each direction. + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + phi(i,j,k) = amplitude + $ * sin(kx * (x(i,j,k) - 0.5d0) * cpi) + $ * sin(ky * (y(i,j,k) - 0.5d0) * cpi) + $ * sin(kz * (z(i,j,k) - 0.5d0) * cpi) + + phi_p(i,j,k) = phi(i,j,k) + $ * cos(omega * (cctk_time - dt) * cpi) + + phi_p_p(i,j,k) = phi(i,j,k) + $ * cos(omega * (cctk_time - 2*dt) * cpi) + + phi(i,j,k) = phi(i,j,k) * cos(omega * cctk_time * cpi) + + end do + end do + end do + + else if (CCTK_EQUALS(initial_data, "1/r")) then + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + pi(i,j,k) = 0.0 + phi(i,j,k) = 1 / sqrt(x(i,j,k)**2 + y(i,j,k)**2 + & + z(i,j,k)**2) + ri3 = phi(i,j,k)**3 + phix(i,j,k) = - x(i,j,k) * ri3 + phiy(i,j,k) = - y(i,j,k) * ri3 + phiz(i,j,k) = - z(i,j,k) * ri3 + + pi_p(i,j,k) = pi(i,j,k) + pi_p_p(i,j,k) = pi(i,j,k) + phix_p(i,j,k) = phix(i,j,k) + phix_p_p(i,j,k) = phix(i,j,k) + phiy_p(i,j,k) = phiy(i,j,k) + phiy_p_p(i,j,k) = phiy(i,j,k) + phiz_p(i,j,k) = phiz(i,j,k) + phiz_p_p(i,j,k) = phiz(i,j,k) + phi_p(i,j,k) = phi(i,j,k) + phi_p_p(i,j,k) = phi(i,j,k) + + end do + end do + end do + + else + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + + phi(i,j,k) = 0.0d0 + phi_p(i,j,k) = 0.0d0 + phi_p_p(i,j,k) = 0.0d0 + + end do + end do + end do + + end if + + end -- cgit v1.2.3