From ff1ce7cfc67cea59c735c97d9ac98c44b90274b4 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Wed, 5 Jun 2024 14:17:29 -0800 Subject: [PATCH 1/3] add fgout output_style parameter and support for array of output_times --- src/2d/shallow/fgout_module.f90 | 105 +++++++++++++++++------------- src/python/geoclaw/fgout_tools.py | 30 +++++++-- 2 files changed, 83 insertions(+), 52 deletions(-) diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90 index bd9effda3..8f422730e 100644 --- a/src/2d/shallow/fgout_module.f90 +++ b/src/2d/shallow/fgout_module.f90 @@ -10,7 +10,7 @@ module fgout_module real(kind=8), pointer :: late(:,:,:) ! Geometry - integer :: num_vars(2),mx,my,point_style,fgno,output_format + integer :: num_vars(2),mx,my,point_style,fgno,output_format,output_style real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi ! Time Tracking and output types @@ -84,18 +84,27 @@ subroutine set_fgout(rest,fname) 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%output_style read(unit,*) fg%num_output + allocate(fg%output_times(fg%num_output)) + allocate(fg%output_frames(fg%num_output)) + + if (fg%output_style == 1) then + read(unit,*) fg%start_time + read(unit,*) fg%end_time + else if (fg%output_style == 2) then + read(unit,*) (fg%output_times(k), k=1,fg%num_output) + fg%start_time = fg%output_times(1) + fg%end_time = fg%output_times(fg%num_output) + endif + 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 @@ -104,49 +113,55 @@ subroutine set_fgout(rest,fname) 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 + + if (fg%output_style == 1) then + ! Setup data for this grid + ! Set fg%dt (the timestep length between outputs) + 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 + ! only a single fgout time: + 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 + enddo + endif + endif + endif + + do k=1,fg%num_output + if (rest) then + ! don't write initial time or earlier + ts = tstart_thisrun+FGOUT_ttol else - fg%dt = 0.d0 + ! do write initial time + ts = tstart_thisrun-FGOUT_ttol 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 + + 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 - 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 + ! will be reset to frameno when this is written + fg%output_frames(k) = -1 endif - endif + enddo diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py index 11494b124..b140c98e2 100644 --- a/src/python/geoclaw/fgout_tools.py +++ b/src/python/geoclaw/fgout_tools.py @@ -178,6 +178,7 @@ def __init__(self,fgno=None,outdir=None,output_format=None): self.npts = None self.nx = None self.ny = None + self.output_style = 1 self.tstart = None self.tend = None self.nout = None @@ -372,18 +373,33 @@ def write_to_fgout_data(self, fid): errmsg = "fgout output_format must be ascii, binary32, or binary64" raise NotImplementedError(errmsg) - assert self.tstart is not None, 'Need to set tstart' - assert self.tend is not None, 'Need to set tend' - assert self.nout is not None, 'Need to set nout' - assert self.point_style is not None, 'Need to set point_style' + # write header, independent of point_style: #fid = open(self.input_file_name,'w') fid.write("\n") fid.write("%i # fgno\n" % self.fgno) - fid.write("%16.10e # tstart\n" % self.tstart) - fid.write("%16.10e # tend\n" % self.tend) - fid.write("%i %s # nout\n" % (self.nout, 11*" ")) + + fid.write("%i # output_style\n" \ + % self.output_style) + + if self.output_style == 1: + assert self.tstart is not None, 'Need to set tstart' + assert self.tend is not None, 'Need to set tend' + assert self.nout is not None, 'Need to set nout' + fid.write("%i %s # nout\n" % (self.nout, 11*" ")) + fid.write("%16.10e # tstart\n" % self.tstart) + fid.write("%16.10e # tend\n" % self.tend) + elif self.output_style == 2: + self.nout = len(self.output_times) + fid.write("%i %s # nout\n" % (self.nout, 11*" ")) + + # remove [] and , from list of times: + output_times_str = repr(list(self.output_times))[1:-1] + output_times_str = output_times_str.replace(',','') + fid.write("%s # output_times\n" % output_times_str) + else: + raise ValueError('fgout output_style must be 1 or 2') fid.write("%i %s # point_style\n" \ % (self.point_style,12*" ")) fid.write("%i %s # output_format\n" \ From 564cbbc69f0214bc9e39e6f1e56fa658d5a01c19 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Wed, 5 Jun 2024 14:31:21 -0800 Subject: [PATCH 2/3] remove trailing whitespaces in fgout_module.f990 --- src/2d/shallow/fgout_module.f90 | 258 ++++++++++++++++---------------- 1 file changed, 129 insertions(+), 129 deletions(-) diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90 index 8f422730e..69cce221b 100644 --- a/src/2d/shallow/fgout_module.f90 +++ b/src/2d/shallow/fgout_module.f90 @@ -8,18 +8,18 @@ module fgout_module ! Grid data real(kind=8), pointer :: early(:,:,:) real(kind=8), pointer :: late(:,:,:) - + ! Geometry integer :: num_vars(2),mx,my,point_style,fgno,output_format,output_style 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 + end type fgout_grid logical, private :: module_setup = .false. @@ -32,27 +32,27 @@ module fgout_module 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 @@ -75,7 +75,7 @@ subroutine set_fgout(rest,fname) write(parmunit,*) ' No fixed grids specified for output' return endif - + ! Allocate fixed grids (not the data yet though) allocate(FGOUT_fgrids(FGOUT_num_grids)) @@ -97,28 +97,28 @@ subroutine set_fgout(rest,fname) fg%start_time = fg%output_times(1) fg%end_time = fg%output_times(fg%num_output) endif - + 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 - - + + ! 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 - - if (fg%output_style == 1) then + + if (fg%output_style == 1) then ! Setup data for this grid ! Set fg%dt (the timestep length between outputs) if (fg%end_time <= fg%start_time) then - if (fg%num_output > 1) 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' @@ -142,7 +142,7 @@ subroutine set_fgout(rest,fname) endif endif endif - + do k=1,fg%num_output if (rest) then ! don't write initial time or earlier @@ -151,7 +151,7 @@ subroutine set_fgout(rest,fname) ! 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) @@ -184,30 +184,30 @@ subroutine set_fgout(rest,fname) else print *,'set_fgout: ERROR for fixed grid', i print *,'y grid points my <= 0, set my >= 1' - endif - + 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) = 6 - + ! 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 @@ -216,10 +216,10 @@ end subroutine set_fgout subroutine fgout_interp(fgrid_type,fgrid, & t,q,meqn,mxc,myc,mbc,dxc,dyc,xlowc,ylowc, & maux,aux) - - use geoclaw_module, only: dry_tolerance + + use geoclaw_module, only: dry_tolerance implicit none - + ! Subroutine arguments integer, intent(in) :: fgrid_type type(fgout_grid), intent(inout) :: fgrid @@ -227,9 +227,9 @@ subroutine fgout_interp(fgrid_type,fgrid, & 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 @@ -240,18 +240,18 @@ subroutine fgout_interp(fgrid_type,fgrid, & ! 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) @@ -264,103 +264,103 @@ subroutine fgout_interp(fgrid_type,fgrid, & stop ! fgrid_type==3 is deprecated, use fgmax grids instead endif - - xhic = xlowc + dxc*mxc - yhic = ylowc + dyc*myc - + + xhic = xlowc + dxc*mxc + yhic = ylowc + dyc*myc + ! Find indices of various quantities in the fgrid arrays - bathy_index = 4 ! works for both shallow and bouss + bathy_index = 4 ! works for both shallow and bouss eta_index = 5 ! works for both shallow and bouss - + !write(59,*) '+++ ifg,jfg,eta,geometry at t = ',t - - ! Primary interpolation loops + + ! 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 + + ! 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 + 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) + interpolate(q(m,ic1:ic2,jc1:jc2), geometry) end forall - fg_data(bathy_index,ifg,jfg) = & + 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 + ! 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 @@ -372,12 +372,12 @@ 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 @@ -385,14 +385,14 @@ subroutine fgout_write(fgrid,out_time,out_index) 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 + + ! 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',/," // & @@ -402,7 +402,7 @@ subroutine fgout_write(fgrid,out_time,out_index) "e26.16,' ylow', /," // & "e26.16,' dx', /," // & "e26.16,' dy',/)" - + character(len=*), parameter :: t_file_format = "(e18.8,' time', /," // & "i6,' meqn'/," // & "i6,' ngrids'/," // & @@ -410,64 +410,64 @@ subroutine fgout_write(fgrid,out_time,out_index) "i6,' ndim'/," // & "i6,' nghost'/," // & "a10,' format'/,/)" - + ! Other locals integer :: i,j,m real(kind=8) :: t0,tf,tau, qaug(6) real(kind=8), allocatable :: qeta(:,:,:) real(kind=4), allocatable :: qeta4(:,:,:) real(kind=8) :: h_early,h_late,topo_early,topo_late - + allocate(qeta(4, fgrid%mx, fgrid%my)) ! to store h,hu,hv,eta - - - ! Interpolate the grid in time, to the output time, using - ! the solution in fgrid1 and fgrid2, which represent the + + + ! 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 + + ! 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(4,i,j) topo_late = fgrid%late(4,i,j) @@ -485,13 +485,13 @@ subroutine fgout_write(fgrid,out_time,out_index) endif endif endif - + ! Output the conserved quantities and eta value 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(5) ! eta - + enddo enddo @@ -514,8 +514,8 @@ subroutine fgout_write(fgrid,out_time,out_index) cframeno(ipos:ipos) = char(ichar('0') + idigit) iframe1 = iframe1/10 enddo - - fg_filename = 'fgout' // cfgno // '.q' // cframeno + + fg_filename = 'fgout' // cfgno // '.q' // cframeno open(unit,file=fg_filename,status='unknown',form='formatted') @@ -524,17 +524,17 @@ subroutine fgout_write(fgrid,out_time,out_index) 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 @@ -544,20 +544,20 @@ subroutine fgout_write(fgrid,out_time,out_index) enddo write(unit,*) ' ' ! blank line required between rows enddo - endif - + 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 + 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 + fg_filename = 'fgout' // cfgno // '.b' // cframeno open(unit=unit, file=fg_filename, status="unknown", & access='stream') allocate(qeta4(4, fgrid%mx, fgrid%my)) ! for 4-byte binary output @@ -566,44 +566,44 @@ subroutine fgout_write(fgrid,out_time,out_index) 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 + 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 + + 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, 4, 1, 0, 2, 0, file_format close(unit) - + 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) @@ -612,21 +612,21 @@ real(kind=8) function nan() 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 + ! 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 @@ -634,8 +634,8 @@ real(kind=8) pure function interpolate(points,geometry) result(interpolant) 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 From bcc4288b9a740cbe8e5c7d92a4eac447cb2f0fcd Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Wed, 5 Jun 2024 18:41:00 -0800 Subject: [PATCH 3/3] Refactor fgout_module.f90 so it works for either GeoClaw or D-Claw and support in fgout_tools.py for new dclaw attribute dclaw to set in setrun.py to indicate D-Claw, in which case 7 components of q are output instead of 4. Support for eventually indicating fewer components to output. --- src/2d/shallow/fgout_module.f90 | 181 ++++++++++++++++++++++-------- src/python/geoclaw/fgout_tools.py | 2 + 2 files changed, 136 insertions(+), 47 deletions(-) diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90 index 69cce221b..981d11ac3 100644 --- a/src/2d/shallow/fgout_module.f90 +++ b/src/2d/shallow/fgout_module.f90 @@ -10,7 +10,7 @@ module fgout_module real(kind=8), pointer :: late(:,:,:) ! Geometry - integer :: num_vars(2),mx,my,point_style,fgno,output_format,output_style + integer :: num_vars,mx,my,point_style,fgno,output_format,output_style real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi ! Time Tracking and output types @@ -19,6 +19,12 @@ module fgout_module integer, allocatable :: output_frames(:) real(kind=8), allocatable :: output_times(:) + + integer :: nqout ! number of q components to output (+1 for eta) + logical, allocatable :: iqout(:) ! which components to output + integer :: bathy_index,eta_index + logical :: dclaw ! False for GeoClaw + end type fgout_grid @@ -103,6 +109,7 @@ subroutine set_fgout(rest,fname) read(unit,*) fg%mx, fg%my read(unit,*) fg%x_low, fg%y_low read(unit,*) fg%x_hi, fg%y_hi + read(unit,*) fg%dclaw ! Initialize next_output_index @@ -186,15 +193,72 @@ subroutine set_fgout(rest,fname) 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) = 6 + + ! For now, hard-wire with defaults for either GeoClaw or D-Claw + ! need to save q plus topo, eta, t for interp in space-time + + if (fg%dclaw) then + ! For D-Claw: + fg%num_vars = 10 + ! for h, hu, hv, hm, pb, hchi, b_eroded, bathy, eta, time + else + ! GeoClaw: + fg%num_vars = 6 + ! for h, hu, hv, bathy, eta, time + endif + + ! specify which components of q (plus eta?) to output: + ! (eventually this should be set from user input) + + if (fg%num_vars == 6) then + ! GeoClaw + ! indexes used in early and late arrays: + ! 1:3 are q variables, 6 is time + fg%bathy_index = 4 + fg%eta_index = 5 + + allocate(fg%iqout(4)) + fg%iqout(1) = .true. ! output h? + fg%iqout(2) = .true. ! output hu? + fg%iqout(3) = .true. ! output hv? + fg%iqout(4) = .true. ! output eta? + fg%nqout = 0 + do k=1,4 + if (fg%iqout(k)) fg%nqout = fg%nqout + 1 + enddo + endif + + if (fg%num_vars == 10) then + ! D-Claw: + ! indexes used in early and late arrays: + ! 1:7 are q variables, 10 is time + fg%bathy_index = 8 + fg%eta_index = 9 + + allocate(fg%iqout(8)) + fg%iqout(1) = .true. ! output h? + fg%iqout(2) = .true. ! output hu? + fg%iqout(3) = .true. ! output hv? + fg%iqout(4) = .true. ! output hm? + fg%iqout(5) = .true. ! output pb? + fg%iqout(6) = .true. ! output hchi? + fg%iqout(7) = .true. ! output beroded? + fg%iqout(8) = .true. ! output eta? + fg%nqout = 0 + do k=1,8 + if (fg%iqout(k)) fg%nqout = fg%nqout + 1 + enddo + endif + + write(6,*) '+++ nqout = ',fg%nqout - ! Allocate new fixed grid data array - allocate(fg%early(fg%num_vars(1), fg%mx,fg%my)) + ! Allocate new fixed grid data arrays at early, late time: + ! dimension (10,:,:) to work for either GeoClaw or D-Claw + + allocate(fg%early(10, fg%mx,fg%my)) fg%early = nan() - allocate(fg%late(fg%num_vars(1), fg%mx,fg%my)) + allocate(fg%late(10, fg%mx,fg%my)) fg%late = nan() enddo @@ -232,7 +296,6 @@ subroutine fgout_interp(fgrid_type,fgrid, & ! Indices integer :: ifg,jfg,m,ic1,ic2,jc1,jc2 - integer :: bathy_index,eta_index ! Tolerances real(kind=8) :: total_depth,depth_indicator,nan_check @@ -254,10 +317,10 @@ subroutine fgout_interp(fgrid_type,fgrid, & ! Setup aliases for specific fixed grid if (fgrid_type == 1) then - num_vars = fgrid%num_vars(1) + num_vars = fgrid%num_vars fg_data => fgrid%early else if (fgrid_type == 2) then - num_vars = fgrid%num_vars(1) + num_vars = fgrid%num_vars fg_data => fgrid%late else write(6,*) '*** Unexpected fgrid_type = ', fgrid_type @@ -268,9 +331,6 @@ subroutine fgout_interp(fgrid_type,fgrid, & xhic = xlowc + dxc*mxc yhic = ylowc + dyc*myc - ! Find indices of various quantities in the fgrid arrays - bathy_index = 4 ! works for both shallow and bouss - eta_index = 5 ! works for both shallow and bouss !write(59,*) '+++ ifg,jfg,eta,geometry at t = ',t @@ -299,12 +359,12 @@ subroutine fgout_interp(fgrid_type,fgrid, & fg_data(m,ifg,jfg) = q(m,ic1,jc1) end forall - fg_data(bathy_index,ifg,jfg) = aux(1,ic1,jc1) + fg_data(fgrid%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) + fg_data(fgrid%eta_index,ifg,jfg) = fg_data(1,ifg,jfg) & + + fg_data(fgrid%bathy_index,ifg,jfg) else if (method == 1) then @@ -335,7 +395,7 @@ subroutine fgout_interp(fgrid_type,fgrid, & interpolate(q(m,ic1:ic2,jc1:jc2), geometry) end forall - fg_data(bathy_index,ifg,jfg) = & + fg_data(fgrid%bathy_index,ifg,jfg) = & interpolate(aux(1,ic1:ic2,jc1:jc2),geometry) @@ -346,7 +406,7 @@ subroutine fgout_interp(fgrid_type,fgrid, & ! 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) + fg_data(fgrid%eta_index,ifg,jfg) = interpolate(eta,geometry) ! NEED TO FIX endif @@ -412,14 +472,14 @@ subroutine fgout_write(fgrid,out_time,out_index) "a10,' format'/,/)" ! Other locals - integer :: i,j,m - real(kind=8) :: t0,tf,tau, qaug(6) + integer :: i,j,m,iq,k + real(kind=8) :: t0,tf,tau, qaug(10) real(kind=8), allocatable :: qeta(:,:,:) real(kind=4), allocatable :: qeta4(:,:,:) real(kind=8) :: h_early,h_late,topo_early,topo_late - allocate(qeta(4, fgrid%mx, fgrid%my)) ! to store h,hu,hv,eta - + allocate(qeta(fgrid%nqout, fgrid%mx, fgrid%my)) ! to store h,hu,hv,eta + ! or subset ! Interpolate the grid in time, to the output time, using ! the solution in fgrid1 and fgrid2, which represent the @@ -428,7 +488,7 @@ subroutine fgout_write(fgrid,out_time,out_index) 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) + forall(m=1:fgrid%num_vars-1,abs(fgrid%early(m,i,j)) < 1d-90) fgrid%early(m,i,j) = 0.d0 end forall @@ -449,12 +509,12 @@ subroutine fgout_write(fgrid,out_time,out_index) ! 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) + t0 = fgrid%early(fgrid%num_vars,i,j) + tf = fgrid%late(fgrid%num_vars,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) + forall(m=1:fgrid%num_vars-1,abs(fgrid%late(m,i,j)) < 1d-90) fgrid%late(m,i,j) = 0.d0 end forall @@ -487,11 +547,51 @@ subroutine fgout_write(fgrid,out_time,out_index) endif ! Output the conserved quantities and eta value - 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(5) ! eta - + iq = 1 + ! qaug(1:3) are h,hu,hv for both GeoClaw and D-Claw: + if (fgrid%iqout(1)) then + qeta(iq,i,j) = qaug(1) ! h + iq = iq+1 + endif + if (fgrid%iqout(2)) then + qeta(iq,i,j) = qaug(2) ! hu + iq = iq+1 + endif + if (fgrid%iqout(3)) then + qeta(iq,i,j) = qaug(3) ! hv + iq = iq+1 + endif + + if (fgrid%num_vars == 6) then + ! GeoClaw: + if (fgrid%iqout(4)) then + qeta(iq,i,j) = qaug(5) ! eta since qaug(4)=topo + iq = iq+1 + endif + + else if (fgrid%num_vars == 10) then + ! D-Claw: + if (fgrid%iqout(4)) then + qeta(iq,i,j) = qaug(4) ! hm + iq = iq+1 + endif + if (fgrid%iqout(5)) then + qeta(iq,i,j) = qaug(5) ! pb + iq = iq+1 + endif + if (fgrid%iqout(6)) then + qeta(iq,i,j) = qaug(6) ! hchi + iq = iq+1 + endif + if (fgrid%iqout(7)) then + qeta(iq,i,j) = qaug(7) ! b_eroded + iq = iq+1 + endif + if (fgrid%iqout(8)) then + qeta(iq,i,j) = qaug(9) ! eta since qaug(8)=topo + iq = iq+1 + endif + endif enddo enddo @@ -519,18 +619,6 @@ subroutine fgout_write(fgrid,out_time,out_index) 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 @@ -539,8 +627,7 @@ subroutine fgout_write(fgrid,out_time,out_index) ! 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) + write(unit, "(50e26.16)") (qeta(k,i,j), k=1,fgrid%nqout) enddo write(unit,*) ' ' ! blank line required between rows enddo @@ -560,7 +647,7 @@ subroutine fgout_write(fgrid,out_time,out_index) fg_filename = 'fgout' // cfgno // '.b' // cframeno open(unit=unit, file=fg_filename, status="unknown", & access='stream') - allocate(qeta4(4, fgrid%mx, fgrid%my)) ! for 4-byte binary output + allocate(qeta4(fgrid%nqout, fgrid%mx, fgrid%my)) qeta4 = real(qeta, kind=4) write(unit) qeta4 deallocate(qeta4) @@ -587,7 +674,7 @@ subroutine fgout_write(fgrid,out_time,out_index) 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, 4, 1, 0, 2, 0, file_format + write(unit, t_file_format) out_time, fgrid%nqout, 1, 0, 2, 0,file_format close(unit) print "(a,i4,a,i4,a,e18.8)",'Writing fgout grid #',fgrid%fgno, & diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py index b140c98e2..87ae51aa3 100644 --- a/src/python/geoclaw/fgout_tools.py +++ b/src/python/geoclaw/fgout_tools.py @@ -185,6 +185,7 @@ def __init__(self,fgno=None,outdir=None,output_format=None): self.fgno = fgno self.outdir = outdir self.output_format = output_format + self.dclaw = False self.drytol = 1e-3 # used for computing u,v from hu,hv @@ -430,6 +431,7 @@ def write_to_fgout_data(self, fid): print(" upper right = (%15.10f,%15.10f)" % (x2,y2)) print(" dx = %15.10e, dy = %15.10e" % (dx,dy)) + fid.write("%s # dclaw" % self.dclaw) def read_frame(self, frameno):