aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--param.ccl5
-rw-r--r--src/Dissipation_2_1.F9025
-rw-r--r--src/Dissipation_4_2.F9025
-rw-r--r--src/Dissipation_4_3_min_err_coeff.F9025
-rw-r--r--src/Dissipation_6_3.F9025
-rw-r--r--src/Dissipation_6_5_min_err_coeff.F9022
-rw-r--r--src/Dissipation_8_4.F9025
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)
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-02 22:07:41 +0000