From e324b28ed17601a378e050d7a54867226851d981 Mon Sep 17 00:00:00 2001 From: allen Date: Sun, 14 Feb 1999 20:35:42 +0000 Subject: *** empty log message *** git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/CartGrid3D/trunk@2 c78560ca-4b45-4335-b268-5f3340f3cb52 --- README | 20 ++++++ interface.ccl | 23 +++++++ param.ccl | 52 +++++++++++++++ schedule.ccl | 11 +++ src/CartGrid3D.F | 173 +++++++++++++++++++++++++++++++++++++++++++++++ src/WaveToy.F | 192 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 9 +++ 7 files changed, 480 insertions(+) create mode 100644 README create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/CartGrid3D.F create mode 100644 src/WaveToy.F create mode 100644 src/make.code.defn diff --git a/README b/README new file mode 100644 index 0000000..625d155 --- /dev/null +++ b/README @@ -0,0 +1,20 @@ +Cactus Code Thorn CartGrid3D +Authors : ... +Managed by : ... <...@...........> +Version : ... +CVS info : $Header$ +-------------------------------------------------------------------------- + +1. Purpose of the thorn + +This thorn does ... + +2. Dependencies of the thorn + +This thorn additionally requires thorns ... + +3. Thorn distribution + +This thorn is available to ... + +4. Additional information diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..100bdb9 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,23 @@ +# Interface definition for thorn CartGrid3D +# $Header$ + +implements: grid +friend: driver + +public: + +real gridspacings type=SCALAR +{ + coarse_dx, + coarse_dy, + coarse_dz +} "3D Cartesian grid spacings" + +real coordinates type = GF +{ + x, + y, + z, + r +} "3D Cartesian grid coordinates" + diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..4832cfa --- /dev/null +++ b/param.ccl @@ -0,0 +1,52 @@ +# Parameter definitions for thorn CartGrid3D +# $Header$ + +public: + +REAL grid_dx "Grid spacing in x-direction" +{ + 0: :: +} 0.3 +REAL grid_dy "Grid spacing in y-direction" +{ + 0: :: +} 0.3 +REAL grid_dz "Grid spacing in z-direction" +{ + 0: :: +} 0.3 + +REAL grid_xmin "Coordinate minimum in x-direction" +{ + : :: +} -1.0 +REAL grid_ymin "Coordinate minimum in y-direction" +{ + : :: +} -1.0 +REAL grid_zmin "Coordinate minimum in z-direction" +{ + : :: +} -1.0 +REAL grid_xmax "Coordinate maximum in x-direction" +{ + : :: +} 1.0 +REAL grid_ymax "Coordinate maximum in y-direction" +{ + : :: +} 1.0 +REAL grid_zmax "Coordinate maximum in z-direction" +{ + : :: +} 1.0 + + +KEYWORD grid "Grid type" +{ + "box" :: "Box" + "octant" :: "Octant" + "full" :: "Full" + "minmax" :: "Minmax" +} "box" + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..964b37b --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,11 @@ +# Schedule definitions for thorn CartGrid3D +# $Header$ + +STORAGE: grid::coordinates + +schedule CartGrid3D at CACTUS_BASEGRID +{ + LANG:Fortran +} "Set up spatial 3D Cartesian coordinates on the GH" + + diff --git a/src/CartGrid3D.F b/src/CartGrid3D.F new file mode 100644 index 0000000..5a9332d --- /dev/null +++ b/src/CartGrid3D.F @@ -0,0 +1,173 @@ +#define THORN_IS_grid + +#include "cctk.h" +#include "declare_arguments.h" +#include "declare_parameters.h" + + subroutine CartGrid3D(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + + integer :: iflag, iconv, i, j, k + REAL :: x0,y0,z0, dx, dy, dz + + integer :: CCTK_Equals + + iflag = 0 + + iconv = 2**convlevel + + print *,"grid_xmin = ",grid_xmin + +c Defining all the parameters +c --------------------------- + grid_xmin = -1.0 + grid_xmax = 1.0 + grid_ymin = -1.0 + grid_ymax = 1.0 + grid_zmin = -1.0 + grid_zmax = 1.0 + +c First, grids that ignore dx,dy,dz +c --------------------------------- + +c if (CCTK_Equals(grid,'minmax')) then + + iflag = iflag +1 + +c Set minimum values of coordinates + x0 = grid_xmin + y0 = grid_ymin + z0 = grid_zmin + +c dx,dy,dz on the coarsest grid of each GH + coarse_dx = (grid_xmax-grid_xmin)/max(global_sh(1)-1,1) + coarse_dy = (grid_ymax-grid_ymin)/max(global_sh(2)-1,1) + coarse_dz = (grid_zmax-grid_zmin)/max(global_sh(3)-1,1) + +c dx,dy,dz on the grid we are on + dx = coarse_dx/levfac + dy = coarse_dy/levfac + dz = coarse_dz/levfac + + write(*,'(1X,A,1X,3(A,G12.7,2X))') + & 'Minmax grid: ','dx=>',dx,'dy=>',dy,'dz=>',dz + +c elseif (CCTK_Equals(grid,'box')) then + +c iflag = iflag + 1 + +c x0 = -0.5 +c y0 = -0.5 +c z0 = -0.5 +c dx = 1.d0/max(global_nx-1,1) +c dy = 1.d0/max(global_ny-1,1) +c dz = 1.d0/max(global_nz-1,1) +c if (global_nx == 1) x0 = 0.0D0 +c if (global_ny == 1) y0 = 0.0D0 +c if (global_nz == 1) z0 = 0.0D0 + +c write(*,'(1X,A,1X,3(A,G12.7,2X))') +c & 'Box grid: Resetting ','dx=>',dx,'dy=>',dy,'dz=>',dz + +c Now, grids that use dx,dy,dz +c ---------------------------- + +c else + +c dx = grid_dx*iconv +c dy = grid_dy*iconv +c dz = grid_dz*iconv + +c if (CCTK_Equals(grid,'bitant')) then + +c iflag = iflag + 1 +c x0 = (0.5 - global_nx/2)*dx +c y0 = (0.5 - global_ny/2)*dy + +c if(nghostzones == 2) then +c z0 = -1.5*dz +c else +c z0 = -0.5*dz +c end if + +c end if + +c if (CCTK_Equals(grid,'quadrant')) then + +c iflag = iflag + 1 + +c if(nghostzones == 2) then +c x0 = -1.5D0*dx +c y0 = -1.5D0*dy +c else +c x0 = -0.5D0*dx +c y0 = -0.5D0*dy +c end if + +c z0 = (0.5D0 - global_nz/2)*dz + +c end if + +c if (CCTK_Equals(grid,'octant')) then + +c if(nghostzones == 2) then + +c iflag = iflag + 1 +c x0 = -1.5*dx +c y0 = -1.5*dy +c z0 = -1.5*dz + +c else + +c iflag = iflag + 1 +c x0 = -0.5*dx +c y0 = -0.5*dy +c z0 = -0.5*dz + +c end if + +c end if + +c if (CCTK_Equals(grid,'full')) then + +c iflag = iflag + 1 +c x0 = (0.5 - global_nx/2)*dx +c y0 = (0.5 - global_ny/2)*dy +c z0 = (0.5 - global_nz/2)*dz + +c end if + +c end if + + +c No grid was set up +c ------------------ +c if (iflag.ne.1) then +c write(*,*) 'setbasegrid error: check your grid keyword' +c STOP +c endif + + +c Set spatial coordinates +c ----------------------- + do i=1,sh(1) + do j=1,sh(2) + do k=1,sh(3) + x(i,j,k) = dx*(i+lb(1)-1) + x0 + y(i,j,k) = dy*(j+lb(2)-1) + y0 + z(i,j,k) = dz*(k+lb(3)-1) + z0 + end do + end do + end do + + r = sqrt(x**2 + y**2 + z**2) + + return + end subroutine CartGrid3D + + + diff --git a/src/WaveToy.F b/src/WaveToy.F new file mode 100644 index 0000000..f23ec94 --- /dev/null +++ b/src/WaveToy.F @@ -0,0 +1,192 @@ +#include "cctk.h" +#include "declare_parameters.h" +#include "declare_arguments.h" + + subroutine WaveToy_Boundaries(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + + INTEGER CCTK_Equals + + WaveToy_bound = "none" + + if(.NOT.CCTK_Equals(WaveToy_Bound, "periodic") .and. + & .NOT.CCTK_Equals(WaveToy_Bound,"none")) then + scalar_field(1,:,:) = 0.0 + scalar_field(:,1,:) = 0.0 + scalar_field(:,:,1) = 0.0 + + scalar_field(sh(1),:,:) = 0.0 + scalar_field(:,sh(2),:) = 0.0 + scalar_field(:,:,sh(3)) = 0.0 + end if + +c Should have a call to symmetry bounds here + + end subroutine WaveToy_boundaries + + + + subroutine WaveToy_evolution(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + +c Declare local variables + INTEGER :: i,j,k + INTEGER :: istart, jstart, kstart, iend, jend, kend + REAL :: dx,dy,dz + + dx = coarse_dx/levfac + dy = coarse_dy/levfac + dz = coarse_dz/levfac + + istart = 2 + jstart = 2 + kstart = 2 + + iend = sh(1)-1 + jend = sh(2)-1 + kend = sh(3)-1 + + do k = kstart, kend + do j = jstart, jend + do i = istart, iend + + scalar_field_tmp(i,j,k) = + 1 2.0*(1.0 - (dt**2)*(1.0/dx**2 + + 2 1.0/dy**2 +1.0/dz**2))*scalar_field(i,j,k) - + 3 scalar_field_old(i,j,k) + + 4 (dt**2) * + 5 ((scalar_field(i+1,j,k)+scalar_field(i-1,j,k))/dx**2 + 6 +(scalar_field(i,j+1,k)+scalar_field(i,j-1,k))/dy**2 + 7 +(scalar_field(i,j,k+1)+scalar_field(i,j,k-1))/dz**2) + + end do + end do + end do + + scalar_field_old = scalar_field + scalar_field = scalar_field_tmp + + call WaveToy_boundaries(CCTK_FARGUMENTS) + +c call CCTK_SyncGroup("scalarfields") + + end subroutine WaveToy_evolution + + + + subroutine WaveToy_init(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_PARAMETERS + + REAL :: omega, pi, kx, ky, kz, radius, amplitude + + integer CCTK_Equals + + pi = 4.0*atan(1.0) + +c Set parameters manually for now +c ------------------------------- + WaveToy_dtfac = 0.5 + WaveToy_kx = 2 + WaveToy_ky = 2 + WaveToy_kz = 2 + WaveToy_radius = 0.25 + WaveToy_amplitude = 1.0 + WaveToy_initialdata = "plane" + +c Calculate timestep + dt = WaveToy_dtfac*coarse_dx/levfac + time = 0 + +c Shorthand + kx = WaveToy_kx + ky = WaveToy_ky + kz = WaveToy_kz + radius = WaveToy_radius + amplitude = WaveToy_amplitude + + call WaveToy_boundaries(CCTK_FARGUMENTS) + +c call CCTK_SyncGroup("scalarfields") + + omega= sqrt(kx**2+ky**2+kz**2) + +c Initialise the scalar field. Can be plane waves, spherical waves +c or waves in a box. + + if (CCTK_Equals(WaveToy_InitialData,"plane")) then + scalar_field = amplitude* + & cos(kx*x+ky*y+kz*z+omega*time) + scalar_field_old = amplitude* + & cos(kx*x+ky*y+kz*z+omega*(time-dt)) + + else if (CCTK_Equals(WaveToy_initialdata,"spherical")) then + + scalar_field = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*time) + scalar_field_old = amplitude*cos(kx*sqrt(x**2+y**2+z**2)+omega*(time-dt)) + + else if (CCTK_Equals(WaveToy_InitialData, "box")) then + +c Only currently implemented for the "box" grid type. + +c if (CCTK_Equals(grid, "box")) then + +c Must have all the wavenumbers non-zero for this to work. + +c if (kx == 0) then +c write(*,*) "Box mode selected. Setting kx to 1" +c kx = 1 +c end if + +c if (ky == 0) then +c write(*,*) "Box mode selected. Setting ky to 1" +c ky = 1 +c end if + +c if (kz == 0) then +c write(*,*) "Box mode selected. Setting kz to 1" +c kz = 1 +c end if + +c Set up the initial data. +c Use kx,ky,kz as number of modes in each direction. + +c scalar_field = amplitude*sin(kx*(x-0.5)*pi)* +c $ sin(ky*(y-0.5)*pi)* +c $ sin(kz*(z-0.5)*pi)* +c $ cos(omega*time*pi) + +c scalar_field_old = amplitude*sin(kx*(x-0.5)*pi)* +c $ sin(ky*(y-0.5)*pi)* +c $ sin(kz*(z-0.5)*pi)* +c $ cos(omega*(time-dt)*pi) +c else +c write(*,*) "WaveToy box mode selected, but grid not set to box" +c STOP +c end if + + else if (CCTK_Equals(WaveToy_InitialData,"tsgaussian")) then + + scalar_field = exp(-amplitude*(sqrt(x**2+y**2+z**2)-radius)**2) + scalar_field_old = scalar_field + + else + + print *,"Unknown initialisation mode in WaveToy" +c CCTK_Warning(GH,0,"Unknown initialisation mode in WaveToy") +c CCTK_Stop(GH,0) + + end if + + end subroutine WaveToy_Init diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..bf2ffff --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn CartGrid3D +# $Header$ + +# Source files in this directory +SRCS = CartGrid3D.F + +# Subdirectories containing source files +SUBDIRS = + -- cgit v1.2.3