Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow specifying which components of q to output on each fgout grid #624

Merged
merged 11 commits into from
Jul 30, 2024
Merged

Allow specifying which components of q to output on each fgout grid #624

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
556c769
Allow specifying which components of q to output on each fgout grid
rjleveque Jul 14, 2024
b5c898a
working on dclaw variables
rjleveque Jul 17, 2024
70ee91e
Fix fgout_tools.py to support fgout.q_out_vars also for Boussinesq or…
rjleveque Jul 18, 2024
981ab8c
fix data_file -> data_path
rjleveque Jul 24, 2024
9378f12
add FGoutGrid.read_fgout_grids_data_pre511 for reading legacy fgout_g…
rjleveque Jul 29, 2024
665bbda
add to doc string for pre511
rjleveque Jul 29, 2024
f373e23
print more helpful error msg if fgout_grids.data seems to be pre511
rjleveque Jul 29, 2024
6c849f2
force user to call read_fgout_grids_data before read_frame
rjleveque Jul 29, 2024
178403e
Fix chiles2010_fgmax-fgout scripts to make fgout animations
rjleveque Jul 29, 2024
1ec2512
Leave output_format as an argument to FGoutGrids.__init__ but not that
rjleveque Jul 29, 2024
6e67973
fix qmap=='geoclaw-bouss' mapping
rjleveque Jul 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,6 @@
fgno = 1 # which fgout grid

outdir = '_output'
format = 'binary' # format of fgout grid output

if 1:
# use all fgout frames in outdir:
Expand All @@ -48,7 +47,8 @@
fgframes = range(1,26) # frames of fgout solution to use in animation

# Instantiate object for reading fgout frames:
fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format)
fgout_grid = fgout_tools.FGoutGrid(fgno, outdir)
fgout_grid.read_fgout_grids_data()


# Plot one frame of fgout data and define the Artists that will need to
Expand Down
4 changes: 2 additions & 2 deletions examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,6 @@
fgno = 1 # which fgout grid

outdir = '_output'
format = 'binary' # format of fgout grid output

if 1:
# use all fgout frames in outdir:
Expand All @@ -52,7 +51,8 @@
fgframes = range(1,26) # frames of fgout solution to use in animation

# Instantiate object for reading fgout frames:
fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format)
fgout_grid = fgout_tools.FGoutGrid(fgno, outdir)
fgout_grid.read_fgout_grids_data()


# Plot one frame of fgout data and define the Artists that will need to
Expand Down
4 changes: 2 additions & 2 deletions src/2d/shallow/amr2.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_qinit() ! specifies file with dh if this used instead
call set_fgout(rest) ! Fixed grid settings
call set_fgout(rest,nvar) ! Fixed grid settings
call setup_variable_friction() ! Variable friction parameter
!call set_multilayer() ! Set multilayer SWE parameters
call set_storm() ! Set storm parameters
Expand Down Expand Up @@ -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_qinit() ! specifies file with dh if this used instead
call set_fgout(rest) ! Fixed grid settings
call set_fgout(rest,nvar) ! Fixed grid settings
call setup_variable_friction() ! Variable friction parameter
call set_multilayer() ! Set multilayer SWE parameters
call set_storm() ! Set storm parameters
Expand Down
198 changes: 62 additions & 136 deletions src/2d/shallow/fgout_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ module fgout_module
real(kind=8), pointer :: late(:,:,:)

! Geometry
integer :: num_vars,mx,my,point_style,fgno,output_format,output_style
integer :: 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
Expand All @@ -19,11 +19,13 @@ module fgout_module

integer, allocatable :: output_frames(:)
real(kind=8), allocatable :: output_times(:)


integer :: num_vars ! number of variables to interpolate (num_eqn+3)
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
logical, allocatable :: q_out_vars(:) ! which components to print

!logical, allocatable :: iqout(:) ! which components to output
integer :: bathy_index,eta_index,time_index

end type fgout_grid

Expand All @@ -43,14 +45,16 @@ module fgout_module
! Setup routine that reads in the fixed grids data file and sets up the
! appropriate data structures

subroutine set_fgout(rest,fname)
subroutine set_fgout(rest,num_eqn,fname)

use amr_module, only: parmunit, tstart_thisrun
use utility_module, only: parse_values

implicit none

! Subroutine arguments
logical :: rest ! restart?
logical, intent(in) :: rest ! restart?
integer, intent(in) :: num_eqn
character(len=*), optional, intent(in) :: fname

! Local storage
Expand All @@ -59,6 +63,9 @@ subroutine set_fgout(rest,fname)
type(fgout_grid), pointer :: fg
real(kind=8) :: ts

integer :: n, iq
real(kind=8) :: values(num_eqn+2)
character(len=80) :: str

if (.not.module_setup) then

Expand Down Expand Up @@ -109,7 +116,25 @@ 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

allocate(fg%q_out_vars(num_eqn+2)) ! for q + eta, topo
do iq=1,num_eqn+2
fg%q_out_vars(iq) = .false.
enddo
read(unit,'(a)') str
call parse_values(str, n, values)
do k=1,n
iq = nint(values(k))
fg%q_out_vars(iq) = .true.
enddo
write(6,*) '+++ q_out_vars:', fg%q_out_vars

fg%num_vars = num_eqn + 3 ! also interp topo,eta,time
fg%nqout = 0 ! count how many are to be output
do k=1,num_eqn+2
if (fg%q_out_vars(k)) fg%nqout = fg%nqout + 1
enddo
!fg%nqout = fg%nqout + 1 ! for eta


! Initialize next_output_index
Expand Down Expand Up @@ -193,68 +218,15 @@ subroutine set_fgout(rest,fname)
print *,'y grid points my <= 0, set my >= 1'
endif


! 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
fg%bathy_index = num_eqn + 2
fg%eta_index = num_eqn + 1
fg%time_index = num_eqn + 3

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 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()

Expand Down Expand Up @@ -412,7 +384,7 @@ subroutine fgout_interp(fgrid_type,fgrid, &


! save the time this fgout point was computed:
fg_data(num_vars,ifg,jfg) = t
fg_data(fgrid%time_index,ifg,jfg) = t


endif ! if fgout point is on this grid
Expand Down Expand Up @@ -441,7 +413,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
! I/O
integer, parameter :: unit = 87
character(len=15) :: fg_filename
character(len=4) :: cfgno, cframeno
character(len=8) :: cfgno, cframeno
character(len=8) :: file_format
integer :: grid_number,ipos,idigit,out_number,columns
integer :: ifg,ifg1, iframe,iframe1
Expand Down Expand Up @@ -472,7 +444,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
"a10,' format'/,/)"

! Other locals
integer :: i,j,m,iq,k
integer :: i,j,m,iq,k,jq,num_eqn
real(kind=8) :: t0,tf,tau, qaug(10)
real(kind=8), allocatable :: qeta(:,:,:)
real(kind=4), allocatable :: qeta4(:,:,:)
Expand Down Expand Up @@ -546,76 +518,34 @@ subroutine fgout_write(fgrid,out_time,out_index)
endif
endif

! Output the conserved quantities and eta value
! Output the requested quantities and eta value:

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
num_eqn = size(fgrid%q_out_vars) - 2 ! last components eta,B
do jq=1,num_eqn+2
if (fgrid%q_out_vars(jq)) then
qeta(iq,i,j) = qaug(jq)
iq = iq+1
endif
endif
enddo

! now append eta value: (this is now included in loop above)
!qeta(iq,i,j) = qaug(fgrid%eta_index)
! NOTE: qaug(fgrid%bathy_index) contains topo B value

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
! convert fgrid%fgno to a string of length at least 4 with leading 0's
! e.g. '0001' or '12345':
write(cfgno, '(I0.4)') fgrid%fgno

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
! convert out_index to a string of length at least 4 with leading 0's
write(cframeno, '(I0.4)') out_index

fg_filename = 'fgout' // cfgno // '.q' // cframeno

fg_filename = 'fgout' // trim(cfgno) // '.q' // trim(cframeno)

open(unit,file=fg_filename,status='unknown',form='formatted')

Expand All @@ -637,14 +567,14 @@ subroutine fgout_write(fgrid,out_time,out_index)

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' // trim(cfgno) // '.b' // trim(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' // trim(cfgno) // '.b' // trim(cframeno)
open(unit=unit, file=fg_filename, status="unknown", &
access='stream')
allocate(qeta4(fgrid%nqout, fgrid%mx, fgrid%my))
Expand All @@ -671,7 +601,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
write(6,*) '*** should be ascii, binary32, or binary64'
endif

fg_filename = 'fgout' // cfgno // '.t' // cframeno
fg_filename = 'fgout' // trim(cfgno) // '.t' // trim(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, fgrid%nqout, 1, 0, 2, 0,file_format
Expand All @@ -680,10 +610,6 @@ subroutine fgout_write(fgrid,out_time,out_index)
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


Expand Down
2 changes: 1 addition & 1 deletion src/python/geoclaw/data.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
- GeoClawData
- RefinementData
- TopographyData
- FixedGridData
- FGoutData
- FGmaxData
- DTopoData
- QinitData
Expand Down
Loading
Loading