aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl2
-rw-r--r--src/EHFinder_Constants.F902
-rw-r--r--src/EHFinder_FindSurface.F9050
-rw-r--r--src/EHFinder_fuzzy.F9083
-rw-r--r--src/make.code.defn3
-rw-r--r--src/make.code.deps3
6 files changed, 118 insertions, 25 deletions
diff --git a/interface.ccl b/interface.ccl
index ae4d0e7..e63187f 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -75,7 +75,7 @@ CCTK_INT niter_reinit TYPE=SCALAR
#CCTK_INT carpet_re_init_control TYPE=SCALAR
# Grid function used in counting surfaces.
-CCTK_INT surface_index TYPE=GF TIMELEVELS=1 tags='tensortypealias="Scalar" Prolongation="None"'
+CCTK_REAL surface_index TYPE=GF TIMELEVELS=1 tags='tensortypealias="Scalar"
{
sc
}
diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90
index 2162631..c2c316e 100644
--- a/src/EHFinder_Constants.F90
+++ b/src/EHFinder_Constants.F90
@@ -16,4 +16,6 @@ module EHFinder_Constants
CCTK_REAL, parameter :: fourthirds = four / three
CCTK_REAL, parameter :: quarter = 0.25d0
CCTK_REAL, parameter :: pi = 3.1415926535897932385d0
+ CCTK_REAL, parameter :: small = 1.0d-10
+ integer, parameter :: wp = kind(zero)
end module EHFinder_Constants
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 63981bf..4662360 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -9,6 +9,7 @@
subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
use EHFinder_mod
+ use EHFinder_fuzzy
implicit none
@@ -115,13 +116,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
( f(i,j,k,l) * f(i,j+1,k,l) .lt. zero ) .or. &
( f(i,j,k,l) * f(i,j,k-1,l) .lt. zero ) .or. &
( f(i,j,k,l) * f(i,j,k+1,l) .lt. zero ) ) .and. &
- ( ( sc(i,j,k) .eq. sn ) .or. &
- ( sc(i-1,j,k) .eq. sn ) .or. &
- ( sc(i+1,j,k) .eq. sn ) .or. &
- ( sc(i,j-1,k) .eq. sn ) .or. &
- ( sc(i,j+1,k) .eq. sn ) .or. &
- ( sc(i,j,k-1) .eq. sn ) .or. &
- ( sc(i,j,k+1) .eq. sn ) ) ) then
+ ( fuzzy_eq ( sc(i,j,k), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i-1,j,k), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i+1,j,k), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i,j-1,k), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i,j+1,k), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i,j,k-1), real(sn,wp), small ) .or. &
+ fuzzy_eq ( sc(i,j,k+1), real(sn,wp), small ) ) ) then
nc_loc = nc_loc + 1
center_loc(1) = center_loc(1) + x(i,j,k)
center_loc(2) = center_loc(2) + y(i,j,k)
@@ -144,12 +145,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
! Try to figure out the symmetries. First find the range in coordinates
! for the points inside the current surface.
- crange_min_loc(1) = minval ( x, mask = sc .eq. sn )
- crange_min_loc(2) = minval ( y, mask = sc .eq. sn )
- crange_min_loc(3) = minval ( z, mask = sc .eq. sn )
- crange_max_loc(1) = maxval ( x, mask = sc .eq. sn )
- crange_max_loc(2) = maxval ( y, mask = sc .eq. sn )
- crange_max_loc(3) = maxval ( z, mask = sc .eq. sn )
+ crange_min_loc(1) = minval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+ crange_min_loc(2) = minval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+ crange_min_loc(3) = minval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+ crange_max_loc(1) = maxval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+ crange_max_loc(2) = maxval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
+ crange_max_loc(3) = maxval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) )
! Reduce over all processors.
call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
@@ -741,7 +742,7 @@ subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS)
more_surfaces = 1
surface_counter = 0
- sc = 0
+ sc = zero
end subroutine EHFinder_CountSurfacesInit
@@ -749,6 +750,7 @@ end subroutine EHFinder_CountSurfacesInit
subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
use EHFinder_mod
+ use EHFinder_fuzzy
implicit none
@@ -767,8 +769,8 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
! Find the location of the minimum among active points not already found
! to be inside a surface.
minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter), &
- ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. &
- ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
+ fuzzy_eq ( sc(ixl:ixr,jyl:jyr,kzl:kzr), zero, small ) .and. &
+ ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
! Find the value of f at that location.
minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter)
@@ -804,7 +806,7 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
points_counter = 0
more_points = 1
if ( procpos_loc .eq. procpos_glob ) then
- sc(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) = surface_counter
+ sc(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) = real(surface_counter,wp)
end if
else
more_surfaces = 0
@@ -818,6 +820,7 @@ end subroutine EHFinder_CountSurfaces
subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
use EHFinder_mod
+ use EHFinder_fuzzy
implicit none
@@ -846,7 +849,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
marked = .false.
if ( ( eh_mask(i,j,k,levelset_counter) .ge. 0 ) .and. &
( f(i,j,k,levelset_counter) .lt. 0 ) .and. &
- ( sc(i,j,k) .ne. surface_counter ) ) then
+ fuzzy_ne ( sc(i,j,k), real(surface_counter,wp), small ) ) then
! Find out if we are on the boundary
if ( i + cctk_lbnd(1) .eq. 1 ) then
@@ -881,18 +884,21 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
end if
! Then check if any valid neighbours are marked
- if ( any ( sc(i+i1:i+i2,j+j1:j+j2,k) .eq. surface_counter ) ) then
+ if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j+j1:j+j2,k), &
+ real(surface_counter,wp), small ) ) ) then
marked = .true.
end if
- if ( any ( sc(i+i1:i+i2,j,k+k1:k+k2) .eq. surface_counter ) ) then
+ if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j,k+k1:k+k2), &
+ real(surface_counter,wp), small ) ) ) then
marked = .true.
end if
- if ( any ( sc(i,j+j1:j+j2,k+k1:k+k2) .eq. surface_counter ) ) then
+ if ( any ( fuzzy_eq ( sc(i,j+j1:j+j2,k+k1:k+k2), &
+ real(surface_counter,wp), small ) ) ) then
marked = .true.
end if
if ( marked ) then
- sc(i,j,k) = surface_counter
+ sc(i,j,k) = real(surface_counter,wp)
n_loc = n_loc + 1
end if
diff --git a/src/EHFinder_fuzzy.F90 b/src/EHFinder_fuzzy.F90
new file mode 100644
index 0000000..a856fef
--- /dev/null
+++ b/src/EHFinder_fuzzy.F90
@@ -0,0 +1,83 @@
+! Module with routines for fuzzy comparisons.
+! $Header$
+#include "cctk.h"
+
+module EHFinder_fuzzy
+
+ use EHFinder_Constants
+
+ interface fuzzy_eq
+ module procedure fuzzy_eq_0d, fuzzy_eq_1d, fuzzy_eq_2d, fuzzy_eq_3d
+ end interface
+
+ interface fuzzy_ne
+ module procedure fuzzy_ne_0d, fuzzy_ne_1d, fuzzy_ne_2d, fuzzy_ne_3d
+ end interface
+
+contains
+
+ function fuzzy_eq_0d ( x, y, eps )
+ CCTK_REAL, intent(IN) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical :: fuzzy_eq_0d
+
+ fuzzy_eq_0d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+ end function fuzzy_eq_0d
+
+ function fuzzy_eq_1d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1)) :: fuzzy_eq_1d
+
+ fuzzy_eq_1d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+ end function fuzzy_eq_1d
+
+ function fuzzy_eq_2d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:,:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1),size(x,2)) :: fuzzy_eq_2d
+
+ fuzzy_eq_2d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+ end function fuzzy_eq_2d
+
+ function fuzzy_eq_3d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:,:,:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1),size(x,2),size(x,3)) :: fuzzy_eq_3d
+
+ fuzzy_eq_3d = ( ( x .gt. y-eps ) .and. ( x .lt. y+eps ) )
+ end function fuzzy_eq_3d
+
+ function fuzzy_ne_0d ( x, y, eps )
+ CCTK_REAL, intent(IN) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical :: fuzzy_ne_0d
+
+ fuzzy_ne_0d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+ end function fuzzy_ne_0d
+
+ function fuzzy_ne_1d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1)) :: fuzzy_ne_1d
+
+ fuzzy_ne_1d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+ end function fuzzy_ne_1d
+
+ function fuzzy_ne_2d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:,:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1),size(x,2)) :: fuzzy_ne_2d
+
+ fuzzy_ne_2d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+ end function fuzzy_ne_2d
+
+ function fuzzy_ne_3d ( x, y, eps )
+ CCTK_REAL, intent(IN), dimension(:,:,:) :: x
+ CCTK_REAL, intent(IN) :: y, eps
+ logical, dimension(size(x,1),size(x,2),size(x,3)) :: fuzzy_ne_3d
+
+ fuzzy_ne_3d = ( ( x .lt. y-eps ) .or. ( x .gt. y+eps ) )
+ end function fuzzy_ne_3d
+
+end module EHFinder_fuzzy
diff --git a/src/make.code.defn b/src/make.code.defn
index 701a619..d24f363 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -18,7 +18,8 @@ SRCS = EHFinder_mod.F90 \
EHFinder_Generator_Sources.F90 \
EHFinder_Generator_Sources2.F90 \
EHFinder_ReInitialize.F90 \
- EHFinder_CarpetSupport.cc
+ EHFinder_CarpetSupport.cc \
+ EHFinder_fuzzy.F90
# Subdirectories containing source files
SUBDIRS =
diff --git a/src/make.code.deps b/src/make.code.deps
index 0b34034..21d5cc2 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -3,6 +3,7 @@
EHFinder_mod.F90.o: EHFinder_Constants.F90.o
EHFinder_IsoSurface_mod.F90.o: EHFinder_Constants.F90.o
+EHFinder_fuzzy.F90.o: EHFinder_Constants.F90.o
EHFinder_Init.F90.o: EHFinder_mod.F90.o
EHFinder_Sources.F90.o: EHFinder_mod.F90.o
EHFinder_Generator_Sources.F90.o: EHFinder_mod.F90.o
@@ -14,5 +15,5 @@ EHFinder_SetSym.F90.o: EHFinder_mod.F90.o
EHFinder_ReadData.F90.o: EHFinder_mod.F90.o
EHFinder_Check.F90.o: EHFinder_mod.F90.o
EHFinder_Integrate.F90.o: EHFinder_mod.F90.o
-EHFinder_FindSurface.F90.o: EHFinder_mod.F90.o
+EHFinder_FindSurface.F90.o: EHFinder_mod.F90.o EHFinder_fuzzy.F90.o
EHFinder_IsoSurface.F90.o: EHFinder_mod.F90.o IsoSurface_mod.F90.o
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-21 14:11:50 +0000