From 629afedf629e04240b8b2e6b3d648c4f5135bdde Mon Sep 17 00:00:00 2001 From: shawley <> Date: Tue, 26 Feb 2002 18:43:00 +0000 Subject: added/fixed gaussian initial data. uses "cartesian style" instead of "radial style" darcs-hash:20020226184332-e415b-abf73446aecba961351dc139eb34a3d1ff8d1cf8.gz --- CarpetExtra/IDFOScalarWave/src/InitialData.F77 | 174 +++++++++++++++---------- 1 file changed, 108 insertions(+), 66 deletions(-) (limited to 'CarpetExtra/IDFOScalarWave/src') diff --git a/CarpetExtra/IDFOScalarWave/src/InitialData.F77 b/CarpetExtra/IDFOScalarWave/src/InitialData.F77 index b996fdf3f..f5fb023b0 100644 --- a/CarpetExtra/IDFOScalarWave/src/InitialData.F77 +++ b/CarpetExtra/IDFOScalarWave/src/InitialData.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/InitialData.F77,v 1.3 2002/02/26 16:44:10 shawley Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/InitialData.F77,v 1.4 2002/02/26 19:43:32 shawley Exp $ /*@@ @file InitialData.F77 @@ -42,6 +42,7 @@ c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDFOScalarWave/src/Ini CCTK_REAL dt,omega, cpi CCTK_REAL x,y,z, r, ri3 CCTK_REAL dxi,dyi,dzi + CCTK_REAL tmp c call CCTK_INFO ("IDFOScalarWave_InitialData") @@ -75,81 +76,122 @@ c call CCTK_INFO ("IDFOScalarWave_InitialData") end do else if (CCTK_EQUALS(initial_data,"gaussian")) then -c -c A shell with gaussian cross-section, time reflection symmetry -c moving at unit velocity -c phi(r,t) = (r-t)/r * amp/2 * exp(-(r-r0-t)**2/sigma**2) -c + (r+t)/r * amp/2 * exp(-(r-r0+t)**2/sigma**2) do k=1, cctk_lsh(3) do j=1, cctk_lsh(2) do i=1, cctk_lsh(1) - r = spher3d_r(i,j,k) + x = cart3d_x(i,j,k) + y = cart3d_y(i,j,k) + z = cart3d_z(i,j,k) phi(i,j,k) = amplitude - $ * exp(- (r - radius)**2 / sigma**2) + $ * exp(- (x - radius)**2 / sigma**2) + $ * exp(- (y - radius)**2 / sigma**2) + $ * exp(- (z - radius)**2 / sigma**2) pi(i,j,k) = 0.0 - phix(i,j,k) = phi(i,j,k)* (-2) * (r - radius) * - $ cart3d_x(i,j,k) / r / sigma**2 - phiy(i,j,k) = phi(i,j,k)* (-2) * (r - radius) * - $ cart3d_y(i,j,k) / r / sigma**2 - phiz(i,j,k) = phi(i,j,k)* (-2) * (r - radius) * - $ cart3d_z(i,j,k) / r / sigma**2 - - - pi_p(i,j,k) = amplitude * - & exp(- (r - radius - dt)**2 / sigma**2) * - & ( (r - dt)/r * (r - radius - dt) / sigma**2 - 1/2/r ) - & + amplitude * - & exp(- (r - radius + dt)**2 / sigma**2) * - & ( -(r + dt)/r * (r - radius + dt) / sigma**2 + 1/2/r) - - pi_p_p(i,j,k) = amplitude * - & exp(- (r - radius - 2*dt)**2 / sigma**2) * - & ( (r - 2*dt)/r * (r - radius - 2*dt) / sigma**2 - & - 1/2/r ) + amplitude * - & exp(- (r - radius + 2*dt)**2 / sigma**2) * - & ( -(r + 2*dt)/r * (r - radius + 2*dt) / sigma**2 - & + 1/2/r ) - -c first store partial phi / partial r in phix as a temporary variable, -c then come back and multiply by x/r - phix_p(i,j,k) = amplitude * - & exp(- (r - radius - dt)**2 / sigma**2) * - & (dt/r**2 - (r-dt)/r*(r - radius - dt) / sigma**2 ) - & - amplitude * - & exp(- (r - radius + dt)**2 / sigma**2) * - & (dt/r**2 + (r+dt)/r*(r - radius + dt) / sigma**2 ) -c do not mix up the order of these next three lines! - phiy_p(i,j,k) = phix_p(i,j,k) * cart3d_y(i,j,k) / r - phiz_p(i,j,k) = phix_p(i,j,k) * cart3d_z(i,j,k) / r - phix_p(i,j,k) = phix_p(i,j,k) * cart3d_x(i,j,k) / r - - phix_p(i,j,k) = amplitude * - & exp(- (r - radius - 2*dt)**2 / sigma**2) * - & (2*dt/r**2 - (r-2*dt)/r*(r - radius - 2*dt) - & / sigma**2 ) - & - amplitude * - & exp(- (r - radius + 2*dt)**2 / sigma**2) * - & (2*dt/r**2 + (r+2*dt)/r*(r - radius + 2*dt) - & / sigma**2 ) -c do not mix up the order of these next three lines! - phiy_p_p(i,j,k) = phix_p_p(i,j,k) * cart3d_y(i,j,k) / r - phiz_p_p(i,j,k) = phix_p_p(i,j,k) * cart3d_z(i,j,k) / r - phix_p_p(i,j,k) = phix_p_p(i,j,k) * cart3d_x(i,j,k) / r - - phi_p(i,j,k) = amplitude/2 * (r - dt) / r - $ * exp(- (r - radius - dt)**2 / sigma**2) - $ + amplitude/2 * (r + dt) / r - $ * exp(- (r - radius + dt)**2 / sigma**2) - - phi_p_p(i,j,k) = amplitude/2 * (r - 2*dt) / r - $ * exp(- (r - radius - 2*dt)**2 / sigma**2) - $ + amplitude/2 * (r + 2*dt) / r - $ * exp(- (r - radius + 2*dt)**2 / sigma**2) + phix(i,j,k) = phi(i,j,k)* (-2) * (x - radius) / sigma**2 + phiy(i,j,k) = phi(i,j,k)* (-2) * (y - radius) / sigma**2 + phiz(i,j,k) = phi(i,j,k)* (-2) * (z - radius) / sigma**2 + resid(i,j,k) = 0.0 + pi_p(i,j,k) = amplitude * + & exp( -(x - radius - dt)**2 / sigma**2 )* + & exp( -(y - radius - dt)**2 / sigma**2 )* + & exp( -(z - radius - dt)**2 / sigma**2 )* + & ((x-radius-dt) + (y-radius-dt) + (z-radius-dt))/sigma**2 + & - amplitude * + & exp( -(x - radius + dt)**2 / sigma**2 )* + & exp( -(y - radius + dt)**2 / sigma**2 )* + & exp( -(z - radius + dt)**2 / sigma**2 )* + & ((x-radius+dt) + (y-radius+dt) + (z-radius+dt))/sigma**2 + + pi_p_p(i,j,k) = amplitude * + & exp( -(x - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z - radius - 2*dt)**2 / sigma**2 )* + & ((x-radius-2*dt) + (y-radius-2*dt) + (z-radius-2*dt))/sigma**2 + & - amplitude * + & exp( -(x - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z - radius + 2*dt)**2 / sigma**2 )* + & ((x-radius+2*dt) + (y-radius+2*dt) + (z-radius+2*dt))/sigma**2 + + phix_p(i,j,k) = - amplitude * (x - radius - dt) / sigma**2 + & * exp( -(x - radius - dt)**2 / sigma**2 )* + & exp( -(y - radius - dt)**2 / sigma**2 )* + & exp( -(z - radius - dt)**2 / sigma**2 ) + & - amplitude * (x - radius + dt) / sigma**2 + & * exp( -(x - radius + dt)**2 / sigma**2 )* + & exp( -(y - radius + dt)**2 / sigma**2 )* + & exp( -(z - radius + dt)**2 / sigma**2 ) + + phix_p_p(i,j,k) = - amplitude * (x - radius - 2*dt) / sigma**2 + & * exp( -(x - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (x - radius + 2*dt) / sigma**2 + & * exp( -(x - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z - radius + 2*dt)**2 / sigma**2 ) + + phiy_p(i,j,k) = - amplitude * (y - radius - dt) / sigma**2 + & * exp( -(x - radius - dt)**2 / sigma**2 )* + & exp( -(y - radius - dt)**2 / sigma**2 )* + & exp( -(z - radius - dt)**2 / sigma**2 ) + & - amplitude * (y - radius + dt) / sigma**2 + & * exp( -(x - radius + dt)**2 / sigma**2 )* + & exp( -(y - radius + dt)**2 / sigma**2 )* + & exp( -(z - radius + dt)**2 / sigma**2 ) + + phiy_p_p(i,j,k) = - amplitude * (y - radius - 2*dt) / sigma**2 + & * exp( -(x - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (y - radius + 2*dt) / sigma**2 + & * exp( -(x - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z - radius + 2*dt)**2 / sigma**2 ) + + phiz_p(i,j,k) = - amplitude * (z - radius - dt) / sigma**2 + & * exp( -(x - radius - dt)**2 / sigma**2 )* + & exp( -(y - radius - dt)**2 / sigma**2 )* + & exp( -(z - radius - dt)**2 / sigma**2 ) + & - amplitude * (z - radius + dt) / sigma**2 + & * exp( -(x - radius + dt)**2 / sigma**2 )* + & exp( -(y - radius + dt)**2 / sigma**2 )* + & exp( -(z - radius + dt)**2 / sigma**2 ) + + phiz_p_p(i,j,k) = - amplitude * (z - radius - 2*dt) / sigma**2 + & * exp( -(x - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z - radius - 2*dt)**2 / sigma**2 ) + & - amplitude * (z - radius + 2*dt) / sigma**2 + & * exp( -(x - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z - radius + 2*dt)**2 / sigma**2 ) + + phi_p(i,j,k) = amplitude / 2 * + & exp( -(x - radius - dt)**2 / sigma**2 )* + & exp( -(y - radius - dt)**2 / sigma**2 )* + & exp( -(z - radius - dt)**2 / sigma**2 ) + & + amplitude / 2 * + & exp( -(x - radius + dt)**2 / sigma**2 )* + & exp( -(y - radius + dt)**2 / sigma**2 )* + & exp( -(z - radius + dt)**2 / sigma**2 ) + + phi_p_p(i,j,k) = amplitude / 2 * + & exp( -(x - radius - 2*dt)**2 / sigma**2 )* + & exp( -(y - radius - 2*dt)**2 / sigma**2 )* + & exp( -(z - radius - 2*dt)**2 / sigma**2 ) + & + amplitude / 2 * + & exp( -(x - radius + 2*dt)**2 / sigma**2 )* + & exp( -(y - radius + 2*dt)**2 / sigma**2 )* + & exp( -(z - radius + 2*dt)**2 / sigma**2 ) + + resid_p(i,j,k) = 0.0 + resid_p_p(i,j,k) = 0.0 end do end do -- cgit v1.2.3