From 1587fd0060f15d41893f30e592b1de7c2ce94032 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 8 Apr 2020 13:17:04 -0400 Subject: [PATCH] Updated profile code to better connect the Pike (1977) and Fassett & Thomson (2014) models. Also improved the complex rim code to better control the shape of the scallops --- src/Makefile.am | 1 + src/crater/crater_dimensions.f90 | 17 ++- src/crater/crater_emplace.f90 | 2 - src/crater/crater_form_interior.f90 | 33 +++--- src/crater/crater_profile.f90 | 121 +++++++++++++++++++-- src/crater/crater_realistic_topography.f90 | 48 ++++---- src/crater/module_crater.f90 | 40 +++++-- src/ejecta/ejecta_table_define.f90 | 17 +-- 8 files changed, 198 insertions(+), 81 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index 18b3ae2a..1e0e80b9 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -86,6 +86,7 @@ crater/crater_critical_slope.f90\ crater/crater_degradation_function.f90\ crater/crater_visibility.f90\ crater/crater_get_degradation_state.f90\ +crater/crater_profile.f90\ init/init_domain.f90\ init/init_dist.f90\ init/init_surf.f90\ diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90 index 6ea7302a..9495a9d3 100644 --- a/src/crater/crater_dimensions.f90 +++ b/src/crater/crater_dimensions.f90 @@ -39,11 +39,14 @@ subroutine crater_dimensions(user,crater,domain) ! Set the crater morphology select case(crater%morphtype) - case("SIMPLE","TRANSITION") ! Following Pike (1977) - crater%rimheight = 0.036_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat - crater%rimwidth = 0.257_DP * (crater%fcrat * 1e-3_DP)**(1.011_DP) * 1e3_DP !* crater%fcrat - crater%floordepth = 0.196_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat - crater%floordiam = 0.031_DP * (crater%fcrat * 1e-3_DP)**(1.765_DP) * 1e3_DP !* crater%fcrat + case("SIMPLE","TRANSITION") ! A hybrid model between Pike (1977) and Fassett & Thomson (2014) + !crater%rimheight = 0.036_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat ! Pike model + crater%rimheight = 0.043_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat ! Closer to Fassett & Thomson + crater%rimwidth = 0.257_DP * (crater%fcrat * 1e-3_DP)**(1.011_DP) * 1e3_DP !* crater%fcrat ! Pike model + !crater%floordepth = 0.196_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat ! Pike model + crater%floordepth = 0.224_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat ! Closer to Fassett & Thomson + !crater%floordiam = 0.031_DP * (crater%fcrat * 1e-3_DP)**(1.765_DP) * 1e3_DP !* crater%fcrat ! Pike model + crater%floordiam = 0.200_DP * (crater%fcrat * 1e-3_DP)**(1.143_DP) * 1e3_DP !* crater%fcrat ! Fassett & Thomson for D~1km, Pike for D~20km case("COMPLEX","PEAKRING","MULTIRING") ! Following Pike (1977) crater%rimheight = 0.236_DP * (crater%fcrat * 1e-3_DP)**(0.399_DP) * 1e3_DP !* crater%fcrat crater%rimwidth = 0.467_DP * (crater%fcrat * 1e-3_DP)**(0.836_DP) * 1e3_DP !* crater%fcrat @@ -63,6 +66,10 @@ subroutine crater_dimensions(user,crater,domain) ! Adjust the floor depth to measure from the pre-existing level surface, rather than the rim crater%floordepth = crater%floordepth - crater%rimheight + ! 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 + !find rim for counting purposes crater%frim = RIMFAC * crater%frad diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index b67c0983..c9544dcb 100644 --- a/src/crater/crater_emplace.f90 +++ b/src/crater/crater_emplace.f90 @@ -80,8 +80,6 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot) deltaMtot = 0.0_DP !ejbmass incsq = inc**2 - ! TEMP UNTIL BETTER FORMULA - crater%ejrad = crater%frad ! This loop may not be parallelizable because of the linked list operation inside crater_form_interior do j=-inc,inc ! Do the loop in pixel space do i=-inc,inc diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index a1ea3572..369ba7a9 100644 --- a/src/crater/crater_form_interior.f90 +++ b/src/crater/crater_form_interior.f90 @@ -54,30 +54,29 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele real(DP),parameter :: outer_c2 = 0.018_DP real(DP),parameter :: outer_c3 = 0.015_DP - real(DP) :: c0,c1,c2,c3,flrad,rh,fld + !real(DP) :: c0,c1,c2,c3,flrad,rh,fld ! Executable code r = sqrt(x_relative**2+y_relative**2) / crater%frad - rh = crater%rimheight - fld = -crater%floordepth - flrad = 0.5_DP * crater%floordiam / crater%frad + !rh = crater%rimheight + !fld = -crater%floordepth + !flrad = 0.5_DP * crater%floordiam / crater%frad ! Use polynomial crater profile similar to that of Fassett et al. (2014), but the parameters are set by the crater dimensions - c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP) - c0 = rh - (7._DP / 6._DP) * c1 - c2 = c1 / 3._DP - c3 = -c2 / 2._DP - - if (r < flrad) then - cform = fld - else - cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 - end if - ! TEMP UNTIL I WRITE THE SOLUTION PROPERLY - if (cform > 0.0_DP .and. r * crater%frad < crater%ejrad) crater%ejrad = r * crater%frad + !c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP) + !c0 = rh - (7._DP / 6._DP) * c1 + !c2 = c1 / 3._DP + !c3 = -c2 / 2._DP +! +! if (r < flrad) then +! cform = fld +! else +! cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 +! end if + !if (r < r_floor) then ! cform = -simple_depth_diam * crater%fcrat @@ -86,6 +85,8 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele !else ! cform = (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat !end if + + cform = crater_profile(user,crater,r) newdem = newelev + cform !if (crater%fcrat > crater%cxtran * 2) then diff --git a/src/crater/crater_profile.f90 b/src/crater/crater_profile.f90 index f6dce262..ce033da3 100644 --- a/src/crater/crater_profile.f90 +++ b/src/crater/crater_profile.f90 @@ -27,29 +27,134 @@ function crater_profile(user,crater,r) result(h) type(usertype),intent(in) :: user type(cratertype),intent(in) :: crater real(DP),intent(in) :: r + + ! Result variable real(DP) :: h ! Internal variables real(DP) :: flrad,c0,c1,c2,c3 + real(DP),parameter :: A = 4._DP / 11._DP + real(DP),parameter :: B = -32._DP / 187._DP ! Executable code - flrad = 0.5_DP * crater%floordiam / crater%frad + flrad = crater%floordiam / crater%fcrat ! Use polynomial crater profile similar to that of Fassett et al. (2014), but the parameters are set by the crater dimensions - c1 = (-crater%floordepth - crater%rimheight) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP) - c0 = crater%rimheight - (7._DP / 6._DP) * c1 - c2 = c1 / 3._DP - c3 = -c2 / 2._DP + c1 = (-crater%floordepth - crater%rimheight) / (flrad - 1._DP + A * (flrad**2 - 1._DP) + B * (flrad**3 - 1._DP)) + c0 = crater%rimheight - c1 * (1._DP + A + B) + c2 = A * c1 + c3 = B * c1 if (r < flrad) then h = -crater%floordepth - if (r > crater%frad) - h = crater%rimheight * ((crater%frad / lrad)**RIMDROP) + else if (r > 1.0_DP) then + h = crater%rimheight * (r**(-RIMDROP)) else h = c0 + c1 * r + c2 * r**2 + c3 * r**3 end if + return +end function crater_profile + + +function crater_profile_find_r_inner_wall(user,crater) result(r_inner_wall) + !Uses Brent's method to solve for the point in the crater profile where the inner wall of the crater profile meets the original surface + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_profile_find_r_inner_wall + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + + ! Result variable + real(DP) :: r_inner_wall + + integer(I4B),parameter :: maxIterations=100 + real(DP) :: valueAtRoot + + real(DP),parameter :: Tolerance = 1e-8_DP + real(DP),parameter :: FPP = 1.e-11_DP + real(DP),parameter :: nearzero = 1.e-20_DP + + integer :: niter + real(DP) :: resultat,AA,BB,CC,DD,EE,FA,FB,FC,Tol1,PP,QQ,RR,SS,xm,x1,x2 + integer :: i, done + + i = 0 + done = 0 + + x1 = 0.0_DP + x2 = 1.0_DP + + AA = x1 + BB = x2 + FA = crater_profile(user,crater,AA) + FB = crater_profile(user,crater,BB) + + FC = FB; + do while (done.eq.0.and.i < maxIterations) + if (.not.util_rootbracketed(FC,FB)) then + CC = AA; FC = FA; DD = BB - AA; EE = DD + endif + if (abs(FC) < abs(FB)) then + AA = BB; BB = CC; CC = AA + FA = FB; FB = FC; FC = FA + endif + Tol1 = 2 * FPP * abs(BB) + 0.5_DP * Tolerance + xm = 0.5_DP * (CC-BB) + if ((abs(xm) <= Tol1).or.(abs(FA) < nearzero)) then + ! A root has been found + resultat = BB; + done = 1 + valueAtRoot = crater_profile(user,crater,resultat) + else + if ((abs(EE) >= Tol1).and.(abs(FA) > abs(FB))) then + SS = FB/ FA; + if (abs(AA - CC) < nearzero) then + PP = 2 * xm * SS; + QQ = 1._DP - SS; + else + QQ = FA/FC; + RR = FB /FC; + PP = SS * (2 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1._DP)); + QQ = (QQ - 1._DP) * (RR - 1._DP) * (SS - 1._DP); + endif + if (PP > nearzero) QQ = -QQ; + PP = abs(PP); + if ((2 * PP) < min(3 * xm * QQ - abs(Tol1 * QQ),abs( EE * QQ))) then + EE = DD; DD = PP/QQ; + else + DD = xm; EE = DD; + endif + else + DD = xm; + EE = DD; + endif + AA = BB; + FA = FB; + if (abs(DD) > Tol1) then + BB = BB + DD; + else + if (xm > 0._DP) then + BB = BB + abs(Tol1) + else + BB = BB - abs(Tol1) + endif + endif + FB = crater_profile(user,crater,BB) + i=i+1 + endif + end do + if (i >= maxIterations) then + write(*,*) 'Too many iterations in profile solver.' + r_inner_wall = 1._DP + return + end if + niter = i + r_inner_wall = resultat return -end subroutine crater_profile +end function crater_profile_find_r_inner_wall diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index b2e96fe7..207c6940 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -106,7 +106,7 @@ end subroutine complex_rim call complex_rim(user,surf,crater,rn,deltaMtot) !call complex_terrace(user,surf,crater,rn,deltaMtot) - call crater_slope_collapse(user,surf,crater,domain,(0.35_DP * user%pix)**2,deltaMtot) + !call crater_slope_collapse(user,surf,crater,domain,(0.35_DP * user%pix)**2,deltaMtot) !call complex_wall_texture(user,surf,crater,rn,deltaMtot) call complex_floor(user,surf,crater,rn,deltaMtot) call complex_peak(user,surf,crater,rn,deltaMtot) @@ -343,22 +343,22 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) real(DP) :: xynoise, znoise real(DP) :: noise,dnoise - ! Complex crater wall - integer(I4B) :: offset ! Scales the random xy-offset so that each crater's random noise is unique - real(DP) :: xy_noise_fac ! Spatial "size" of noise features at the first octave - real(DP) :: freq ! Spatial size scale factor multiplier at each octave level - real(DP) :: pers ! The relative size scaling at each octave level - real(DP) :: noise_height ! Vertical height of noise features as a function of crater radius at the first octave + ! Complex crater rim + integer(I4B) :: nscallops ! Approximate number of scallop features on edge of crater + integer(I4B),parameter :: offset = 3000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP),parameter :: noise_height = 0.15_DP ! Vertical height of noise features as a function of crater radius at the first octave + ! Higher values make the scallop features broader + real(DP),parameter :: noise_slope = 0.6_DP ! The slope of the power law function that defines the noise shape for scallop features + ! Higher values makes the inner sides of the scallop features more rounded, lower values + ! make them have sharper points at the inflections + + nscallops = 16 ! determine area to effect ! First make the interior of the crater rad = 1.25_DP * crater%frad ! Create diffusion noise for scalloped rim - offset = 3000 - xy_noise_fac = 5.00_DP - noise_height = 1.00_DP - inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1 crater%maxinc = max(crater%maxinc,inc) @@ -378,28 +378,22 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) r = sqrt(xbar**2 + ybar**2) / crater%frad ! Make scalloped rim - znoise = noise_height * crater%fcrat - xynoise = xy_noise_fac / crater%fcrat + znoise = noise_height**(1._DP / (2 * noise_slope)) + xynoise = (nscallops / PI) / crater%fcrat dnoise = util_perlin_noise(xynoise * xbar + offset * rn(1), & xynoise * ybar + offset * rn(2)) * znoise - noise = sqrt(sqrt(dnoise**2)) + noise = (dnoise**2)**noise_slope + hprof = r**(-1) + tprof = 0.5_DP * crater%ejrim + crater%melev - if (r >= 1.0) then - hprof = crater%ejrim * r**(-3) + crater%melev - tprof = 0.5_DP * crater%ejrim + crater%melev !* ((2.0_DP - r)**(-4) - 1.0_DP) + crater%melev - else - hprof = crater%ejrim * (2.0_DP - r)**(-2) + crater%melev - tprof = 0.5_DP * crater%ejrim + crater%melev ! (r**(-4) - 1.0_DP) + crater%melev - end if - - if (crater%melev + crater%ejrim - noise < hprof) then - newdem = tprof - noise + if (1.0_DP - noise < hprof) then + newdem = min(tprof,newdem) + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem end if - elchange = newdem - surf(xpi,ypi)%dem - deltaMtot = deltaMtot + elchange - surf(xpi,ypi)%dem = newdem end do end do diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index a3b37775..3db27529 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -17,15 +17,15 @@ module module_crater public save -interface - subroutine crater_mass_conservation(user,surf,crater) - use module_globals - implicit none - type(usertype),intent(in) :: user - type(surftype),dimension(:,:),intent(inout) :: surf - type(cratertype),intent(in) :: crater - end subroutine crater_mass_conservation -end interface + interface + subroutine crater_mass_conservation(user,surf,crater) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(in) :: crater + end subroutine crater_mass_conservation + end interface interface @@ -313,4 +313,26 @@ end function crater_visibility end interface + interface + function crater_profile(user,crater,r) result(h) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r + real(DP) :: h + end function crater_profile + end interface + + + interface + function crater_profile_find_r_inner_wall(user,crater) result(r_inner_wall) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP) :: r_inner_wall + end function crater_profile_find_r_inner_wall + end interface + end module diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index 5d5dd1ca..8f4aa22e 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -20,6 +20,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) use module_globals use module_util use module_regolith + use module_crater use module_ejecta, EXCEPT_THIS_ONE => ejecta_table_define implicit none @@ -33,23 +34,12 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! Internal variables integer(I4B) :: k - real(DP) :: erad,eradold,thick,vejsq,ejang,lrad + real(DP) :: erad,eradold,thick,vejsq,ejang,lrad,r logical :: firstrun ! Regotrack internal variables real(DP) :: rmelt,depthb,dimp,vimp - ! This will be replaced by its own function - real(DP) :: c0,c1,c2,c3,flrad,rh,fld,r - rh = crater%rimheight - fld = -crater%floordepth - flrad = 0.5_DP * crater%floordiam / crater%frad - c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP) - c0 = rh - (7._DP / 6._DP) * c1 - c2 = c1 / 3._DP - c3 = -c2 / 2._DP - !^^^^^^^^^^^^^^^ - ! Executable code ! Get estimate of size of ejb table @@ -90,8 +80,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) if (lrad >= crater%frad) then thick = crater%rimheight * r**(-3.0_DP) else - !thick = 0.14_DP * crater%frad**(0.74_DP) / (crater%frad - crater%ejrad) * (lrad - crater%ejrad) - thick = c0 + c1 * r + c2 * r**2 + c3 * r**3 + thick = max(crater_profile(user,crater,r),VSMALL) end if ejb(k)%thick = log(thick) ejb(k)%vesq = vejsq