diff options
Diffstat (limited to 'src/InitialData.cc')
-rw-r--r-- | src/InitialData.cc | 104 |
1 files changed, 76 insertions, 28 deletions
diff --git a/src/InitialData.cc b/src/InitialData.cc index c8d790d..8b770dc 100644 --- a/src/InitialData.cc +++ b/src/InitialData.cc @@ -8,23 +8,24 @@ @enddesc @@*/ +#include <math.h> + #include "cctk.h" #include "cctk_Flesh.h" #include "cctk_parameters.h" #include "cctk_Groups.h" #include "cctk_arguments.h" -#include "cctk_Comm.h" +#include "cctk_Misc.h" -#include "CactusBase/CartGrid3D/src/Symmetry.h" +static char *rcsid = "$Header$"; -#include <math.h> - -inline double sqr(double val) +inline CCTK_REAL sqr(CCTK_REAL val) { - return val*val; + return val*val; } + /*@@ @routine IDScalarWave_InitialData @date @@ -35,39 +36,86 @@ inline double sqr(double val) @calls @calledby @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. @endhistory @@*/ extern "C" void IDScalarWaveCXX_InitialData(CCTK_CARGUMENTS) { - DECLARE_CCTK_CARGUMENTS - DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_CARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL dt = CCTK_DELTA_TIME; + + if(CCTK_Equals(initial_data, "plane")) + { + CCTK_REAL omega = sqrt(sqr(kx)+sqr(ky)+sqr(kz)); + + for(int k=0; k<cctk_lsh[2]; k++) + { + for(int j=0; j<cctk_lsh[1]; j++) + { + for(int i=0; i<cctk_lsh[0]; i++) + { + + int index = CCTK_GFINDEX3D(cctkGH,i,j,k); + + phi[index] = amplitude*cos(kx*x[index]+ky*y[index]+kz*z[index]+omega*cctk_time); + phi_old[index] = amplitude*cos(kx*x[index]+ky*y[index]+kz*z[index]+omega*(cctk_time-dt)); + } + } + } + } + else if(CCTK_Equals(initial_data, "gaussian")) + { + for(int k=0; k<cctk_lsh[2]; k++) + { + for(int j=0; j<cctk_lsh[1]; j++) + { + for(int i=0; i<cctk_lsh[0]; i++) + { + int index = CCTK_GFINDEX3D(cctkGH,i,j,k); - for(int k=0; k<cctk_lsh[2]; k++) - for(int j=0; j<cctk_lsh[1]; j++) - for(int i=0; i<cctk_lsh[0]; i++) - { - int index = CCTK_GFINDEX3D(cctkGH,i,j,k); + CCTK_REAL X = x[index], Y = y[index], Z = z[index]; + CCTK_REAL R = sqrt(X*X + Y*Y + Z*Z); - double X = x[index], Y = y[index], Z = z[index]; - double R = sqrt(X*X + Y*Y + Z*Z); + phi[index] = amplitude*exp( - sqr( (R - radius) / sigma ) ); + phi_old[index] = amplitude*exp( - sqr( (R - radius - dt) / sigma ) ); - phi_old[index] = phi[index] = - amplitude*exp( - sqr( (R - radius) / sigma ) ); - } + } + } + } + } + else if(CCTK_Equals(initial_data, "box")) + { + CCTK_REAL pi = 4.0*atan(1.0); + CCTK_REAL omega = sqrt(sqr(kx)+sqr(ky)+sqr(kz)); - // - // Apply symmetry boundary conditions - // - ApplySymmetry(cctkGH,"wavetoy::scalarevolve"); + for(int k=0; k<cctk_lsh[2]; k++) + { + for(int j=0; j<cctk_lsh[1]; j++) + { + for(int i=0; i<cctk_lsh[0]; i++) + { + int index = CCTK_GFINDEX3D(cctkGH,i,j,k); - // - // Synchronise - // - CCTK_SyncGroup(cctkGH,"wavetoy::scalarevolve"); + phi[index] = amplitude*sin(kx*(x[index]-0.5)*pi)* + sin(ky*(y[index]-0.5)*pi)* + sin(kz*(z[index]-0.5)*pi)* + cos(omega*cctk_time*pi); - return; + phi_old[index] = amplitude*sin(kx*(x[index]-0.5)*pi)* + sin(ky*(y[index]-0.5)*pi)* + sin(kz*(z[index]-0.5)*pi)* + cos(omega*(cctk_time-dt)*pi); + } + } + } + } + return; } |