Changes related to the mask, and to dynamical excision. Everything about

excision here is protected by ifdef's.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@53 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 89 insertions, 10 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index 8918fdd..7ccc12b 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -9,8 +9,10 @@
 @@*/
 
 #include "cctk.h"
+
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
+#include "cctk_DefineThorn.h"
 #include "CactusEinstein/Einstein/src/Einstein.h"
 
       subroutine AHFinder(CCTK_FARGUMENTS)
@@ -41,6 +43,7 @@
       CCTK_REAL aux,expin,expout,small1,small2
       CCTK_REAL c0_temp(18)
       CCTK_REAL maxval(3),minval(3)
+      CCTK_REAL rhor
 
       character*200 almf,logf,areaf,massf,circeqf,merip1f,merip2f
 
@@ -117,6 +120,17 @@
       firstint = .true.
 
 
+!     ***************************
+!     ***   INITIALIZE MASK   ***
+!     ***************************
+
+!     Initialize mask to one, but only on first iteration.
+
+      if (cctk_iteration.le.1) then
+         ahmask = one
+      end if
+
+
 !     ***********************************************
 !     ***   CHECK IF WE WANT TO FIND 3 HORIZONS   ***
 !     ***********************************************
@@ -469,9 +483,7 @@
 
       if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3
 
-!     Check if we are using Cartoon_2d, in which case
-!     we override the value of "ahf_phi" since we have
-!     axisymmetry.
+!     Check if we are using cartoon.
 
       cartoon = .false.
 
@@ -1373,14 +1385,79 @@
 !     ***   MASK   ***
 !     ****************
 
-      if (horizon.and.(ahf_mask.ne.0)) then
-         call AHFinder_mask(CCTK_FARGUMENTS)
-         call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask")
-      else if ((mtype.eq.1).and.(ahf_weakmask.ne.0)) then
-         call AHFinder_mask(CCTK_FARGUMENTS)
+!     Check if we are using dynamical excision with the mask off.
+
+#ifdef DEVELOPMENT_SIMPLEEXCISION
+      if (CCTK_EQUALS(excisiontype,'dynamic')) then
+         if ((ahf_mask.eq.0).and.(ahf_weakmask.eq.0)) then
+            call CCTK_WARN(1,'AHFinder: You selected Dynamic excision but the mask is off!')
+         end if
+      end if
+#endif
+
+!     The mask is only different from 1 inside a horizon.
+!     I therefore need to check if I do have a horizon, and
+!     also if the mask is desired.
+
+      if ((horizon.and.(ahf_mask.ne.0)).or.
+     .   ((mtype.eq.1).and.(ahf_weakmask.ne.0))) then
+
+         call AHFinder_mask(CCTK_FARGUMENTS,rhor)
          call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask")
-      else if (cctk_iteration.lt.1) then
-         ahmask = one
+
+#ifdef DEVELOPMENT_SIMPLEEXCISION
+
+!        Here I alter the parameters for SimpleExcision.  Notice that
+!        at the moment, I only allow the excised cube to grow, never
+!        to shrink.
+
+         if (find3) then
+
+            if (mfind.eq.0) then
+
+!              First hole.
+
+               if (dhole1_xmin.gt.(xc-0.5D0*rhor)) dhole1_xmin = xc-0.5D0*rhor
+               if (dhole1_ymin.gt.(yc-0.5D0*rhor)) dhole1_ymin = yc-0.5D0*rhor
+               if (dhole1_zmin.gt.(zc-0.5D0*rhor)) dhole1_zmin = zc-0.5D0*rhor
+
+               if (dhole1_xmax.lt.(xc+0.5D0*rhor)) dhole1_xmax = xc+0.5D0*rhor
+               if (dhole1_ymax.lt.(yc+0.5D0*rhor)) dhole1_ymax = yc+0.5D0*rhor
+               if (dhole1_zmax.lt.(zc+0.5D0*rhor)) dhole1_zmax = zc+0.5D0*rhor
+
+            else if (mfind.eq.1) then
+
+!              Second hole.
+
+               if (dhole2_xmin.gt.(xc-0.5D0*rhor)) dhole2_xmin = xc-0.5D0*rhor
+               if (dhole2_ymin.gt.(yc-0.5D0*rhor)) dhole2_ymin = yc-0.5D0*rhor
+               if (dhole2_zmin.gt.(zc-0.5D0*rhor)) dhole2_zmin = zc-0.5D0*rhor
+
+               if (dhole2_xmax.lt.(xc+0.5D0*rhor)) dhole2_xmax = xc+0.5D0*rhor
+               if (dhole2_ymax.lt.(yc+0.5D0*rhor)) dhole2_ymax = yc+0.5D0*rhor
+               if (dhole2_zmax.lt.(zc+0.5D0*rhor)) dhole2_zmax = zc+0.5D0*rhor
+
+            else
+
+!              Third hole.  This is not yet here since the parameters
+!              do not exist yet.
+
+            end if
+
+         else
+
+            if (dhole1_xmin.gt.(xc-0.5D0*rhor)) dhole1_xmin = xc-0.5D0*rhor
+            if (dhole1_ymin.gt.(yc-0.5D0*rhor)) dhole1_ymin = yc-0.5D0*rhor
+            if (dhole1_zmin.gt.(zc-0.5D0*rhor)) dhole1_zmin = zc-0.5D0*rhor
+
+            if (dhole1_xmax.lt.(xc+0.5D0*rhor)) dhole1_xmax = xc+0.5D0*rhor
+            if (dhole1_ymax.lt.(yc+0.5D0*rhor)) dhole1_ymax = yc+0.5D0*rhor
+            if (dhole1_zmax.lt.(zc+0.5D0*rhor)) dhole1_zmax = zc+0.5D0*rhor
+
+         end if
+
+#endif
+
       end if
 
 
@@ -1398,6 +1475,7 @@
          gzz=gzz/psi**4
       end if
 
+
 !     **********************************************
 !     ***   UPDATING COUNTER AND LOOP IF FIND3   ***
 !     **********************************************
@@ -1411,6 +1489,7 @@
          end if
       end if
 
+
 !     ***************
 !     ***   END   ***
 !     ***************