Skip to content

Fix missing actualstep in b4step calls #670

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

Merged
merged 5 commits into from
Sep 7, 2025
Merged

Fix missing actualstep in b4step calls #670

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
25852e4
Modify source to add and subract
mandli Sep 4, 2025
007e095
Update remaining source
mandli Sep 5, 2025
514fca5
Add fix for undefined tend_force_dry
mandli Sep 5, 2025
2a56549
Bug fixes to gauge code in multilayer
mandli Sep 5, 2025
ee024d7
Remove actualstep from stepgrid signatures
mandli Sep 5, 2025
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
41 changes: 20 additions & 21 deletions examples/multi-layer/bowl-radial/setrun.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,8 +87,8 @@ def setrun(claw_pkg='geoclaw'):
# Index of aux array corresponding to capacity function, if there is one:
clawdata.capa_index = 0



# -------------
# Initial time:
# -------------
Expand All @@ -98,7 +98,7 @@ def setrun(claw_pkg='geoclaw'):

# 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
# restart_file 'fort.chkNNNNN' specified below should be in
# the OUTDIR indicated in Makefile.

clawdata.restart = False # True to restart from prior results
Expand Down Expand Up @@ -129,7 +129,7 @@ def setrun(claw_pkg='geoclaw'):
clawdata.output_step_interval = 1
clawdata.total_steps = 500
clawdata.output_t0 = True


clawdata.output_format = 'ascii' # 'ascii', 'binary32' or 'binary64'

Expand Down Expand Up @@ -185,20 +185,20 @@ def setrun(claw_pkg='geoclaw'):

# 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

# 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 = 6
# List of limiters to use for each wave family:

# List of limiters to use for each wave family:
# Required: len(limiter) == num_waves
# Some options:
# 0 or 'none' ==> no limiter (Lax-Wendroff)
Expand All @@ -209,10 +209,10 @@ def setrun(claw_pkg='geoclaw'):
clawdata.limiter = ['mc', 'mc', 'mc', 'mc', 'mc', 'mc']

clawdata.use_fwaves = True # True ==> use f-wave version of algorithms

# 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 == 1 or 'godunov' ==> Godunov (1st order) splitting used,
# src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended.
clawdata.source_split = 'godunov'

Expand Down Expand Up @@ -250,7 +250,7 @@ def setrun(claw_pkg='geoclaw'):
pass

elif np.abs(clawdata.checkpt_style) == 2:
# Specify a list of checkpoint times.
# Specify a list of checkpoint times.
clawdata.checkpt_times = [0.1,0.15]

elif np.abs(clawdata.checkpt_style) == 3:
Expand Down Expand Up @@ -295,10 +295,10 @@ def setrun(claw_pkg='geoclaw'):
amrdata.clustering_cutoff = 0.700000

# print info about each regridding up to this level:
amrdata.verbosity_regrid = 0
amrdata.verbosity_regrid = 0


# ----- For developers -----
# ----- For developers -----
# Toggle debugging print statements:
amrdata.dprint = False # print domain flags
amrdata.eprint = False # print err est flags
Expand All @@ -310,7 +310,7 @@ def setrun(claw_pkg='geoclaw'):
amrdata.sprint = False # space/memory output
amrdata.tprint = False # time step reporting each level
amrdata.uprint = False # update/upbnd reporting

# More AMR parameters can be set -- see the defaults in pyclaw/data.py

# == setregions.data values ==
Expand Down Expand Up @@ -343,7 +343,7 @@ def setrun(claw_pkg='geoclaw'):
x = (r + .001) / np.sqrt(2.)
y = (r + .001) / np.sqrt(2.)
rundata.gaugedata.gauges.append([gaugeno, x, y, 0., 1e10])


return rundata
# end of function setrun
Expand All @@ -364,7 +364,7 @@ def setgeo(rundata):
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
Expand Down Expand Up @@ -406,7 +406,7 @@ def setgeo(rundata):

# == fgout grids ==
# new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated)
# set rundata.fgout_data.fgout_grids to be a
# set rundata.fgout_data.fgout_grids to be a
# list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid:
#rundata.fgout_data.fgout_grids = []

Expand All @@ -422,11 +422,11 @@ def set_multilayer(rundata):
data.num_layers = 2
data.eta = [0.0, -20]
data.layer_index = 1

# Algorithm parameters
data.eigen_method = 2
data.inundation_method = 2
data.richardson_tolerance = np.infty
data.richardson_tolerance = np.inf
data.wave_tolerance = [1e-3,1e-2]

rundata.qinit_data.qinit_type = 4
Expand All @@ -441,4 +441,3 @@ def set_multilayer(rundata):
import sys
rundata = setrun(*sys.argv[1:])
rundata.write()

58 changes: 29 additions & 29 deletions src/1d_classic/shallow/claw1.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c ==============================================================
c
c Solves a hyperbolic system of conservation laws in one space dimension
c of the general form
c of the general form
c
c capa * q_t + A q_x = psi
c
Expand All @@ -24,7 +24,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c
c The user must supply the following subroutines:
c
c bc1, rp1 subroutines specifying the boundary conditions and
c bc1, rp1 subroutines specifying the boundary conditions and
c Riemann solver.
c These are described in greater detail below.
c
Expand All @@ -44,7 +44,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c These routines must be declared EXTERNAL in the main program.
c For description of the calling sequences, see below.
c
c Dummy routines b4step1.f and src1.f are available in
c Dummy routines b4step1.f and src1.f are available in
c claw/clawpack/1d/lib
c
c
Expand Down Expand Up @@ -74,12 +74,12 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c physical domain. In addition there are mbc grid cells
c along each edge of the grid that are used for boundary
c conditions.
c
c q(meqn, 1-mbc:mx+mbc)
c
c q(meqn, 1-mbc:mx+mbc)
c On input: initial data at time tstart.
c On output: final solution at time tend.
c q(m,i) = value of mth component in the i'th cell.
c Values within the physical domain are in q(m,i)
c Values within the physical domain are in q(m,i)
c for i = 1,2,...,mx
c mbc extra cells on each end are needed for boundary conditions
c as specified in the routine bc1.
Expand All @@ -102,15 +102,15 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c stored in aux(mcapa,i).
c In this case we require method(7).ge.mcapa.
c
c dx = grid spacing in x.
c dx = grid spacing in x.
c (for a computation in ax <= x <= bx, set dx = (bx-ax)/mx.)
c
c tstart = initial time.
c
c tend = Desired final time (on input).
c If tend<tstart, then claw1 returns after a single successful
c time step has been taken (single-step mode).
c Otherwise, as many steps are taken as needed to reach tend,
c Otherwise, as many steps are taken as needed to reach tend,
c up to a maximum of nv(1).
c = Actual time reached (on output).
c
Expand Down Expand Up @@ -152,13 +152,13 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c comes from the previous step, the Courant number will
c not in general be exactly equal to the desired value
c If the actual Courant number in the next step is
c greater than 1, then this step is redone with a
c greater than 1, then this step is redone with a
c smaller dt.
c
c method(2) = 1 if Godunov's method is to be used, with no 2nd order
c corrections.
c = 2 if second order correction terms are to be added, with
c a flux limiter as specified by mthlim.
c a flux limiter as specified by mthlim.
c
c method(3) is not used in one-dimension.
c
Expand Down Expand Up @@ -223,7 +223,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c If mwork is too small then the program returns with info = 4
c and prints the necessary value of mwork to unit 6.
c
c
c
c info = output value yielding error information:
c = 0 if normal return.
c = 1 if mbc.lt.2
Expand All @@ -241,7 +241,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c User-supplied subroutines
c -------------------------
c
c bc1 = subroutine that specifies the boundary conditions.
c bc1 = subroutine that specifies the boundary conditions.
c This subroutine should extend the values of q from cells
c 1:mx to the mbc ghost cells along each edge of the domain.
c
Expand Down Expand Up @@ -290,7 +290,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c are passed in using auxl and auxr. Again, in the standard routines
c auxl=auxr=aux in the call to rp1.
c
c On output,
c On output,
c wave(m,mw,i) is the m'th component of the jump across
c wave number mw in the ith Riemann problem.
c s(mw,i) is the wave speed of wave number mw in the
Expand All @@ -307,7 +307,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c into the specification of amdq and apdq.
c
c src1 = subroutine for the source terms that solves the equation
c capa * q_t = psi
c capa * q_t = psi
c over time dt.
c
c If method(5)=0 then the equation does not contain a source
Expand Down Expand Up @@ -338,7 +338,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c other tasks which must be done every time step.
c
c The form of this subroutine is
c
c
c -------------------------------------------------
c subroutine b4step1(mbc,mx,meqn,q,xlower,dx,time,dt,maux,aux)
c implicit double precision (a-h,o-z)
Expand All @@ -353,12 +353,12 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c
c Copyright 1994 -- 2002 R. J. LeVeque
c
c This software is made available for research and instructional use only.
c This software is made available for research and instructional use only.
c You may copy and use this software without charge for these non-commercial
c purposes, provided that the copyright notice and associated text is
c reproduced on all copies. For all other uses (including distribution of
c modified versions), please contact the author at the address given below.
c
c modified versions), please contact the author at the address given below.
c
c *** This software is made available "as is" without any assurance that it
c *** will work for your purposes. The software may in fact have defects, so
c *** use the software at your own risk.
Expand All @@ -370,18 +370,18 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c Author: Randall J. LeVeque
c Applied Mathematics
c Box 352420
c University of Washington,
c University of Washington,
c Seattle, WA 98195-2420
c [email protected]
c =========================================================================
c
c
c
c
c
c ======================================================================
c Beginning of claw1 code
c ======================================================================
c
c

use gauges_module, only: num_gauges, update_gauges,
& print_gauges_and_reset_nextLoc
Expand Down Expand Up @@ -435,7 +435,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
i0wave = i0f + (mx + 2*mbc) * meqn
i0s = i0wave + (mx + 2*mbc) * meqn * mwaves
i0dtdx = i0s + (mx + 2*mbc) * mwaves
i0qwork = i0dtdx + (mx + 2*mbc)
i0qwork = i0dtdx + (mx + 2*mbc)
i0amdq = i0qwork + (mx + 2*mbc) * meqn
i0apdq = i0amdq + (mx + 2*mbc) * meqn
i0dtdx = i0apdq + (mx + 2*mbc) * meqn
Expand Down Expand Up @@ -463,7 +463,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
c # save old q in case we need to retake step with smaller dt:
call copyq1(meqn,mbc,mx,q,work(i0qwork))
endif
c
c
40 continue
dt2 = dt / 2.d0
thalf = t + dt2 !# midpoint in time for Strang splitting
Expand Down Expand Up @@ -532,8 +532,8 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
endif
dtmin = dmin1(dt,dtmin)
dtmax = dmax1(dt,dtmax)


else
dt = dtv(2)
endif
Expand Down Expand Up @@ -563,13 +563,13 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
go to 900
endif
endif

if (.not. topo_finalized) then
! modify topo using dtopo arrays:
call topo_update(t)
call setaux(mbc,mx,xlower,dx,maux,aux)
endif

c
c # see if we are done:
nv(2) = nv(2) + 1
Expand All @@ -578,7 +578,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
100 continue
c
900 continue
c
c
c # return information
c
if (method(1).eq.1 .and. t.lt.tend .and. nv(2) .eq. maxn) then
Expand All @@ -602,5 +602,5 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
call print_gauges_and_reset_nextLoc(ii,meqn)
end do

return
return
end
Loading
Loading