Event detection is now modified to check for sign changes between the…

… initial condition and the first accepted step
egiovan · May 6, 2020 · 03ea266 · 03ea266
1 parent e8052d2
commit 03ea266
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 8 deletions.
diff --git a/src/ERKIntegrate.f90 b/src/ERKIntegrate.f90
@@ -240,6 +240,10 @@ module subroutine erk_int(me, X0, Y0, Xf, Yf, StepSz, UseConstStepSz, IntStepsOn
             !me%Yint(:,1) = Y0
         end if
 
+        ! Evaluate the event function at the initial condition to check for
+        ! events between X0 and first integrated state
+        if (EventsOn) call pDiffEqSys%G(X0, Y0, EV0, EventDir, EventTerm)
+
         ! Initialize the number of steps taken and other counters
         X = X0
         Y1 = Y0
@@ -298,8 +302,9 @@ module subroutine erk_int(me, X0, Y0, Xf, Yf, StepSz, UseConstStepSz, IntStepsOn
                     end if
 
                     ! Events checking start here, ignore event checking at IC
+
                     BipComputed = .FALSE.
-                    if (EventsOn .AND. AcceptedSteps > 1) then
+                    if (EventsOn) then
                         ! generate X values at which event triggers are checked
                         if (EventStepSz == 0.0 .OR. EventStepSz >= abs(h)) then
                             Eventh = h
@@ -323,9 +328,6 @@ module subroutine erk_int(me, X0, Y0, Xf, Yf, StepSz, UseConstStepSz, IntStepsOn
                         EventX0 = X
                         EventX1 = X + Eventh
 
-                        ! The event value at the first integrated state right after the initial condition
-                        if (AcceptedSteps == 2) call pDiffEqSys%G(X, Y1, EV0, EventDir, EventTerm)
-
                         ! check for event triggers between X and X+h
              EventLoop: do
                             block

diff --git a/src/tests/DiffEq.f90 b/src/tests/DiffEq.f90
@@ -155,15 +155,25 @@ subroutine SampleEventCR3BP(me, X, Y, Value, Direction, Terminal)
         logical, dimension(:), intent(out) :: Terminal
 
         Value = [Y(1),Y(2)]
-
+
+        ! Detect discontinuous events
+        if (X > 0 .AND. X < 1) then
+            Value(1) = 1        
+        else
+            Value(1) = -1
+        end if
+
         ! only detect Y-crossings in the region X<0
         if (Y(1) > 0.0)   Value(2) = IEEE_VALUE(1.0,IEEE_QUIET_NAN) 
+
 
-        Direction = [-1,1]
+
+        Direction = [0,1]
         Terminal = [.FALSE.,.FALSE.]
 
     end subroutine SampleEventCR3BP          
 
+
     function TBIOM(mu, X)
         real(wp), intent(in) :: mu
         real(wp), dimension(:), intent(in) :: X

diff --git a/src/tests/test.f90 b/src/tests/test.f90
@@ -39,7 +39,7 @@ program TestFLINT
     integer, parameter :: nloops = 1
 
     ! Turn on the events
-    logical, parameter :: EventsEnabled = .FALSE.
+    logical, parameter :: EventsEnabled = .TRUE.
 
     ! scalar tolerance
     real(wp),parameter :: tol   = 1.0e-11_WP  
@@ -89,7 +89,7 @@ program TestFLINT
     Xarr = [(x0+ipdx*itr,itr=0,(nIp-1))]
     allocate(Yarr(6,size(Xarr)))
 
-    stepsz0 = 1.0E-3_WP    ! let FLINT compute the initial step-size    
+    stepsz0 = 0.0E-3_WP    ! let FLINT compute the initial step-size    
     stifftestval = 1  ! check for stiffness and stop integration if stiff
 
     ! create results file