diff options
| author | schnetter <schnetter@850bcc8b-0e4f-0410-8c26-8d28fbf1eda9> | 2004-05-14 13:39:52 +0000 |
|---|---|---|
| committer | schnetter <schnetter@850bcc8b-0e4f-0410-8c26-8d28fbf1eda9> | 2004-05-14 13:39:52 +0000 |
| commit | 3c4c7191cc89fc255f5ac66156eec06a84c77f9d (patch) | |
| tree | 1557b393d7056ebd93ebadbe1361eb515423ec4d /src/apply_dissipation.F77 | |
| parent | a5dfe864857e7252dab39257cec3a75565305455 (diff) | |
Fix a severe bug: forgot to divide by the grid spacing.
Remove spurious factor of 1/16.
Add second order dissipation, which is first order correct.
Checked stability interval. It should now be as described.
git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Dissipation/trunk@4 850bcc8b-0e4f-0410-8c26-8d28fbf1eda9
Diffstat (limited to 'src/apply_dissipation.F77')
| -rw-r--r-- | src/apply_dissipation.F77 | 47 |
1 files changed, 36 insertions, 11 deletions
diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77 index 9363824..ab78713 100644 --- a/src/apply_dissipation.F77 +++ b/src/apply_dissipation.F77 @@ -2,26 +2,51 @@ c $Header$ #include "cctk.h" - subroutine apply_dissipation (var, rhs, ni, nj, nk, epsdis) + subroutine apply_dissipation (var, rhs, ni, nj, nk, dx, order, epsdis) implicit none integer ni, nj, nk CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk) + CCTK_REAL dx(3) + CCTK_INT order CCTK_REAL epsdis integer i, j, k - do k = 3, nk-2 - do j = 3, nj-2 - do i = 3, ni-2 - - rhs(i,j,k) = rhs(i,j,k) - epsdis / 16 - $ * ( var(i-2,j,k) - 4*var(i-1,j,k) + 6*var(i,j,k) - 4*var(i+1,j,k) + var(i+2,j,k) - $ + var(i,j-2,k) - 4*var(i,j-1,k) + 6*var(i,j,k) - 4*var(i,j+1,k) + var(i,j+2,k) - $ + var(i,j,k-2) - 4*var(i,j,k-1) + 6*var(i,j,k) - 4*var(i,j,k+1) + var(i,j,k+2)) - + if (order .eq. 2) then + + do k = 2, nk-1 + do j = 2, nj-1 + do i = 2, ni-1 + + rhs(i,j,k) = rhs(i,j,k) - epsdis + $ * (+ (var(i-1,j,k) - 2*var(i,j,k) - var(i+1,j,k)) / dx(1) + $ + (var(i,j-1,k) - 2*var(i,j,k) - var(i,j+1,k)) / dx(2) + $ + (var(i,j,k-1) - 2*var(i,j,k) - var(i,j,k+1)) / dx(3)) + + end do end do end do - end do + + else if (order .eq. 4) then + + do k = 3, nk-2 + do j = 3, nj-2 + do i = 3, ni-2 + + rhs(i,j,k) = rhs(i,j,k) - epsdis + $ * (+ (var(i-2,j,k) - 4*var(i-1,j,k) + 6*var(i,j,k) - 4*var(i+1,j,k) + var(i+2,j,k)) / dx(1) + $ + (var(i,j-2,k) - 4*var(i,j-1,k) + 6*var(i,j,k) - 4*var(i,j+1,k) + var(i,j+2,k)) / dx(2) + $ + (var(i,j,k-2) - 4*var(i,j,k-1) + 6*var(i,j,k) - 4*var(i,j,k+1) + var(i,j,k+2)) / dx(3)) + + end do + end do + end do + + else + + call CCTK_WARN (0, "internal error") + + end if end |