From 3b7b0ad6e22d49d9d09a194c70297c3be1ff4dd0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 12 Mar 2020 16:10:32 -0400 Subject: [PATCH] Added template code to generate realistic topography. --- src/Makefile.am | 5 +- src/crater/crater_emplace.f90 | 1 + src/crater/crater_form_interior.f90 | 58 ++++--- src/crater/crater_populate.f90 | 4 + src/crater/crater_realistic_topography.f90 | 166 +++++++++++++++++++++ src/crater/module_crater.f90 | 12 ++ src/globals/module_globals.f90 | 37 ++--- src/io/io_input.f90 | 7 + src/util/util_perlin_noise.f90 | 26 ++-- 9 files changed, 265 insertions(+), 51 deletions(-) create mode 100644 src/crater/crater_realistic_topography.f90 diff --git a/src/Makefile.am b/src/Makefile.am index 6508f5c7..de726747 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -6,9 +6,9 @@ bin_PROGRAMS = CTEM #AM_FCFLAGS = -O3 -p -g -openmp -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR #gfortran optimized flags -#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace +AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace #gfortran debug flags -AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form +#AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form CTEM_SOURCES = globals/module_globals.f90\ util/module_util.f90\ @@ -71,6 +71,7 @@ crater/crater_generate.f90\ crater/crater_find_visible.f90\ crater/crater_averages.f90\ crater/crater_emplace.f90\ +crater/crater_realistic_topography.f90\ crater/crater_form_interior.f90\ crater/crater_form_exterior.f90\ crater/crater_record.f90\ diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index 79b91a44..a05c5938 100644 --- a/src/crater/crater_emplace.f90 +++ b/src/crater/crater_emplace.f90 @@ -71,6 +71,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot) ! Executable code + ! determine area to effect ! First make the interior of the crater inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1 diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index b9f6d4ae..889b8d01 100644 --- a/src/crater/crater_form_interior.f90 +++ b/src/crater/crater_form_interior.f90 @@ -33,40 +33,60 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele real(DP),intent(out) :: deltaMi ! Internal variables - real(DP) :: cform,newdem,elchange,pikeD,r + real(DP) :: cform,newdem,elchange,pikeD,r, Hpeak integer(I4B) :: layer - ! A list for poped data + ! A list for popped data type(regolisttype),pointer :: poppedlist + + ! Empirical crater shape parameters from Fassett et al. (2014) + real(DP),parameter :: r_floor = 0.2_DP + real(DP),parameter :: simple_depth_diam = 0.181_DP + real(DP),parameter :: r_rim = 0.98_DP + real(DP),parameter :: inner_c0 = -0.229_DP + real(DP),parameter :: inner_c1 = 0.228_DP + real(DP),parameter :: inner_c2 = 0.083_DP + real(DP),parameter :: inner_c3 = -0.039_DP + + real(DP),parameter :: outer_c0 = 0.188_DP + real(DP),parameter :: outer_c1 = -0.187_DP + real(DP),parameter :: outer_c2 = 0.018_DP + real(DP),parameter :: outer_c3 = 0.015_DP + + ! Executable code !change digital elevation map r = sqrt(x_relative**2+y_relative**2) / crater%frad ! Use empirical crater form from Fassett et al. 2014 - if (r < 0.2_DP) then - cform = -0.181_DP * crater%fcrat - else if (r < 0.98_DP) then - cform = (-0.229_DP + 0.228_DP * r + 0.083_DP * r**2 - 0.039_DP * r**3) * crater%fcrat + if (r < r_floor) then + cform = -simple_depth_diam * crater%fcrat + else if (r < r_rim) then + cform = (inner_c0 + inner_c1 * r + inner_c2 * r**2 + inner_c3 * r**3) * crater%fcrat else - cform = (0.188_DP - 0.187_DP * r + 0.018_DP * r**2 + 0.015_DP * r**3) * crater%fcrat + cform = (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat end if newdem = newelev + cform - pikeD = 1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP) ! Pike (1977) - if ((crater%fcrat > crater%cxtran * 2) .and. newdem < (crater%melev - pikeD)) then - newdem = crater%melev - pikeD ! Flatten out the bottom of the crater - !write(*,*) x_relative,y_relative,util_perlin_noise(x_relative,y_relative)*1e3_DP - !newdem = newdem + max(util_perlin_noise(x_relative/crater%fcrat,y_relative/crater%fcrat)*1e3_DP,0.0_DP) - !newdem = newdem + abs(util_perlin_noise(x_relative,y_relative)*1e3_DP) - end if - if (newdem < (crater%melev - user%deplimit)) then - newdem = crater%melev - user%deplimit ! Flatten out the bottom of the crater - do layer = 1,user%numlayers ! Remove all pre-existing craters from this current pixel - call util_remove_from_layer(surfi,layer) - end do + pikeD = min(1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP), user%deplimit) ! Pike (1977) depth/diameter ratio of complex craters + if (crater%fcrat > crater%cxtran * 2) then + ! Make this crater complex + + !if (r < r_rim) then + if (newdem < (crater%melev - pikeD)) then ! Flatten out the bottom of the crater + do layer = 1,user%numlayers ! Remove all pre-existing craters from this current pixel + call util_remove_from_layer(surfi,layer) + end do + newdem = crater%melev - pikeD + + + end if + !end if + end if + newdem = min(newdem,surfi%dem) ! Only allow excavation, no deposition elchange = newdem - surfi%dem deltaMi = elchange diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index e0f60bb9..969ecd0e 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -195,6 +195,10 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Place crater onto the surface call crater_emplace(user,surf,crater,domain,ejbmass) + + + if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejbmass) + call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 new file mode 100644 index 00000000..9fd12825 --- /dev/null +++ b/src/crater/crater_realistic_topography.f90 @@ -0,0 +1,166 @@ +!***** crater/crater_realistic_topography +! Name +! crater_realistic_topography +! SYNOPSIS +! This uses +! * module_globals +! * module_util +! * module_crater +! +! call crater_realistic_topography(user,surf,crater,domain,deltaMtot) +! +! DESCRIPTION +! Creates realistic topography for each fresh crater using Perlin noise. +! +! +! ARGUMENTS +! Input +! * user -- User input parameters +! * surf -- Surface grid +! * crater -- Crater dimension container +! * domain -- Simulation domain variable container +! * deltaMtot -- Total displaced mass used to calculate mass-conserving ejecta blanket +! +! Output +! * surf -- Outputs the new ejecta blanket onto the grid +! * crater -- May affects the value of the maximum affected distance +! +! Notes +! +!********************************************************************************************************************************** +subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(out) :: deltaMtot + + ! Internal variables + real(DP) :: lradsq,newelev, x_relative, y_relative + integer(I4B) :: xpi,ypi,i,j,inc,incsq,iradsq + real(DP) :: xp,yp,fradsq,deltaMi,rimheight + logical :: lastloop + real(DP), dimension(2) :: rn + + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + + + ! Complex crater floor + integer(I4B), parameter :: complex_floor_num_octaves = 8 ! Number of Perlin noise octaves + integer(I4B), parameter :: complex_floor_offset = 10000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP), parameter :: complex_floor_xy_noise_fac = 2.0_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: complex_floor_noise_height = 1e-3_DP ! Vertical height of noise features as a function of crater radius at the first octave + real(DP), parameter :: complex_floor_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level + real(DP), parameter :: complex_floor_pers = 0.80_DP ! The relative size scaling at each octave level + + ! Complex crater peak + integer(I4B), parameter :: complex_peak_num_octaves = 8 ! Number of Perlin noise octaves + integer(I4B), parameter :: complex_peak_offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP), parameter :: complex_peak_rad = 0.20_DP ! Radial size of central peak relative to frad + real(DP), parameter :: complex_peak_xy_noise_fac = 16.0_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: complex_peak_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level + real(DP), parameter :: complex_peak_pers = 0.80_DP ! The relative size scaling at each octave level + + ! Complex crater wall + integer(I4B), parameter :: complex_wall_num_octaves = 3 ! Number of Perlin noise octaves + integer(I4B), parameter :: complex_wall_offset = 2000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP), parameter :: complex_wall_xy_noise_fac = 2.0_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: complex_wall_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level + real(DP), parameter :: complex_wall_pers = 0.80_DP ! The relative size scaling at each octave level + real(DP), parameter :: complex_wall_noise_height = 0.10_DP ! Vertical height of noise features as a function of crater radius at the first octave + real(DP), parameter :: complex_wall_slope = 25._DP * DEG2RAD + integer(I4B), parameter :: complex_wall_terracefac = 16 ! Spatial size factor for terraces relative to crater radius + integer(I4B) :: complex_wall_terrace_num + + + ! Executable code + + call random_number(rn) + + ! determine area to effect + ! First make the interior of the crater + inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1 + crater%maxinc = max(crater%maxinc,inc) + fradsq = crater%frad**2 + deltaMtot = 0.0_DP + incsq = inc**2 + + do j=-inc,inc ! Do the loop in pixel space + do i=-inc,inc + iradsq = i**2 + j**2 + ! find distance from crater center + + ! find elevation and grid point + newelev = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + xp = xpi * user%pix + yp = ypi * user%pix + + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + x_relative = (crater%xl - xp) + y_relative = (crater%yl - yp) + + lradsq = x_relative**2 + y_relative**2 + + if (lradsq > crater%frad**2) cycle + + + + + deltaMtot = deltaMtot + deltaMi + + end do + end do + + +! ! Make the floor topographically "noisy" +! noise = 0.0_DP +! do octave = 1, complex_floor_num_octaves +! xynoise = complex_floor_xy_noise_fac * complex_floor_freq ** octave / crater%fcrat +! znoise = complex_floor_noise_height * complex_floor_pers ** (octave - 1) * crater%fcrat +! noise = noise + util_perlin_noise(xynoise * x_relative + complex_floor_offset * rn(1), & +! xynoise * y_relative + complex_floor_offset * rn(2)) * znoise +! end do +! newdem = newdem + max(noise,0.0_DP) +! +! ! Add in "noisy" central peak +! if (r < complex_peak_rad) then +! Hpeak = 0.032e3_DP * (crater%fcrat * 1e-3_DP)**(0.900_DP) ! Pike (1977) +! noise = 0.0_DP +! do octave = 1, complex_peak_num_octaves +! xynoise = complex_peak_xy_noise_fac * complex_peak_freq ** octave / crater%fcrat +! znoise = Hpeak * (1.0_DP - r / complex_peak_rad) * complex_peak_pers ** (octave - 1) +! noise = noise + util_perlin_noise(xynoise * x_relative + complex_peak_offset * rn(1), & +! xynoise * y_relative + complex_peak_offset * rn(2)) * znoise +! end do +! newdem = newdem + max(noise,0.0_DP) +! end if +! else ! Terraced walls +! complex_wall_terrace_num = int(complex_wall_terracefac * r) + 1 +! noise = 0.0_DP +! do octave = 1, complex_wall_num_octaves +! xynoise = complex_wall_xy_noise_fac * complex_wall_freq ** octave / crater%fcrat +! znoise = complex_wall_noise_height * complex_wall_pers ** (octave - 1) * crater%fcrat +! noise = noise + util_perlin_noise(xynoise * x_relative + complex_wall_terrace_num * complex_wall_offset * rn(1), & +! xynoise * y_relative + complex_wall_terrace_num * complex_wall_offset * rn(2)) & +! * znoise +! end do +! newdem = max(newdem + noise,crater%melev - pikeD) !/ cos(complex_wall_slope) + + + return +end subroutine crater_realistic_topography + diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 7e071b99..e006b621 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -107,6 +107,18 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot) end subroutine crater_emplace end interface + interface + subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(out) :: deltaMtot + end subroutine crater_realistic_topography + end interface + interface subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative, newelev,deltaMi) use module_globals diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index da48fb27..ae3571f0 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -150,10 +150,10 @@ module module_globals integer(I4B) :: gridsize ! Resolution integer(I4B) :: numlayers ! Number of perched layers real(DP) :: pix ! Pixel size (m) - real(DP) :: mu_b ! Crater scaling exponential constant (ignored for basins) - real(DP) :: kv_b ! Crater scaling linear constant - real(DP) :: mu_r ! Crater scaling exponential constant (ignored for basins) - real(DP) :: kv_r ! Crater scaling linear constant + real(DP) :: mu_b ! Crater scaling exponential constant (ignored for basins) + real(DP) :: kv_b ! Crater scaling linear constant + real(DP) :: mu_r ! Crater scaling exponential constant (ignored for basins) + real(DP) :: kv_r ! Crater scaling linear constant integer(I4B) :: seed ! Random number generator seed (only used in non-IDL driven mode) real(DP) :: trho_b ! Target bedrock density real(DP) :: trho_r ! Target surface regolith layer density @@ -165,23 +165,24 @@ module module_globals real(DP) :: prho ! Projectile density character(STRMAX) :: sfdfile ! Name of size distribution file character(STRMAX) :: velfile ! Name of velocity distribution file - real(DP) :: seisk,cohaccel ! seismic keff, cohesion breaking acceleration + real(DP) :: seisk,cohaccel ! seismic keff, cohesion breaking acceleration ! Optional input variables - logical :: docollapse ! Set T to use the slope collapse model (turning off speeds up the code for testing) - logical :: doangle ! Set to F to only do vertical impacts, otherwise do range of angles (default is T) - logical :: doporosity ! Porosity on/off flg. Set to F to turn the model off. Default F. - real(DP) :: basinimp ! Impactor size to switch to multiring basin - real(DP) :: maxcrat ! fraction that maximum crater can be relative to grid - real(DP) :: deplimit ! complex crater depth limit + logical :: docollapse ! Set T to use the slope collapse model (turning off speeds up the code for testing) + logical :: doangle ! Set to F to only do vertical impacts, otherwise do range of angles (default is T) + logical :: doporosity ! Porosity on/off flg. Set to F to turn the model off. Default F. + real(DP) :: basinimp ! Impactor size to switch to multiring basin + real(DP) :: maxcrat ! fraction that maximum crater can be relative to grid + real(DP) :: deplimit ! complex crater depth limit + logical :: dorealistic ! Set to T to enable realistic crater morphology. Default is F. ! Seismic input variables - logical :: doseismic ! Set to T if you want to do the seismic shaking model - real(DP) :: seisq ! Seismic energy attenuation quality factor (Q) - real(DP) :: neff ! impact seismic energy efficiency factor - real(DP) :: tvel ! target P-wave (body wave) speed (m/s) - real(DP) :: tfrac ! mean free path for seismic wave scattering in medium - real(DP) :: regcoh ! target surface regolith layer cohesion + logical :: doseismic ! Set to T if you want to do the seismic shaking model + real(DP) :: seisq ! Seismic energy attenuation quality factor (Q) + real(DP) :: neff ! impact seismic energy efficiency factor + real(DP) :: tvel ! target P-wave (body wave) speed (m/s) + real(DP) :: tfrac ! mean free path for seismic wave scattering in medium + real(DP) :: regcoh ! target surface regolith layer cohesion ! Crater diffusion input parameters real(DP) :: Kd1 ! Degradation function coefficient (from Minton et al. (2019)) @@ -230,6 +231,8 @@ module module_globals real(DP) :: shadedminh ! Minimum height for shaded relief map (m) real(DP) :: shadedmaxh ! Maximum height for shaded relief map (m) character(STRMAX) :: sfdcompare ! Type of run: 0 for normal, 1 for statistical (domain is reset between intervals) + + end type usertype ! Derived data type for the ejecta blanket table elements diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index c8fa61a6..b090134b 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -92,6 +92,7 @@ subroutine io_input(infile,user) user%rbreak = 0.1e3_DP user%fe = 5.0_DP user%ejecta_truncation = 5.0_DP + user%dorealistic = .false. open(unit=LUN,file=infile,status="old",iostat=ierr) if (ierr /= 0) then @@ -466,6 +467,12 @@ subroutine io_input(infile,user) token = line(ifirst:ilast) read(token, *) user%fe + case ("DOREALISTIC") + ifirst = ilast + 1 + call io_get_token(line, ilength, ifirst, ilast, ierr) + token = line(ifirst:ilast) + read(token, *) user%dorealistic + !************************************************************************** ! The following is for backwards compatibility with older style input files !************************************************************************** diff --git a/src/util/util_perlin_noise.f90 b/src/util/util_perlin_noise.f90 index b260d9a1..797458c5 100644 --- a/src/util/util_perlin_noise.f90 +++ b/src/util/util_perlin_noise.f90 @@ -79,22 +79,22 @@ function util_perlin_noise(xx,yy,zz) result (noise) w = fade(z) A = p(xi) + yi - B = p(xi+1) + yi + B = p(xi + 1) + yi AA = p(A) + zi BA = p(B) + zi - AB = p(A+1) + zi - BB = p(B+1) + zi + AB = p(A + 1) + zi + BB = p(B + 1) + zi - noise = lerp(w, lerp(v, lerp(u, grad(p(AA ), x , y, z),& - grad(p(BA ), x-1, y, z)),& - lerp(u, grad(p(AB ), x , y-1, z),& - grad(p(BB ), x-1, y-1, z))),& - lerp(v, lerp(u, grad(p(AA+1), x , y, z-1),& - grad(p(BA+1), x-1, y, z-1)),& - lerp(u, grad(p(AB+1), x , y-1, z-1),& - grad(p(BB+1), x-1, y-1, z-1)))) + noise = lerp(w, lerp(v, lerp(u, grad(p(AA ), x , y, z ),& + grad(p(BA ), x - 1, y, z )),& + lerp(u, grad(p(AB ), x , y - 1, z ),& + grad(p(BB ), x - 1, y - 1, z ))),& + lerp(v, lerp(u, grad(p(AA + 1), x , y , z - 1),& + grad(p(BA + 1), x - 1, y , z - 1)),& + lerp(u, grad(p(AB + 1), x , y - 1, z - 1),& + grad(p(BB + 1), x - 1, y - 1, z - 1)))) return contains function fade(t) result(f) @@ -135,12 +135,12 @@ function grad(hash,x,y,z) result(g) else v = z end if - if (iand(h,1) == 0) then + if (iand(h,1 ) == 0) then g = u else g = -u end if - if (iand(h,2) == 0) then + if (iand(h, 2) == 0) then g = g + v else g = g - v