diff options
-rw-r--r-- | CarpetExtra/IDScalarWaveMoL/param.ccl | 20 | ||||
-rw-r--r-- | CarpetExtra/IDScalarWaveMoL/src/initialdata.F77 | 27 |
2 files changed, 47 insertions, 0 deletions
diff --git a/CarpetExtra/IDScalarWaveMoL/param.ccl b/CarpetExtra/IDScalarWaveMoL/param.ccl index 548d55551..c06c04c20 100644 --- a/CarpetExtra/IDScalarWaveMoL/param.ccl +++ b/CarpetExtra/IDScalarWaveMoL/param.ccl @@ -4,6 +4,7 @@ KEYWORD initial_data "Initial data" { "plane wave" :: "Plane wave" "Gaussian pulse" :: "Gaussian pulse" + "Gaussian" :: "Gaussian" } "plane wave" @@ -36,6 +37,25 @@ CCTK_REAL pulse_offset[3] "Pulse offset" +# Gaussian specification + +REAL radius "Radius" +{ + *:* :: "" +} 0.0 + +REAL sigma "Width" +{ + (0:* :: "" +} 1.0 + +REAL origin[3] "Origin" +{ + *:* :: "" +} 0.0 + + + # Generic specifications CCTK_REAL time_offset "Time offset" diff --git a/CarpetExtra/IDScalarWaveMoL/src/initialdata.F77 b/CarpetExtra/IDScalarWaveMoL/src/initialdata.F77 index 17eed204b..dc5ee0369 100644 --- a/CarpetExtra/IDScalarWaveMoL/src/initialdata.F77 +++ b/CarpetExtra/IDScalarWaveMoL/src/initialdata.F77 @@ -43,5 +43,32 @@ 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) + tmp = sqrt( (x(i,j,k)-origin(1))**2 + $ + (y(i,j,k)-origin(2))**2 + $ + (z(i,j,k)-origin(3))**2) + if (tmp .gt. 1.0d-6 * sigma) then + phi(i,j,k) = amplitude/2 + $ * ( (cctk_time + tmp)/tmp * exp (-0.5d0 * ((cctk_time + tmp - radius) / sigma)**2) + $ - (cctk_time - tmp)/tmp * exp (-0.5d0 * ((cctk_time - tmp - radius) / sigma)**2)) + psi(i,j,k) = amplitude/2 + $ * (+ (1 + (cctk_time + tmp) * (cctk_time + tmp - radius) / sigma**2) / tmp * exp (-0.5d0 * ((cctk_time + tmp - radius) / sigma)**2) + $ - (1 + (cctk_time - tmp) * (cctk_time - tmp - radius) / sigma**2) / tmp * exp (-0.5d0 * ((cctk_time - tmp - radius) / sigma)**2)) + else + phi(i,j,k) = amplitude + $ * ( (radius * cctk_time - cctk_time**2 + sigma**2) / sigma**2 + $ + (radius**3 * cctk_time - cctk_time**4 + 6 * cctk_time**2 * sigma**2 - 3 * sigma**4 - 3 * radius**2 * (cctk_time - sigma)**2 + 3 * radius * (cctk_time**3 - 3 * cctk_time * sigma**2)) / (6 * sigma**6) * tmp**2) + $ * exp (-0.5d0 * (radius - cctk_time)**2 / sigma**2) + psi(i,j,k) = amplitude + $ * ( (radius**2 * cctk_time - 2 * radius * cctk_time**2 + cctk_time**3 + 2 * radius * sigma**2 - 3 * cctk_time * sigma**2) / sigma**4 + $ + (radius**4 * cctk_time + cctk_time**5 - 10 * cctk_time**3 * sigma**2 + 15 * cctk_time * sigma**4 - 4 * radius**4 * (cctk_time**2 - sigma**2) + 6 * radius**2 * (cctk_time**3 - 3 * cctk_time * sigma**2) - 4 * radius * (cctk_time**4 - 6 * cctk_time**2 * sigma**2 + 3 * sigma**4)) / (6 * sigma**8) * tmp**2) + $ * exp (-0.5d0 * (radius - cctk_time)**2 / sigma**2) + end if + end do + end do + end do end if end |