diff --git a/examples/radial_slide/pyvista_plotting.py b/examples/radial_slide/pyvista_plotting.py new file mode 100644 index 0000000..d9895fb --- /dev/null +++ b/examples/radial_slide/pyvista_plotting.py @@ -0,0 +1,205 @@ +from pylab import * +import glob, os +import pyvista as pv +from clawpack.geoclaw import topotools +from clawpack.visclaw import animation_tools +from clawpack.geoclaw import fgout_tools +from clawpack.visclaw import colormaps + +warpfactor = 5 # vertical amplification of elevations in 3D plots + +outdir = '_output' +fgno = 1 # fgout grid number + +if 1: + # determine how many fgout frames are in outdir: + fgout_frames = glob.glob(os.path.join(outdir, \ + 'fgout%s.t*' % str(fgno).zfill(4))) + nfgout = len(fgout_frames) + fgframenos = array(range(1, nfgout+1)) + print('Found %i fgout frames' % nfgout) +else: + # specify which frames to use: + fgframenos = array(range(1,50,2)) + + +# where to save a png file for each frame, for constructing an animation: +framedir = '_frames' +os.system('mkdir -p %s' % framedir) +os.system('rm %s/*' % framedir) # remove frames from previous version + +window_size = (1200,1200) + +tan = [0.8,0.5,0.2] +light_blue = [0.0,1.0,1.0] +cmap_massfrac = colormaps.make_colormap({0:light_blue, 0.1:light_blue, + 1:tan}) + +def plot_topo(): + """ + plot topo from topofile + """ + topo = topotools.Topography('basal_topo.tt3') + z = array([0.]) + x = topo.x + y = -topo.y + X,Y,Z = meshgrid(x, y, z, indexing='ij') + topoxyz = pv.StructuredGrid(X,Y,Z) + + B = flipud(topo.Z) + B = fliplr(B) + topoxyz.point_data['B'] = B.flatten(order='C') + + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) + + p = pv.Plotter() + p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) + p.add_title('Topography') + p.show(window_size=window_size) + +def plot_topo_fgout(): + """ + plot topo from an fgout file + """ + fgno = 1 # which fgout grid + + fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw') + + fgout = fgout_grid.read_frame(1) + + x = fgout.x + y = fgout.y + z = array([0.]) + X,Y,Z = meshgrid(x, y, z, indexing='ij') + topoxyz = pv.StructuredGrid(X,Y,Z) + + B = fgout.B + topoxyz.point_data['B'] = B.flatten(order='C') + warpfactor = 5 # amplification of elevations + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) + + h = fgout.h + topoxyz.point_data['h'] = h.flatten(order='C') + warpfactor = 5 # amplification of elevations + hwarp = topoxyz.warp_by_scalar('h', factor=warpfactor) + + eta = fgout.eta + topoxyz.point_data['eta'] = B.flatten(order='C') + warpfactor = 5 # amplification of elevations + etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor) + + p = pv.Plotter() + p.window_size = window_size + p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) + #p.add_mesh(topowarp,color='g') + p.add_title('Topography fgout.B') + p.show() + + +def plot_fgout(make_animation=False): + """ + If make_animation == False, plot one fgout frame with a slider bar to + change frame number. Also prints the camera position after changing frame, + which you can copy and paste into this function in order to set the + camera position for making an animation. + + If make_animation == True, make an mp4 animation of all frames + specified by fgframenos set above. + """ + + global etamesh + + fgno = 1 # which fgout grid + warpfactor = 5 # amplification of elevations + + # Instantiate object for reading fgout frames: + fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw') + + fgout = fgout_grid.read_frame(1) + + x = fgout.x + y = fgout.y + z = array([0.]) + X,Y,Z = meshgrid(x, y, z, indexing='ij') + topoxyz = pv.StructuredGrid(X,Y,Z) + + B = fgout.B + topoxyz.point_data['B'] = B.flatten(order='C') + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) + + p = pv.Plotter(off_screen=make_animation, lighting='three lights') + p.window_size = window_size + #p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) + p.add_mesh(topowarp,color='lightgreen') + etamesh = None + + + def set_frameno(fgframeno): + global etamesh + fgframeno = int(round(fgframeno)) + fgout = fgout_grid.read_frame(fgframeno) + tsec = fgout.t + print('Frame %i, t = %.1f seconds' % (fgframeno, fgout.t)) + + # replace land surface by nan in eta so it only shows water: + # (big tolerance for landslide, would normally be smaller) + eta = where(fgout.h>1, fgout.eta, nan) + + topoxyz.point_data['eta'] = eta.flatten(order='C') + etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor) + if etamesh: + p.remove_actor(etamesh) + + # color the mesh based on mass fraction, using cmap_massfrac + h = fgout.h + massfrac = divide(fgout.hm, h, where=h>0, out=nan*ones(h.shape)) + etamesh = p.add_mesh(etawarp,scalars=massfrac, + colormap=cmap_massfrac, clim=(0,0.6)) + + p.add_title('Frame %i at time %.1f seconds' % (fgframeno,tsec)) + + if not make_animation: + # print camera position so that this can be copied and pasted + # into this script after adjusting (and then sliding frameno) + print('p.camera_position = ', p.camera_position) + + + # initial camera position: + p.camera_position = [ + (9034.796206317897, -2484.6171987620633, 3703.6128928929047), + (1500.0, 1500.0, 950.89111328125), + (-0.29984440650740857, 0.08916816605502331, 0.949811755059182)] + + + if not make_animation: + fgfr1 = fgframenos[0] + fgfr2 = fgframenos[-1] + p.add_slider_widget(set_frameno, [fgfr1,fgfr2], + value=fgfr1, title='Frame', + pointa=(0.4,0.85), pointb=(0.9,0.85), color='blue', + slider_width=0.02, tube_width=0.005) + + p.show() + + else: + + # make a png file for each frame: + for fgframeno in fgframenos: + set_frameno(fgframeno) + fname_png = '%s/PyVistaFrame%s.png' \ + % (framedir, str(fgframeno).zfill(4)) + p.screenshot(fname_png) + print('Created ',fname_png) + + p.close() + + # combine png files into mp4 animation: + anim = animation_tools.make_anim(framedir, + fname_pattern='PyVistaFrame*.png') + fname_mp4 = 'radial_slide_animation.mp4' + animation_tools.make_mp4(anim, fname_mp4) + print('Created ',fname_mp4) + +if __name__=='__main__': + + plot_fgout(make_animation=False) diff --git a/examples/radial_slide/setrun.py b/examples/radial_slide/setrun.py index ffeb955..fdfe319 100644 --- a/examples/radial_slide/setrun.py +++ b/examples/radial_slide/setrun.py @@ -16,6 +16,8 @@ raise Exception("*** Must first set CLAW environment variable") from clawpack.amrclaw.data import FlagRegion +from clawpack.geoclaw import fgout_tools + #------------------------------ @@ -124,8 +126,8 @@ def setrun(claw_pkg='dclaw'): if clawdata.output_style==1: # Output nout frames at equally spaced times up to tfinal: - clawdata.num_output_times = 30 #240 - clawdata.tfinal = 120. #240. + clawdata.num_output_times = 10 #240 + clawdata.tfinal = 100. #240. clawdata.output_t0 = True # output at initial (or restart) time? elif clawdata.output_style == 2: @@ -280,7 +282,7 @@ def setrun(claw_pkg='dclaw'): amrdata = rundata.amrdata # max number of refinement levels: - amrdata.amr_levels_max = 1 + amrdata.amr_levels_max = 2 # List of refinement ratios at each level (length at least mxnest-1) # dx = dy = 2', 10", 2", 1/3": @@ -394,8 +396,8 @@ def setrun(claw_pkg='dclaw'): etafile = 'surface_topo.tt3' qinitdclaw_data.qinitfiles.append([3, 8, 1, 2, etafile]) - #mfile = 'mass_frac.tt3' - mfile = 'mass_frac.tt3' # with m0 = 0 below + mfile = 'mass_frac.tt3' + #mfile = 'mass_frac0.tt3' # with m0 = 0 below qinitdclaw_data.qinitfiles.append([3, 4, 1, 2, mfile]) #hfile = 'landslide_depth.tt3' @@ -460,7 +462,32 @@ def setrun(claw_pkg='dclaw'): #flowgrades_data.flowgrades.append([1.0e-6, 2, 1, 1]) #flowgrades_data.flowgrades.append([1.0e-6, 1, 1, 1]) - # ----- For developers ----- + + # == fgout_grids.data values == + # NEW IN v5.9.0 + # Set rundata.fgout_data.fgout_grids to be a list of + # objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + fgout_grids = rundata.fgout_data.fgout_grids # empty list initially + + fgout = fgout_tools.FGoutGrid() + fgout.fgno = 1 + fgout.point_style = 2 # will specify a 2d grid of points + #fgout.output_format = 'binary32' # 4-byte, float32 + fgout.output_format = 'ascii' # 4-byte, float32 + fgout.nx = 300 + fgout.ny = 300 + fgout.x1 = 0. # specify edges (fgout pts will be cell centers) + fgout.x2 = 3e3 + fgout.y1 = 0. + fgout.y2 = 3e3 + fgout.tstart = 0. + fgout.tend = 100. + fgout.nout = 101 + fgout.q_out_vars = [1,4,8] + fgout_grids.append(fgout) # written to fgout_grids.data + + + # ----- For developers ----- # Toggle debugging print statements: amrdata.dprint = False # print domain flags amrdata.eprint = False # print err est flags diff --git a/src/2d/dig/Makefile.dclaw b/src/2d/dig/Makefile.dclaw index a36a65a..1b13fc4 100644 --- a/src/2d/dig/Makefile.dclaw +++ b/src/2d/dig/Makefile.dclaw @@ -23,14 +23,14 @@ COMMON_MODULES += \ $(GEOLIB)/multilayer/multilayer_module.f90 \ $(GEOLIB)/friction_module.f90 \ $(DIGLIB)/refinement_module.f90 \ - $(GEOLIB)/adjointsup_module.f90 + $(GEOLIB)/adjointsup_module.f90 \ + $(GEOLIB)/fgout_module.f90 \ #list of common modules needed for dclaw codes COMMON_MODULES += \ $(DIGLIB)/digclaw_module.f90 \ $(DIGLIB)/qinit_module.f90 \ $(DIGLIB)/auxinit_module.f90 \ - $(DIGLIB)/fgout_module.f90 \ # list of source files from AMR library. diff --git a/src/2d/dig/amr2.f90 b/src/2d/dig/amr2.f90 index d253892..d70d23d 100644 --- a/src/2d/dig/amr2.f90 +++ b/src/2d/dig/amr2.f90 @@ -491,7 +491,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_auxinit() ! specifies file(s) for auxiliary variables - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call set_regions() ! Set refinement regions call set_dig(naux) call set_pinit() @@ -526,7 +526,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_auxinit() ! specifies file(s) for auxiliary variables - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call set_regions() ! Set refinement regions call set_gauges(rest, nvar, naux) ! Set gauge output call set_fgmax() diff --git a/src/2d/dig/fgout_module.f90 b/src/2d/dig/fgout_module.f90 deleted file mode 100644 index 0b1a733..0000000 --- a/src/2d/dig/fgout_module.f90 +++ /dev/null @@ -1,647 +0,0 @@ -module fgout_module - - implicit none - save - - ! Container for fixed grid data, geometry and output settings - type fgout_grid - ! Grid data - real(kind=8), pointer :: early(:,:,:) - real(kind=8), pointer :: late(:,:,:) - - ! Geometry - integer :: num_vars(2),mx,my,point_style,fgno,output_format - real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi - - ! Time Tracking and output types - integer :: num_output,next_output_index - real(kind=8) :: start_time,end_time,dt - - integer, allocatable :: output_frames(:) - real(kind=8), allocatable :: output_times(:) - end type fgout_grid - - - logical, private :: module_setup = .false. - - ! Fixed grid arrays and sizes - integer :: FGOUT_num_grids - type(fgout_grid), target, allocatable :: FGOUT_fgrids(:) - real(kind=8) :: FGOUT_tcfmax - real(kind=8), parameter :: FGOUT_ttol = 1.d-6 ! tolerance for times - - -contains - - - ! Setup routine that reads in the fixed grids data file and sets up the - ! appropriate data structures - - subroutine set_fgout(rest,fname) - - use amr_module, only: parmunit, tstart_thisrun - - implicit none - - ! Subroutine arguments - logical :: rest ! restart? - character(len=*), optional, intent(in) :: fname - - ! Local storage - integer, parameter :: unit = 7 - integer :: i,k - type(fgout_grid), pointer :: fg - real(kind=8) :: ts - - - if (.not.module_setup) then - - write(parmunit,*) ' ' - write(parmunit,*) '--------------------------------------------' - write(parmunit,*) 'SETFGOUT:' - write(parmunit,*) '-----------' - - ! Open data file - if (present(fname)) then - call opendatafile(unit,fname) - else - call opendatafile(unit,'fgout_grids.data') - endif - - ! Read in data - read(unit,'(i2)') FGOUT_num_grids - write(parmunit,*) ' mfgrids = ',FGOUT_num_grids - if (FGOUT_num_grids == 0) then - write(parmunit,*) ' No fixed grids specified for output' - return - endif - - ! Allocate fixed grids (not the data yet though) - allocate(FGOUT_fgrids(FGOUT_num_grids)) - - ! Read in data for each fixed grid - do i=1,FGOUT_num_grids - fg => FGOUT_fgrids(i) - ! Read in this grid's data - read(unit,*) fg%fgno - read(unit,*) fg%start_time - read(unit,*) fg%end_time - read(unit,*) fg%num_output - read(unit,*) fg%point_style - read(unit,*) fg%output_format - read(unit,*) fg%mx, fg%my - read(unit,*) fg%x_low, fg%y_low - read(unit,*) fg%x_hi, fg%y_hi - - allocate(fg%output_times(fg%num_output)) - allocate(fg%output_frames(fg%num_output)) - - ! Initialize next_output_index - ! (might be reset below in case of a restart) - fg%next_output_index = 1 - - if (fg%point_style .ne. 2) then - print *, 'set_fgout: ERROR, unrecognized point_style = ',\ - fg%point_style - endif - - ! Setup data for this grid - ! Set dtfg (the timestep length between outputs) for each grid - if (fg%end_time <= fg%start_time) then - if (fg%num_output > 1) then - print *,'set_fgout: ERROR for fixed grid', i - print *,'start_time <= end_time yet num_output > 1' - print *,'set end_time > start_time or set num_output = 1' - stop - else - fg%dt = 0.d0 - endif - else - if (fg%num_output < 2) then - print *,'set_fgout: ERROR for fixed grid', i - print *,'end_time > start_time, yet num_output = 1' - print *,'set num_output > 2' - stop - else - fg%dt = (fg%end_time - fg%start_time) & - / (fg%num_output - 1) - do k=1,fg%num_output - fg%output_times(k) = fg%start_time + (k-1)*fg%dt - if (rest) then - ! don't write initial time or earlier - ts = tstart_thisrun+FGOUT_ttol - else - ! do write initial time - ts = tstart_thisrun-FGOUT_ttol - endif - - if (fg%output_times(k) < ts) then - ! will not output this time in this run - ! (might have already be done when restarting) - fg%output_frames(k) = -2 - fg%next_output_index = k+1 - else - ! will be reset to frameno when this is written - fg%output_frames(k) = -1 - endif - enddo - endif - endif - - - - ! Set spatial intervals dx and dy on each grid - if (fg%mx > 1) then - !fg%dx = (fg%x_hi - fg%x_low) / (fg%mx - 1) ! points - fg%dx = (fg%x_hi - fg%x_low) / fg%mx ! cells - else if (fg%mx == 1) then - fg%dx = 0.d0 - else - print *,'set_fgout: ERROR for fixed grid', i - print *,'x grid points mx <= 0, set mx >= 1' - endif - - if (fg%my > 1) then - !fg%dy = (fg%y_hi - fg%y_low) / (fg%my - 1) ! points - fg%dy = (fg%y_hi - fg%y_low) / fg%my ! cells - else if (fg%my == 1) then - fg%dy = 0.d0 - else - print *,'set_fgout: ERROR for fixed grid', i - print *,'y grid points my <= 0, set my >= 1' - endif - - ! set the number of variables stored for each grid - ! this should be (the number of variables you want to write out + 1) - fg%num_vars(1) = 10 ! DIG different number - - ! for geoclaw num_vars is 6 - ! h, hu, hv, bathy, eta, time - ! so that interpolation can be done in time and space - - ! for dclaw num_vars is 10 - ! h, hu, hv, hm, pb, hchi, b_eroded, bathy, eta, time - - ! Allocate new fixed grid data array - allocate(fg%early(fg%num_vars(1), fg%mx,fg%my)) - fg%early = nan() - - allocate(fg%late(fg%num_vars(1), fg%mx,fg%my)) - fg%late = nan() - - enddo - close(unit) - - FGOUT_tcfmax=-1.d16 - - module_setup = .true. - end if - - end subroutine set_fgout - - - !=====================FGOUT_INTERP======================================= - ! This routine interpolates q and aux on a computational grid - ! to an fgout grid not necessarily aligned with the computational grid - ! using bilinear interpolation defined on computational grid - !======================================================================= - subroutine fgout_interp(fgrid_type,fgrid, & - t,q,meqn,mxc,myc,mbc,dxc,dyc,xlowc,ylowc, & - maux,aux) - - use geoclaw_module, only: dry_tolerance - implicit none - - ! Subroutine arguments - integer, intent(in) :: fgrid_type - type(fgout_grid), intent(inout) :: fgrid - integer, intent(in) :: meqn,mxc,myc,mbc,maux - real(kind=8), intent(in) :: t,dxc,dyc,xlowc,ylowc - real(kind=8), intent(in) :: q(meqn,1-mbc:mxc+mbc,1-mbc:myc+mbc) - real(kind=8), intent(in) :: aux(maux,1-mbc:mxc+mbc,1-mbc:myc+mbc) - - integer, parameter :: method = 0 ! interpolate in space? - - ! Indices - integer :: ifg,jfg,m,ic1,ic2,jc1,jc2 - integer :: bathy_index,eta_index - - ! Tolerances - real(kind=8) :: total_depth,depth_indicator,nan_check - - ! Geometry - real(kind=8) :: xfg,yfg,xc1,xc2,yc1,yc2,xhic,yhic - real(kind=8) :: geometry(4) - - real(kind=8) :: points(2,2), eta_tmp - - ! Work arrays for eta interpolation - real(kind=8) :: eta(2,2),h(2,2) - - - ! Alias to data in fixed grid - integer :: num_vars - real(kind=8), pointer :: fg_data(:,:,:) - - - ! Setup aliases for specific fixed grid - if (fgrid_type == 1) then - num_vars = fgrid%num_vars(1) - fg_data => fgrid%early - else if (fgrid_type == 2) then - num_vars = fgrid%num_vars(1) - fg_data => fgrid%late - else - write(6,*) '*** Unexpected fgrid_type = ', fgrid_type - stop - ! fgrid_type==3 is deprecated, use fgmax grids instead - endif - - xhic = xlowc + dxc*mxc - yhic = ylowc + dyc*myc - - ! Find indices of various quantities in the fgrid arrays - bathy_index = meqn + 1 - eta_index = meqn + 2 - - !write(59,*) '+++ ifg,jfg,eta,geometry at t = ',t - - ! Primary interpolation loops - do ifg=1,fgrid%mx - xfg=fgrid%x_low + (ifg-0.5d0)*fgrid%dx ! cell centers - do jfg=1,fgrid%my - yfg=fgrid%y_low + (jfg-0.5d0)*fgrid%dy ! cell centers - - ! Check to see if this coordinate is inside of this grid - if (.not.((xfg < xlowc.or.xfg > xhic) & - .or.(yfg < ylowc.or.yfg > yhic))) then - - ! find where xfg,yfg is in the computational grid and - ! compute the indices - ! (Note: may be subject to rounding error if fgout point - ! is right on a cell edge!) - ic1 = int((xfg - xlowc + dxc)/dxc) - jc1 = int((yfg - ylowc + dyc)/dyc) - - if (method == 0) then - - ! piecewise constant: take values from cell (ic1,jc1): - - forall (m=1:meqn) - fg_data(m,ifg,jfg) = q(m,ic1,jc1) - end forall - - fg_data(bathy_index,ifg,jfg) = aux(1,ic1,jc1) - - ! for pw constant we take B, h, eta from same cell, - ! so setting eta = h+B should be fine even near shore: - fg_data(eta_index,ifg,jfg) = fg_data(1,ifg,jfg) & - + fg_data(bathy_index,ifg,jfg) - - - else if (method == 1) then - - ! bilinear used to interpolate to xfg,yfg - ! (not recommended) - - ! define constant parts of bilinear - if (ic1.eq.mxc) ic1=mxc-1 - if (jc1.eq.myc) jc1=myc-1 - ic2 = ic1 + 1 - jc2 = jc1 + 1 - - xc1 = xlowc + dxc * (ic1 - 0.5d0) - yc1 = ylowc + dyc * (jc1 - 0.5d0) - xc2 = xlowc + dxc * (ic2 - 0.5d0) - yc2 = ylowc + dyc * (jc2 - 0.5d0) - - geometry = [(xfg - xc1) / dxc, & - (yfg - yc1) / dyc, & - (xfg - xc1) * (yfg - yc1) / (dxc*dyc), & - 1.d0] - - - ! Interpolate all conserved quantities and bathymetry - forall (m=1:meqn) - fg_data(m,ifg,jfg) = & - interpolate(q(m,ic1:ic2,jc1:jc2), geometry) - end forall - - fg_data(bathy_index,ifg,jfg) = & - interpolate(aux(1,ic1:ic2,jc1:jc2),geometry) - - - ! surface eta = h + B: - - ! Note that for pw bilinear interp there may - ! be a problem interpolating each separately since - ! interpolated h + interpolated B may be much larger - ! than eta should be offshore. - eta = q(1,ic1:ic2,jc1:jc2) + aux(1,ic1:ic2,jc1:jc2) - fg_data(eta_index,ifg,jfg) = interpolate(eta,geometry) - ! NEED TO FIX - endif - - - ! save the time this fgout point was computed: - fg_data(num_vars,ifg,jfg) = t - - - endif ! if fgout point is on this grid - enddo ! fgout y-coordinate loop - enddo ! fgout x-coordinte loop - - end subroutine fgout_interp - - - !================ fgout_write ========================================== - ! This routine interpolates in time and then outputs a grid at - ! time=out_time - ! - ! files now have the same format as frames produced by outval - !======================================================================= - subroutine fgout_write(fgrid,out_time,out_index) - - use geoclaw_module, only: dry_tolerance - implicit none - - ! Subroutine arguments - type(fgout_grid), intent(inout) :: fgrid - real(kind=8), intent(in) :: out_time - integer, intent(in) :: out_index - - ! I/O - integer, parameter :: unit = 87 - character(len=15) :: fg_filename - character(len=4) :: cfgno, cframeno - character(len=8) :: file_format - integer :: grid_number,ipos,idigit,out_number,columns - integer :: ifg,ifg1, iframe,iframe1 - - integer, parameter :: method = 0 ! interpolate in time? - - ! Output format strings - ! These are now the same as in outval for frame data, for compatibility - ! For fgout grids there is only a single grid (ngrids=1) - ! and we set AMR_level=0, naux=0, nghost=0 (so no extra cells in binary) - - character(len=*), parameter :: header_format = & - "(i6,' grid_number',/," // & - "i6,' AMR_level',/," // & - "i6,' mx',/," // & - "i6,' my',/" // & - "e26.16,' xlow', /, " // & - "e26.16,' ylow', /," // & - "e26.16,' dx', /," // & - "e26.16,' dy',/)" - - character(len=*), parameter :: t_file_format = "(e18.8,' time', /," // & - "i6,' meqn'/," // & - "i6,' ngrids'/," // & - "i6,' naux'/," // & - "i6,' ndim'/," // & - "i6,' nghost'/," // & - "a10,' format'/,/)" - - ! Other locals - integer :: i,j,m - real(kind=8) :: t0,tf,tau, qaug(10) ! DIG switch from 6 to 10 - real(kind=8), allocatable :: qeta(:,:,:) - real(kind=4), allocatable :: qeta4(:,:,:) - real(kind=8) :: h_early,h_late,topo_early,topo_late - - allocate(qeta(8, fgrid%mx, fgrid%my)) ! to store h,hu,hv,hm,pb,hchi,b_eroded,eta # DIG - - - ! Interpolate the grid in time, to the output time, using - ! the solution in fgrid1 and fgrid2, which represent the - ! solution on the fixed grid at the two nearest computational times - do j=1,fgrid%my - do i=1,fgrid%mx - - ! Check for small numbers - forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%early(m,i,j)) < 1d-90) - fgrid%early(m,i,j) = 0.d0 - end forall - - if (method == 0) then - - ! no interpolation in time, use solution from full step: - qaug = fgrid%early(:,i,j) - - ! note that CFL condition ==> waves can't move more than 1 - ! cell per time step on each level, so solution from nearest - ! full step should be correct to within a cell width - ! Better to use early than late since for refinement tracking - ! wave moving out into still water. - - else if (method == 1) then - - ! interpolate in time. May have problems near shore? - - ! Fetch times for interpolation, this is done per grid point - ! since each grid point may come from a different source - t0 = fgrid%early(fgrid%num_vars(1),i,j) - tf = fgrid%late(fgrid%num_vars(1),i,j) - tau = (out_time - t0) / (tf - t0) - - ! check for small values: - forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%late(m,i,j)) < 1d-90) - fgrid%late(m,i,j) = 0.d0 - end forall - - ! linear interpolation: - qaug = (1.d0-tau)*fgrid%early(:,i,j) + tau*fgrid%late(:,i,j) - - ! If resolution changed between early and late time, may be - ! problems near shore when interpolating B, h, eta - ! separately (at least in case when B changed and point - ! was dry at one time and wet the other). - ! Switch back to fgrid%early values, only in this case. - ! This is implemented below but not extensively tested. - - if (qaug(1) > 0.d0) then - topo_early = fgrid%early(8,i,j) ! DIG num_var - 2 (was 4 is now 8) - topo_late = fgrid%late(8,i,j) ! DIG num_var - 2 (was 4 is now 8) - if (topo_early .ne. topo_late) then - ! AMR resolution changed between times - h_early = fgrid%early(1,i,j) - h_late = fgrid%late(1,i,j) - if (((h_early < dry_tolerance) & - .and. (h_late >= dry_tolerance)) & - .or. ((h_late < dry_tolerance) & - .and. (h_early >= dry_tolerance))) then - ! point changed between wet and dry - qaug = fgrid%early(:,i,j) ! don't interpolate - endif - endif - endif - endif - - ! Output the conserved quantities and eta value - ! DIG change assignments. - qeta(1,i,j) = qaug(1) ! h - qeta(2,i,j) = qaug(2) ! hu - qeta(3,i,j) = qaug(3) ! hv - qeta(4,i,j) = qaug(4) ! hm - qeta(5,i,j) = qaug(5) ! pb - qeta(6,i,j) = qaug(6) ! hchi - qeta(7,i,j) = qaug(7) ! b_eroded - qeta(8,i,j) = qaug(9) ! eta - ! bathy is qaug(8), time is qaug(10) - - - enddo - enddo - - - ! Make the file names and open output files - cfgno = '0000' - ifg = fgrid%fgno - ifg1 = ifg - do ipos=4,1,-1 - idigit = mod(ifg1,10) - cfgno(ipos:ipos) = char(ichar('0') + idigit) - ifg1 = ifg1/10 - enddo - - cframeno = '0000' - iframe = out_index - iframe1 = iframe - do ipos=4,1,-1 - idigit = mod(iframe1,10) - cframeno(ipos:ipos) = char(ichar('0') + idigit) - iframe1 = iframe1/10 - enddo - - fg_filename = 'fgout' // cfgno // '.q' // cframeno - - open(unit,file=fg_filename,status='unknown',form='formatted') - - ! Determine number of columns that will be written out - columns = fgrid%num_vars(1) - 1 - if (fgrid%num_vars(2) > 1) then - columns = columns + 2 - endif - - !write(6,*) '+++ fgout out_time = ',out_time - !write(6,*) '+++ fgrid%num_vars: ',fgrid%num_vars(1),fgrid%num_vars(2) - - ! Write out header in .q file: - !write(unit,header_format) out_time,fgrid%mx,fgrid%my, & - ! fgrid%x_low,fgrid%y_low, fgrid%x_hi,fgrid%y_hi, columns - - write(unit,header_format) fgrid%fgno, 0, fgrid%mx,fgrid%my, & - fgrid%x_low,fgrid%y_low, fgrid%dx, fgrid%dy - - if (fgrid%output_format == 1) then - ! ascii output added to .q file: - do j=1,fgrid%my - do i=1,fgrid%mx - write(unit, "(50e26.16)") qeta(1,i,j),qeta(2,i,j), & - qeta(3,i,j),qeta(4,i,j), & - qeta(5,i,j),qeta(6,i,j), & - qeta(7,i,j),qeta(8,i,j) ! DIG add additional elements of qeta - enddo - write(unit,*) ' ' ! blank line required between rows - enddo - endif - - close(unit) - - if (fgrid%output_format == 3) then - ! binary output goes in .b file as full 8-byte (float64): - fg_filename = 'fgout' // cfgno // '.b' // cframeno - open(unit=unit, file=fg_filename, status="unknown", & - access='stream') - write(unit) qeta - close(unit) - else if (fgrid%output_format == 2) then - ! binary output goes in .b file as 4-byte (float32): - fg_filename = 'fgout' // cfgno // '.b' // cframeno - open(unit=unit, file=fg_filename, status="unknown", & - access='stream') - allocate(qeta4(8, fgrid%mx, fgrid%my)) ! for 4-byte binary output ! DIG fix to size 8*mx*my - qeta4 = real(qeta, kind=4) - write(unit) qeta4 - deallocate(qeta4) - close(unit) - endif - - deallocate(qeta) - - ! time info .t file: - - - if (fgrid%output_format == 1) then - file_format = 'ascii' - else if (fgrid%output_format == 2) then - file_format = 'binary32' - else if (fgrid%output_format == 3) then - file_format = 'binary64' - else - write(6,*) '*** unrecognized fgrid%output_format = ', & - fgrid%output_format - write(6,*) '*** should be ascii, binary32, or binary64' - endif - - fg_filename = 'fgout' // cfgno // '.t' // cframeno - open(unit=unit, file=fg_filename, status='unknown', form='formatted') - ! time, num_eqn+1, num_grids, num_aux, num_dim, num_ghost: - write(unit, t_file_format) out_time, 8, 1, 0, 2, 0, file_format ! DIG changed num_eqn+1 - close(unit) - - ! DIG. might want to eventually write out aux for entrainment? - ! This will *break* for entrainment in that our B is aux(1)-q(7) - ! Fixed grids not yet supported for entrainment. - - - print "(a,i4,a,i4,a,e18.8)",'Writing fgout grid #',fgrid%fgno, & - ' frame ',out_index,' at time =',out_time - - ! Index into qeta for binary output - ! Note that this implicitly assumes that we are outputting only h, hu, hv - ! and will not output more (change num_eqn parameter above) - - end subroutine fgout_write - - - ! ========================================================================= - ! Utility functions for this module - ! Returns back a NaN - - real(kind=8) function nan() - real(kind=8) dnan - integer inan(2) - equivalence (dnan,inan) - inan(1)=2147483647 - inan(2)=2147483647 - nan=dnan - end function nan - - ! Interpolation function (in space) - ! Given 4 points (points) and geometry from x,y,and cross terms - - real(kind=8) pure function interpolate(points,geometry) result(interpolant) - - implicit none - - ! Function signature - real(kind=8), intent(in) :: points(2,2) - real(kind=8), intent(in) :: geometry(4) - integer :: icell, jcell - - ! pw bilinear - ! This is set up as a dot product between the approrpriate terms in - ! the input data. This routine could be vectorized or a BLAS routine - ! used instead of the intrinsics to ensure that the fastest routine - ! possible is being used - interpolant = sum([points(2,1)-points(1,1), & - points(1,2)-points(1,1), & - points(1,1) + points(2,2) - (points(2,1) + points(1,2)), & - points(1,1)] * geometry) - - end function interpolate - - -end module fgout_module diff --git a/src/python/dclaw/data.py b/src/python/dclaw/data.py index 44df7e0..9207e95 100644 --- a/src/python/dclaw/data.py +++ b/src/python/dclaw/data.py @@ -156,10 +156,9 @@ class QinitDClawData(clawpack.clawutil.data.ClawData): depth * chi, the fraction of species 1 If specified this way, provide chi rather than h*chi TODO: should be chi_init_val rather than 0.5 by default - - - q7, b_eroded: depth of material removed by - entrainment (cannot be set in this way as it - must start as zero). + - q7, bdif: depth of material removed by + entrainment (cannot be set in this way as it + must start as zero). - q8, eta: h+b diff --git a/tests/geoclaw_bouss_radial_slide/Makefile b/tests/geoclaw_bouss_radial_slide/Makefile new file mode 100644 index 0000000..bc0f196 --- /dev/null +++ b/tests/geoclaw_bouss_radial_slide/Makefile @@ -0,0 +1,126 @@ +# Makefile for Clawpack code in this directory. + +# Version for Boussinesq solvers in GeoClaw, which requires PETSc and MPI +# See https://www.clawpack.org/bouss2d.html +# +# Execute this command to check that enviornment variables set properly: +# make check +# before compiling or executing the code! +# Then do `make new` if using the Boussinesq version for the first time. + +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use + +# These environment variables need to be set properly, usually in your +# shell, but they could be explicitly set here: +#PETSC_DIR=/full/path/to/petsc +#PETSC_ARCH=arch-darwin-c-opt + +# The file petscMPIoptions sets parameters for MPI and the iterative solver: +# PETSC_OPTIONS must be set as an environment variable! +# e.g. in bash: +# export PETSC_OPTIONS="-options_file $CLAW/geoclaw/examples/bouss/petscMPIoptions" + +ifndef PETSC_DIR + $(error PETSC_DIR not set) +endif + +ifndef PETSC_ARCH + $(error PETSC_ARCH not set) +endif + +ifndef PETSC_OPTIONS + PETSC_OPTIONS=MISSING + $(error PETSC_OPTIONS must be declared as environment variable) +endif + +# How many MPI processes to use: +BOUSS_MPI_PROCS ?= 6 + +EXE = $(PWD)/xgeoclaw +RUNEXE="${PETSC_DIR}/${PETSC_ARCH}/bin/mpiexec -n ${BOUSS_MPI_PROCS}" +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +FC = gfortran + +# Some compiler flags below are needed for PETSc +PETSC_INCLUDE = $(PETSC_DIR)/include $(PETSC_DIR)/$(PETSC_ARCH)/include +INCLUDE += $(PETSC_INCLUDE) +PETSC_LFLAGS = $(shell PKG_CONFIG_PATH=$(PETSC_DIR)/$(PETSC_ARCH)/lib/pkgconfig pkg-config --libs-only-L --libs-only-l PETSc) + +FFLAGS ?= -O -gno-strict-dwarf -fbounds-check -fopenmp \ + -std=legacy -ffpe-trap='invalid,overflow,zero' +FFLAGS += -DHAVE_PETSC -ffree-line-length-none +LFLAGS += $(PETSC_LFLAGS) -fopenmp + + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +# BOUSSLIB contains library code for Boussinesq solvers, +# AMRLIB and GEOLIB are standard libraries, defined in case you need to +# exclude some modele or source file from one of these. +BOUSSLIB = $(CLAW)/geoclaw/src/2d/bouss +AMRLIB = $(CLAW)/amrclaw/src/2d +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow + +# See this Makefile for the list of library routines used: +include $(BOUSSLIB)/Makefile.bouss + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + + +MODULES = \ + +SOURCES = \ + ./set_eta_init.f90 \ + + + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +.PHONY: check + +check: + @echo =================== + @echo CLAW = $(CLAW) + @echo OMP_NUM_THREADS = $(OMP_NUM_THREADS) + @echo BOUSS_MPI_PROCS = $(BOUSS_MPI_PROCS) + @env | grep PETSC_OPTIONS + @echo PETSC_DIR = $(PETSC_DIR) + @echo PETSC_ARCH = $(PETSC_ARCH) + @echo RUNEXE = $(RUNEXE) + @echo FFLAGS = $(FFLAGS) + @echo =================== diff --git a/tests/geoclaw_bouss_radial_slide/maketopo_radial.py b/tests/geoclaw_bouss_radial_slide/maketopo_radial.py new file mode 100644 index 0000000..2311c93 --- /dev/null +++ b/tests/geoclaw_bouss_radial_slide/maketopo_radial.py @@ -0,0 +1,161 @@ + +from pylab import * +from clawpack.geoclaw import topotools +from clawpack.visclaw import colormaps +from clawpack.geoclaw.data import Rearth # radius of earth +from clawpack.clawutil.data import ClawData +from scipy.interpolate import interp1d + + +# # Radially symmetric landslide (and tsunami) + +m0 = 0.63 # mass fraction in landslide + +r0 = 0. # left boundary (center of radially symmetric) +r1 = 1100. # top of landslide +r2 = 1300. +r3 = 1600. # toe of landslide +r4 = 2000. # start of flat +r5 = 5000.*sqrt(2) # right boundary for radial coordinate + +B4 = B5 = 0 # flat region +Bslope = 25 # degrees +B3 = (r4 - r3) * sin(Bslope * pi/180) +B1 = (r4 - r1) * sin(Bslope * pi/180) +B2 = (r4 - r2) * sin(Bslope * pi/180) +B0 = B1 +etaslope = 40 # degrees +eta2 = B1 + +rr = array([r0,r1,r2,r3,r4,r5]) +Br = array([B0,B1,B2,B3,B4,B5]) +etar = Br.copy() +etar[2] = eta2 + +def make_plots(): + figure(figsize=(10,3)) + fill_between(rr,Br,etar,color='brown') + plot(rr,Br,'g') + grid(True) + #axis('equal') + plot([r1,r1],[0,500],'k--') + text(r1,-20,'r1\n%.0f' % r1,va='top',ha='center') + plot([r2,r2],[0,500],'k--') + text(r2,-20,'r2\n%.0f' % r2,va='top',ha='center') + plot([r3,r3],[0,500],'k--') + text(r3,-20,'r3\n%.0f' % r3,va='top',ha='center') + plot([r4,r4],[0,500],'k--') + text(r4,-20,'r4\n%.0f' % r4,va='top',ha='center') + axis([0,3000,-100,500]) + fname = '1d_topo.png' + savefig(fname) + print('Created ',fname) + + basal = topotools.Topography('basal_topo.tt3',3) + basal.plot() + title('Basal topo') + fname = 'basal_topo.png' + savefig(fname) + print('Created ',fname) + + eta = topotools.Topography('surface_topo.tt3',3) + eta.plot() + title('Surface topo eta') + fname = 'surface_topo.png' + savefig(fname) + print('Created ',fname) + + h = eta.Z - basal.Z + figure() + pcolormesh(eta.X,eta.Y,h,cmap=colormaps.white_red) + axis('equal') + colorbar() + title('Landslide depth') + fname = 'landslide_depth.png' + savefig(fname) + print('Created ',fname) + + +landslide_depth = eta2 - B2 +print('Maximum landslide depth: %.2f m' % landslide_depth) + +#x1d, z1d = loadtxt('1d_radial/celledges.data',skiprows=1,unpack=True) +B1d_func = interp1d(rr, Br, bounds_error=False, fill_value = Br[-1]) + +def basal(x,y): + """ + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + z = B1d_func(d) + return z + +eta1d_func = interp1d(rr, etar, bounds_error=False, fill_value = etar[-1]) + +def eta(x,y): + """ + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + z = eta1d_func(d) + return z + + +m1d_func = interp1d(rr, etar, bounds_error=False, fill_value = 0.) + +def mfrac(x,y): + """ + mass fraction + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + eta = eta1d_func(d) + B = B1d_func(d) + mfrac = where(eta-B > 0, m0, 0.) + return mfrac + + +xylim2d = 4e3 + +def maketopo(): + """ + Output topography file for the entire domain + """ + nxpoints = 501 + nypoints = 501 + xlower= -xylim2d + xupper= xylim2d + ylower= -xylim2d + yupper= xylim2d + outfile= "basal_topo.tt3" + topotools.topo3writer(outfile,basal,xlower,xupper,ylower,yupper,nxpoints,nypoints) + outfile= "mass_frac.tt3" + topotools.topo3writer(outfile,mfrac,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +def make_surface(): + """ + Output surface topography file for the entire domain + (Could be for smaller region) + """ + nxpoints = 501 + nypoints = 501 + xlower= -xylim2d + xupper= xylim2d + ylower= -xylim2d + yupper= xylim2d + outfile= "surface_topo.tt3" + topotools.topo3writer(outfile,eta,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +if __name__=='__main__': + maketopo() + make_surface() + make_plots() diff --git a/tests/geoclaw_bouss_radial_slide/set_eta_init.f90 b/tests/geoclaw_bouss_radial_slide/set_eta_init.f90 new file mode 100644 index 0000000..9e2a2aa --- /dev/null +++ b/tests/geoclaw_bouss_radial_slide/set_eta_init.f90 @@ -0,0 +1,153 @@ + +subroutine set_eta_init(mbc,mx,my,xlow,ylow,dx,dy,t,veta) + + ! This routine is called only if variable_eta_init = .true. + ! The input is a single grid patch at time t, as specified by + ! mbc,mx,my,xlow,ylow,dx,dy,t + ! The output is an array + ! veta(1-mbc:mx+mbc,1-mbc:my+mbc) + ! with the desired initial surface elevation at all points on this patch. + + ! This is called by qinit and also by filpatch and filval when refining. + + ! This default version sets veta to the value sea_level everywhere and + ! then adjusts this based on any dtopo deformation that has occured + ! up to this time t, as determined by interpolation from the + ! (possibly time-dependent) ! dtopo files. + + ! There is also a commented-out section below indicating how you might set + ! a higher value in some region that contains an onshore lake, for example. + + + use topo_module + use geoclaw_module, only: sea_level + + implicit none + + ! Arguments + integer, intent(in) :: mbc,mx,my + real(kind=8), intent(in) :: xlow,ylow,dx,dy,t + real(kind=8), intent(inout) :: veta(1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Local + integer :: i,j,i1,i2,j1,j2,idtopo,jdtopo,kdtopo,m + integer(kind=8) :: index0_dtopowork, ij + real(kind=8) :: x,y,r2,r,r3,eta2,eta3,eta_slope + real(kind=8) :: x1_lake,x2_lake,y1_lake,y2_lake, lake_level + + veta = sea_level ! initialize to sea_level, update below at some (i,j) + lake_level = 30.d0 + + !r2 = 2100.d0**2 + r2 = 1300.d0 + eta2 = 380.35d0 + r3 = 1600.d0 + eta3 = 169.047d0 + eta_slope = (eta3 - eta2) / (r3 - r2) + + do j=1,my + y = ylow + (j-0.5d0)*dy + do i=1,mx + x = xlow + (i-0.5d0)*dx + r = sqrt(x**2 + y**2) + !veta(i,j) = 200.d0 * min(1.d0, (1.d0 - (r - 2.d3)/10.d0)) + if (r > 1100.d0 .and. r < r3) then + if (r < r2) then + veta(i,j) = eta2 + else + veta(i,j) = eta2 + (r-r2)*eta_slope + endif + endif + + !veta(i,j) = 200.d0 * min(1.d0, (1.d0 - (r - 2.d3)/10.d0)) + !if (x**2 + y**2 < r2) then + ! veta(i,j) = lake_level + !endif + enddo + enddo + + if (.false.) then + ! This illustrates how you might set a higher value than sea_level + ! in a rectangle containing an onshore lake, as an example. + ! Depending on the geometry, a simple rectangle might not suffice. + ! What's below has been used for Lake Ozette on the Washington coast. + lake_level = 11.d0 ! meters relative to 0 datum of topofiles + x1_lake = -124.67d0 + x2_lake = -124.58d0 + y1_lake = 48.03d0 + y2_lake = 48.16d0 + i1 = max(int(floor((x1_lake - xlow)/dx)), 1) + i2 = min(int(floor((x2_lake - xlow)/dx)), mx) + j1 = max(int(floor((y1_lake - ylow)/dy)), 1) + j2 = min(int(floor((y2_lake - ylow)/dy)), my) + + forall(i=i1:i2, j=j1:j2) + veta(i,j) = lake_level + end forall + + endif + + ! The code below adjusts veta values set above for any uplift or + ! subsidence found in this region at the current time t, which assumes + ! that the water surface moves with the grount initially. + + if (num_dtopo == 0) return + + do m=1,num_dtopo + if (t < t0dtopo(m)) cycle ! this dtopo isn't moving yet + + ! compute indices in patch that overlap with this dtopo, + ! cycling to the next dtopo array if we discover there's no overlap: + + i1 = max(int(floor((xlowdtopo(m) - xlow)/dx)), 1) + if (i1 >= mx) cycle + + i2 = min(int(floor((xhidtopo(m) - xlow)/dx)), mx) + if (i2 < 1) cycle + + j1 = max(int(floor((ylowdtopo(m) - ylow)/dy)), 1) + if (j1 >= my) cycle + + j2 = min(int(floor((yhidtopo(m) - ylow)/dy)), my) + if (j2 < 1) cycle + + ! There is some overlap of dtopo with this patch + ! Next figure out index into time-dependent dtopo based on t: + + if (mtdtopo(m) == 1) then + ! Special case: instantaneous displacement at one instant in time + kdtopo = 1 + else + kdtopo = int(floor((t-t0dtopo(m))/dtdtopo(m)))+1 + kdtopo = min(kdtopo,mtdtopo(m)) + kdtopo = max(kdtopo,1) + endif + + index0_dtopowork = i0dtopo(m) + int(kdtopo-1, 8)*int(mxdtopo(m), 8)*int(mydtopo(m), 8) + !write(6,*) '+++ index0_dtopowork = ',index0_dtopowork + + ! Adjust eta_init by dtopo on part of patch that overlaps dtopo. + ! Code below assumes dtopo is smooth enough that we can just + ! evaluate at one point in space and time, not doing interpolation. + ! This could be improved as in topo_update.f, but probably not + ! necessary (??). + + do i=i1,i2 + x = xlow + (i-0.5d0)*dx + do j=j1,j2 + y = ylow + (j-0.5d0)*dy + idtopo = int(floor((x-xlowdtopo(m))/dxdtopo(m))) + 1 + idtopo = max(1, min(mxdtopo(m)-1, idtopo)) + jdtopo = int(floor((yhidtopo(m)-y)/dydtopo(m))) + 1 + jdtopo = max(1, min(mydtopo(m)-1, jdtopo)) + ij = index0_dtopowork + int(jdtopo-1, 8)*int(mxdtopo(m) + idtopo-1, 8) + + veta(i,j) = veta(i,j) + dtopowork(ij) + + enddo + enddo + + enddo + return + +end subroutine set_eta_init diff --git a/tests/geoclaw_bouss_radial_slide/setplot.py b/tests/geoclaw_bouss_radial_slide/setplot.py new file mode 100644 index 0000000..4635623 --- /dev/null +++ b/tests/geoclaw_bouss_radial_slide/setplot.py @@ -0,0 +1,473 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +import numpy as np +import matplotlib.pyplot as plt +import cmocean +import matplotlib as mpl + +from clawpack.visclaw import geoplot +from clawpack.visclaw import colormaps +from clawpack.visclaw import gridtools + + +import os,sys + +sea_level = 200. + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + import clawpack.dclaw.plot as dplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + #plotdata.format = 'binary' + + + + def timeformat(t): + from numpy import mod + hours = int(t/3600.) + tmin = mod(t,3600.) + min = int(tmin/60.) + sec = int(mod(tmin,60.)) + timestr = '%s:%s:%s' % (hours,str(min).zfill(2),str(sec).zfill(2)) + return timestr + + def title_hours(current_data): + from pylab import title + t = current_data.t + timestr = timeformat(t) + title('t = %s' % timestr) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Computational domain', figno=0) + plotfigure.kwargs = {'figsize':(8,7)} + plotfigure.show = False + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + + def aa(current_data): + from pylab import ticklabel_format, xticks, gca, cos, pi, savefig + gca().set_aspect(1.) + title_hours(current_data) + ticklabel_format(useOffset=False) + xticks(rotation=20) + + + plotaxes.afteraxes = aa + + + # Hillshade + plotitem = plotaxes.new_plotitem(plot_type="2d_hillshade") + plotitem.show = False + plotitem.plot_var = dplot.eta + plotitem.add_colorbar = False + + # Surface + plotitem = plotaxes.new_plotitem(plot_type="2d_imshow") + plotitem.plot_var = dplot.surface_solid_frac_lt03 + plotitem.add_colorbar = True + plotitem.colorbar_kwargs = { + "shrink": 0.9, + "location": "bottom", + "orientation": "horizontal", + } + plotitem.colorbar_label = "Surface (m)" + plotitem.imshow_cmap = cmocean.cm.curl + plotitem.imshow_cmin = 90 + plotitem.imshow_cmax = 110 + # Debris + plotitem = plotaxes.new_plotitem(plot_type="2d_imshow") + plotitem.plot_var = dplot.solid_frac_gt03 + plotitem.imshow_cmap = cmocean.cm.turbid + plotitem.imshow_cmin = 0.3 + plotitem.imshow_cmax = 1 + plotitem.add_colorbar = False + + #----------------------------------------- + # Figure for water / landslide separately + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Landslide/water surface', figno=1) + plotfigure.figsize=(8,8) + plotfigure.facecolor='w' + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.15,.5,.7,.45])' + plotaxes.title = 'Landslide / Water' + plotaxes.scaled = True + plotaxes.xlimits = [-3e3,3e3] + plotaxes.ylimits = [-3e3,3e3] + + def pure_water(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + eta = dplot.eta(current_data) + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + water = np.where(np.logical_and(h>1e-3, m<0.1), + eta, np.nan) + #import pdb; pdb.set_trace() + return water + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + #plotitem.show = False + #plotitem.plot_var = pure_water + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = sea_level - 20. + plotitem.pcolor_cmax = sea_level + 20. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def landslide(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + eta = dplot.eta(current_data) + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + landslide = np.where(np.logical_and(h>1e-3, m>0.1), + h, np.nan) + #import pdb; pdb.set_trace() + return landslide + + # Landslide + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + plotitem.show = False + plotitem.plot_var = landslide + cmap_mass = colormaps.make_colormap({0.:'w', 1.:'brown'}) + plotitem.pcolor_cmap = cmap_mass + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 80. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def land(current_data): + """ + Return a masked array containing the surface elevation only in dry cells. + """ + drytol = 1e-3 + q = current_data.q + h = q[0,...] + eta = q[-1,...] + land = np.ma.masked_where(h>drytol, eta) + #import pdb; pdb.set_trace() + #land = eta - h # everywhere + return land + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figure for cross section compared to 1d_radial + #----------------------------------------- + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('radial slice') + plotaxes.axescmd = 'axes([.1,.1,.8,.3])' + plotaxes.title = 'Diagonal Transect' + #plotaxes.scaled = True + + def plot_xsec(current_data): + from pylab import plot,linspace,zeros,ones,legend,xlabel,\ + sqrt,grid,xlim,fill_between,logical_and + from pylab import nan,where,ylim,loadtxt,arange + from clawpack.pyclaw import Solution + pd = current_data.plotdata + frameno = current_data.frameno + framesoln = Solution(frameno, path=pd.outdir, file_format=pd.format) + + #xout = linspace(-3e3,3e3,1000) + #yout = ones(xout.shape) # near x-axis + #rout = xout + + xout = linspace(-3e3,3e3,1000) + yout = xout + rout = xout * sqrt(2) + + etaout = gridtools.grid_output_2d(framesoln, -1, xout, yout) + hout = gridtools.grid_output_2d(framesoln, 0, xout, yout) + zetaout = where(hout>0.001, etaout, nan) + Bout = etaout - hout + if 0: + hmout = gridtools.grid_output_2d(framesoln, 3, xout, yout) + with np.errstate(divide="ignore", invalid="ignore"): + mout = hmout / hout + water = where(mout<0.1, etaout, nan) + landslide = where(mout>0.5, etaout, nan) + mixed = where(logical_and(0.11e-3, m<0.1), + h, np.nan) + import pdb; pdb.set_trace() + return water + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + #plotitem.show = False + plotitem.plot_var = 0 + plotitem.pcolor_cmap = colormaps.white_red + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 100. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def land(current_data): + """ + Return a masked array containing the surface elevation only in dry cells. + """ + drytol = 1e-3 + q = current_data.q + h = q[0,...] + eta = q[-1,...] + land = np.ma.masked_where(h>drytol, eta) + #import pdb; pdb.set_trace() + #land = eta - h # everywhere + return land + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figure for mass fraction + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Mass Fraction', figno=6) + plotfigure.kwargs = {'figsize':(8,7)} + plotfigure.show = False + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'mass fraction' + plotaxes.scaled = True + + plotaxes.scaled = True + plotaxes.xlimits = [-3e3,3e3] + plotaxes.ylimits = [-3e3,3e3] + + def mass_frac(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + mwet = np.where(h > 0.01, m, np.nan) + return mwet + + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = mass_frac + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 0.65 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + + # Figure for scatter plot + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter', figno=3) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0, 3e3*np.sqrt(2)] + plotaxes.ylimits = 'auto' + plotaxes.title = 'Scatter plot' + plotaxes.grid = True + + # Set up for item on these axes: scatter of 2d data + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + + def h_vs_r(current_data): + # Return radius of each grid cell and p value in the cell + from pylab import sqrt + x = current_data.x + y = current_data.y + r = np.sqrt(x**2 + y**2) + q = current_data.q + h = q[0,:,:] + return r,h + + plotitem.map_2d_to_1d = h_vs_r + plotitem.plot_var = 0 + plotitem.plotstyle = '.' + plotitem.color = 'b' + plotitem.kwargs = {'markersize':1} + plotitem.show = True # show on plot? + + + # ----------------------- + # Figure for scatter plot of speed + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter_speed', figno=9) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0, 3e3*np.sqrt(2)] + plotaxes.ylimits = 'auto' + plotaxes.title = 'speed' + plotaxes.grid = True + + + def s_vs_r(current_data): + # Return radius of each grid cell and p value in the cell + from pylab import sqrt,nan,where + x = current_data.x + y = current_data.y + r = np.sqrt(x**2 + y**2) + q = current_data.q + s = sqrt(q[1,:,:]**2 + q[2,:,:]**2) + s = sqrt(q[1,:,:]**2 + q[2,:,:]**2) + h = where(q[0,:,:]>1e-2, q[0,:,:], nan) + s = s/h + return r,s + + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + plotitem.map_2d_to_1d = s_vs_r + #plotitem.plot_var = 0 + plotitem.plotstyle = '+' + plotitem.color = 'r' + #plotitem.kwargs = {'markersize':1} + plotitem.show = True # show on plot? + + #------------------------------------- + + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + import os + from clawpack.visclaw import plot_timing_stats + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + units = {'comptime':'minutes', 'simtime':'minutes', 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, make_pngs=True, + plotdir=timing_plotdir, units=units) + os.system('cp %s/timing.* %s' % (plotdata.outdir, timing_plotdir)) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata diff --git a/tests/geoclaw_bouss_radial_slide/setrun.py b/tests/geoclaw_bouss_radial_slide/setrun.py new file mode 100644 index 0000000..fd3dbc5 --- /dev/null +++ b/tests/geoclaw_bouss_radial_slide/setrun.py @@ -0,0 +1,502 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os, sys +import numpy as np + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW environment variable") + +from clawpack.amrclaw.data import FlagRegion + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + #probdata.add_param('variable_eta_init', True) # now in qinit info + + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -3e3 + clawdata.upper[0] = 3e3 + + clawdata.lower[1] = -3e3 + clawdata.upper[1] = 3e3 + + + # Number of grid cells: Coarsest grid + clawdata.num_cells[0] = 240 + clawdata.num_cells[1] = 240 + + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 5 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 1 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 40 #240 + clawdata.tfinal = 160. #240. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [0.5, 1.0] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 10 + clawdata.output_t0 = True + + + clawdata.output_format = 'ascii' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = True # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.0001 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + # D-Claw requires CFL<0.5 + clawdata.cfl_desired = 0.85 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = [4, 4, 4] # TODO VERIFY THAT 4 in old and new are the same + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + # TODO This is not in old setrun.py + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + amrdata.max1d = 300 + + # max number of refinement levels: + amrdata.amr_levels_max = 1 + + # List of refinement ratios at each level (length at least mxnest-1) + # dx = dy = 2', 10", 2", 1/3": + amrdata.refinement_ratios_x = [2,2] + amrdata.refinement_ratios_y = [2,2] + amrdata.refinement_ratios_t = [2,2] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = [ + "center"] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + # NO OLD STYLE REGIONS USED HERE + + + # --------------- + # NEW flagregions + # --------------- + + + flagregions = rundata.flagregiondata.flagregions # initialized to [] + + # now append as many flagregions as desired to this list: + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges = [] + + # Set GeoClaw specific runtime parameters. + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 1 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 200.0 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True # TODO change? + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.01 + + # == settopo.data values == + topofiles = rundata.topo_data.topofiles + # for topography, append lines of the form + # [topotype, fname] + topofiles.append([3, 'basal_topo.tt3']) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + + rundata.qinit_data.variable_eta_init = True + + if 0: + # == setqinit.data values == + qinitdclaw_data = rundata.qinitdclaw_data # initialized when rundata instantiated + + etafile = 'surface_topo.tt3' + qinitdclaw_data.qinitfiles.append([3, 8, 1, 2, etafile]) + + mfile = 'mass_frac.tt3' + qinitdclaw_data.qinitfiles.append([3, 4, 1, 2, mfile]) + + #hfile = 'landslide_depth.tt3' + #qinitdclaw_data.qinitfiles.append([3, 1, 1, 2, hfile]) + + # == setauxinit.data values == + #auxinitdclaw_data = rundata.auxinitdclaw_data # initialized when rundata instantiated + + # == fgmax.data values == + #fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + + # == setdclaw.data values == + dclaw_data = rundata.dclaw_data # initialized when rundata instantiated + + dclaw_data.c1 = 1.0 # do we want to remove this? + dclaw_data.rho_f = 1000.0 + dclaw_data.rho_s = 2700.0 + dclaw_data.phi_bed = 32.0 + dclaw_data.theta_input = 0.0 + dclaw_data.mu = 0.005 + dclaw_data.m0 = 0.63 + dclaw_data.m_crit = 0.64 + dclaw_data.kappita = 1.e-8 + #dclaw_data.kappita_diff = 1 + #dclaw_data.chi_init_val=0.5 # not currently used. + dclaw_data.alpha_c = 0.05 + dclaw_data.alpha_seg = 0.0 + #dclaw_data.phi_seg_coeff = 0.0 + dclaw_data.delta = 0.001 + dclaw_data.bed_normal = 0 + dclaw_data.entrainment = 0 + dclaw_data.entrainment_rate = 0.0 + dclaw_data.sigma_0 = 1.0e3 + #dclaw_data.mom_autostop = True + #dclaw_data.momlevel = 1 + #dclaw_data.mom_perc = 0.0 + + # == pinitdclaw.data values == + pinitdclaw_data = rundata.pinitdclaw_data # initialized when rundata instantiated + + pinitdclaw_data.init_ptype = 0 # hydrostatic (-1 ==> zero everywhere) + pinitdclaw_data.init_pmax_ratio = 0.00e0 + pinitdclaw_data.init_ptf = 0.0 + pinitdclaw_data.init_ptf2 = 0.0 + + # == flowgrades.data values == + flowgrades_data = rundata.flowgrades_data # initialized when rundata instantiated + + flowgrades_data.flowgrades = [] + # for using flowgrades for refinement append lines of the form + # [flowgradevalue, flowgradevariable, flowgradetype, flowgrademinlevel] + # where: + # flowgradevalue: floating point relevant flowgrade value for following measure: + # flowgradevariable: 1=depth, 2= momentum, 3 = sign(depth)*(depth+topo) (0 at sealevel or dry land). + # flowgradetype: 1 = norm(flowgradevariable), 2 = norm(grad(flowgradevariable)) + # flowgrademinlevel: refine to at least this level if flowgradevalue is exceeded. + + + #flowgrades_data.keep_fine = True + #flowgrades_data.flowgrades.append([1.0e-6, 2, 1, 1]) + #flowgrades_data.flowgrades.append([1.0e-6, 1, 1, 1]) + + + # To use Boussinesq solver, add bouss_data parameters here + # Also make sure to use the correct Makefile pointing to bouss version + # and set clawdata.num_eqn = 5 + + from clawpack.geoclaw.data import BoussData + rundata.add_data(BoussData(),'bouss_data') + + rundata.bouss_data.bouss_equations = 2 # 0=SWE, 1=MS, 2=SGN + rundata.bouss_data.bouss_min_level = 1 # coarsest level to apply bouss + rundata.bouss_data.bouss_max_level = 10 # finest level to apply bouss + rundata.bouss_data.bouss_min_depth = -100. # depth to switch to SWE + rundata.bouss_data.bouss_solver = 3 # 1=GMRES, 2=Pardiso, 3=PETSc + rundata.bouss_data.bouss_tstart = 0. # time to switch from SWE + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + + return rundata + + + + + # end of function setrun + # ---------------------- + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + diff --git a/tests/geoclaw_radial_slide/Makefile b/tests/geoclaw_radial_slide/Makefile new file mode 100644 index 0000000..15fb09c --- /dev/null +++ b/tests/geoclaw_radial_slide/Makefile @@ -0,0 +1,69 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + + +MODULES = \ + +SOURCES = \ + set_eta_init.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + $(CLAW_PYTHON) maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls + diff --git a/tests/geoclaw_radial_slide/maketopo_radial.py b/tests/geoclaw_radial_slide/maketopo_radial.py new file mode 100644 index 0000000..2311c93 --- /dev/null +++ b/tests/geoclaw_radial_slide/maketopo_radial.py @@ -0,0 +1,161 @@ + +from pylab import * +from clawpack.geoclaw import topotools +from clawpack.visclaw import colormaps +from clawpack.geoclaw.data import Rearth # radius of earth +from clawpack.clawutil.data import ClawData +from scipy.interpolate import interp1d + + +# # Radially symmetric landslide (and tsunami) + +m0 = 0.63 # mass fraction in landslide + +r0 = 0. # left boundary (center of radially symmetric) +r1 = 1100. # top of landslide +r2 = 1300. +r3 = 1600. # toe of landslide +r4 = 2000. # start of flat +r5 = 5000.*sqrt(2) # right boundary for radial coordinate + +B4 = B5 = 0 # flat region +Bslope = 25 # degrees +B3 = (r4 - r3) * sin(Bslope * pi/180) +B1 = (r4 - r1) * sin(Bslope * pi/180) +B2 = (r4 - r2) * sin(Bslope * pi/180) +B0 = B1 +etaslope = 40 # degrees +eta2 = B1 + +rr = array([r0,r1,r2,r3,r4,r5]) +Br = array([B0,B1,B2,B3,B4,B5]) +etar = Br.copy() +etar[2] = eta2 + +def make_plots(): + figure(figsize=(10,3)) + fill_between(rr,Br,etar,color='brown') + plot(rr,Br,'g') + grid(True) + #axis('equal') + plot([r1,r1],[0,500],'k--') + text(r1,-20,'r1\n%.0f' % r1,va='top',ha='center') + plot([r2,r2],[0,500],'k--') + text(r2,-20,'r2\n%.0f' % r2,va='top',ha='center') + plot([r3,r3],[0,500],'k--') + text(r3,-20,'r3\n%.0f' % r3,va='top',ha='center') + plot([r4,r4],[0,500],'k--') + text(r4,-20,'r4\n%.0f' % r4,va='top',ha='center') + axis([0,3000,-100,500]) + fname = '1d_topo.png' + savefig(fname) + print('Created ',fname) + + basal = topotools.Topography('basal_topo.tt3',3) + basal.plot() + title('Basal topo') + fname = 'basal_topo.png' + savefig(fname) + print('Created ',fname) + + eta = topotools.Topography('surface_topo.tt3',3) + eta.plot() + title('Surface topo eta') + fname = 'surface_topo.png' + savefig(fname) + print('Created ',fname) + + h = eta.Z - basal.Z + figure() + pcolormesh(eta.X,eta.Y,h,cmap=colormaps.white_red) + axis('equal') + colorbar() + title('Landslide depth') + fname = 'landslide_depth.png' + savefig(fname) + print('Created ',fname) + + +landslide_depth = eta2 - B2 +print('Maximum landslide depth: %.2f m' % landslide_depth) + +#x1d, z1d = loadtxt('1d_radial/celledges.data',skiprows=1,unpack=True) +B1d_func = interp1d(rr, Br, bounds_error=False, fill_value = Br[-1]) + +def basal(x,y): + """ + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + z = B1d_func(d) + return z + +eta1d_func = interp1d(rr, etar, bounds_error=False, fill_value = etar[-1]) + +def eta(x,y): + """ + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + z = eta1d_func(d) + return z + + +m1d_func = interp1d(rr, etar, bounds_error=False, fill_value = 0.) + +def mfrac(x,y): + """ + mass fraction + Cartesian: x,y in meters + """ + import numpy as np + x0 = 0. + y0 = 0. + d = np.sqrt((x-x0)**2 + (y-y0)**2) + eta = eta1d_func(d) + B = B1d_func(d) + mfrac = where(eta-B > 0, m0, 0.) + return mfrac + + +xylim2d = 4e3 + +def maketopo(): + """ + Output topography file for the entire domain + """ + nxpoints = 501 + nypoints = 501 + xlower= -xylim2d + xupper= xylim2d + ylower= -xylim2d + yupper= xylim2d + outfile= "basal_topo.tt3" + topotools.topo3writer(outfile,basal,xlower,xupper,ylower,yupper,nxpoints,nypoints) + outfile= "mass_frac.tt3" + topotools.topo3writer(outfile,mfrac,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +def make_surface(): + """ + Output surface topography file for the entire domain + (Could be for smaller region) + """ + nxpoints = 501 + nypoints = 501 + xlower= -xylim2d + xupper= xylim2d + ylower= -xylim2d + yupper= xylim2d + outfile= "surface_topo.tt3" + topotools.topo3writer(outfile,eta,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +if __name__=='__main__': + maketopo() + make_surface() + make_plots() diff --git a/tests/geoclaw_radial_slide/set_eta_init.f90 b/tests/geoclaw_radial_slide/set_eta_init.f90 new file mode 100644 index 0000000..9e2a2aa --- /dev/null +++ b/tests/geoclaw_radial_slide/set_eta_init.f90 @@ -0,0 +1,153 @@ + +subroutine set_eta_init(mbc,mx,my,xlow,ylow,dx,dy,t,veta) + + ! This routine is called only if variable_eta_init = .true. + ! The input is a single grid patch at time t, as specified by + ! mbc,mx,my,xlow,ylow,dx,dy,t + ! The output is an array + ! veta(1-mbc:mx+mbc,1-mbc:my+mbc) + ! with the desired initial surface elevation at all points on this patch. + + ! This is called by qinit and also by filpatch and filval when refining. + + ! This default version sets veta to the value sea_level everywhere and + ! then adjusts this based on any dtopo deformation that has occured + ! up to this time t, as determined by interpolation from the + ! (possibly time-dependent) ! dtopo files. + + ! There is also a commented-out section below indicating how you might set + ! a higher value in some region that contains an onshore lake, for example. + + + use topo_module + use geoclaw_module, only: sea_level + + implicit none + + ! Arguments + integer, intent(in) :: mbc,mx,my + real(kind=8), intent(in) :: xlow,ylow,dx,dy,t + real(kind=8), intent(inout) :: veta(1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Local + integer :: i,j,i1,i2,j1,j2,idtopo,jdtopo,kdtopo,m + integer(kind=8) :: index0_dtopowork, ij + real(kind=8) :: x,y,r2,r,r3,eta2,eta3,eta_slope + real(kind=8) :: x1_lake,x2_lake,y1_lake,y2_lake, lake_level + + veta = sea_level ! initialize to sea_level, update below at some (i,j) + lake_level = 30.d0 + + !r2 = 2100.d0**2 + r2 = 1300.d0 + eta2 = 380.35d0 + r3 = 1600.d0 + eta3 = 169.047d0 + eta_slope = (eta3 - eta2) / (r3 - r2) + + do j=1,my + y = ylow + (j-0.5d0)*dy + do i=1,mx + x = xlow + (i-0.5d0)*dx + r = sqrt(x**2 + y**2) + !veta(i,j) = 200.d0 * min(1.d0, (1.d0 - (r - 2.d3)/10.d0)) + if (r > 1100.d0 .and. r < r3) then + if (r < r2) then + veta(i,j) = eta2 + else + veta(i,j) = eta2 + (r-r2)*eta_slope + endif + endif + + !veta(i,j) = 200.d0 * min(1.d0, (1.d0 - (r - 2.d3)/10.d0)) + !if (x**2 + y**2 < r2) then + ! veta(i,j) = lake_level + !endif + enddo + enddo + + if (.false.) then + ! This illustrates how you might set a higher value than sea_level + ! in a rectangle containing an onshore lake, as an example. + ! Depending on the geometry, a simple rectangle might not suffice. + ! What's below has been used for Lake Ozette on the Washington coast. + lake_level = 11.d0 ! meters relative to 0 datum of topofiles + x1_lake = -124.67d0 + x2_lake = -124.58d0 + y1_lake = 48.03d0 + y2_lake = 48.16d0 + i1 = max(int(floor((x1_lake - xlow)/dx)), 1) + i2 = min(int(floor((x2_lake - xlow)/dx)), mx) + j1 = max(int(floor((y1_lake - ylow)/dy)), 1) + j2 = min(int(floor((y2_lake - ylow)/dy)), my) + + forall(i=i1:i2, j=j1:j2) + veta(i,j) = lake_level + end forall + + endif + + ! The code below adjusts veta values set above for any uplift or + ! subsidence found in this region at the current time t, which assumes + ! that the water surface moves with the grount initially. + + if (num_dtopo == 0) return + + do m=1,num_dtopo + if (t < t0dtopo(m)) cycle ! this dtopo isn't moving yet + + ! compute indices in patch that overlap with this dtopo, + ! cycling to the next dtopo array if we discover there's no overlap: + + i1 = max(int(floor((xlowdtopo(m) - xlow)/dx)), 1) + if (i1 >= mx) cycle + + i2 = min(int(floor((xhidtopo(m) - xlow)/dx)), mx) + if (i2 < 1) cycle + + j1 = max(int(floor((ylowdtopo(m) - ylow)/dy)), 1) + if (j1 >= my) cycle + + j2 = min(int(floor((yhidtopo(m) - ylow)/dy)), my) + if (j2 < 1) cycle + + ! There is some overlap of dtopo with this patch + ! Next figure out index into time-dependent dtopo based on t: + + if (mtdtopo(m) == 1) then + ! Special case: instantaneous displacement at one instant in time + kdtopo = 1 + else + kdtopo = int(floor((t-t0dtopo(m))/dtdtopo(m)))+1 + kdtopo = min(kdtopo,mtdtopo(m)) + kdtopo = max(kdtopo,1) + endif + + index0_dtopowork = i0dtopo(m) + int(kdtopo-1, 8)*int(mxdtopo(m), 8)*int(mydtopo(m), 8) + !write(6,*) '+++ index0_dtopowork = ',index0_dtopowork + + ! Adjust eta_init by dtopo on part of patch that overlaps dtopo. + ! Code below assumes dtopo is smooth enough that we can just + ! evaluate at one point in space and time, not doing interpolation. + ! This could be improved as in topo_update.f, but probably not + ! necessary (??). + + do i=i1,i2 + x = xlow + (i-0.5d0)*dx + do j=j1,j2 + y = ylow + (j-0.5d0)*dy + idtopo = int(floor((x-xlowdtopo(m))/dxdtopo(m))) + 1 + idtopo = max(1, min(mxdtopo(m)-1, idtopo)) + jdtopo = int(floor((yhidtopo(m)-y)/dydtopo(m))) + 1 + jdtopo = max(1, min(mydtopo(m)-1, jdtopo)) + ij = index0_dtopowork + int(jdtopo-1, 8)*int(mxdtopo(m) + idtopo-1, 8) + + veta(i,j) = veta(i,j) + dtopowork(ij) + + enddo + enddo + + enddo + return + +end subroutine set_eta_init diff --git a/tests/geoclaw_radial_slide/setplot.py b/tests/geoclaw_radial_slide/setplot.py new file mode 100644 index 0000000..4635623 --- /dev/null +++ b/tests/geoclaw_radial_slide/setplot.py @@ -0,0 +1,473 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +import numpy as np +import matplotlib.pyplot as plt +import cmocean +import matplotlib as mpl + +from clawpack.visclaw import geoplot +from clawpack.visclaw import colormaps +from clawpack.visclaw import gridtools + + +import os,sys + +sea_level = 200. + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + import clawpack.dclaw.plot as dplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + #plotdata.format = 'binary' + + + + def timeformat(t): + from numpy import mod + hours = int(t/3600.) + tmin = mod(t,3600.) + min = int(tmin/60.) + sec = int(mod(tmin,60.)) + timestr = '%s:%s:%s' % (hours,str(min).zfill(2),str(sec).zfill(2)) + return timestr + + def title_hours(current_data): + from pylab import title + t = current_data.t + timestr = timeformat(t) + title('t = %s' % timestr) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Computational domain', figno=0) + plotfigure.kwargs = {'figsize':(8,7)} + plotfigure.show = False + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + + def aa(current_data): + from pylab import ticklabel_format, xticks, gca, cos, pi, savefig + gca().set_aspect(1.) + title_hours(current_data) + ticklabel_format(useOffset=False) + xticks(rotation=20) + + + plotaxes.afteraxes = aa + + + # Hillshade + plotitem = plotaxes.new_plotitem(plot_type="2d_hillshade") + plotitem.show = False + plotitem.plot_var = dplot.eta + plotitem.add_colorbar = False + + # Surface + plotitem = plotaxes.new_plotitem(plot_type="2d_imshow") + plotitem.plot_var = dplot.surface_solid_frac_lt03 + plotitem.add_colorbar = True + plotitem.colorbar_kwargs = { + "shrink": 0.9, + "location": "bottom", + "orientation": "horizontal", + } + plotitem.colorbar_label = "Surface (m)" + plotitem.imshow_cmap = cmocean.cm.curl + plotitem.imshow_cmin = 90 + plotitem.imshow_cmax = 110 + # Debris + plotitem = plotaxes.new_plotitem(plot_type="2d_imshow") + plotitem.plot_var = dplot.solid_frac_gt03 + plotitem.imshow_cmap = cmocean.cm.turbid + plotitem.imshow_cmin = 0.3 + plotitem.imshow_cmax = 1 + plotitem.add_colorbar = False + + #----------------------------------------- + # Figure for water / landslide separately + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Landslide/water surface', figno=1) + plotfigure.figsize=(8,8) + plotfigure.facecolor='w' + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.15,.5,.7,.45])' + plotaxes.title = 'Landslide / Water' + plotaxes.scaled = True + plotaxes.xlimits = [-3e3,3e3] + plotaxes.ylimits = [-3e3,3e3] + + def pure_water(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + eta = dplot.eta(current_data) + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + water = np.where(np.logical_and(h>1e-3, m<0.1), + eta, np.nan) + #import pdb; pdb.set_trace() + return water + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + #plotitem.show = False + #plotitem.plot_var = pure_water + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = sea_level - 20. + plotitem.pcolor_cmax = sea_level + 20. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def landslide(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + eta = dplot.eta(current_data) + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + landslide = np.where(np.logical_and(h>1e-3, m>0.1), + h, np.nan) + #import pdb; pdb.set_trace() + return landslide + + # Landslide + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + plotitem.show = False + plotitem.plot_var = landslide + cmap_mass = colormaps.make_colormap({0.:'w', 1.:'brown'}) + plotitem.pcolor_cmap = cmap_mass + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 80. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def land(current_data): + """ + Return a masked array containing the surface elevation only in dry cells. + """ + drytol = 1e-3 + q = current_data.q + h = q[0,...] + eta = q[-1,...] + land = np.ma.masked_where(h>drytol, eta) + #import pdb; pdb.set_trace() + #land = eta - h # everywhere + return land + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figure for cross section compared to 1d_radial + #----------------------------------------- + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('radial slice') + plotaxes.axescmd = 'axes([.1,.1,.8,.3])' + plotaxes.title = 'Diagonal Transect' + #plotaxes.scaled = True + + def plot_xsec(current_data): + from pylab import plot,linspace,zeros,ones,legend,xlabel,\ + sqrt,grid,xlim,fill_between,logical_and + from pylab import nan,where,ylim,loadtxt,arange + from clawpack.pyclaw import Solution + pd = current_data.plotdata + frameno = current_data.frameno + framesoln = Solution(frameno, path=pd.outdir, file_format=pd.format) + + #xout = linspace(-3e3,3e3,1000) + #yout = ones(xout.shape) # near x-axis + #rout = xout + + xout = linspace(-3e3,3e3,1000) + yout = xout + rout = xout * sqrt(2) + + etaout = gridtools.grid_output_2d(framesoln, -1, xout, yout) + hout = gridtools.grid_output_2d(framesoln, 0, xout, yout) + zetaout = where(hout>0.001, etaout, nan) + Bout = etaout - hout + if 0: + hmout = gridtools.grid_output_2d(framesoln, 3, xout, yout) + with np.errstate(divide="ignore", invalid="ignore"): + mout = hmout / hout + water = where(mout<0.1, etaout, nan) + landslide = where(mout>0.5, etaout, nan) + mixed = where(logical_and(0.11e-3, m<0.1), + h, np.nan) + import pdb; pdb.set_trace() + return water + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + #plotitem.plot_var = geoplot.surface_or_depth + #plotitem.show = False + plotitem.plot_var = 0 + plotitem.pcolor_cmap = colormaps.white_red + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 100. + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + def land(current_data): + """ + Return a masked array containing the surface elevation only in dry cells. + """ + drytol = 1e-3 + q = current_data.q + h = q[0,...] + eta = q[-1,...] + land = np.ma.masked_where(h>drytol, eta) + #import pdb; pdb.set_trace() + #land = eta - h # everywhere + return land + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figure for mass fraction + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Mass Fraction', figno=6) + plotfigure.kwargs = {'figsize':(8,7)} + plotfigure.show = False + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'mass fraction' + plotaxes.scaled = True + + plotaxes.scaled = True + plotaxes.xlimits = [-3e3,3e3] + plotaxes.ylimits = [-3e3,3e3] + + def mass_frac(current_data): + q = current_data.q + h = q[0,:,:] + hm = q[3,:,:] + with np.errstate(divide="ignore", invalid="ignore"): + m = hm / h + mwet = np.where(h > 0.01, m, np.nan) + return mwet + + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = mass_frac + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = 0. + plotitem.pcolor_cmax = 0.65 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.show = False + plotitem.plot_var = land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 400.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 0 + + + + # Figure for scatter plot + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter', figno=3) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0, 3e3*np.sqrt(2)] + plotaxes.ylimits = 'auto' + plotaxes.title = 'Scatter plot' + plotaxes.grid = True + + # Set up for item on these axes: scatter of 2d data + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + + def h_vs_r(current_data): + # Return radius of each grid cell and p value in the cell + from pylab import sqrt + x = current_data.x + y = current_data.y + r = np.sqrt(x**2 + y**2) + q = current_data.q + h = q[0,:,:] + return r,h + + plotitem.map_2d_to_1d = h_vs_r + plotitem.plot_var = 0 + plotitem.plotstyle = '.' + plotitem.color = 'b' + plotitem.kwargs = {'markersize':1} + plotitem.show = True # show on plot? + + + # ----------------------- + # Figure for scatter plot of speed + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter_speed', figno=9) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0, 3e3*np.sqrt(2)] + plotaxes.ylimits = 'auto' + plotaxes.title = 'speed' + plotaxes.grid = True + + + def s_vs_r(current_data): + # Return radius of each grid cell and p value in the cell + from pylab import sqrt,nan,where + x = current_data.x + y = current_data.y + r = np.sqrt(x**2 + y**2) + q = current_data.q + s = sqrt(q[1,:,:]**2 + q[2,:,:]**2) + s = sqrt(q[1,:,:]**2 + q[2,:,:]**2) + h = where(q[0,:,:]>1e-2, q[0,:,:], nan) + s = s/h + return r,s + + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + plotitem.map_2d_to_1d = s_vs_r + #plotitem.plot_var = 0 + plotitem.plotstyle = '+' + plotitem.color = 'r' + #plotitem.kwargs = {'markersize':1} + plotitem.show = True # show on plot? + + #------------------------------------- + + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + import os + from clawpack.visclaw import plot_timing_stats + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + units = {'comptime':'minutes', 'simtime':'minutes', 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, make_pngs=True, + plotdir=timing_plotdir, units=units) + os.system('cp %s/timing.* %s' % (plotdata.outdir, timing_plotdir)) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata diff --git a/tests/geoclaw_radial_slide/setrun.py b/tests/geoclaw_radial_slide/setrun.py new file mode 100644 index 0000000..b1d7f5c --- /dev/null +++ b/tests/geoclaw_radial_slide/setrun.py @@ -0,0 +1,487 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os, sys +import numpy as np + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW environment variable") + +from clawpack.amrclaw.data import FlagRegion + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + #probdata.add_param('variable_eta_init', True) # now in qinit info + + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -3e3 + clawdata.upper[0] = 3e3 + + clawdata.lower[1] = -3e3 + clawdata.upper[1] = 3e3 + + + # Number of grid cells: Coarsest grid + clawdata.num_cells[0] = 240 + clawdata.num_cells[1] = 240 + + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 1 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 40 #240 + clawdata.tfinal = 160. #240. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [0.5, 1.0] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 10 + clawdata.output_t0 = True + + + clawdata.output_format = 'ascii' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = True # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.0001 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + # D-Claw requires CFL<0.5 + clawdata.cfl_desired = 0.85 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = [4, 4, 4] # TODO VERIFY THAT 4 in old and new are the same + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + # TODO This is not in old setrun.py + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 1 + + # List of refinement ratios at each level (length at least mxnest-1) + # dx = dy = 2', 10", 2", 1/3": + amrdata.refinement_ratios_x = [2,2] + amrdata.refinement_ratios_y = [2,2] + amrdata.refinement_ratios_t = [2,2] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = [ + "center"] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + # NO OLD STYLE REGIONS USED HERE + + + # --------------- + # NEW flagregions + # --------------- + + + flagregions = rundata.flagregiondata.flagregions # initialized to [] + + # now append as many flagregions as desired to this list: + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges = [] + + # Set GeoClaw specific runtime parameters. + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 1 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 200.0 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True # TODO change? + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.01 + + # == settopo.data values == + topofiles = rundata.topo_data.topofiles + # for topography, append lines of the form + # [topotype, fname] + topofiles.append([3, 'basal_topo.tt3']) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + + rundata.qinit_data.variable_eta_init = True + + if 0: + # == setqinit.data values == + qinitdclaw_data = rundata.qinitdclaw_data # initialized when rundata instantiated + + etafile = 'surface_topo.tt3' + qinitdclaw_data.qinitfiles.append([3, 8, 1, 2, etafile]) + + mfile = 'mass_frac.tt3' + qinitdclaw_data.qinitfiles.append([3, 4, 1, 2, mfile]) + + #hfile = 'landslide_depth.tt3' + #qinitdclaw_data.qinitfiles.append([3, 1, 1, 2, hfile]) + + # == setauxinit.data values == + #auxinitdclaw_data = rundata.auxinitdclaw_data # initialized when rundata instantiated + + # == fgmax.data values == + #fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + + # == setdclaw.data values == + dclaw_data = rundata.dclaw_data # initialized when rundata instantiated + + dclaw_data.c1 = 1.0 # do we want to remove this? + dclaw_data.rho_f = 1000.0 + dclaw_data.rho_s = 2700.0 + dclaw_data.phi_bed = 32.0 + dclaw_data.theta_input = 0.0 + dclaw_data.mu = 0.005 + dclaw_data.m0 = 0.63 + dclaw_data.m_crit = 0.64 + dclaw_data.kappita = 1.e-8 + #dclaw_data.kappita_diff = 1 + #dclaw_data.chi_init_val=0.5 # not currently used. + dclaw_data.alpha_c = 0.05 + dclaw_data.alpha_seg = 0.0 + #dclaw_data.phi_seg_coeff = 0.0 + dclaw_data.delta = 0.001 + dclaw_data.bed_normal = 0 + dclaw_data.entrainment = 0 + dclaw_data.entrainment_rate = 0.0 + dclaw_data.sigma_0 = 1.0e3 + #dclaw_data.mom_autostop = True + #dclaw_data.momlevel = 1 + #dclaw_data.mom_perc = 0.0 + + # == pinitdclaw.data values == + pinitdclaw_data = rundata.pinitdclaw_data # initialized when rundata instantiated + + pinitdclaw_data.init_ptype = 0 # hydrostatic (-1 ==> zero everywhere) + pinitdclaw_data.init_pmax_ratio = 0.00e0 + pinitdclaw_data.init_ptf = 0.0 + pinitdclaw_data.init_ptf2 = 0.0 + + # == flowgrades.data values == + flowgrades_data = rundata.flowgrades_data # initialized when rundata instantiated + + flowgrades_data.flowgrades = [] + # for using flowgrades for refinement append lines of the form + # [flowgradevalue, flowgradevariable, flowgradetype, flowgrademinlevel] + # where: + # flowgradevalue: floating point relevant flowgrade value for following measure: + # flowgradevariable: 1=depth, 2= momentum, 3 = sign(depth)*(depth+topo) (0 at sealevel or dry land). + # flowgradetype: 1 = norm(flowgradevariable), 2 = norm(grad(flowgradevariable)) + # flowgrademinlevel: refine to at least this level if flowgradevalue is exceeded. + + + #flowgrades_data.keep_fine = True + #flowgrades_data.flowgrades.append([1.0e-6, 2, 1, 1]) + #flowgrades_data.flowgrades.append([1.0e-6, 1, 1, 1]) + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + amrdata.max1d = 300 + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + return rundata + + + + + # end of function setrun + # ---------------------- + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() +