Merge branch 'debugGlassMerge'

MintonGroup · Jan 28, 2022 · a26b05e · a26b05e
2 parents c4c5cae + 988a780
commit a26b05e
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 14 deletions.
diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90
@@ -68,7 +68,7 @@ subroutine crater_dimensions(user,crater,domain)
 
    ! Calculate the radius where the inner wall meets the original pre-existing surface
    ! This is used to demark the location where excavation transitions to deposition
-   crater%ejrad = crater_profile_find_r_inner_wall(user,crater) * crater%frad
+   crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad)
 
    !find rim for counting purposes
    crater%frim = RIMFAC * crater%frad

diff --git a/src/ejecta/ejecta_rootfind.f90 b/src/ejecta/ejecta_rootfind.f90
@@ -43,10 +43,10 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
    real(DP) :: resultat,AA,BB,CC,DD,EE,FA,FB,FC,Tol1,PP,QQ,RR,SS,xm
    integer(I4B) :: i,j, ev,br
    integer(I4B),parameter :: NTRY = 50
-   integer(I4B),parameter :: NBRACKET = 5
+   integer(I4B),parameter :: NBRACKET =20 
    real(DP),parameter :: FIRSTFACTOR = 1.6_DP
    real(DP) :: factor
-   real(DP) :: f1,f2
+   real(DP) :: f1, f2, fmin
    integer(I4B) :: numev,numbr
 
    ! Executable code
@@ -64,22 +64,32 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
          numbr=numbr+1
          erad = starterad
          x2 = starterad/factor
-
          ! First bracket the root
          f1 = ejecta_blanket_func(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
          f2 = ejecta_blanket_func(user,crater,domain,x2,lrad,vejsq,ejang,firstrun)
+         fmin = abs(f1)
          do j=1,NTRY
-            if (f1*f2 < 0._DP) exit
+            if ((f1 * f2 < 0._DP)) exit bracket
             if (abs(f1) < abs(f2)) then
-               erad=erad+factor*(erad-x2)
-               f1=ejecta_blanket_func(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
+               erad = erad + factor * (erad - x2)
+               if (erad < 0.0_DP) exit
+               f1 = ejecta_blanket_func(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
+               if (abs(f1) < fmin) then
+                  fmin = abs(f1)
+                  starterad = erad
+               end if
             else
-               x2=x2+factor*(x2-erad)
-               f2=ejecta_blanket_func(user,crater,domain,x2,lrad,vejsq,ejang,firstrun)
+               x2 = x2 + factor * (x2 - erad)
+               if (x2 < 0.0_DP) exit
+               f2 = ejecta_blanket_func(user,crater,domain,x2,lrad,vejsq,ejang,firstrun)
+               if (abs(f2) < fmin) then
+                  fmin = abs(f2)
+                  starterad = x2
+               end if
             end if
          end do
-         if ((erad<=crater%ejrad).and.(erad>0._DP)) exit
-         factor=0.5_DP*(factor+1._DP)
+         erad = x2
+         factor = 0.5_DP * (factor + 1._DP)
       end do bracket
 
       ! Now do a Brent's method to find the root
@@ -153,7 +163,7 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
       end if
       niter = i
       erad = resultat
-      if ((erad <= crater%ejrad).and.(lrad >= crater%ejrad).and.(erad >= 0._DP)) exit
+      if (erad >= 0._DP) exit
       factor = 0.5_DP  * (factor + 1.0_DP) ! Failed. Try again with a new factor
    end do everything
 

diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90
@@ -47,8 +47,9 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
    crater%ejdis = DISEJB * crater%continuous
                                                         ! We go out a factor of 3 to get the discontinuous ejecta thickness 
    domain%ejbres = (log(crater%ejdis) - log(crater%ejrad)) / EJBTABSIZE
-   lrad = crater%ejrad 
-   erad = crater%rad 
+   lrad = crater%ejrad
+   lrad = exp(log(lrad) + domain%ejbres)
+   erad = crater%ejrad 
    ejtble = EJBTABSIZE
    firstrun = .true.
    thick = 0._DP