diff options
Diffstat (limited to 'src/InitialData.c')
-rw-r--r-- | src/InitialData.c | 94 |
1 files changed, 42 insertions, 52 deletions
diff --git a/src/InitialData.c b/src/InitialData.c index 2261b97..18cee00 100644 --- a/src/InitialData.c +++ b/src/InitialData.c @@ -1,17 +1,17 @@ /*@@ - @file InitialData.c - @date - @author Werner Benger - @desc - Initial data for the 3D Wave Equation - Derived from Tom Goodale - @enddesc - @version $Header$ + @file InitialData.c + @date + @author Werner Benger + @desc + Initial data for the 3D Wave Equation + Derived from Tom Goodale + @enddesc + @version $Id$ @@*/ #include <math.h> -#include "cctk.h" +#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" @@ -19,31 +19,26 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusWave_IDScalarWaveC_InitialData_c) -static CCTK_REAL sqr(CCTK_REAL val) -{ - return val*val; -} +#define SQR(val) ((val) * (val)) void IDScalarWaveC_InitialData(CCTK_ARGUMENTS); /*@@ @routine IDScalarWaveC_InitialData - @date + @date @author Tom Goodale - @desc + @desc Set up initial data for the wave equation - @enddesc - @calls - @calledby - @history + @enddesc + @history @hdate Mon Oct 11 11:48:03 1999 @hauthor Werner Benger @hdesc Converted to C++ @hdate Mon Oct 11 11:48:20 1999 @hauthor Tom Goodale @hdesc Added the rest of the initial data. @hdate Thu Feb 17 09:22:01 2000 @hauthor Tom Goodale @hdesc Converted to C - @endhistory + @endhistory @@*/ @@ -51,12 +46,12 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - + int i,j,k; CCTK_REAL dt; CCTK_REAL omega; - int index; + int vindex; CCTK_REAL X, Y, Z, R; CCTK_REAL pi; @@ -64,7 +59,7 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) if(CCTK_Equals(initial_data, "plane")) { - omega = sqrt(sqr(kx)+sqr(ky)+sqr(kz)); + omega = sqrt(SQR(kx)+SQR(ky)+SQR(kz)); for(k=0; k<cctk_lsh[2]; k++) { @@ -72,10 +67,10 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) { for(i=0; i<cctk_lsh[0]; i++) { - index = CCTK_GFINDEX3D(cctkGH,i,j,k); + vindex = CCTK_GFINDEX3D(cctkGH,i,j,k); - phi[index] = amplitude*cos(kx*x[index]+ky*y[index]+kz*z[index]+omega*cctk_time); - phi_p[index] = amplitude*cos(kx*x[index]+ky*y[index]+kz*z[index]+omega*(cctk_time-dt)); + phi[vindex] = amplitude*cos(kx*x[vindex]+ky*y[vindex]+kz*z[vindex]+omega*cctk_time); + phi_p[vindex] = amplitude*cos(kx*x[vindex]+ky*y[vindex]+kz*z[vindex]+omega*(cctk_time-dt)); } } } @@ -88,35 +83,35 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) { for(i=0; i<cctk_lsh[0]; i++) { - index = CCTK_GFINDEX3D(cctkGH,i,j,k); + vindex = CCTK_GFINDEX3D(cctkGH,i,j,k); - X = x[index]; - Y = y[index]; - Z = z[index]; + X = x[vindex]; + Y = y[vindex]; + Z = z[vindex]; R = sqrt(X*X + Y*Y + Z*Z); - phi[index] = amplitude*exp( - sqr( (R - radius) / sigma ) ); + phi[vindex] = amplitude*exp( - SQR( (R - radius) / sigma ) ); - if (R == 0.0) + if (R == 0.0) { - phi_p[index] = amplitude*(1.0 - 2.0*dt*dt/sigma)*exp(-dt*dt/sigma); + phi_p[vindex] = amplitude*(1.0 - 2.0*dt*dt/sigma)*exp(-dt*dt/sigma); } else { - phi_p[index] = amplitude/2.0*(R-dt)/R* - exp( - sqr( (R - radius - dt)/ sigma ) ) + phi_p[vindex] = amplitude/2.0*(R-dt)/R* + exp( - SQR( (R - radius - dt)/ sigma ) ) + amplitude/2.0*(R+dt)/R* - exp( - sqr( (R - radius + dt)/ sigma ) ); + exp( - SQR( (R - radius + dt)/ sigma ) ); } - } + } } } } else if(CCTK_Equals(initial_data, "box")) { pi = 4.0*atan(1.0); - omega = sqrt(sqr(kx)+sqr(ky)+sqr(kz)); + omega = sqrt(SQR(kx)+SQR(ky)+SQR(kz)); for(k=0; k<cctk_lsh[2]; k++) { @@ -124,16 +119,16 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) { for(i=0; i<cctk_lsh[0]; i++) { - index = CCTK_GFINDEX3D(cctkGH,i,j,k); + vindex = CCTK_GFINDEX3D(cctkGH,i,j,k); - phi[index] = amplitude*sin(kx*(x[index]-0.5)*pi)* - sin(ky*(y[index]-0.5)*pi)* - sin(kz*(z[index]-0.5)*pi)* + phi[vindex] = amplitude*sin(kx*(x[vindex]-0.5)*pi)* + sin(ky*(y[vindex]-0.5)*pi)* + sin(kz*(z[vindex]-0.5)*pi)* cos(omega*cctk_time*pi); - phi_p[index] = amplitude*sin(kx*(x[index]-0.5)*pi)* - sin(ky*(y[index]-0.5)*pi)* - sin(kz*(z[index]-0.5)*pi)* + phi_p[vindex] = amplitude*sin(kx*(x[vindex]-0.5)*pi)* + sin(ky*(y[vindex]-0.5)*pi)* + sin(kz*(z[vindex]-0.5)*pi)* cos(omega*(cctk_time-dt)*pi); } } @@ -147,16 +142,11 @@ void IDScalarWaveC_InitialData(CCTK_ARGUMENTS) { for(i=0; i<cctk_lsh[0]; i++) { - index = CCTK_GFINDEX3D(cctkGH,i,j,k); + vindex = CCTK_GFINDEX3D(cctkGH,i,j,k); - phi[index] = 0.0; - - phi_p[index] = 0.0; + phi[vindex] = phi_p[vindex] = 0.0; } } } } - - return; } - |