From 4ea97d9b69c26ce524be1d403fbb960eb2e7528e Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 25 May 2005 01:24:48 +0000 Subject: Dissipation operator for the optimized 4-3 restricted full norm case. Testing still in progress, but it looks like the convergence order is right. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@33 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- param.ccl | 5 + src/Dissipation_4_3_min_err_coeff.F90 | 497 ++++++++++++++++++++++++++++++++++ src/dissipation.c | 16 ++ src/make.code.defn | 1 + 4 files changed, 519 insertions(+) create mode 100644 src/Dissipation_4_3_min_err_coeff.F90 diff --git a/param.ccl b/param.ccl index 0866282..ebbe4d4 100644 --- a/param.ccl +++ b/param.ccl @@ -27,6 +27,11 @@ REAL epsdis "Dissipation strength" STEERABLE=always *:* :: "Values typical between 0 and 1" } 0.2 +REAL diss_fraction "Fractional size of the transition region for the full restricted dissipation operator" +{ + 0:0.5 :: "" +} 0.2 + STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always { .* :: "Must be a valid list of grid functions" diff --git a/src/Dissipation_4_3_min_err_coeff.F90 b/src/Dissipation_4_3_min_err_coeff.F90 new file mode 100644 index 0000000..eefbecf --- /dev/null +++ b/src/Dissipation_4_3_min_err_coeff.F90 @@ -0,0 +1,497 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + + +subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, dfl, rhs) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_INT, dimension(3), intent(in) :: lsh, gsh, lbnd + CCTK_INT :: ni, nj, nk + CCTK_REAL, dimension(lsh(1),lsh(2),lsh(3)), intent(in) :: var + CCTK_REAL, dimension(lsh(1),lsh(2),lsh(3)), intent(inout) :: rhs + CCTK_INT, dimension(6), intent(in) :: bb + CCTK_INT, dimension(3), intent(in) :: gsize + CCTK_REAL, dimension(3), intent(in) :: delta + CCTK_REAL, intent(in) :: epsdisl, dfl + + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + integer :: i, j, k + CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + ni = lsh(1); nj = lsh(2); nk = lsh(3) + allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) ) + a = zero; d = zero; b = zero; h = zero; + + call set_dmatrix ( d, bb(1:2) ) + + call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl ) + + call set_hmatrix ( h, bb(1:2) ) + + a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + idel = epsdisl / 16.0_wp + + if ( bb(1) == 0 ) then + il = 1 + gsize(1) + else + rhs(1,:,:) = rhs(1,:,:) + & + ( a(1,1) * var(1,:,:) + a(2,1) * var(2,:,:) + & + a(3,1) * var(3,:,:) ) * idel + rhs(2,:,:) = rhs(2,:,:) + & + ( a(1,2) * var(1,:,:) + a(2,2) * var(2,:,:) + & + a(3,2) * var(3,:,:) + a(4,2) * var(4,:,:) + & + a(5,2) * var(5,:,:) + a(6,2) * var(6,:,:) + & + a(7,2) * var(7,:,:) ) * idel + rhs(3,:,:) = rhs(3,:,:) + & + ( a(1,3) * var(1,:,:) + a(2,3) * var(2,:,:) + & + a(3,3) * var(3,:,:) + a(4,3) * var(4,:,:) + & + a(5,3) * var(5,:,:) + a(6,3) * var(6,:,:) + & + a(7,3) * var(7,:,:) ) * idel + rhs(4,:,:) = rhs(4,:,:) + & + ( a(1,4) * var(1,:,:) + a(2,4) * var(2,:,:) + & + a(3,4) * var(3,:,:) + a(4,4) * var(4,:,:) + & + a(5,4) * var(5,:,:) + a(6,4) * var(6,:,:) + & + a(7,4) * var(7,:,:) ) * idel + rhs(5,:,:) = rhs(5,:,:) + & + ( a(1,5) * var(1,:,:) + a(2,5) * var(2,:,:) + & + a(3,5) * var(3,:,:) + a(4,5) * var(4,:,:) + & + a(5,5) * var(5,:,:) + a(6,5) * var(6,:,:) + & + a(7,5) * var(7,:,:) ) * idel + il = 6 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize(1) + else + rhs(ni-4,:,:) = rhs(ni-4,:,:) + & + ( a(ni-6,ni-4) * var(ni-6,:,:) + & + a(ni-5,ni-4) * var(ni-5,:,:) + & + a(ni-4,ni-4) * var(ni-4,:,:) + & + a(ni-3,ni-4) * var(ni-3,:,:) + & + a(ni-2,ni-4) * var(ni-2,:,:) + & + a(ni-1,ni-4) * var(ni-1,:,:) + & + a(ni,ni-4) * var(ni,:,:) ) * idel + rhs(ni-3,:,:) = rhs(ni-3,:,:) + & + ( a(ni-6,ni-3) * var(ni-6,:,:) + & + a(ni-5,ni-3) * var(ni-5,:,:) + & + a(ni-4,ni-3) * var(ni-4,:,:) + & + a(ni-3,ni-3) * var(ni-3,:,:) + & + a(ni-2,ni-3) * var(ni-2,:,:) + & + a(ni-1,ni-3) * var(ni-1,:,:) + & + a(ni,ni-3) * var(ni,:,:) ) * idel + rhs(ni-2,:,:) = rhs(ni-2,:,:) + & + ( a(ni-6,ni-2) * var(ni-6,:,:) + & + a(ni-5,ni-2) * var(ni-5,:,:) + & + a(ni-4,ni-2) * var(ni-4,:,:) + & + a(ni-3,ni-2) * var(ni-3,:,:) + & + a(ni-2,ni-2) * var(ni-2,:,:) + & + a(ni-1,ni-2) * var(ni-1,:,:) + & + a(ni,ni-2) * var(ni,:,:) ) * idel + rhs(ni-1,:,:) = rhs(ni-1,:,:) + & + ( a(ni-6,ni-1) * var(ni-6,:,:) + & + a(ni-5,ni-1) * var(ni-5,:,:) + & + a(ni-4,ni-1) * var(ni-4,:,:) + & + a(ni-3,ni-1) * var(ni-3,:,:) + & + a(ni-2,ni-1) * var(ni-2,:,:) + & + a(ni-1,ni-1) * var(ni-1,:,:) + & + a(ni,ni-1) * var(ni,:,:) ) * idel + rhs(ni,:,:) = rhs(ni,:,:) + & + ( a(ni-2,ni) * var(ni-2,:,:) + & + a(ni-1,ni) * var(ni-1,:,:) + & + a(ni,ni) * var(ni,:,:) ) * idel + ir = ni - 5 + end if + do i = il, ir + rhs(i,:,:) = rhs(i,:,:) + & + ( a(i-2,i) * var(i-2,:,:) + & + a(i-1,i) * var(i-1,:,:) + & + a(i,i) * var(i,:,:) + & + a(i+1,i) * var(i+1,:,:) + & + a(i+2,i) * var(i+2,:,:) ) * idel + end do + + deallocate ( a, d, b, h ) + + allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj) ) + a = zero; d = zero; b = zero; h = zero; + + call set_dmatrix ( d, bb(3:4) ) + + call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl ) + + call set_hmatrix ( h, bb(3:4) ) + + a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + idel = epsdisl / 16.0_wp + + if ( bb(3) == 0 ) then + jl = 1 + gsize(2) + else + rhs(:,1,:) = rhs(:,1,:) + & + ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + & + a(3,1) * var(:,3,:) ) * idel + rhs(:,2,:) = rhs(:,2,:) + & + ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + & + a(3,2) * var(:,3,:) + a(4,2) * var(:,4,:) + & + a(5,2) * var(:,5,:) + a(6,2) * var(:,6,:) + & + a(7,2) * var(:,7,:) ) * idel + rhs(:,3,:) = rhs(:,3,:) + & + ( a(1,3) * var(:,1,:) + a(2,3) * var(:,2,:) + & + a(3,3) * var(:,3,:) + a(4,3) * var(:,4,:) + & + a(5,3) * var(:,5,:) + a(6,3) * var(:,6,:) + & + a(7,3) * var(:,7,:) ) * idel + rhs(:,4,:) = rhs(:,4,:) + & + ( a(1,4) * var(:,1,:) + a(2,4) * var(:,2,:) + & + a(3,4) * var(:,3,:) + a(4,4) * var(:,4,:) + & + a(5,4) * var(:,5,:) + a(6,4) * var(:,6,:) + & + a(7,4) * var(:,7,:) ) * idel + rhs(:,5,:) = rhs(:,5,:) + & + ( a(1,5) * var(:,1,:) + a(2,5) * var(:,2,:) + & + a(3,5) * var(:,3,:) + a(4,5) * var(:,4,:) + & + a(5,5) * var(:,5,:) + a(6,5) * var(:,6,:) + & + a(7,5) * var(:,7,:) ) * idel + jl = 6 + end if + if ( bb(4) == 0 ) then + jr = nj - gsize(2) + else + rhs(:,nj-4,:) = rhs(:,nj-4,:) + & + ( a(nj-6,nj-4) * var(:,nj-6,:) + & + a(nj-5,nj-4) * var(:,nj-5,:) + & + a(nj-4,nj-4) * var(:,nj-4,:) + & + a(nj-3,nj-4) * var(:,nj-3,:) + & + a(nj-2,nj-4) * var(:,nj-2,:) + & + a(nj-1,nj-4) * var(:,nj-1,:) + & + a(nj,nj-4) * var(:,nj,:) ) * idel + rhs(:,nj-3,:) = rhs(:,nj-3,:) + & + ( a(nj-6,nj-3) * var(:,nj-6,:) + & + a(nj-5,nj-3) * var(:,nj-5,:) + & + a(nj-4,nj-3) * var(:,nj-4,:) + & + a(nj-3,nj-3) * var(:,nj-3,:) + & + a(nj-2,nj-3) * var(:,nj-2,:) + & + a(nj-1,nj-3) * var(:,nj-1,:) + & + a(nj,nj-3) * var(:,nj,:) ) * idel + rhs(:,nj-2,:) = rhs(:,nj-2,:) + & + ( a(nj-6,nj-2) * var(:,nj-6,:) + & + a(nj-5,nj-2) * var(:,nj-5,:) + & + a(nj-4,nj-2) * var(:,nj-4,:) + & + a(nj-3,nj-2) * var(:,nj-3,:) + & + a(nj-2,nj-2) * var(:,nj-2,:) + & + a(nj-1,nj-2) * var(:,nj-1,:) + & + a(nj,nj-2) * var(:,nj,:) ) * idel + rhs(:,nj-1,:) = rhs(:,nj-1,:) + & + ( a(nj-6,nj-1) * var(:,nj-6,:) + & + a(nj-5,nj-1) * var(:,nj-5,:) + & + a(nj-4,nj-1) * var(:,nj-4,:) + & + a(nj-3,nj-1) * var(:,nj-3,:) + & + a(nj-2,nj-1) * var(:,nj-2,:) + & + a(nj-1,nj-1) * var(:,nj-1,:) + & + a(nj,nj-1) * var(:,nj,:) ) * idel + rhs(:,nj,:) = rhs(:,nj,:) + & + ( a(nj-2,nj) * var(:,nj-2,:) + & + a(nj-1,nj) * var(:,nj-1,:) + & + a(nj,nj) * var(:,nj,:) ) * idel + jr = nj - 5 + end if + + do j = jl, jr + rhs(:,j,:) = rhs(:,j,:) + & + ( a(j-2,j) * var(:,j-2,:) + & + a(j-1,j) * var(:,j-1,:) + & + a(j,j) * var(:,j,:) + & + a(j+1,j) * var(:,j+1,:) + & + a(j+2,j) * var(:,j+2,:) ) * idel + end do + + deallocate ( a, d, b, h ) + + allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) ) + a = zero; d = zero; b = zero; h = zero; + + call set_dmatrix ( d, bb(5:6) ) + + call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl ) + + call set_hmatrix ( h, bb(5:6) ) + + a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + idel = epsdisl / 16.0_wp + + + if ( bb(5) == 0 ) then + kl = 1 + gsize(3) + else + rhs(:,:,1) = rhs(:,:,1) + & + ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) + & + a(3,1) * var(:,:,3) ) * idel + rhs(:,:,2) = rhs(:,:,2) + & + ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + & + a(3,2) * var(:,:,3) + a(4,2) * var(:,:,4) + & + a(5,2) * var(:,:,5) + a(6,2) * var(:,:,6) + & + a(7,2) * var(:,:,7) ) * idel + rhs(:,:,3) = rhs(:,:,3) + & + ( a(1,3) * var(:,:,1) + a(2,3) * var(:,:,2) + & + a(3,3) * var(:,:,3) + a(4,3) * var(:,:,4) + & + a(5,3) * var(:,:,5) + a(6,3) * var(:,:,6) + & + a(7,3) * var(:,:,7) ) * idel + rhs(:,:,4) = rhs(:,:,4) + & + ( a(1,4) * var(:,:,1) + a(2,4) * var(:,:,2) + & + a(3,4) * var(:,:,3) + a(4,4) * var(:,:,4) + & + a(5,4) * var(:,:,5) + a(6,4) * var(:,:,6) + & + a(7,4) * var(:,:,7) ) * idel + rhs(:,:,5) = rhs(:,:,5) + & + ( a(1,5) * var(:,:,1) + a(2,5) * var(:,:,2) + & + a(3,5) * var(:,:,3) + a(4,5) * var(:,:,4) + & + a(5,5) * var(:,:,5) + a(6,5) * var(:,:,6) + & + a(7,5) * var(:,:,7) ) * idel + kl = 6 + end if + if ( bb(6) == 0 ) then + kr = nk - gsize(3) + else + rhs(:,:,nk-4) = rhs(:,:,nk-4) + & + ( a(nk-6,nk-4) * var(:,:,nk-6) + & + a(nk-5,nk-4) * var(:,:,nk-5) + & + a(nk-4,nk-4) * var(:,:,nk-4) + & + a(nk-3,nk-4) * var(:,:,nk-3) + & + a(nk-2,nk-4) * var(:,:,nk-2) + & + a(nk-1,nk-4) * var(:,:,nk-1) + & + a(nk,nk-4) * var(:,:,nk) ) * idel + rhs(:,:,nk-3) = rhs(:,:,nk-3) + & + ( a(nk-6,nk-3) * var(:,:,nk-6) + & + a(nk-5,nk-3) * var(:,:,nk-5) + & + a(nk-4,nk-3) * var(:,:,nk-4) + & + a(nk-3,nk-3) * var(:,:,nk-3) + & + a(nk-2,nk-3) * var(:,:,nk-2) + & + a(nk-1,nk-3) * var(:,:,nk-1) + & + a(nk,nk-3) * var(:,:,nk) ) * idel + rhs(:,:,nk-2) = rhs(:,:,nk-2) + & + ( a(nk-6,nk-2) * var(:,:,nk-6) + & + a(nk-5,nk-2) * var(:,:,nk-5) + & + a(nk-4,nk-2) * var(:,:,nk-4) + & + a(nk-3,nk-2) * var(:,:,nk-3) + & + a(nk-2,nk-2) * var(:,:,nk-2) + & + a(nk-1,nk-2) * var(:,:,nk-1) + & + a(nk,nk-2) * var(:,:,nk) ) * idel + rhs(:,:,nk-1) = rhs(:,:,nk-1) + & + ( a(nk-6,nk-1) * var(:,:,nk-6) + & + a(nk-5,nk-1) * var(:,:,nk-5) + & + a(nk-4,nk-1) * var(:,:,nk-4) + & + a(nk-3,nk-1) * var(:,:,nk-3) + & + a(nk-2,nk-1) * var(:,:,nk-2) + & + a(nk-1,nk-1) * var(:,:,nk-1) + & + a(nk,nk-1) * var(:,:,nk) ) * idel + rhs(:,:,nk) = rhs(:,:,nk) + & + ( a(nk-2,nk) * var(:,:,nk-2) + & + a(nk-1,nk) * var(:,:,nk-1) + & + a(nk,nk) * var(:,:,nk) ) * idel + kr = nk - 5 + end if + + do k = kl, kr + rhs(:,:,k) = rhs(:,:,k) + & + ( a(k-2,k) * var(:,:,k-2) + & + a(k-1,k) * var(:,:,k-1) + & + a(k,k) * var(:,:,k) + & + a(k+1,k) * var(:,:,k+1) + & + a(k+2,k) * var(:,:,k+2) ) * idel + end do + + deallocate ( a, d, b, h ) + +contains + + subroutine set_dmatrix ( d, bb ) + + implicit none + + CCTK_REAL, dimension(:,:), intent(out) :: d + CCTK_INT, dimension(2), intent(in) :: bb + CCTK_INT :: n + CCTK_REAL, dimension(3), save :: ac = (/ 1.0_wp, -2.0_wp, 1.0_wp /) + CCTK_INT :: i + + n = size(d,1) + if ( bb(1) == 0 ) then + d(1,1:2) = ac(2:3) + else + d(1,1:3) = ac + end if + if ( bb(2) == 0 ) then + d(n,n-1:n) = ac(1:2) + else + d(n,n-2:n) = ac + end if + do i = 2, n-1 + d(i,i-1:i+1) = ac + end do + end subroutine set_dmatrix + + subroutine set_bmatrix ( b, bb, lsh, gsh, lbnd, h, dfl ) + + implicit none + + CCTK_REAL, dimension(:,:), intent(out) :: b + CCTK_INT, dimension(2), intent(in) :: bb + CCTK_INT, intent(in) :: lsh, gsh, lbnd + CCTK_REAL, intent(in) :: h, dfl + CCTK_INT :: n, nwt + CCTK_INT :: i + CCTK_REAL :: xp, xmax + + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + + n = size(d,1) + xp = ( gsh - 1 ) * dfl - 1 + do i = 1, n + b(i,i) = bf ( real(i+lbnd-1,wp), xp, real(gsh-1,wp), h ) + end do + + end subroutine set_bmatrix + + subroutine set_hmatrix ( h, bb ) + + implicit none + + CCTK_REAL, dimension(:,:), intent(out) :: h + CCTK_INT, dimension(2), intent(in) :: bb + CCTK_INT :: n, i + + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + + n = size(h,1) + do i = 1, n + h(i,i) = 1.0_wp + end do + if ( bb(1) /= 0 ) then + h(1:5,1:5) = sigma ( ) + end if + if ( bb(2) /= 0 ) then + h(n:n-4:-1,n:n-4:-1) = sigma ( ) + end if + + end subroutine set_hmatrix + + function bf ( x, xp, xmax, h ) + + implicit none + + CCTK_REAL :: bf + CCTK_REAL, intent(in) :: x, xp, xmax, h + + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + + if ( x < xp ) then + bf = ( 1.0_wp - h ) / xp * x + h + else if ( x > xmax - xp ) then + bf = ( h - 1.0_wp ) / xp * ( x + xp - xmax ) + 1 + else + bf = 1.0_wp + end if + + end function bf + + function sigma ( ) + + implicit none + + CCTK_REAL, dimension(5,5) :: sigma + + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + + sigma(1,1) = 4.186595370392226897362216859769846226369_wp + sigma(2,1) = 0_wp + sigma(3,1) = 0_wp + sigma(4,1) = 0_wp + sigma(5,1) = 0_wp + sigma(1,2) = 0_wp + sigma(2,2) = 0.6725191921225620731888714836983116420871_wp + sigma(3,2) = 0.3613418181134949259370502966736306984367_wp + sigma(4,2) = -0.2021316117293899791481674539631879662707_wp + sigma(5,2) = 0.03455320708729270824077678274955265350304_wp + sigma(1,3) = 0_wp + sigma(2,3) = 0.3613418181134949259370502966736306984367_wp + sigma(3,3) = 0.7206133711630147057720442098623847362950_wp + sigma(4,3) = 0.1376472340546569368321616389764958792591_wp + sigma(5,3) = -0.04136405531324488624637892257286207044784_wp + sigma(1,4) = 0_wp + sigma(2,4) = -0.2021316117293899791481674539631879662707_wp + sigma(3,4) = 0.1376472340546569368321616389764958792591_wp + sigma(4,4) = 0.9578653607931026822074133441449909689509_wp + sigma(5,4) = 0.02069353627247161734563597102894256809696_wp + sigma(1,5) = 0_wp + sigma(2,5) = 0.03455320708729270824077678274955265350304_wp + sigma(3,5) = -0.04136405531324488624637892257286207044784_wp + sigma(4,5) = 0.02069353627247161734563597102894256809696_wp + sigma(5,5) = 0.9908272703370861473007798925906968380654_wp + + end function sigma + + subroutine set_coeff ( a, dx ) + + implicit none + + CCTK_REAL, dimension(7,5), intent(OUT) :: a + CCTK_REAL, intent(in) :: dx + CCTK_REAL :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_REAL, dimension(4) :: b + + b(1) = dx + b(2) = ( 1 + 3*b(1) ) / 4 + b(3) = ( 1 + b(1) ) / 2 + b(4) = ( 3 + b(1) )/4 + + a(1,1) = -2.001729494283060564327579028280434295231079652820_wp*b(1)-2.001729494283060564327579028280434295231079652820_wp*b(2) + a(2,1) = 4.0034589885661211286551580565608685904621593056401_wp*b(1)+4.0034589885661211286551580565608685904621593056401_wp*b(2) + a(3,1) = -2.001729494283060564327579028280434295231079652820_wp*b(1)-2.001729494283060564327579028280434295231079652820_wp*b(2) + a(4,1) = 0.0_wp + a(5,1) = 0.0_wp + a(6,1) = 0.0_wp + a(7,1) = 0.0_wp + a(1,2) = 33.349009885955840248954469595024528895341083236112_wp*b(1)+33.349009885955840248954469595024528895341083236112_wp*b(2) + a(2,2) = -66.698019771911680497908939190049057790682166472223_wp*b(1)-66.698019771911680497908939190049057790682166472223_wp*b(2)-14.236744515267393430983451831106076148131079582685_wp*b(3) + a(3,2) = 33.349009885955840248954469595024528895341083236112_wp*b(1)+33.349009885955840248954469595024528895341083236112_wp*b(2)+28.473489030534786861966903662212152296262159165371_wp*b(3)-5.0440331538485788446043928053469767557227163958272_wp*b(4) + a(4,2) = -1.3951356569032137488502820136180622005437534045293_wp-14.236744515267393430983451831106076148131079582685_wp*b(3)+10.088066307697157689208785610693953511445432791654_wp*b(4) + a(5,2) = 5.9576955605539351261459121347762848951403545560302_wp-5.0440331538485788446043928053469767557227163958272_wp*b(4) + a(6,2) = -7.7299841503982290057409782286983831886494488984724_wp + a(7,2) = 3.1674242467475076284453481075401604940528477469716_wp + a(1,3) = -3.9903407976170909924931618645027940552738274990905_wp*b(1)-3.9903407976170909924931618645027940552738274990905_wp*b(2) + a(2,3) = 7.9806815952341819849863237290055881105476549981809_wp*b(1)+7.9806815952341819849863237290055881105476549981809_wp*b(2)+2.0281925060753259562659269077917853627407035645458_wp*b(3) + a(3,3) = -3.9903407976170909924931618645027940552738274990905_wp*b(1)-3.9903407976170909924931618645027940552738274990905_wp*b(2)-4.0563850121506519125318538155835707254814071290916_wp*b(3)+0.85749863949531149345224149788443686605331530370993_wp*b(4) + a(4,3) = -0.17805296968581960088505371131995522597503148483988_wp+2.0281925060753259562659269077917853627407035645458_wp*b(3)-1.7149972789906229869044829957688737321066306074199_wp*b(4) + a(5,3) = -0.029159711483604197600932897250942820943416950674515_wp+0.85749863949531149345224149788443686605331530370993_wp*b(4) + a(6,3) = 0.59247833202466719785702692846175131981192835586867_wp + a(7,3) = -0.38526565085524339937104031989085327289347992035428_wp + a(1,4) = -10.828009944579909415707962754026179346885875504376_wp*b(1)-10.828009944579909415707962754026179346885875504376_wp*b(2) + a(2,4) = 21.656019889159818831415925508052358693771751008752_wp*b(1)+21.656019889159818831415925508052358693771751008752_wp*b(2)+4.8030111910833762941985823652201696605573253856589_wp*b(3) + a(3,4) = -10.828009944579909415707962754026179346885875504376_wp*b(1)-10.828009944579909415707962754026179346885875504376_wp*b(2)-9.6060223821667525883971647304403393211146507713177_wp*b(3)+1.8413742894937777907392755981492717093567250921496_wp*b(4) + a(4,4) = 0.52169344312294565855548171683437104934836522915934_wp+4.8030111910833762941985823652201696605573253856589_wp*b(3)-3.6827485789875555814785511962985434187134501842992_wp*b(4) + a(5,4) = -2.3211687019578831529367545613086339586352127503243_wp+1.8413742894937777907392755981492717093567250921496_wp*b(4) + a(6,4) = 3.0772570745469293302070639721141547692253298131705_wp + a(7,4) = -1.2777818157119918358257911276398918599384822920056_wp + a(1,5) = -6.7201141443502586562617365349711742609991754142974_wp*b(1)-6.7201141443502586562617365349711742609991754142974_wp*b(2) + a(2,5) = 13.440228288700517312523473069942348521998350828595_wp*b(1)+13.440228288700517312523473069942348521998350828595_wp*b(2)+2.6601737327460025913616376196819751799013252956745_wp*b(3) + a(3,5) = -6.7201141443502586562617365349711742609991754142974_wp*b(1)-6.7201141443502586562617365349711742609991754142974_wp*b(2)-5.3203474654920051827232752393639503598026505913491_wp*b(3)+0.68697583481584777698277591803888262651185956794776_wp*b(4) + a(4,5) = 1.6888624757937931547697409070602037810047678994127_wp+2.6601737327460025913616376196819751799013252956745_wp*b(3)-1.3739516696316955539655518360777652530237191358955_wp*b(4) + a(5,5) = -4.8610470973404788048372478314704553824811608945344_wp+0.68697583481584777698277591803888262651185956794776_wp*b(4) + a(6,5) = 4.6555067672995781453652729417602994219480180908309_wp + a(7,5) = -1.4833221457528924952977660173500478204716250957091_wp + + end subroutine set_coeff +end subroutine dissipation_4_3_opt diff --git a/src/dissipation.c b/src/dissipation.c index 638a47a..59daa23 100644 --- a/src/dissipation.c +++ b/src/dissipation.c @@ -7,6 +7,16 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" +void CCTK_FCALL CCTK_FNAME(dissipation_4_3_opt) (const CCTK_REAL *var, + const CCTK_INT *lsh, + const CCTK_INT *gsh, + const CCTK_INT *lbnd, + const CCTK_INT *bbox, + const CCTK_INT *gsize, + const CCTK_REAL *dx, + const CCTK_REAL *epsdis, + const CCTK_REAL *diss_fraction, + CCTK_REAL *rhs); void CCTK_FCALL CCTK_FNAME(dissipation_6_5) (const CCTK_REAL *var, const CCTK_INT *ni, const CCTK_INT *nj, @@ -124,6 +134,12 @@ apply (int const varindex, char const * const optstring, void * const arg) } } else { switch(order) { + case 4: { + CCTK_FNAME(dissipation_4_3_opt) + (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, + bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr); + break; + } case 6: { CCTK_FNAME(dissipation_6_5) (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], diff --git a/src/make.code.defn b/src/make.code.defn index 040c89d..171c48f 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -17,6 +17,7 @@ SRCS = call_derivs.c \ set_norm_mask.F90 \ CheckGridSizes.F90 \ dissipation.c \ + Dissipation_4_3_min_err_coeff.F90 \ Dissipation_8_4.F90 \ Dissipation_6_5.F90 \ get_coeffs.c \ -- cgit v1.2.3