Re: [Anuga-user] Modeled discharge doesn't match input discharge?

[Kyle A W. <kyl...@ut...>]

Hi Gareth,

Thanks for the prompt reply, this is helpful information. I think my
estimates of discharge from the raster reflects your second point, as it
not only differs from the input discharge but also varies quite a bit
spatially around the ~780 m3/s mean value, depending on where across the
channel I measure it. These fluctuations look to be a bit larger than I
would expect from typical interpolation errors, so that helps clear up the
cause of that a bit.

To expand on the forcing issue, my reason for wanting to use a BC rather
than an inlet is that I'm trying to reproduce the conditions at a specific
gauge station co-located with my boundary, which gives me a time-series for
both discharge and stage. Ideally, I would like to force this boundary to
behave like that gauge station, which I'm not sure how to achieve if I
can't also prescribe a stage. Perhaps more importantly, the model will
eventually include the effects of tides, and I worry that having a
reflective boundary near the inlet will cause reflections of tidal
constituents to bounce back into the domain. I could maybe create an
artificial upstream buffer for the tides, but I'm not sure how best to do
that.

Best regards,
- Kyle

On Mon, Apr 1, 2019 at 7:07 PM Gareth Davies <gar...@gm...>
wrote:

>
> Hi Kyle,
>
> There are 2 issues to consider: How to force the boundary, and how to
> back-calculate the fluxes.
>
> ## Forcing the boundary ##
>
> Can you explain this a bit more?:
> **For physical reasons, I would prefer to prescribe this discharge using a
> Dirichlet boundary condition (rather than an inlet + reflective BC).**
>
> Without some more context, I'd think the best approach would be "inlet +
> reflective". This has the advantage of ensuring the mass inflow is exactly
> as desired.
>
> In ANUGA, it is hard to control the boundary fluxes using boundary
> conditions, because the boundary condition is used to compute the "boundary
> edge flow state as viewed from outside the domain". The latter is combined
> with the interior flow state to compute the flux (with the approximate
> Riemann solver). Thus with boundary conditions, the flux depends on the
> flow inside the model, as well as the boundary condition, and so is hard to
> control.
>
>
> ## Back-calculate fluxes ##
>
> A discharge measurement based on the raster (like in your snippet below)
> will be inaccurate to some extent. ANUGA outputs flow quantities at
> centroids, but it is difficult to directly access the fluxes. The variables
> like "uh/vh" at centroids will not be the same as the mass fluxes at edges,
> even at steady state, because the edge fluxes are computed with the
> approximate Riemann solver. This is separate to any "raster interpolation
> error", which might also be substantial if the channel is small.
>
> To directly measure the flux when debugging the model, one approach is to
> cut a "pit" somewhere in your domain just downstream of the channel, to
> catch all the flow. Then check how quickly the pit fills.
>
> Cheers,
> Gareth.
>
>
> On 2/4/19 9:39 am, Kyle A Wright wrote:
>
> Hi,
>
> I am currently using Anuga to model a river delta, and need a discharge
> inflow at the upstream boundary. For physical reasons, I would prefer to
> prescribe this discharge using a Dirichlet boundary condition (rather than
> an inlet + reflective BC). I have tried multiple boundary conditions for
> this input, including: the default Dirichlet, the
> Dirichlet_discharge_boundary, the Inflow_boundary, and one additional BC
> module I wrote that uses discharge and stage as inputs. I have tested all
> of these on a simplified case, and with minor differences, all appear to
> work fine.
>
> Here is my problem: when I use these BCs on my full domain, the discharge
> I am measuring across the upstream channel (~780 m3/s) is much lower than
> the discharge I am prescribing at the boundary (1620 m3/s). I do not know
> if this discrepancy is due to the BC, to the model evolution, or the data
> post-processing. Has anybody else encountered this problem, and does
> anybody know what the issue is?
>
> For context, I have even tried multiplying the inflow by the fractional
> difference in the model results (2.09x), and the resulting discharge was
> still too low (~1300 m3/s). Aside from the discharge, everything about the
> model evolution appears normal.
>
> A simplified version of the code/workflow, if needed:
>
> ### Compute inflow:
> Q = 1620 # [m3/s]
> upstream_stage = 0.74 # [m]
> z = pd.read_csv('channel_profile.csv', names=['elev']) # reads in
> 1-meter-res csv of elevations along the boundary
> z = np.array(z['elev'])
> profile = -z[z <= upstream_stage] + upstream_stage # to find channel area
> mean_depth = np.mean(profile)
> channel_area = np.trapz(profile)
> U = Q / channel_area
> momentum_in = mean_depth*U
>
> ### [establish domain, mesh, topography, etc. Manning's is n=0.03 in
> channel]
>
> ### Inflow BC is one of the following:
> Bin = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain,
> stage0 = upstream_stage, wh0 = momentum_in) # my default choice
> Bin = anuga.Dirichlet_boundary([upstream_stage, x_momentum_in,
> y_momentum_in]) # equivalent to above, requires computing x,y components
> Bin = anuga.shallow_water.boundaries.Inflow_boundary(domain, rate = Q) #
> requires changing the mesh boundary & gives unrealistic depths, but still
> results in a discharge
> Bin = anuga.shallow_water.boundaries.Dirichlet_Q_stage_boundary(domain,
> rate = Q, stage = upstream_stage)  # I wrote this to work like an
> Inflow_boundary, but uses a prescribed stage instead of a calculated one,
> and only applies it over BC edges with elevation < upstream_stage
>
> domain.set_boundary({'upstream: Bin, 'bay': Bout, 'sides': Br})
>
> ### [run model, save SWW]
>
> ### Post-processing:
> anuga.utilities.plot_utils.Make_Geotif(
>      swwFile = swwname,
>      output_quantities = ['depth', 'velocity', 'depthIntegratedVelocity'],
>      myTimeStep = 100,
>      CellSize = 10.0, # cell size of raster in meters
>      proj4string = '+proj=utm +zone=15 +datum=NAVD88 +units=m +no_defs',
>      k_nearest_neighbours = 3) # number of nearest neighbors for
> interpolation (I've tried = 1)
>
> ### [load in raster, sum depthIntegratedVelocity across transects of the
> upstream channel, multiply by cell size. Resulting Q is too low]
>
> Any help or insight would be greatly appreciated.
> -- Kyle
>
>
> *----------------------------------------------------------------------------------------*
> *Kyle Wright*
> The University of Texas at Austin
> Department of Environmental & Water Resources Engineering
>
>
> _______________________________________________
> Anuga-user mailing lis...@li...://lists.sourceforge.net/lists/listinfo/anuga-user
>
>
>

-- 
*----------------------------------------------------------------------------------------*
*Kyle Wright, EIT*
The University of Texas at Austin
Department of Environmental & Water Resources Engineering
(512) 712 - 8688