From 74de5846f7ddeb3bf181b3205f6539b4f0081705 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 6 Dec 2005 17:20:11 +0000 Subject: Added a new boolean parameter "scale_with_h" that allows the user to choose to scale the dissipation with the grid spacing h. This results formally in the loss of one order of convergence. The default value is "no". git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@52 f69c4107-0314-4c4f-9ad4-17e986b73f4a --- param.ccl | 5 +++++ src/Dissipation_2_1.F90 | 25 ++++++++++++++++++++----- src/Dissipation_4_2.F90 | 25 ++++++++++++++++++++----- src/Dissipation_4_3_min_err_coeff.F90 | 25 +++++++++++++++++++------ src/Dissipation_6_3.F90 | 25 ++++++++++++++++++++----- src/Dissipation_6_5_min_err_coeff.F90 | 22 +++++++++++++++++----- src/Dissipation_8_4.F90 | 25 ++++++++++++++++++++----- 7 files changed, 121 insertions(+), 31 deletions(-) diff --git a/param.ccl b/param.ccl index d91cd2a..4e3fda5 100644 --- a/param.ccl +++ b/param.ccl @@ -22,6 +22,11 @@ BOOLEAN use_dissipation "Should we add dissipation" { } "no" +# Note: scaling the dissipation operators with h reduces the order by one. +BOOLEAN scale_with_h "Should we scale the dissipation with the grid spacing h" +{ +} "no" + REAL epsdis "Dissipation strength" STEERABLE=always { *:* :: "Values typical between 0 and 1" diff --git a/src/Dissipation_2_1.F90 b/src/Dissipation_2_1.F90 index 62ba030..0904b33 100644 --- a/src/Dissipation_2_1.F90 +++ b/src/Dissipation_2_1.F90 @@ -1,18 +1,21 @@ ! $Header$ #include "cctk.h" +#include "cctk_Parameters.h" -subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) +subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs) implicit none + DECLARE_CCTK_PARAMETERS + integer :: ni, nj, nk CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var CCTK_REAL, dimension(ni,nj,nk), 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) :: epsdis + CCTK_REAL, intent(in) :: epsilon CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -24,7 +27,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 4 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 4 * delta(1) ) + else + idel = epsilon / 4 + end if if ( bb(1) == 0 ) then il = 1 + gsize(1) @@ -55,7 +62,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 4 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 4 * delta(2) ) + else + idel = epsilon / 4 + end if if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -86,7 +97,11 @@ subroutine dissipation_2_1 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 4 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 4 * delta(3) ) + else + idel = epsilon / 4 + end if if ( bb(5) == 0 ) then kl = 1 + gsize(3) diff --git a/src/Dissipation_4_2.F90 b/src/Dissipation_4_2.F90 index 16a870d..f43073e 100644 --- a/src/Dissipation_4_2.F90 +++ b/src/Dissipation_4_2.F90 @@ -1,18 +1,21 @@ ! $Header$ #include "cctk.h" +#include "cctk_Parameters.h" -subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) +subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs) implicit none + DECLARE_CCTK_PARAMETERS + integer :: ni, nj, nk CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var CCTK_REAL, dimension(ni,nj,nk), 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) :: epsdis + CCTK_REAL, intent(in) :: epsilon CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -24,7 +27,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 16 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(1) ) + else + idel = epsilon / 16 + end if if ( bb(1) == 0 ) then il = 1 + gsize(1) @@ -75,7 +82,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 16 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(2) ) + else + idel = epsilon / 16 + end if if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -126,7 +137,11 @@ subroutine dissipation_4_2 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 16 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(3) ) + else + idel = epsilon / 16 + end if if ( bb(5) == 0 ) then kl = 1 + gsize(3) diff --git a/src/Dissipation_4_3_min_err_coeff.F90 b/src/Dissipation_4_3_min_err_coeff.F90 index b4b6ba3..dc4178d 100644 --- a/src/Dissipation_4_3_min_err_coeff.F90 +++ b/src/Dissipation_4_3_min_err_coeff.F90 @@ -4,7 +4,8 @@ #include "cctk_Parameters.h" -subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, dfl, rhs) +subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, & + delta, epsilon, dfl, rhs) implicit none @@ -17,7 +18,7 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, 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, intent(in) :: epsilon, dfl CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -28,6 +29,7 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, 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; @@ -39,7 +41,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - idel = epsdisl / 16.0_wp + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(1) ) + else + idel = epsilon / 16 + end if if ( bb(1) == 0 ) then il = 1 + gsize(1) @@ -132,7 +138,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - idel = epsdisl / 16.0_wp + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(2) ) + else + idel = epsilon / 16 + end if if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -226,8 +236,11 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, delta, epsdisl, a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - idel = epsdisl / 16.0_wp - + if ( scale_with_h > 0 ) then + idel = epsilon / ( 16 * delta(3) ) + else + idel = epsilon / 16 + end if if ( bb(5) == 0 ) then kl = 1 + gsize(3) diff --git a/src/Dissipation_6_3.F90 b/src/Dissipation_6_3.F90 index 4f35eec..1f0c74b 100644 --- a/src/Dissipation_6_3.F90 +++ b/src/Dissipation_6_3.F90 @@ -1,18 +1,21 @@ ! $Header$ #include "cctk.h" +#include "cctk_Parameters.h" -subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) +subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs) implicit none + DECLARE_CCTK_PARAMETERS + integer :: ni, nj, nk CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var CCTK_REAL, dimension(ni,nj,nk), 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) :: epsdis + CCTK_REAL, intent(in) :: epsilon CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -24,7 +27,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 64 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(1) ) + else + idel = epsilon / 64 + end if if ( bb(1) == 0 ) then il = 1 + gsize(1) @@ -101,7 +108,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 64 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(2) ) + else + idel = epsilon / 64 + end if if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -178,7 +189,11 @@ subroutine dissipation_6_3 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 64 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(3) ) + else + idel = epsilon / 64 + end if if ( bb(5) == 0 ) then kl = 1 + gsize(3) diff --git a/src/Dissipation_6_5_min_err_coeff.F90 b/src/Dissipation_6_5_min_err_coeff.F90 index a88d3ad..312ecef 100644 --- a/src/Dissipation_6_5_min_err_coeff.F90 +++ b/src/Dissipation_6_5_min_err_coeff.F90 @@ -6,7 +6,7 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & - delta, epsdisl, dfl, rhs) + delta, epsilon, dfl, rhs) implicit none @@ -20,7 +20,7 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & 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, intent(in) :: epsilon, dfl CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -79,7 +79,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & end if - idel = epsdisl / 64.0_wp + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(1) ) + else + idel = epsilon / 64 + end if if ( bb(1) == 0 ) then @@ -275,7 +279,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & end if - idel = epsdisl / 64.0_wp + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(2) ) + else + idel = epsilon / 64 + end if if ( bb(3) == 0 ) then @@ -471,7 +479,11 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & end if - idel = epsdisl / 64.0_wp + if ( scale_with_h > 0 ) then + idel = epsilon / ( 64 * delta(3) ) + else + idel = epsilon / 64 + end if if ( bb(5) == 0 ) then diff --git a/src/Dissipation_8_4.F90 b/src/Dissipation_8_4.F90 index bc225af..ec98407 100644 --- a/src/Dissipation_8_4.F90 +++ b/src/Dissipation_8_4.F90 @@ -1,18 +1,21 @@ ! $Header$ #include "cctk.h" +#include "cctk_Parameters.h" -subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) +subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs) implicit none + DECLARE_CCTK_PARAMETERS + integer :: ni, nj, nk CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var CCTK_REAL, dimension(ni,nj,nk), 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) :: epsdis + CCTK_REAL, intent(in) :: epsilon CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) @@ -23,7 +26,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 256 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 256 * delta(1) ) + else + idel = epsilon / 256 + end if if ( bb(1) == 0 ) then il = 1 + gsize(1) @@ -134,7 +141,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 256 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 256 * delta(2) ) + else + idel = epsilon / 256 + end if if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -245,7 +256,11 @@ subroutine dissipation_8_4 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs) call set_coeff ( a ) - idel = epsdis / 256 + if ( scale_with_h > 0 ) then + idel = epsilon / ( 256 * delta(3) ) + else + idel = epsilon / 256 + end if if ( bb(5) == 0 ) then kl = 1 + gsize(3) -- cgit v1.2.3