aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl96
-rw-r--r--param.ccl90
-rw-r--r--schedule.ccl145
-rw-r--r--src/NoExcision_mod.F90266
-rw-r--r--src/cg.F90634
-rw-r--r--src/getlevelinfo.cc45
-rw-r--r--src/make.code.defn2
-rw-r--r--src/make.code.deps1
-rw-r--r--src/overwrite.F90206
-rw-r--r--src/overwriteBSSN.F90248
-rw-r--r--src/reduce.F9097
-rw-r--r--src/reduce.c34
-rw-r--r--src/smooth.F90150
13 files changed, 1926 insertions, 88 deletions
diff --git a/interface.ccl b/interface.ccl
index 0aa4916..968f4cf 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -3,4 +3,98 @@
IMPLEMENTS: NoExcision
-INHERITS: ADMBase StaticConformal grid
+INHERITS: ADMBase StaticConformal grid SphericalSurface Boundary
+
+USES INCLUDE: Boundary.h
+
+USES INCLUDE: carpet.hh
+
+CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
+USES FUNCTION MoLQueryEvolvedRHS
+
+CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \
+ CCTK_INT IN faces, \
+ CCTK_INT IN boundary_width, \
+ CCTK_INT IN table_handle, \
+ CCTK_STRING IN var_name, \
+ CCTK_STRING IN bc_name)
+
+USES FUNCTION Boundary_SelectGroupForBC
+
+CCTK_INT smask type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+ nes_mask
+} "mask for smoothing"
+
+CCTK_REAL cg_res_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz
+} "Conjugate Gradient residual for the metric"
+
+CCTK_REAL cg_res_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz
+} "Conjugate Gradient residual for the extrinsic curvature"
+
+CCTK_REAL cg_res_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+ resx, resy, resz
+} "Conjugate Gradient residual for the shift"
+
+CCTK_REAL cg_res_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+ res
+} "Conjugate Gradient residual for the lapse"
+
+
+CCTK_REAL cg_d_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ dgxx, dgxy, dgxz, dgyy, dgyz, dgzz
+} "Conjugate Gradient d for the metric"
+
+CCTK_REAL cg_d_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz
+} "Conjugate Gradient d for the extrinsic curvature"
+
+CCTK_REAL cg_d_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+ dx, dy, dz
+} "Conjugate Gradient d for the shift"
+
+CCTK_REAL cg_d_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+ d
+} "Conjugate Gradient d for the lapse"
+
+
+CCTK_REAL cg_q_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ qgxx, qgxy, qgxz, qgyy, qgyz, qgzz
+} "Conjugate Gradient q for the metric"
+
+CCTK_REAL cg_q_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+ qkxx, qkxy, qkxz, qkyy, qkyz, qkzz
+} "Conjugate Gradient q for the extrinsic curvature"
+
+CCTK_REAL cg_q_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+ qx, qy, qz
+} "Conjugate Gradient q for the shift"
+
+CCTK_REAL cg_q_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+ q
+} "Conjugate Gradient q for the lapse"
+
+
+CCTK_REAL cg_red_all type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+ redgxx, redgxy, redgxz, redgyy, redgyz, redgzz
+ redkxx, redkxy, redkxz, redkyy, redkyz, redkzz
+ red, redx, redy, redz
+} "Conjugate Gradient red for all variables"
+
+# Control variable for smoothing loop
+CCTK_INT loop_control TYPE=SCALAR
diff --git a/param.ccl b/param.ccl
index 1fe7342..6d93bd3 100644
--- a/param.ccl
+++ b/param.ccl
@@ -5,13 +5,18 @@ BOOLEAN verbose "Produce some screen output"
{
} "no"
+KEYWORD method "Method to use"
+{
+ "old" :: "Use old method"
+ "new" :: "Use new method"
+} "old"
+
+
INT num_regions "Number of no-excision regions"
{
0:10 :: ""
} 0
-
-
BOOLEAN overwrite_geometry[10] "Set the geometry to Minkowski"
{
} "yes"
@@ -40,27 +45,100 @@ BOOLEAN overwrite_shift[10] "Set the shift to zero"
-REAL centre_x[10] "X-coordinate of the centre of the region"
+REAL centre_x[10] "x-coordinate of the centre of the region"
{
(*:*) :: ""
} 0.0
-REAL centre_y[10] "Y-coordinate of the centre of the region"
+REAL centre_y[10] "y-coordinate of the centre of the region"
{
(*:*) :: ""
} 0.0
-REAL centre_z[10] "Z-coordinate of the centre of the region"
+REAL centre_z[10] "z-coordinate of the centre of the region"
{
(*:*) :: ""
} 0.0
+KEYWORD region_shape[10] "Shape of the region"
+{
+ "sphere" :: "use radius"
+ "ellipsoid" :: "use radius_x, radius_y, and radius_z"
+ "surface" :: "use a spherical surface shape"
+} "sphere"
+
REAL radius[10] "Radius of the region"
{
0.0:*) :: ""
} 1.0
-REAL smoothing_zone_width[10] "Width of smoothing zone inside the region"
+REAL radius_x[10] "x-radius of the region"
+{
+ 0.0:*) :: ""
+} 1.0
+
+REAL radius_y[10] "y-radius of the region"
{
0.0:*) :: ""
+} 1.0
+
+REAL radius_z[10] "z-radius of the region"
+{
+ 0.0:*) :: ""
+} 1.0
+
+INT surface_index[10] "Spherical surface index"
+{
+ 0:* :: "must be an index of a spherical surface"
+} 0
+
+
+
+BOOLEAN reduce_rhs[10] "Reduce RHS"
+{
+} "no"
+
+REAL reduction_factor[10] "Reduction factor for RHS (0=complete, 1=no reduction)"
+{
+ *:* :: ""
+} 0.0
+
+
+
+KEYWORD smoothing_function[10] "Smoothing function"
+{
+ "linear" :: "linear ramp"
+ "spline" :: "cubic spline ramp"
+ "cosine" :: "cosine ramp"
+} "linear"
+
+REAL smoothing_zone_width[10] "Relative width of smoothing zone inside the region"
+{
+ 0.0:1.0 :: ""
} 0.0
+
+INT smoothing_iterations "Smoothing iterations"
+{
+ 0:* :: ""
+} 10
+
+REAL smoothing_factor "Initial moothing factor"
+{
+ (0:2) :: ""
+} 1.2
+
+
+
+BOOLEAN smooth_regions "Smooth overwritten regions?"
+{
+} "no"
+
+INT smoothing_order "Order of the derivatives used for CG smoothing"
+{
+ 2:6:2 :: ""
+} 6
+
+REAL smoothing_eps "CG smoothing stop criteria"
+{
+ (0.0:* :: ""
+} 1e-6
diff --git a/schedule.ccl b/schedule.ccl
index ac7b7b1..54d402b 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,7 +1,150 @@
# Schedule definitions for thorn NoExcision
# $Header$
-SCHEDULE NoExcision_Overwrite AT initial AFTER (ADMBase_InitialData ADMBase_InitialGauge)
+SCHEDULE NoExcision_Overwrite IN ADMBase_PostInitial
{
LANG: Fortran
} "Overwrite regions with Minkowski"
+
+if (CCTK_Equals(method,"new")) {
+ schedule NoExcision_SetSym at CCTK_BASEGRID
+ {
+ LANG: Fortran
+ } "Register the symmetries for the conjugate gradient functions."
+}
+
+if (smooth_regions) {
+ if (CCTK_Equals(method,"old")) {
+ SCHEDULE NoExcision_Smooth IN ADMBase_PostInitial AFTER NoExcision_Overwrite
+ {
+ LANG: Fortran
+ } "Smooth regions"
+ } else {
+
+ SCHEDULE GROUP NoExcision_CGSmoothing IN ADMBase_PostInitial
+ {
+ STORAGE: cg_res_metric, cg_res_curv, cg_res_shift, cg_res_lapse
+ STORAGE: cg_d_metric, cg_d_curv, cg_d_shift, cg_d_lapse
+ STORAGE: cg_q_metric, cg_q_curv, cg_q_shift, cg_q_lapse
+ STORAGE: cg_red_all
+ STORAGE: smask
+ STORAGE: loop_control
+
+ } "Conjugate Gradient smoothing"
+
+ SCHEDULE NoExcision_CGInit_1 IN NoExcision_CGSmoothing
+ {
+ LANG: Fortran
+ } "Initialise the conjugate gradient method 1"
+
+ SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p1 IN NoExcision_CGSmoothing AFTER NoExcision_CGInit_1
+ {
+ LANG: Fortran
+ SYNC: cg_d_metric
+ SYNC: cg_d_curv
+ SYNC: cg_d_shift
+ SYNC: cg_d_lapse
+ SYNC: cg_res_metric
+ SYNC: cg_res_curv
+ SYNC: cg_res_shift
+ SYNC: cg_res_lapse
+ SYNC: cg_red_all
+ } "Select variables for boundary conditions 1"
+
+ SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p1 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p1
+ {
+ } "Apply boundary conditions (Symmetries) 1"
+
+ SCHEDULE NoExcision_CGInit_2 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p1
+ {
+ LANG: Fortran
+ OPTIONS: level
+ } "Initialise the conjugate gradient method 2"
+
+ SCHEDULE GROUP NoExcision_Smoothing IN NoExcision_CGSmoothing AFTER NoExcision_CGInit_2 while NoExcision::loop_control
+ {
+ } "Smooth regions"
+
+ SCHEDULE NoExcision_CG_1 IN NoExcision_Smoothing
+ {
+ LANG: Fortran
+ } "Conjugate gradients step 1"
+
+ SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p2 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+ {
+ LANG: Fortran
+ SYNC: cg_q_metric
+ SYNC: cg_q_curv
+ SYNC: cg_q_shift
+ SYNC: cg_q_lapse
+ SYNC: cg_red_all
+ } "Select variables for boundary conditions 2"
+
+ SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p2 IN NoExcision_Smoothing AFTER NoExcision_CGApplySym_p2
+ {
+ } "Apply boundary conditions (Symmetries) 2"
+
+ SCHEDULE NoExcision_CG_2 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+ {
+ LANG: Fortran
+ OPTIONS: level
+ } "Conjugate gradients step 2"
+
+ SCHEDULE NoExcision_CG_3 IN NoExcision_Smoothing
+ {
+ LANG: Fortran
+ } "Conjugate gradients step 3"
+
+ SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p3 IN NoExcision_Smoothing AFTER NoExcision_CG_3
+ {
+ LANG: Fortran
+ SYNC: cg_res_metric
+ SYNC: cg_res_curv
+ SYNC: cg_res_shift
+ SYNC: cg_res_lapse
+ SYNC: metric
+ SYNC: curv
+ SYNC: shift
+ SYNC: lapse
+ SYNC: cg_red_all
+ } "Select variables for boundary conditions 3"
+
+ SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p3 IN NoExcision_Smoothing AFTER NoExcision_CGApplySym_p3
+ {
+ } "Apply boundary conditions (Symmetries) 3"
+
+ SCHEDULE NoExcision_CG_4 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+ {
+ LANG: Fortran
+ OPTIONS: level
+ } "Conjugate gradients step 4"
+
+ SCHEDULE NoExcision_CG_5 IN NoExcision_Smoothing
+ {
+ LANG: Fortran
+ } "Conjugate gradients step 5"
+
+ SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p4 IN NoExcision_CGSmoothing AFTER NoExcision_CG_5
+ {
+ LANG: Fortran
+ SYNC: cg_d_metric
+ SYNC: cg_d_curv
+ SYNC: cg_d_shift
+ SYNC: cg_d_lapse
+ } "Select variables for boundary conditions 4"
+
+ SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p4 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p4
+ {
+ } "Apply boundary conditions (Symmetries) 4"
+ }
+}
+
+#SCHEDULE NoExcision_OverwriteBSSN AT poststep
+#{
+# LANG: Fortran
+#} "Overwrite regions with Minkowski"
+
+SCHEDULE NoExcision_Reduce IN MoL_PostRHS
+{
+ LANG: C
+} "Reduce RHS"
diff --git a/src/NoExcision_mod.F90 b/src/NoExcision_mod.F90
new file mode 100644
index 0000000..646d335
--- /dev/null
+++ b/src/NoExcision_mod.F90
@@ -0,0 +1,266 @@
+! $Header$
+
+#include "cctk.h"
+
+module NoExcision_mod
+
+ implicit none
+
+ CCTK_REAL, parameter :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_REAL, parameter :: one = 1.0_wp
+ CCTK_REAL, parameter :: two = 2.0_wp
+ CCTK_INT :: infnorm_handle, sum_handle
+ CCTK_INT :: sym_selector
+ CCTK_INT :: loop_counter
+ CCTK_REAL, dimension(16) :: infnormresid
+ CCTK_REAL, dimension(16) :: delta_new, delta_old, alpha, beta, absres
+ character(len=18), dimension(16) :: red_names = (/ 'NoExcision::redgxx', &
+ 'NoExcision::redgxy', &
+ 'NoExcision::redgxz', &
+ 'NoExcision::redgyy', &
+ 'NoExcision::redgyz', &
+ 'NoExcision::redgzz', &
+ 'NoExcision::redkxx', &
+ 'NoExcision::redkxy', &
+ 'NoExcision::redkxz', &
+ 'NoExcision::redkyy', &
+ 'NoExcision::redkyz', &
+ 'NoExcision::redkzz', &
+ 'NoExcision::red', &
+ 'NoExcision::redx', &
+ 'NoExcision::redy', &
+ 'NoExcision::redz' /)
+ character(len=5), dimension(16) :: var_names = (/ ' gxx', ' gxy', &
+ ' gxz', ' gyy', &
+ ' gyz', ' gzz', &
+ ' kxx', ' kxy', &
+ ' kxz', ' kyy', &
+ ' kyz', ' kzz', &
+ ' alp', 'betax', &
+ 'betay', 'betaz' /)
+ integer :: nx, ny, nz
+ integer :: ii
+ logical, dimension(16) :: cont = (/ (.true., ii=1, 16) /)
+
+ contains
+
+ ! This is the basic residual calculation routine...
+
+ subroutine residual (val, mask, res, si, order )
+ CCTK_REAL, dimension(:,:,:), intent(in) :: val
+ CCTK_INT, dimension(:,:,:), intent(in) :: mask
+ CCTK_REAL, dimension(:,:,:), intent(out) :: res
+ CCTK_REAL, intent(in) :: si
+ CCTK_INT, intent(in) :: order
+ CCTK_REAL, parameter :: six = 6.0_wp
+ CCTK_REAL, parameter :: sixth = 1.0_wp / six
+ CCTK_REAL, parameter :: four = 4.0_wp
+ CCTK_REAL, parameter :: eighteen = 18.0_wp
+ CCTK_REAL, parameter :: eighteenth = 1.0_wp / eighteen
+ CCTK_REAL, parameter :: fifteen = 15.0_wp
+ CCTK_REAL, parameter :: sixty = 60.0_wp
+ CCTK_REAL, parameter :: sixtieth = 1.0_wp / sixty
+ integer :: i, j, k
+
+ select case (order)
+ case (2)
+
+ do k = 2, nz-1
+ do j = 2, ny-1
+ do i = 2, nx-1
+
+ if (mask(i,j,k) > 0) then
+
+ res(i,j,k) = val(i+1,j,k) + val(i-1,j,k) + val(i,j+1,k) + &
+ val(i,j-1,k) + val(i,j,k+1) + val(i,j,k-1) - &
+ six * val(i,j,k)
+ res(i,j,k) = sign(one,si) * res(i,j,k)
+
+ end if
+
+ end do
+ end do
+ end do
+
+ case (4)
+
+ do k = 3, nz-2
+ do j = 3, ny-2
+ do i = 3, nx-2
+
+ if (mask(i,j,k) > 0) then
+
+ res(i,j,k) = - ( val(i+2,j,k) + val(i-2,j,k) + &
+ val(i,j+2,k) + val(i,j-2,k) + &
+ val(i,j,k+2) + val(i,j,k-2) ) + &
+ four * ( val(i+1,j,k) + val(i-1,j,k) + &
+ val(i,j+1,k) + val(i,j-1,k) + &
+ val(i,j,k+1) + val(i,j,k-1) ) - &
+ eighteen * val(i,j,k)
+ res(i,j,k) = sign(one,si) * res(i,j,k)
+
+ end if
+
+ end do
+ end do
+ end do
+
+ case (6)
+
+ do k = 4, nz-3
+ do j = 4, ny-3
+ do i = 4, nx-3
+
+ if (mask(i,j,k) > 0) then
+
+ res(i,j,k) = ( val(i+3,j,k) + val(i-3,j,k) + &
+ val(i,j+3,k) + val(i,j-3,k) + &
+ val(i,j,k+3) + val(i,j,k-3) ) - &
+ six * ( val(i+2,j,k) + val(i-2,j,k) + &
+ val(i,j+2,k) + val(i,j-2,k) + &
+ val(i,j,k+2) + val(i,j,k-2) ) + &
+ fifteen * ( val(i+1,j,k) + val(i-1,j,k) + &
+ val(i,j+1,k) + val(i,j-1,k) + &
+ val(i,j,k+1) + val(i,j,k-1) ) - &
+ sixty * val(i,j,k)
+ res(i,j,k) = sign(one,si) * res(i,j,k)
+
+ end if
+
+ end do
+ end do
+ end do
+
+ case default
+
+ call CCTK_WARN ( 0, 'Internal error: Order out of range' )
+
+ end select
+
+ end subroutine residual
+
+
+ ! This is the routine calculating residuals for all variables, that
+ ! have not yet converged.
+
+ subroutine residual_all ( v1, v2, v3, v4, v5, v6, v7, v8, &
+ v9, v10, v11, v12, v13, v14, v15, v16, &
+ r1, r2, r3, r4, r5, r6, r7, r8, &
+ r9, r10, r11, r12, r13, r14, r15, r16, &
+ mask, si, order )
+
+ CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+ v7, v8, v9, v10, v11, v12, &
+ v13, v14, v15, v16
+ CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+ r7, r8, r9, r10, r11, r12, &
+ r13, r14, r15, r16
+ CCTK_INT, dimension(:,:,:), intent(in) :: mask
+ CCTK_REAL, intent(in) :: si
+ CCTK_INT, intent(in) :: order
+
+ if ( cont(1) ) call residual ( v1, mask, r1, si, order )
+ if ( cont(2) ) call residual ( v2, mask, r2, si, order )
+ if ( cont(3) ) call residual ( v3, mask, r3, si, order )
+ if ( cont(4) ) call residual ( v4, mask, r4, si, order )
+ if ( cont(5) ) call residual ( v5, mask, r5, si, order )
+ if ( cont(6) ) call residual ( v6, mask, r6, si, order )
+ if ( cont(7) ) call residual ( v7, mask, r7, si, order )
+ if ( cont(8) ) call residual ( v8, mask, r8, si, order )
+ if ( cont(9) ) call residual ( v9, mask, r9, si, order )
+ if ( cont(10) ) call residual ( v10, mask, r10, si, order )
+ if ( cont(11) ) call residual ( v11, mask, r11, si, order )
+ if ( cont(12) ) call residual ( v12, mask, r12, si, order )
+ if ( cont(13) ) call residual ( v13, mask, r13, si, order )
+ if ( cont(14) ) call residual ( v14, mask, r14, si, order )
+ if ( cont(15) ) call residual ( v15, mask, r15, si, order )
+ if ( cont(16) ) call residual ( v16, mask, r16, si, order )
+
+ end subroutine residual_all
+
+
+ ! This routine multiplies two variables. Only does it according
+ ! to the mask and only for variables that have not converged yet.
+
+ subroutine multiply ( u1, u2, u3, u4, u5, u6, u7, u8, &
+ u9, u10, u11, u12, u13, u14, u15, u16, &
+ v1, v2, v3, v4, v5, v6, v7, v8, &
+ v9, v10, v11, v12, v13, v14, v15, v16, &
+ r1, r2, r3, r4, r5, r6, r7, r8, &
+ r9, r10, r11, r12, r13, r14, r15, r16, &
+ mask )
+
+ CCTK_REAL, dimension(:,:,:), intent(in) :: u1, u2, u3, u4, u5, u6, &
+ u7, u8, u9, u10, u11, u12, &
+ u13, u14, u15, u16
+ CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+ v7, v8, v9, v10, v11, v12, &
+ v13, v14, v15, v16
+ CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+ r7, r8, r9, r10, r11, r12, &
+ r13, r14, r15, r16
+ CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+ if ( cont(1) ) where ( mask > 0 ) r1 = u1*v1
+ if ( cont(2) ) where ( mask > 0 ) r2 = u2*v2
+ if ( cont(3) ) where ( mask > 0 ) r3 = u3*v3
+ if ( cont(4) ) where ( mask > 0 ) r4 = u4*v4
+ if ( cont(5) ) where ( mask > 0 ) r5 = u5*v5
+ if ( cont(6) ) where ( mask > 0 ) r6 = u6*v6
+ if ( cont(7) ) where ( mask > 0 ) r7 = u7*v7
+ if ( cont(8) ) where ( mask > 0 ) r8 = u8*v8
+ if ( cont(9) ) where ( mask > 0 ) r9 = u9*v9
+ if ( cont(10) ) where ( mask > 0 ) r10 = u10*v10
+ if ( cont(11) ) where ( mask > 0 ) r11 = u11*v11
+ if ( cont(12) ) where ( mask > 0 ) r12 = u12*v12
+ if ( cont(13) ) where ( mask > 0 ) r13 = u13*v13
+ if ( cont(14) ) where ( mask > 0 ) r14 = u14*v14
+ if ( cont(15) ) where ( mask > 0 ) r15 = u15*v15
+ if ( cont(16) ) where ( mask > 0 ) r16 = u16*v16
+
+ end subroutine multiply
+
+
+ ! This routine takes the product of a scalar and a variable and
+ ! adds it to the product of another scalar and variable. Does it
+ ! according to the mask and only for variables that have not yet
+ ! converged.
+
+ subroutine multiply_sum ( u1, u2, u3, u4, u5, u6, u7, u8, &
+ u9, u10, u11, u12, u13, u14, u15, u16, &
+ v1, v2, v3, v4, v5, v6, v7, v8, &
+ v9, v10, v11, v12, v13, v14, v15, v16, &
+ c1, c2, mask )
+
+ CCTK_REAL, dimension(:,:,:), intent(inout) :: u1, u2, u3, u4, u5, u6, &
+ u7, u8, u9, u10, u11, u12, &
+ u13, u14, u15, u16
+ CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+ v7, v8, v9, v10, v11, v12, &
+ v13, v14, v15, v16
+ CCTK_REAL, dimension(16), intent(in) :: c1, c2
+ CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+
+ if ( cont(1) ) where ( mask > 0 ) u1 = c1(1)*u1 + c2(1)*v1
+ if ( cont(2) ) where ( mask > 0 ) u2 = c1(2)*u2 + c2(2)*v2
+ if ( cont(3) ) where ( mask > 0 ) u3 = c1(3)*u3 + c2(3)*v3
+ if ( cont(4) ) where ( mask > 0 ) u4 = c1(4)*u4 + c2(4)*v4
+ if ( cont(5) ) where ( mask > 0 ) u5 = c1(5)*u5 + c2(5)*v5
+ if ( cont(6) ) where ( mask > 0 ) u6 = c1(6)*u6 + c2(6)*v6
+ if ( cont(7) ) where ( mask > 0 ) u7 = c1(7)*u7 + c2(7)*v7
+ if ( cont(8) ) where ( mask > 0 ) u8 = c1(8)*u8 + c2(8)*v8
+ if ( cont(9) ) where ( mask > 0 ) u9 = c1(9)*u9 + c2(9)*v9
+ if ( cont(10) ) where ( mask > 0 ) u10 = c1(10)*u10 + c2(10)*v10
+ if ( cont(11) ) where ( mask > 0 ) u11 = c1(11)*u11 + c2(11)*v11
+ if ( cont(12) ) where ( mask > 0 ) u12 = c1(12)*u12 + c2(12)*v12
+ if ( cont(13) ) where ( mask > 0 ) u13 = c1(13)*u13 + c2(13)*v13
+ if ( cont(14) ) where ( mask > 0 ) u14 = c1(14)*u14 + c2(14)*v14
+ if ( cont(15) ) where ( mask > 0 ) u15 = c1(15)*u15 + c2(15)*v15
+ if ( cont(16) ) where ( mask > 0 ) u16 = c1(16)*u16 + c2(16)*v16
+
+ end subroutine multiply_sum
+
+
+end module NoExcision_mod
diff --git a/src/cg.F90 b/src/cg.F90
new file mode 100644
index 0000000..71f56c0
--- /dev/null
+++ b/src/cg.F90
@@ -0,0 +1,634 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_CGInit_1 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT :: my_level, n_levels
+ CCTK_INT :: ierr
+ CCTK_INT, dimension(3) :: sym
+
+ nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3)
+
+ loop_control = 0
+
+ call NoExcision_levelinfo ( cctkGH, my_level, n_levels )
+
+ if ( my_level == n_levels-1 ) then
+
+ where (abs(gxx) < 1d-8 .and. abs(alp) < 1d-8)
+ nes_mask = 1
+ elsewhere
+ nes_mask = 0
+ end where
+
+ call CCTK_INFO ( 'Starting smoothing procedure' )
+
+ resgxx = zero; resgxy = zero; resgxz = zero
+ resgyy = zero; resgyz = zero; resgzz = zero
+ reskxx = zero; reskxy = zero; reskxz = zero
+ reskyy = zero; reskyz = zero; reskzz = zero
+ res = zero; resx = zero; resy = zero; resz = zero
+
+ dgxx = zero; dgxy = zero; dgxz = zero
+ dgyy = zero; dgyz = zero; dgzz = zero
+ dkxx = zero; dkxy = zero; dkxz = zero
+ dkyy = zero; dkyz = zero; dkzz = zero
+ d = zero; dx = zero; dy = zero; dz = zero
+
+ qgxx = zero; qgxy = zero; qgxz = zero
+ qgyy = zero; qgyz = zero; qgzz = zero
+ qkxx = zero; qkxy = zero; qkxz = zero
+ qkyy = zero; qkyz = zero; qkzz = zero
+ q = zero; qx = zero; qy = zero; qz = zero
+
+ redgxx = zero; redgxy = zero; redgxz = zero
+ redgyy = zero; redgyz = zero; redgzz = zero
+ redkxx = zero; redkxy = zero; redkxz = zero
+ redkyy = zero; redkyz = zero; redkzz = zero
+ red = zero; redx = zero; redy = zero; redz = zero
+
+ ! r = b - A x.
+ ! Since x=0 and we actually use A': b = -A' 0 and r = b = -A' 0.
+
+ call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+ kxx, kxy, kxz, kyy, kyz, kzz, &
+ alp, betax, betay, betaz, &
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ nes_mask, -one, smoothing_order )
+
+ ! d = r = b.
+
+ call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+ kxx, kxy, kxz, kyy, kyz, kzz, &
+ alp, betax, betay, betaz, &
+ dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+ d, dx, dy, dz, &
+ nes_mask, -one, smoothing_order )
+
+ ! red = r*r.
+
+ call multiply ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+ redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+ red, redx, redy, redz, nes_mask )
+
+ call CCTK_ReductionArrayHandle ( sum_handle, 'sum' )
+ if ( sum_handle .lt. 0 ) then
+ call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
+ end if
+
+ call CCTK_ReductionArrayHandle ( infnorm_handle, 'norm_inf' )
+ if ( infnorm_handle .lt. 0 ) then
+ call CCTK_WARN(0,'Could not obtain a handle for norm_inf reduction')
+ end if
+
+ loop_counter = 0
+ loop_control = 1
+
+ sym_selector = 1
+
+ end if
+
+end subroutine NoExcision_CGInit_1
+
+
+subroutine NoExcision_CGInit_2 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT :: res_var, ierr
+ integer :: i
+ character(len=56) :: conv_message
+
+ if ( loop_control == 1 ) then
+ do i = 1, 16
+ if ( cont(i) ) then
+ call CCTK_VarIndex ( res_var, red_names(i) )
+ if ( res_var < 0 ) then
+ call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+ end if
+
+ ! delta_new = r^T r.
+
+ call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+ CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
+
+ call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
+ CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
+
+ ! Since we are reducing r*r, we have to take the square root.
+
+ infnormresid(i) = sqrt(infnormresid(i))
+
+ end if
+ end do
+
+ ! Check if some variables has already converged. This happens when the
+ ! variable is identically zero.
+
+ do i = 1, 16
+ if ( cont(i) ) then
+ if ( infnormresid(i) < smoothing_eps ) then
+ write ( conv_message, '(a23,i8,a20,a5)' ) 'CG method converged in ', &
+ loop_counter, ' steps for variable ', var_names(i)
+ call CCTK_INFO ( conv_message )
+ cont(i) = .false.
+ end if
+ end if
+ end do
+
+ if ( .not. any ( cont ) ) then
+ loop_control = 0
+ end if
+
+ end if
+
+end subroutine NoExcision_CGInit_2
+
+
+subroutine NoExcision_CG_1 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+ loop_counter = loop_counter + 1
+
+ ! Since d is zero outside of the active region we have:
+ ! q = A d = A' d.
+
+ call residual_all ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+ d, dx, dy, dz, &
+ qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+ qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+ q, qx, qy, qz, &
+ nes_mask, one, smoothing_order )
+
+ ! red = d*q = d*(A d)
+
+ call multiply ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+ d, dx, dy, dz, &
+ qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+ qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+ q, qx, qy, qz, &
+ redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+ redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+ red, redx, redy, redz, nes_mask )
+
+ sym_selector = 2
+
+end subroutine NoExcision_CG_1
+
+
+subroutine NoExcision_CG_2 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT res_var, ierr
+
+ integer :: i
+
+ if ( loop_control == 1 ) then
+ do i = 1, 16
+ if ( cont(i) ) then
+ call CCTK_VarIndex ( res_var, red_names(i) )
+ if ( res_var < 0 ) then
+ call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+ end if
+
+ ! alpha = delta_new / ( d^T A d ).
+ call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+ CCTK_VARIABLE_REAL, alpha(i), 1, res_var )
+
+ alpha(i) = delta_new(i) / alpha(i)
+
+ end if
+ end do
+
+ end if
+
+end subroutine NoExcision_CG_2
+
+
+subroutine NoExcision_CG_3 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+ ! x = x + alpha * d.
+
+ call multiply_sum ( gxx, gxy, gxz, gyy, gyz, gzz, &
+ kxx, kxy, kxz, kyy, kyz, kzz, &
+ alp, betax, betay, betaz, &
+ dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+ d, dx, dy, dz, &
+ (/ (one, i=1, 16) /), alpha, nes_mask )
+
+ if ( mod ( loop_counter, 50 ) == 0 ) then
+
+ ! Restart:
+ ! r = b - A x = - A' x
+
+ call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+ kxx, kxy, kxz, kyy, kyz, kzz, &
+ alp, betax, betay, betaz, &
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ nes_mask, -one, smoothing_order )
+ else
+
+ ! r = r - alpha q.
+
+ call multiply_sum ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+ qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+ q, qx, qy, qz, &
+ (/ (one, i=1, 16) /), -alpha, nes_mask )
+
+ end if
+
+ ! delta_old = delta_new.
+
+ where ( cont ) delta_old = delta_new
+
+ ! red = r*r.
+ call multiply ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+ redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+ red, redx, redy, redz, nes_mask )
+
+ sym_selector = 3
+
+end subroutine NoExcision_CG_3
+
+
+subroutine NoExcision_CG_4 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT res_var, ierr
+
+ integer :: i
+
+ if ( loop_control == 1 ) then
+ do i = 1, 16
+ if ( cont(i) ) then
+ call CCTK_VarIndex ( res_var, red_names(i) )
+ if ( res_var < 0 ) then
+ call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+ end if
+
+ ! delta_new = r^T r.
+ call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+ CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
+
+ call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
+ CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
+ infnormresid(i) = sqrt(infnormresid(i))
+ end if
+ end do
+
+ end if
+
+end subroutine NoExcision_CG_4
+
+
+subroutine NoExcision_CG_5 (CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: i
+ character(len=192) :: res_message, var_message
+ character(len=56) :: conv_message
+ character(len=18) :: iter_message
+
+ nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+ ! beta = delta_new / delta_old
+
+ where ( cont ) beta = delta_new / delta_old
+
+ ! d = r + beta d
+
+ call multiply_sum ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+ dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+ d, dx, dy, dz, &
+ resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+ reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+ res, resx, resy, resz, &
+ beta, (/ (one, i=1, 16) /), nes_mask )
+
+ if ( verbose > 0 .and. mod(loop_counter,100) == 0 ) then
+ write (var_message, '(16(a3,a9))' ) (' | ', var_names(i), i=1, 16)
+ call CCTK_INFO ( var_message )
+ end if
+ if ( verbose > 0 .and. mod(loop_counter,10) == 0 ) then
+ write (iter_message, '(a10,i8)' ), 'Iteration ', loop_counter
+ write (res_message, '(16(a3,es9.3))' ) (' | ', infnormresid(i), i=1, 16 )
+ call CCTK_INFO ( iter_message )
+ call CCTK_INFO ( res_message )
+ end if
+
+ ! Check if any variables have converged during this iteration.
+
+ do i = 1, 16
+ if ( cont(i) ) then
+ if ( infnormresid(i) < smoothing_eps ) then
+ write ( conv_message, '(a23,i8,a20,a5)' ) 'CG method converged in ', &
+ loop_counter, ' steps for variable ', var_names(i)
+ call CCTK_INFO ( conv_message )
+ cont(i) = .false.
+ end if
+ end if
+ end do
+
+ ! If all variables have converged we exit.
+
+ if ( .not. any ( cont ) ) then
+ loop_control = 0
+ end if
+
+ sym_selector = 4
+
+end subroutine NoExcision_CG_5
+
+
+! I leave this in here in case somebody uses CartGrid3D symmetries...
+! Not tested, though...
+
+subroutine NoExcision_SetSym(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+
+ CCTK_INT :: ierr
+ CCTK_INT, dimension(3) :: sym
+
+ sym = 1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkyy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkzz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::res' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::d' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::q' )
+
+ call SetCartSymGN ( ierr, cctkGH, sym, 'noexcision::cg_red_all' )
+
+ sym(1) = -1; sym(2) = -1; sym(3) = 1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxy' )
+
+ sym(1) = -1; sym(2) = 1; sym(3) = -1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxz' )
+
+ sym(1) = 1; sym(2) = -1; sym(3) = -1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgyz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskyz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgyz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkyz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgyz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkyz' )
+
+ sym(1) = -1; sym(2) = 1; sym(3) = 1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dx' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qx' )
+
+ sym(1) = 1; sym(2) = -1; sym(3) = 1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dy' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qy' )
+
+ sym(1) = 1; sym(2) = 1; sym(3) = -1
+
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dz' )
+ call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qz' )
+
+end subroutine NoExcision_SetSym
+
+
+subroutine NoExcision_CGApplySym(CCTK_ARGUMENTS)
+
+ use NoExcision_mod
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
+
+ CCTK_INT :: ierr
+
+ if ( loop_control > 0 ) then
+
+ if ( sym_selector <1 .or. sym_selector > 4 ) then
+ call CCTK_WARN ( 0, 'Internal error. Inconsistent symmetry selector' )
+ end if
+
+ select case ( sym_selector )
+ case (1, 4)
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_d_lapse', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_d_lapse for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_d_shift', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_d_shift for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_d_curv', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_d_curv for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_d_metric', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_d_metric for boundary condition' )
+ end if
+ end select
+
+ select case ( sym_selector )
+ case (1, 3)
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_res_lapse', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_res_lapse for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_res_shift', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_res_shift for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_res_curv', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_res_curv for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_res_metric', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_res_metric for boundary condition' )
+ end if
+ end select
+
+ select case ( sym_selector )
+ case (1, 2, 3)
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_red_all', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_red_all for boundary condition' )
+ end if
+ end select
+
+ select case ( sym_selector )
+ case (2)
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_q_lapse', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_q_lapse for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_q_shift', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_q_shift for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_q_curv', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_q_curv for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'NoExcision::cg_q_metric', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select cg_q_metric for boundary condition' )
+ end if
+ end select
+
+ select case ( sym_selector )
+ case (3)
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'ADMBase::lapse', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select lapse for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'ADMBase::shift', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select shift for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'ADMBase::curv', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select curv for boundary condition' )
+ end if
+
+ ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+ 'ADMBase::metric', 'None' )
+ if ( ierr /= 0 ) then
+ call CCTK_WARN ( 0, 'Could not select metric for boundary condition' )
+ end if
+ end select
+
+ end if
+
+end subroutine NoExcision_CGApplySym
diff --git a/src/getlevelinfo.cc b/src/getlevelinfo.cc
new file mode 100644
index 0000000..4b68abb
--- /dev/null
+++ b/src/getlevelinfo.cc
@@ -0,0 +1,45 @@
+/* $Header$ */
+
+#include <cstdio>
+#include <cmath>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_DefineThorn.h"
+
+#include "carpet.hh"
+
+namespace NoExcision {
+
+#ifdef HAVE_CARPET
+ using namespace Carpet;
+#endif
+
+ extern "C"
+ void NoExcision_levelinfo(const cGH * cctkGH,
+ CCTK_INT * my_level,
+ CCTK_INT * n_levels)
+ {
+
+ *my_level = 0;
+ *n_levels = 1;
+
+#ifdef HAVE_CARPET
+
+ if ( CCTK_IsThornActive ( "Carpet" ) ) {
+ *my_level = reflevel;
+ *n_levels = reflevels;
+ }
+
+#endif
+
+ }
+
+ extern "C"
+ CCTK_FCALL void CCTK_FNAME (NoExcision_levelinfo) (CCTK_POINTER_TO_CONST * cctkGH, CCTK_INT * my_level, CCTK_INT * n_levels)
+ {
+ NoExcision_levelinfo((const cGH *)* cctkGH, my_level, n_levels);
+ }
+
+} // namespace NoExcision
diff --git a/src/make.code.defn b/src/make.code.defn
index 7dd625c..a174c46 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
# $Header$
# Source files in this directory
-SRCS = overwrite.F90
+SRCS = NoExcision_mod.F90 overwrite.F90 smooth.F90 reduce.F90 reduce.c getlevelinfo.cc cg.F90
# Subdirectories containing source files
SUBDIRS =
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100644
index 0000000..e7f8124
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1 @@
+cg.F90.o: NoExcision_mod.F90.o
diff --git a/src/overwrite.F90 b/src/overwrite.F90
index 6c6089f..3df9ffb 100644
--- a/src/overwrite.F90
+++ b/src/overwrite.F90
@@ -10,85 +10,114 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
-
+
CCTK_REAL, parameter :: zero=0
CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
- CCTK_REAL :: cx, cy, cz, rad, width
+ CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
integer :: n
+
+ integer, parameter :: sm_linear = 1
+ integer, parameter :: sm_spline = 2
+ integer, parameter :: sm_cosine = 3
+ integer :: sm_type
+
+ if (CCTK_EQUALS(method,"old")) then
+ do n = 1, num_regions
+
+ cx = centre_x(n)
+ cy = centre_y(n)
+ cz = centre_z(n)
+ if (CCTK_EQUALS(region_shape(n), "sphere")) then
+ radx = radius(n)
+ rady = radius(n)
+ radz = radius(n)
+ else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+ radx = radius_x(n)
+ rady = radius_y(n)
+ radz = radius_z(n)
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+
+ if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+ sm_type = sm_linear
+ else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+ sm_type = sm_spline
+ else if (CCTK_EQUALS (smoothing_function(n), "cosine")) then
+ sm_type = sm_cosine
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+ width = smoothing_zone_width(n)
+
+ dist2 = ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+ & + ((z - cz) / radz)**2
+
+ if (overwrite_geometry(n) /= 0) then
+
+ if (conformal_state >= 1) then
+ where (dist2 <= 1)
+ psi = 1
+ end where
+ end if
+ if (conformal_state >= 2) then
+ where (dist2 <= 1)
+ psix = 0
+ psiy = 0
+ psiz = 0
+ end where
+ end if
+ if (conformal_state >= 3) then
+ where (dist2 <= 1)
+ psixx = 0
+ psixy = 0
+ psixz = 0
+ psiyy = 0
+ psiyz = 0
+ psizz = 0
+ end where
+ end if
- do n = 1, num_regions
-
- cx = centre_x(n)
- cy = centre_y(n)
- cz = centre_z(n)
- rad = radius(n)
- width = smoothing_zone_width(n)
-
- dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
-
- if (overwrite_geometry(n) /= 0) then
-
- if (conformal_state >= 1) then
- where (dist2 <= rad**2)
- psi = 1
- end where
- end if
- if (conformal_state >= 2) then
- where (dist2 <= rad**2)
- psix = 0
- psiy = 0
- psiz = 0
- end where
- end if
- if (conformal_state >= 3) then
- where (dist2 <= rad**2)
- psixx = 0
- psixy = 0
- psixz = 0
- psiyy = 0
- psiyz = 0
- psizz = 0
- end where
- end if
-
- where (dist2 <= rad**2)
- gxx = smooth (gxx, Minkowski_scale(n), dist2)
- gxy = smooth (gxy, zero , dist2)
- gxz = smooth (gxz, zero , dist2)
- gyy = smooth (gyy, Minkowski_scale(n), dist2)
- gyz = smooth (gyz, zero , dist2)
- gzz = smooth (gzz, Minkowski_scale(n), dist2)
- kxx = smooth (kxx, zero , dist2)
- kxy = smooth (kxy, zero , dist2)
- kxz = smooth (kxz, zero , dist2)
- kyy = smooth (kyy, zero , dist2)
- kyz = smooth (kyz, zero , dist2)
- kzz = smooth (kzz, zero , dist2)
- end where
-
- end if
-
- if (overwrite_lapse(n) /= 0) then
-
- where (dist2 <= rad**2)
- alp = smooth (alp, lapse_scale(n), dist2)
- end where
-
- end if
-
- if (overwrite_shift(n) /= 0) then
-
- if (shift_state /= 0) then
- where (dist2 <= rad**2)
- betax = smooth (betax, zero, dist2)
- betay = smooth (betay, zero, dist2)
- betaz = smooth (betaz, zero, dist2)
- end where
- end if
-
- end if
-
- end do
+ where (dist2 <= 1)
+ gxx = smooth (gxx, Minkowski_scale(n), dist2)
+ gxy = smooth (gxy, zero , dist2)
+ gxz = smooth (gxz, zero , dist2)
+ gyy = smooth (gyy, Minkowski_scale(n), dist2)
+ gyz = smooth (gyz, zero , dist2)
+ gzz = smooth (gzz, Minkowski_scale(n), dist2)
+ kxx = smooth (kxx, zero , dist2)
+ kxy = smooth (kxy, zero , dist2)
+ kxz = smooth (kxz, zero , dist2)
+ kyy = smooth (kyy, zero , dist2)
+ kyz = smooth (kyz, zero , dist2)
+ kzz = smooth (kzz, zero , dist2)
+ end where
+
+ end if
+
+ if (overwrite_lapse(n) /= 0) then
+
+ where (dist2 <= 1)
+ alp = smooth (alp, lapse_scale(n), dist2)
+ end where
+
+ end if
+
+ if (overwrite_shift(n) /= 0) then
+
+ if (shift_state /= 0) then
+ where (dist2 <= 1)
+ betax = smooth (betax, zero, dist2)
+ betay = smooth (betay, zero, dist2)
+ betaz = smooth (betaz, zero, dist2)
+ end where
+ end if
+
+ end if
+
+ end do
+
+ end if
contains
@@ -98,13 +127,32 @@ contains
CCTK_REAL, intent(in) :: goodval
CCTK_REAL, intent(in) :: dist2
- if (dist2 <= (rad-width)**2) then
+ CCTK_REAL, parameter :: zero=0, two=2, half=1/two
+ integer, parameter :: rk = kind(zero)
+ CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk
+ CCTK_REAL :: dist
+ CCTK_REAL :: x, f
+
+ if (dist2 <= (1 - width)**2) then
newval = goodval
- else if (width == 0 .or. dist2 >= rad**2) then
+ else if (width == 0 .or. dist2 >= 1) then
newval = oldval
else
- newval = goodval &
- & + (oldval - goodval) * (sqrt(dist2) - (rad-width)) / width
+ dist = sqrt(dist2)
+ x = 1 - (1 - dist) / width
+ select case (sm_type)
+ case (sm_linear)
+ f = x
+ case (sm_spline)
+ if (x < half) then
+ f = 2 * x**2
+ else
+ f = 1 - 2 * (1-x)**2
+ end if
+ case (sm_cosine)
+ f = (1 - cos (pi * x)) / 2
+ end select
+ newval = oldval * f + goodval * (1 - f)
end if
end function smooth
diff --git a/src/overwriteBSSN.F90 b/src/overwriteBSSN.F90
new file mode 100644
index 0000000..cef380a
--- /dev/null
+++ b/src/overwriteBSSN.F90
@@ -0,0 +1,248 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_OverwriteBSSN (CCTK_ARGUMENTS)
+ implicit none
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL, parameter :: zero=0
+ CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+ CCTK_REAL :: cx, cy, cz, radx, rady, radz, width, scale
+ integer :: timelevels
+ integer :: n
+
+ integer, parameter :: sm_linear = 1
+ integer, parameter :: sm_spline = 2
+ integer :: sm_type
+
+ do n = 1, num_regions
+
+ if (cctk_iteration == iteration(n)) then
+
+ cx = centre_x(n)
+ cy = centre_y(n)
+ cz = centre_z(n)
+ if (CCTK_EQUALS(region_shape(n), "sphere")) then
+ radx = radius(n)
+ rady = radius(n)
+ radz = radius(n)
+ else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+ radx = radius_x(n)
+ rady = radius_y(n)
+ radz = radius_z(n)
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+
+ if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+ sm_type = sm_linear
+ else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+ sm_type = sm_spline
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+ width = smoothing_zone_width(n)
+
+ dist2 = ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+ & + ((z - cz) / radz)**2
+
+ if (overwrite_geometry(n) /= 0) then
+
+ scale = Minkowski_scale(n)
+
+ if (conformal_state >= 1) then
+ where (dist2 <= 1)
+ psi = 1
+ end where
+ end if
+ if (conformal_state >= 2) then
+ where (dist2 <= 1)
+ psix = 0
+ psiy = 0
+ psiz = 0
+ end where
+ end if
+ if (conformal_state >= 3) then
+ where (dist2 <= 1)
+ psixx = 0
+ psixy = 0
+ psixz = 0
+ psiyy = 0
+ psiyz = 0
+ psizz = 0
+ end where
+ end if
+
+ where (dist2 <= 1)
+ ADM_BS_phi = smooth (ADM_BS_phi, zero , dist2)
+ ADM_BS_gxx = smooth (ADM_BS_gxx, scale, dist2)
+ ADM_BS_gxy = smooth (ADM_BS_gxy, zero , dist2)
+ ADM_BS_gxz = smooth (ADM_BS_gxz, zero , dist2)
+ ADM_BS_gyy = smooth (ADM_BS_gyy, scale, dist2)
+ ADM_BS_gyz = smooth (ADM_BS_gyz, zero , dist2)
+ ADM_BS_gzz = smooth (ADM_BS_gzz, scale, dist2)
+ ADM_BS_Gx = smooth (ADM_BS_Gx , zero , dist2)
+ ADM_BS_Gy = smooth (ADM_BS_Gy , zero , dist2)
+ ADM_BS_Gz = smooth (ADM_BS_Gz , zero , dist2)
+ ADM_BS_K = smooth (ADM_BS_K , zero , dist2)
+ ADM_BS_Axx = smooth (ADM_BS_Axx, zero , dist2)
+ ADM_BS_Axy = smooth (ADM_BS_Axy, zero , dist2)
+ ADM_BS_Axz = smooth (ADM_BS_Axz, zero , dist2)
+ ADM_BS_Ayy = smooth (ADM_BS_Ayy, zero , dist2)
+ ADM_BS_Ayz = smooth (ADM_BS_Ayz, zero , dist2)
+ ADM_BS_Azz = smooth (ADM_BS_Azz, zero , dist2)
+ end where
+ call CCTK_ActiveTimeLevelsVN &
+ (timelevels, cctkGH, "ADM_BSSN::ADM_BS_phi")
+ if (timelevels > 1) then
+ where (dist2 <= 1)
+ ADM_BS_phi_p = smooth (ADM_BS_phi_p, zero , dist2)
+ ADM_BS_gxx_p = smooth (ADM_BS_gxx_p, scale, dist2)
+ ADM_BS_gxy_p = smooth (ADM_BS_gxy_p, zero , dist2)
+ ADM_BS_gxz_p = smooth (ADM_BS_gxz_p, zero , dist2)
+ ADM_BS_gyy_p = smooth (ADM_BS_gyy_p, scale, dist2)
+ ADM_BS_gyz_p = smooth (ADM_BS_gyz_p, zero , dist2)
+ ADM_BS_gzz_p = smooth (ADM_BS_gzz_p, scale, dist2)
+ ADM_BS_Gx_p = smooth (ADM_BS_Gx_p , zero , dist2)
+ ADM_BS_Gy_p = smooth (ADM_BS_Gy_p , zero , dist2)
+ ADM_BS_Gz_p = smooth (ADM_BS_Gz_p , zero , dist2)
+ ADM_BS_K_p = smooth (ADM_BS_K_p , zero , dist2)
+ ADM_BS_Axx_p = smooth (ADM_BS_Axx_p, zero , dist2)
+ ADM_BS_Axy_p = smooth (ADM_BS_Axy_p, zero , dist2)
+ ADM_BS_Axz_p = smooth (ADM_BS_Axz_p, zero , dist2)
+ ADM_BS_Ayy_p = smooth (ADM_BS_Ayy_p, zero , dist2)
+ ADM_BS_Ayz_p = smooth (ADM_BS_Ayz_p, zero , dist2)
+ ADM_BS_Azz_p = smooth (ADM_BS_Azz_p, zero , dist2)
+ end where
+ end if
+ if (timelevels > 2) then
+ where (dist2 <= 1)
+ ADM_BS_phi_p_p = smooth (ADM_BS_phi_p_p, zero , dist2)
+ ADM_BS_gxx_p_p = smooth (ADM_BS_gxx_p_p, scale, dist2)
+ ADM_BS_gxy_p_p = smooth (ADM_BS_gxy_p_p, zero , dist2)
+ ADM_BS_gxz_p_p = smooth (ADM_BS_gxz_p_p, zero , dist2)
+ ADM_BS_gyy_p_p = smooth (ADM_BS_gyy_p_p, scale, dist2)
+ ADM_BS_gyz_p_p = smooth (ADM_BS_gyz_p_p, zero , dist2)
+ ADM_BS_gzz_p_p = smooth (ADM_BS_gzz_p_p, scale, dist2)
+ ADM_BS_Gx_p_p = smooth (ADM_BS_Gx_p_p , zero , dist2)
+ ADM_BS_Gy_p_p = smooth (ADM_BS_Gy_p_p , zero , dist2)
+ ADM_BS_Gz_p_p = smooth (ADM_BS_Gz_p_p , zero , dist2)
+ ADM_BS_K_p_p = smooth (ADM_BS_K_p_p , zero , dist2)
+ ADM_BS_Axx_p_p = smooth (ADM_BS_Axx_p_p, zero , dist2)
+ ADM_BS_Axy_p_p = smooth (ADM_BS_Axy_p_p, zero , dist2)
+ ADM_BS_Axz_p_p = smooth (ADM_BS_Axz_p_p, zero , dist2)
+ ADM_BS_Ayy_p_p = smooth (ADM_BS_Ayy_p_p, zero , dist2)
+ ADM_BS_Ayz_p_p = smooth (ADM_BS_Ayz_p_p, zero , dist2)
+ ADM_BS_Azz_p_p = smooth (ADM_BS_Azz_p_p, zero , dist2)
+ end where
+ end if
+
+ end if
+
+ if (overwrite_lapse(n) /= 0) then
+
+ scale = lapse_scale(n)
+
+ where (dist2 <= 1)
+ alp = smooth (alp , scale, dist2)
+ ADM_BS_dtalp = smooth (ADM_BS_dtalp, zero , dist2)
+ end where
+ call CCTK_ActiveTimeLevelsVN &
+ (timelevels, cctkGH, "ADMBase::alp")
+ if (timelevels > 1) then
+ where (dist2 <= 1)
+ alp_p = smooth (alp_p , scale, dist2)
+ ADM_BS_dtalp_p = smooth (ADM_BS_dtalp_p, zero , dist2)
+ end where
+ end if
+ if (timelevels > 2) then
+ where (dist2 <= 1)
+ alp_p_p = smooth (alp_p_p , scale, dist2)
+ ADM_BS_dtalp_p_p = smooth (ADM_BS_dtalp_p_p, zero , dist2)
+ end where
+ end if
+
+ end if
+
+ if (overwrite_shift(n) /= 0) then
+
+ if (shift_state /= 0) then
+ where (dist2 <= 1)
+ betax = smooth (betax , zero, dist2)
+ betay = smooth (betay , zero, dist2)
+ betaz = smooth (betaz , zero, dist2)
+ ADM_BS_Bx = smooth (ADM_BS_Bx, zero, dist2)
+ ADM_BS_By = smooth (ADM_BS_By, zero, dist2)
+ ADM_BS_Bz = smooth (ADM_BS_Bz, zero, dist2)
+ end where
+ call CCTK_ActiveTimeLevelsVN &
+ (timelevels, cctkGH, "ADMBase::betax")
+ if (timelevels > 1) then
+ where (dist2 <= 1)
+ betax_p = smooth (betax_p , zero, dist2)
+ betay_p = smooth (betay_p , zero, dist2)
+ betaz_p = smooth (betaz_p , zero, dist2)
+ ADM_BS_Bx_p = smooth (ADM_BS_Bx_p, zero, dist2)
+ ADM_BS_By_p = smooth (ADM_BS_By_p, zero, dist2)
+ ADM_BS_Bz_p = smooth (ADM_BS_Bz_p, zero, dist2)
+ end where
+ end if
+ if (timelevels > 2) then
+ where (dist2 <= 1)
+ betax_p_p = smooth (betax_p_p , zero, dist2)
+ betay_p_p = smooth (betay_p_p , zero, dist2)
+ betaz_p_p = smooth (betaz_p_p , zero, dist2)
+ ADM_BS_Bx_p_p = smooth (ADM_BS_Bx_p_p, zero, dist2)
+ ADM_BS_By_p_p = smooth (ADM_BS_By_p_p, zero, dist2)
+ ADM_BS_Bz_p_p = smooth (ADM_BS_Bz_p_p, zero, dist2)
+ end where
+ end if
+ end if
+
+ end if
+
+ end if
+
+ end do
+
+contains
+
+ elemental function smooth (oldval, goodval, dist2) result (newval)
+ CCTK_REAL :: newval
+ CCTK_REAL, intent(in) :: oldval
+ CCTK_REAL, intent(in) :: goodval
+ CCTK_REAL, intent(in) :: dist2
+
+ CCTK_REAL, parameter :: two=2, half=1/two
+ CCTK_REAL :: dist
+ CCTK_REAL :: x, f
+
+ if (dist2 <= (1 - width)**2) then
+ newval = goodval
+ else if (width == 0 .or. dist2 >= 1) then
+ newval = oldval
+ else
+ dist = sqrt(dist2)
+ x = 1 - (1 - dist) / width
+ select case (sm_type)
+ case (sm_linear)
+ f = x
+ case (sm_spline)
+ if (x < half) then
+ f = 2 * x**2
+ else
+ f = 1 - 2 * (1-x)**2
+ end if
+ end select
+ newval = oldval * f + goodval * (1 - f)
+ end if
+
+ end function smooth
+
+end subroutine NoExcision_OverwriteBSSN
diff --git a/src/reduce.F90 b/src/reduce.F90
new file mode 100644
index 0000000..ab7b70e
--- /dev/null
+++ b/src/reduce.F90
@@ -0,0 +1,97 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Reduce (cctk_iteration, cctk_lsh, rhs, x, y, z)
+ implicit none
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: cctk_iteration
+ integer :: cctk_lsh(3)
+ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: rhs, x, y, z
+
+ CCTK_REAL, parameter :: zero=0
+ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: dist2
+ CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
+ integer :: n
+
+ integer, parameter :: sm_linear = 1
+ integer, parameter :: sm_spline = 2
+ integer :: sm_type
+
+ do n = 1, num_regions
+ if (reduce_rhs(n) /= 0) then
+
+ cx = centre_x(n)
+ cy = centre_y(n)
+ cz = centre_z(n)
+ if (CCTK_EQUALS(region_shape(n), "sphere")) then
+ radx = radius(n)
+ rady = radius(n)
+ radz = radius(n)
+ else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+ radx = radius_x(n)
+ rady = radius_y(n)
+ radz = radius_z(n)
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+
+ if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+ sm_type = sm_linear
+ else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+ sm_type = sm_spline
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+ width = smoothing_zone_width(n)
+
+ dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
+
+#warning "TODO: reduce not only the RHS, but also drive the variables towards Minkowski"
+ where (dist2 <= 1)
+ rhs = smooth (rhs, rhs * reduction_factor(n), dist2)
+ end where
+
+ end if
+ end do
+
+contains
+
+ elemental function smooth (oldval, goodval, dist2) result (newval)
+ CCTK_REAL :: newval
+ CCTK_REAL, intent(in) :: oldval
+ CCTK_REAL, intent(in) :: goodval
+ CCTK_REAL, intent(in) :: dist2
+
+ CCTK_REAL, parameter :: two=2, half=1/two
+ CCTK_REAL :: dist
+ CCTK_REAL :: x, f
+
+ if (dist2 <= (1 - width)**2) then
+ newval = goodval
+ else if (width == 0 .or. dist2 >= 1) then
+ newval = oldval
+ else
+ dist = sqrt(dist2)
+ x = 1 - (1 - dist) / width
+ select case (sm_type)
+ case (sm_linear)
+ f = x
+ case (sm_spline)
+ if (x < half) then
+ f = 2 * x**2
+ else
+ f = 1 - 2 * (1-x)**2
+ end if
+ end select
+ newval = oldval * f + goodval * (1 - f)
+ end if
+
+ end function smooth
+
+end subroutine NoExcision_Reduce
diff --git a/src/reduce.c b/src/reduce.c
new file mode 100644
index 0000000..0b44ca2
--- /dev/null
+++ b/src/reduce.c
@@ -0,0 +1,34 @@
+/* $Header$ */
+
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void CCTK_FCALL
+CCTK_FNAME(NoExcision_Reduce) (int const * cctk_iteration,
+ int const * cctk_lsh,
+ CCTK_REAL * rhs,
+ CCTK_REAL const * x,
+ CCTK_REAL const * y,
+ CCTK_REAL const * z);
+
+void
+NoExcision_Reduce (CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ int var, rhs;
+ void * ptr;
+
+ for (var=0; var<CCTK_NumVars(); ++var) {
+ rhs = MoLQueryEvolvedRHS (var);
+ if (rhs >= 0) {
+ ptr = CCTK_VarDataPtrI (cctkGH, 0, rhs);
+ assert (ptr);
+ CCTK_FNAME(NoExcision_Reduce) (& cctk_iteration, cctk_lsh, ptr, x, y, z);
+ }
+ }
+}
diff --git a/src/smooth.F90 b/src/smooth.F90
new file mode 100644
index 0000000..d6df156
--- /dev/null
+++ b/src/smooth.F90
@@ -0,0 +1,150 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Smooth (CCTK_ARGUMENTS)
+ implicit none
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL, parameter :: zero=0
+ CCTK_REAL :: mask (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+ CCTK_REAL :: cx, cy, cz, radx, rady, radz
+ integer :: n
+ integer :: iter
+ integer :: ierr
+
+ do iter = 1, smoothing_iterations
+
+ do n = 1, num_regions
+
+ cx = centre_x(n)
+ cy = centre_y(n)
+ cz = centre_z(n)
+ if (CCTK_EQUALS(region_shape(n), "sphere")) then
+ radx = radius(n)
+ rady = radius(n)
+ radz = radius(n)
+ else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+ radx = radius_x(n)
+ rady = radius_y(n)
+ radz = radius_z(n)
+ else
+ call CCTK_WARN (0, "internal error")
+ end if
+
+ mask = 1 - sqrt (((x - cx) / radx)**2 + &
+ & ((y - cy) / rady)**2 + &
+ & ((z - cz) / radz)**2)
+ where (mask <= 0)
+ mask = 0 ! outside
+ else where (mask >= smoothing_zone_width(n))
+ mask = 1 ! far inside
+ elsewhere
+ mask = mask / smoothing_zone_width(n) ! a bit inside
+ end where
+
+ if (overwrite_geometry(n) /= 0) then
+
+ if (conformal_state >= 1) then
+ call smooth (psi)
+ end if
+ if (conformal_state >= 2) then
+ call smooth (psix)
+ call smooth (psiy)
+ call smooth (psiz)
+ end if
+ if (conformal_state >= 3) then
+ call smooth (psixx)
+ call smooth (psixy)
+ call smooth (psixz)
+ call smooth (psiyy)
+ call smooth (psiyz)
+ call smooth (psizz)
+ end if
+
+ call smooth (gxx)
+ call smooth (gxy)
+ call smooth (gxz)
+ call smooth (gyy)
+ call smooth (gyz)
+ call smooth (gzz)
+ call smooth (kxx)
+ call smooth (kxy)
+ call smooth (kxz)
+ call smooth (kyy)
+ call smooth (kyz)
+ call smooth (kzz)
+
+ end if
+
+ if (overwrite_lapse(n) /= 0) then
+
+ call smooth (alp)
+
+ end if
+
+ if (overwrite_shift(n) /= 0) then
+
+ if (shift_state /= 0) then
+ call smooth (betax)
+ call smooth (betay)
+ call smooth (betaz)
+ end if
+
+ end if
+
+ end do
+
+ if (overwrite_geometry(n) /= 0) then
+
+ if (conformal_state >= 1) then
+ call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac")
+ end if
+ if (conformal_state >= 2) then
+ call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_1derivs")
+ end if
+ if (conformal_state >= 3) then
+ call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_2derivs")
+ end if
+
+ call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::metric")
+ call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::curv")
+
+ end if
+
+ if (overwrite_lapse(n) /= 0) then
+
+ call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::lapse")
+
+ end if
+
+ if (overwrite_shift(n) /= 0) then
+
+ if (shift_state /= 0) then
+ call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::shift")
+ end if
+
+ end if
+
+ end do
+
+contains
+
+ subroutine smooth (val)
+ CCTK_REAL, intent(inout) :: val(:,:,:)
+
+ where (mask > 0)
+ val = (1 - mask * smoothing_factor) * val &
+ & + (mask * smoothing_factor / 6) &
+ & * (+ eoshift(val, +1, dim=1) + eoshift(val, -1, dim=1) &
+ & + eoshift(val, +1, dim=2) + eoshift(val, -1, dim=2) &
+ & + eoshift(val, +1, dim=3) + eoshift(val, -1, dim=3))
+ end where
+ end subroutine smooth
+
+end subroutine NoExcision_Smooth
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-18 02:25:49 +0000